Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

Radiative transport in large arteries

Open Access Open Access

Abstract

A refined model for the photon energy distribution in a living artery is established by solving the radiative transfer equation in a cylindrical geometry, using the Monte Carlo method. Combining this model with the most recent experimental values for the optical properties of flowing blood and the biomechanics of a blood-filled artery subject to a pulsatile pressure, we find that the optical intensity transmitted through large arteries decreases linearly with increasing arterial distension. This finding provides a solid theoretical foundation for measuring photoplethysmograms.

© 2013 Optical Society of America

1. Introduction

Many biomedical applications take advantage of the interaction of light and living tissue, as for example spectroscopy, tomography or endoscopy [1]. Among all these photonic sensing principles, absorption spectroscopy in pulse oximetry is the most widely used concept in daily health care.

Standard pulse oximeters measure a pulsatile signal either in transmission or reflection mode at the finger tip or earlobe. In clinical practice, this pulsatile signal is referred to as a photoplethysmogram (PPG). For measuring PPGs, two physiologically relevant factors could be identified: as a consequence of the pulsatile nature of our cardiovascular system, both (i) blood opacity and (ii) arterial geometry alter during a cardiac cycle [3, 4]. On the one hand, arteries distend when subjected to an intra-arterial pressure, leading to changes in blood volume and blood vessel wall movement and therewith to changes in the overall geometry of the artery [5]. On the other hand, blood is a complex suspension of cells in plasma, capable of compensating pressure gradients by varying the morphology and intercellular organization of the red blood cells, which in turn affects the transparency of blood [69]. However, until today it is not fully understood, how these two physiological factors impact the measured PPG quantitatively [2].

In this theoretical study, we investigate the radiative transport inside arteries and quantify the impact of both of the above factors (changes in blood opacity and arterial geometry) on PPGs. The simulation used throughout this study is based on the work of Welch et al. and Zheng et al., who investigated the interaction of light and tissue by use of Monte Carlo methods for more than two decades [1012]. Their simulation, called MCML (Monte Carlo Modeling of Light Transport in Multi-layered Tissues in Standard C), models the radiative transfer through a stack of infinite planes. By extending their simulation code to multilayered cylindrical geometries and by considering the biophysics of an artery subject to a pulsatile pressure, we complete the picture of photon energy distribution in large, blood-filled, distending arteries. We use the linear viscoelastic model proposed by [13] to describe the biomechanics of the arterial wall, while optical parameters of flowing blood measured by Friebel et al. [8] and Enejder et al. [9] support our theoretical studies.

2. Approach

The first part of this section introduces the MCML simulation developed by Wang et al. and the basic theory necessary to understand light propagation in turbid media. The second part points out how we modified the MCML simulation to describe radiative transfer in arteries in a more realistic scenario and results in an extended source code called extMCML. The last part of this section validates extMCML by comparing with MCML.

2.1. MCML

MCML is a simulation tool to solve the radiative transfer equation (RTE) in a stack of infinite planes by interpreting the underlying physical processes of absorption, scattering, refraction, and reflection as stochastic processes and applying the Monte Carlo method to obtain a solution. The geometrical description of the underlying models is illustrated in Fig. 1.

 figure: Fig. 1

Fig. 1 Comparison of the simulation geometry between a) MCML and b) extMCML. a) In MCML, a stack of infinitely extended planes is illuminated by a Dirac beam in the z-direction. b) In extMCML, the Dirac beam strikes an infinitely long multilayered cylinder in the lateral direction.

Download Full Size | PDF

2.1.1. Light propagation in turbid media

The radiative transfer equation,

(s)Iν(r,s)+(μa+μs)μtIν(r,s)=μt4π4πIν(r,s)ps(s,s))dΩ,
describes how monochromatic radiance Iv gains and loses energy due to boundary effects, scattering and absorption inside an observer volume element dV, while propagating inside turbid media [14].

The 1st term on the left side of Eq. (1) considers the radiance that is redirected from direction s into another direction due to transmission or reflection at the boundaries of dV. The 2nd term on the left side describes the loss of radiance in direction s due to absorption inside the volume element dV or redirection of radiance due to a scattering event inside dV. The extinction coefficient μt is the sum of the absorption coefficient μa and the scattering coefficient μs. The right hand side reflects the gain of radiance by a scattering event inside the volume dV, where a photon with initial direction s′ is scattered into direction s with a probability defined by the scattering phase function ps(s′, s).

Since the simulation developed here uses the optical parameters determined experimentally by other groups [8, 9], we followed their use of the Henyey-Greenstein phase function as a probability density function which describes the redirection of radiance from direction s′ into direction s. By assuming a spherically shaped scatterer, the Henyey-Greenstein scattering phase function is only a function of the angle θ between the directions s and s′, defined as

ps(θ)=14π1g2(1+g22gcosθ)3/2
[15]. The anisotropy factor g parametrizes the scattering direction by a scalar and is calculated as the mean cosine of the scattering phase function
g=0πps(θ)cosθ2πsinθdθ.
For example, if g = 1, there is only forward scattering; in case of g = 0, we have isotropic scattering; and for g = −1, backward scattering only. The absorption coefficient μa, the scattering coefficient μs, the anisotropy factor g, and the refractive index n are referred to as intrinsic optical parameters.

2.2. extMCML

To investigate light propagation in large arteries in a realistic scenario, we developed a simulation tool which models radiative transfer inside a multilayered cylinder. Compared to the simulation geometry used in MCML, extMCML allows the Dirac beam to strike a multilayered cylinder in the lateral direction, as seen in Fig. 1(b). extMCML extends the very comprehensively written MCML simulation to allow multi-layered cylindrical boundaries and the calculation of Fresnel reflections at these cylindrical boundaries.

In accordance with the MCML simulation, photon packet movement is interrupted when it crosses an optical interface. To calculate the Fresnel equations, the step size s is reduced and the photon packet is moved forward with the reduced step size s*, until it reaches the optical interface. The reflection and transmission coefficients are then calculated in the plane of incidence at the optical interface using Snell’s law and the corresponding refractive indices ni and nt. Subsequently, a Monte-Carlo step decides if the photon packet is reflected or transmitted (see [11]). In case of reflection, the new direction is

pr=2(pir)pi+pi,
where pi is the incident photon packet direction and r is the direction of the radius vector normal to the optical interface. In case of transmission, the new direction is given by
pt=nintpi+(1nint1(pir)2nipir)r
[16]. The position of the reflected or transmitted photon packet is calculated by combining the new direction with the remaining step size ss*.

The physical quantities calculated by the simulation are the absorption probability density A(x, y, z) and the radiative flux Φ(x, y, z). The absorption probability describes the photon energy being absorbed in a certain volume element, while the radiative flux is the radiance through a certain volume element integrated over all solid angles.

The MCML simulation exploits a rotational symmetry given by the infinite plane geometry (recall Fig. 1(a)), to calculate the absorption probability density and the radiative flux. The cylindrical geometry used here does not show this symmetry, hence we calculate the absorption probability density

A(x,y,z)=NphΔWiNphdV
in Cartesian coordinates. Here, Nph is the number of photons simulated and ΔWi is the energy deposited by each photon in the volume element dV = dxdydz located at (x,y,z). The radiative flux is then given as
Φ(x,y,z)=A(x,y,z)μa.

The radiative flux describes the distribution of photons inside the simulation geometry. To describe the photon packages that exit the simulation geometry into a specific observer region located on the surface of the cylinder, we introduce the observed diffuse reflectance and the observed total transmittance. The observed diffuse reflectance is calculated by the photons back reflected into a specific observer region at the entrance side directly around the light source, while the observed total transmittance is calculated diametrically opposed to it. The observer region is the detector area in Fig. 2.

 figure: Fig. 2

Fig. 2 Simulations by extMCML and MCML. The left plot shows the radiative flux obtained by MCML, the plot on the right presents the radiative flux determined by extMCML. The two different detectors are sketched in red, both detector areas are 1 mm2.

Download Full Size | PDF

2.3. Validation of the source code amendments

To validate the code amendments, we compared the results of the original MCML simulation with our extMCML simulation. The quantities used to compare the two simulations are the observed diffuse reflectance and observed total transmittance. To make both simulation geometries comparable, the observed quantities have to be integrated over a small detector area relative to the size of the simulation volume. This ensures a planar approximation of the cylindrical curvature in case of extMCML and makes both simulation results comparable. We choose the diameter of the cylinder d in extMCML and the width of the simulation geometry w in case of MCML to be 1 cm. In order to fulfill the requirement of planar approximation, we set the detector area to be 1 mm2.

The chosen simulation parameters are presented in Table 1.

Tables Icon

Table 1. Optical parameters to compare extMCML with the original MCML code.

By introducing a step in the relative refractive index between turbid medium and ambient medium, we validate both code amendments, the change in simulation geometry as well as the Fresnel reflections off the cylindrical boundaries. The small value of the absorption coefficient and the high anisotropy factor ensure that enough photons reach the small detector area across the huge simulation geometry to obtain a statistically significant value for observed total transmittance. To model a turbid medium we choose the scattering coefficient to be approximately two orders of magnitude higher than the absorption coefficient.

A total of 10 simulations were performed for both source codes. Figure 2 shows the radiative flux for MCML and extMCML given by one representative of these 10 simulations. In both cases, 108 photons were simulated to yield a statistically significant result for the observed total transmittance and observed diffuse reflectance. The results of mean observed total transmittance, mean observed diffuse reflectance, and the corresponding standard deviations over all 10 simulations are given in Table 2. In both cases, the input energy is normalized to 1 W.

Tables Icon

Table 2. Verification of the extMCML code by comparison of total transmittance and diffuse reflectance with the MCML code provided by [11]. Input power is 1 W, detector area is 1 mm2.

2.4. Discussion of the code extensions

As is seen from the data of Table 2, the values of observed diffuse reflectance and observed total transmittance are close together for the two simulation codes. The relative error between both simulation source codes in case of observed total transmittance is 4.36 %, while the relative error in case of observed diffuse reflectance is 4.04 %. The small relative errors verify our code amendments. However, the corresponding standard deviations indicate a significant deviation between both values. This deviation can be explained by the different geometry close to the detector area in both source codes. In case of MCML, the region surrounding the detector area in lateral direction is turbid. Hence, there is a likelihood for photons to scatter from the surrounding region onto the detector area. In contrast, in case of extMCML, the detector area is surrounded partially by ambient medium (see black area surrounding the body of the cylinder in Fig. 2). Photons in ambient medium are treated as lost, hence there is no possibility for these photons to scatter from the ambient medium onto the detector area in case of extMCML. This causes the small but significant deviations between the observed quantities calculated from the two source codes.

To prove our hypothesis, we decreased the anisotropy factor from g = 0.9 to g = 0.75 and g = 0.5, while the other simulation parameters have been retained unchanged. The smaller values in the anisotropy factor increase the likelihood of a photon to scatter onto the detector area from the region nearby in case of MCML. Table 3 gives the relative error between MCML and extMCML for the various anisotropy factors.

Tables Icon

Table 3. Relative error between MCML and extMCML for various anisotropy factors.

As expected, the relative error between both simulations increases with decreasing anisotropy factor confirming our hypothesis. The significant deviation between both source codes clearly show that the code modifications we made to MCML are necessary to investigate radiative transfer in cylindrical geometries as for example arteries, esophagus, intestine, and bones in a more realistic scenario. Moreover, as the detector area increases relative to the simulation volume, geometry related effects become even more prominent.

3. Light propagation in arteries

We now apply the validated extMCML simulation to investigate radiative transport in arteries. In a first step, we define the geometry of the model artery. In a second step, we introduce the latest intrinsic optical parameters of arterial tissue and blood and study the radiative transport in a static artery. In the last step, we include the dynamics of the artery over a cardiac cycle and study the physiological origins of photoplethysmograms at large arteries.

3.1. Model artery

3.1.1. Geometry

Figure 3 describes the model artery used throughout our simulation with an inner radius of ri = 1.9 mm and an outer radius of ro = 2.6 mm [17]. The length of the cylinder and hence the investigated arterial segment is l = 5 mm. The length of the arterial segment was chosen as a trade-off between the size of the entire simulation volume and the drop-off in photon energy in the axial direction, as will be seen in Fig. 4(b). The arterial wall of the model artery consists of three layers. The individual layers are named based on their location in the radial direction, namely intima, media and adventitia with thicknesses of ti = 100 μm, tm = 400 μm, and ta = 200 μm, respectively [18].

 figure: Fig. 3

Fig. 3 Simulation geometry used to calculate the photon energy distribution inside a static, large artery. The 3D plot of the simulation geometry introduces the thicknesses of the various layers forming the arterial wall, adventitia ta, media tm and intima ti; the length of the investigated arterial segment l; inner and outer radii of the arterial wall ri and ro; the Cartesian coordinate system x,y,z; and the Dirac beam striking the arterial wall in z direction.

Download Full Size | PDF

 figure: Fig. 4

Fig. 4 Decay of the radiative flux inside the model artery at rest. left: normalized radiative flux at y = 0μm as a function of z calculated by integrating along the x direction. The blood layer is marked between the two arterial walls. The three different arterial layers forming the arterial wall are labeled as adventitia ta, media tm, and intima ti (dashed, black). right: Normalized radiative flux in x direction integrated over z at y = 0μm. Since there are no optical boundaries in direction of x, the plot shows no discontinuities. The insets in both graphs show the location of the plotted radiative flux inside the cylinder.

Download Full Size | PDF

3.1.2. Intrinsic optical parameters

The intrinsic optical parameters, at a wavelength of λ = 633 nm, for the tissue of the three different arterial layers forming the arterial wall are listed in Table 4 [8,19,20]. To complete the picture of a blood-filled artery, the corresponding values for blood at the same wavelength are also given.

Tables Icon

Table 4. Optical properties [19,20] of the three different layers of the arterial wall and blood at a wavelength of 633 nm.

3.2. Radiative transport in a static artery

The left plot in Fig. 4 shows the decay of the normalized radiative flux integrated along the x direction directly at the center (y = 0μm). As can be seen in the plot, the radiative flux does not decay exponentially across the arterial diameter. The steeper decay of the radiative flux inside the arterial wall shows that most photons are absorbed by the arterial wall. Additionally, small fluence dips are visible between the different arterial wall layers that arise from Fresnel reflections at these optical boundaries [11].

The right graph in Fig. 4 presents the decay of the radiative flux averaged along the z direction at y = 0μm. The radiative flux in x decays twice as fast as in the z direction. The reason for this is that biological materials are mostly forward scattering media.

The radiative flux at x = 2.5 mm has decreased by 5 orders of magnitude from its maximum value. Therefore, the probability of a photon at x > 2.5 mm contributing to the transmission mode photoplethysmogram (see Section 3.3.4 below), is negligible. Hence, the chosen length of the simulated cylinder of 5 mm is reasonable.

3.3. Dynamic artery

We now make use of extMCML to investigate the behavior of PPGs in large arteries. By introducing the dynamics of the artery in a cardiac cycle, we separate two effects that cause changes in optical transmittance in an artery during a cardiac cycle: (i) arterial distension and (ii) variations of optical transmittance of blood for various flow conditions. The transmission PPGs were evaluated by integrating the diametrically transmitted photons of an infinitely narrow beam, known as Dirac beam, over a detector area of 1 mm2, as depicted in Fig. 5. Since the dynamics of the cardiovascular system is much slower than the velocity of light inside turbid media, we simulated 11 steady-state light distributions with strains ranging from 0% ≤ ε ≤ 5.14 % with a step of Δε ≈ 0.5 % to account for the distension of the model artery over an entire cardiovascular cycle. The maximum strain of εmax = 5.14 % corresponds to a maximum distension in arterial diameter of approximately 270μm [17]. The upper part of Fig. 5 shows examples of two photon energy distributions in the yz-plane at two points at times t0 and t1 to illustrate the approach for calculating the photoplethysmograms in large distending arteries. The lower part of Fig. 5 shows the exact inner (red, dashed) and outer (red, dotted) arterial wall distensions.

 figure: Fig. 5

Fig. 5 Schematic of the arrangement for calculating the photoplethysmograms in large, distending arteries. The upper part shows two photon energy distributions in a cross-sectional view for two different points in time t0 and t1. The cross-sectional view at t0 is to scale while the cross-sectional view at t1 is magnified to point out the change in diameter. The exact change in diameter is presented in the lower part of this figure. Here, the temporal behavior of the inner (red, dashed) and outer (red, dotted) arterial wall distensions ur are plotted. The dotted vertical lines indicate the two states above, where t0 and t1 represent the arterial state at the diastole and systole, respectively.

Download Full Size | PDF

By simulating all steady state light distributions for different optical parameters of flowing blood, we can separate their influence on transmission photoplethysmograms in large arteries. The parameter which describes these different flow conditions is the wall shear rate WSR, which will be discussed in Subsection 3.3.3 below. To obtain the optical parameters of flowing blood, we linearly interpolated the experimental data in steps of ΔWSR = 100 s−1. The interpolation is based on the work of Friebel et al. [8] and Enejder et al. [9], who measured the intrinsic optical parameters of blood as a function of the applied wall shear rate. The range of the wall shear rates was chosen to be 200 s−1 to 1000 s−1, as they occur in large arteries, as for example in the common carotid artery [21, 22].

For the entire survey we performed 188 simulations. The lateral and axial resolution was set to 10μm to resolve the small arterial distensions of about 270μm. In total, the simulation geometry was built up by about 150 ×106 voxels. For each simulation, we used Nph = 1 × 108 photons to achieve a statistically significant result. Because of the large number of voxels and the fact that Monte-Carlo simulations are computationally intensive, we performed all simulations on the Black Forest Grid at the University of Freiburg.

3.3.1. Biomechanics of an artery

With every heart contraction, a pressure wave starts propagating through the arterial system. This pressure wave induces a pressure gradient which acts over a certain arterial segment. As a consequence, an axial vx(r, x, t) and radial vr(r, x, t) directed flow develops. The radial flow causes the arterial wall to distend in radial ur(r, x, t) direction, while the axial flow leads to an axial ux(r, x, t) movement of the arterial wall.

From the extensive literature on this subject, the authors propose [2325] for a complete understanding of the biomechanics of arteries. Here, we follow the solution given by Cox [13], since his solution pays attention to the major physical principles behind arterial blood flow. Although sophisticated, he derives an analytic expression for the desired waveforms describing the distension of the inner and outer arterial wall.

The three parameter model used by Cox to describe the linear viscoelastic dynamics of the arterial wall is given as

μ*=μ01+iωλ21+iωλ1,
where μ* denotes the complex modulus of rigidity, μ0 is the static modulus of rigidity (Young’s modulus) and λ1 and λ2 represent the relaxation and retardation time constants.

The intra-arterial pressure p(r, x, t) is assumed as a finite superposition of harmonic waves

p(r,x,t)=nNPn(r)exp[i(nωtγnx)],
where n is the harmonic number, ω the angular frequency, γn the propagation constant of the pressure wave and Pn(r) the radially varying amplitude of the nth harmonic. The highest harmonic is denoted as N.

According to Cox, the radial displacement of the arterial wall can be written as

ur=NnγnPnnωρ(Δ13Δ11J1(knr)+Δ14Δ11Y1(knr)+Δ15Δ11J1(iγnr)+Δ16Δ11Y1(iγnr))exp[i(nωtγnx)],
where ρ denotes the density of blood, Pn is the averaged intra-arterial pressure obtained by integrating Pn(r) over the arterial lumen, Δij are cofactors obtained from eliminating the ith row and jth column in the frequency matrix as defined by Cox [13]. J1 and Y1 denote the Bessel functions of the first and second kind, respectively. For the sake of readability, the parameter kn collects several constants and is given as
kn2=n2ω2ρwμ*γn2,
where ρw is the mass density of the arterial wall.

3.3.2. Biomechanical parameters

In vivo experiments show that the maximum arterial strain decreases with age and varies in the range of 5 % ≤ εmax ≤ 10 % [17]. The four parameters defining the biomechanics of our model artery, given in Table 5, were modified to yield outer arterial wall strains in the range of roughly 5 %, which corresponds to an arterial displacement of elderly people. Figure 6 shows the measured intra-arterial pressure waveform (blue) and the calculated inner (red, dashed) and outer (red, dotted) arterial wall displacement, according to Equation 10. The pressure waveform used in this approach was measured in the common carotid artery of a domestic pig by an intra-arterial tip catheter previously [26].

Tables Icon

Table 5. Biomechanical parameters of the arterial wall defining the linear viscoelastic behavior given in Equation 8 optimized to yield arterial strains of roughly 5 %.

 figure: Fig. 6

Fig. 6 Measured pressure waveform p(t) and calculated inner (red, dashed) and outer (red, dotted) arterial displacements. The displacement waveforms are calculated according to Equation 10 with respect to the biomechanical parameters of Table 5.

Download Full Size | PDF

The viscosity of blood was taken to be 3.3 cP and the density of blood to be 1.056 gcm−3 [13].

3.3.3. Intrinsic optical parameters of flowing blood

The intrinsic optical parameters of blood have been investigated by several authors with respect to various physiological conditions. Besides hematocrit and oxygen saturation, blood flow affects the intrinsic optical parameters. The parameter used to characterize the relation with the flow is the wall shear rate WSR, which is given by

WSR=vxr|r=ri.

The wall shear rate describes the shear stress induced in the fluid at the arterial wall due to the radial variation of the axial flow velocity. Blood, as a complex suspension of plasma and erythrocytes, shows the ability to accommodate for shear stresses by changes in erythrocyte cell morphology as well as changes in intercellular organization of erythrocytes inside the plasma. Both effects have an influence on the physical processes of scattering and absorption, and hence the intrinsic optical parameters of flowing blood. In the literature, four different regimes have been identified [6, 8, 9]:

  • At low shear rates WSR = 0 – 60 s−1, the erythrocytes build aggregates, known as rouleaux, that tend to grow into large branched rouleaux networks for blood in stasis. In this intercellular organization, a significant amount of photons is able to travel through the large plasma gaps in between the rouleaux network without interacting with the cellular parts of blood making blood most transparent.
  • At intermediate shear rates WSR ≤ 200 s−1, the aggregates are dispersed into randomly distributed erythrocytes inside the blood plasma. In this configuration, the mean free path length between two consecutive scattering or absorption events is smallest and hence the transparency of blood is lowest.
  • At high shear rates WSR ≥ 600 s−1, the red blood cells are aligned with the flow and deform into prolate ellipsoids. The intercellular organization changes into cellular layers that slide on plasma layers. In this configuration, the transparency of blood increases. Furthermore, this configuration allows blood to be less viscous at higher shear rates which is an essential feature of our cardiovascular system [27].
  • At even higher shear rates WSR ≥ 1000 s−1, the erythrocytes progressively elongate, while the intercellular organization remains in the layered configuration. The layered intercellular organization as well as the increasing elongation of erythrocytes into prolate ellipsoids further increase the transparency of blood.

Table 6 introduces the latest experimental results of the intrinsic optical parameters in relation to the wall shear rate [8, 9]. All intrinsic optical parameters are selected according to the following physiological conditions: hematocrit HCT = 42.1 %, oxygen saturation SpO2 ≥ 98 %, wavelength λ = 635 nm, and temperature T = 20 °C.

Unfortunately, the experiments of Enejder et al. [9] were carried out under non-physiological conditions with respect to oxygen saturation and hematocrit. Therefore, the relative changes of μa, μs and g given by Enejder et al. are combined with the absolute values given by Friebel et al. [8] to yield a second set of absolute values for intrinsic optical parameters of flowing blood. The refractive index of blood used in our simulations was 1.40 [19].

Tables Icon

Table 6. Intrinsic optical parameters of flowing blood for shear rates ranging from 200 s−1 to 1000 s−1. Bold values are taken from literature [8, 9], the other values are linearly interpolated.

3.3.4. Transmission mode photoplethysmograms at large distending arteries

The transmission photoplethysmograms were calculated according to the arrangement introduced in Fig. 3. As outlined in Section 3.3.3, the optical parameters of flowing blood can be separated into four regimes of which only three are relevant in macrocirculation, i.e. large arteries [21]. Wall shear rates below 200 s−1 are not considered in the simulations here, since the corresponding experimentally determined optical parameters are affected by microcirculatory effects, such as rouleaux forming.

Figure 7 shows the simulated transmission mode photoplethysmograms in terms of normalized transmittance versus arterial strain. The arterial strain was varied in 10 steps of about 0.5 % up to 5.14 % strain, yielding a maximum diameter of z = 5.26 mm + 270μm. The blue circles display the normalized mean transmittance averaged over wall shear rates ranging from 200 to 1000 s−1. The red error bars show the standard deviation of the simulated transmittances for the accepted wall shear rates. The left plot of Fig. 7 was simulated for the intrinsic optical parameters of Friebel et al., whereas the right plot shows the transmittance versus arterial strain for the optical parameters given by Enejder et al. Both graphs of Fig. 7 show a linear correlation between the transmission and the arterial strain as illustrated by the black line, representing a linear least-squares fit. Furthermore, both graphs show that the normalized transmission mode photoplethysmograms are primarily influenced by arterial strain and only slightly by variations in optical transmittance of flowing blood. At an arterial distension of 5.14 % which corresponds to an artery of elderly people, the normalized transmittance decays by roughly Tnorm ≈ 30 %, while the maximum error bar of the intrinsic optical parameters of Friebel et al. is ΔTnorm ≈ 5 %, and the maximum error bar of the intrinsic optical parameters of Enejder et al. is ΔTnorm ≈ 4 %.

 figure: Fig. 7

Fig. 7 Normalized transmission mode photoplethysmograms (blue, circle) versus arterial strain ε(%). left: data based on the intrinsic optical parameters of blood given by Friebel et al. [8]; right: data based on the intrinsic optical parameters reported by Enejder et al. [9]. The red error bars show the standard deviations of the normalized transmission photoplethysmograms with respect to variations of the intrinsic optical parameters in the wall shear rate range of 200 s−1WSR ≤ 1000 s−1.

Download Full Size | PDF

3.4. Discussion of photoplethysmograms in macrocirculation

A priori, the authors did not expect the sixfold dominance of arterial distension over changes in blood opacity in PPGs, since both, the arterial distension and the intrinsic optical parameters of blood vary in the single-digit percentage range. The extinction coefficient of Friebel et al. varies by Δμt ≈ 7.1 %, while the arterial strain is only εmax = 5.1 %. Furthermore, the simulation results, displayed in Fig. 7, reveal that the transmission as a function of arterial strain is a linearization of the radiative flux Φ(z) in z-direction. This linearity is substantiated by the weak decay of the radiative flux at distances in the range of the arterial diameter, where the arterial diameter varies only in the single-digit percentage range. Quantitatively, the linear fits in Fig. 7 are characterized by correlation coefficients R2 = 0.9990 in case of Friebel et al. and R2 = 0.9988 in case of Enejder et al.

For both cases, the slope of the linear fit is s = −0.056. Multiplication of s with the distension waveform of the outer arterial wall, shown in Fig. 6, yields the derived photoplethysmogram. The negative values of the PPG waveform in Fig. 8 indicate that the measured pressure pulse and the derived photoplethysmogram are inversely related. This inverse relation is well-known from general clinical practice [2]. The derived PPG waveform in Fig. 8 shows all key characteristics of experimentally obtained PPGs, such as diastole, systole, and pulse wave reflection [4, 28].

 figure: Fig. 8

Fig. 8 Derived photoplethysmogram and the corresponding intra-arterial pressure. The PPG is inverted as per common clinical practice.

Download Full Size | PDF

Finally, the size of the error bars and the influence of changes in blood opacity on PPGs in macrocirculation also have to be discussed. Friebel et al. and Enejder et al. used rigid flow cuvettes in their experiments. These rigid flow cuvettes are not capable of accommodating the pressure gradient by a lateral wall movement, as outlined in Subsection 3.3.1. Therefore, the forces inside blood are higher and may result in more pronounced changes in red blood cell morphology, intracellular organization of red blood cells or even density fluctuations inside blood. If so, the experiments of Friebel et al. and Enejder et al. overestimate the dependency on wall shear rates of the intrinsic optical properties of flowing blood. This would lead to a systematic overestimation of the error bars in Fig. 7. Therefore, we expect the impact of the variations in intrinsic optical parameters with varying flow conditions to be even lower than the values presented in Section 3.3.4.

4. Conclusion

For describing the photon energy distribution in large distending arteries, we extended the MCML code of Wang et al. to support cylindrical geometries as well as Fresnel reflections at the cylindrical boundaries and validated our code amendments. The extMCML code is available at our website, thereby offering the ability to study light transport phenomena inside multilayered cylindrically shaped tissue such as, bones, esophagus, intestine, or as in our case blood vessels. Furthermore, we added a section to the very comprehensively written MCML manual by Jacques et al., to continue their idea of an easy-to-use program to study the propagation of light inside living tissue.

After validating our code amendments, we studied the photon energy distribution inside a static artery, resulting in the most detailed description of radiative transfer in arteries till date [29]. These results could be beneficial for emerging new biophotonic applications such as photodynamic therapy [30], or photoplethysmography at large arteries [26, 31, 32].

By bridging the gap between the biomechanics of the artery and the light transport inside the artery, we could study how arterial distension affects photoplethysmograms at large, distending, blood-filled arteries. Our results predict that transmission mode photoplethysmograms measured in macrocirculation are mainly a function of arterial strain. Microcirculatory effects causing changes in optical transmittance of blood are only of subordinate importance. Furthermore we showed, that the high extinction inside the artery in combination with small arterial strains make photoplethysmography at large arteries a linear measurement principle of arterial strain.

Acknowledgments

The authors thank Asst. Prof. Dr. rer. nat. Martina Meinke for supporting our theoretical investigations with experimental data on the intrinsic optical properties of flowing blood.

References and links

1. V. S. Letokhov, “Laser light in biomedicine and the life sciences: From the present to the future,” in Biomedical Photonics Handbook, T. Vo-Dinh, ed. (CRC Press, 2003), pp. 1–16.

2. J. Allen, “Photoplethysmography and its application in clinical physiological measurement,” Physiol. Meas. 28, R1–R39 (2007). [CrossRef]   [PubMed]  

3. A. Reisner, P. A. Shaltis, D. McCombie, and H. Asada, “Utility of the Photoplethysmogram in circulatory monitoring,” Anesthesiology 108, 950–958 (2008). [CrossRef]   [PubMed]  

4. J. T. B. Moyle, “Optical principles,” in Pulse Oximetry (BMJ Books, 2003), pp. 7–14.

5. A. A. R. Kamal, J. B. Harness, G. Irving, and A. J. Mearns, ”Skin photoplethysmography - a review,” Computer Methods and Programs in Biomedicine , 28, 257–269 (1989). [CrossRef]  

6. H. J. Klose, E. Volger, H. Brechtelsbauer, L. Heinich, and H. Schmid-Schoenbein, “Microrheology and light transmission of blood I. The photometric effects of red cell aggregation and red cell orientation,” Pfluegers Arch. 333, 126–139 (1972). [CrossRef]  

7. L. D. Shvartsman and I. Fine, “Optical transmission of blood: Effect of erythrocyte aggregation,” IEEE Trans. Biomed. Eng. 50(8), 1026–1033 (2003). [CrossRef]   [PubMed]  

8. M. Friebel, J. Helfmann, G. Mueller, and M. Meinke, “Influence of shear rate on the optical properties of human blood in the spectral range 250 to 1100 nm,” J. Biomed. Opt. 12(5), 1–8 (2007). [CrossRef]  

9. A. M. K. Enejder, J. Swartling, P. Aruna, and S. Andersson-Engels, “Influence of cell shape and aggregate formation on the optical properties of flowing whole blood,” Appl. Opt. 42(7), 1384–1394 (2003). [CrossRef]   [PubMed]  

10. S. A. Prahl, M. Keijzer, S. L. Jacques, and A. J. Welch, “A Monte Carlo model of light propagation in tissue,” SPIE Proceedings of Dosimetry of Laser Radiation in Medicine and Biology , 5, 102–111 (1989).

11. L. Wang, S. L. Jacques, and L. Zheng, “MCML - Monte Carlo modeling of light transport in multi-layered tissues,” Computer Meth. Prog. Biomed. 47, 131–146 (1995). [CrossRef]  

12. L. Wang, S. L. Jacques, and L. Zheng, “CONV - convolution for responses to a finite diameter photon beam incident on multi-layered tissues,” Computer Meth. Prog. Biomed. 54, 141–150 (1997). [CrossRef]  

13. R. H. Cox, “Wave propagation through a newtonian fluid contained within a thick-walled, viscoelastic tube,” Biophys. J. 8, 691–709 (1968). [CrossRef]  

14. V. V. Tuchin, “Part I: An introduction to tissue optics,” in Tissue optics: Light Scattering Methods and Instruments for Medical Diagnosis, V. V. Tuchin, ed. (SPIE, 2007), pp. 3–142.

15. L. G. Henyey and J. L. Greenstein, “Diffuse radiation in the galaxy,” Astrophys. J. 93, 70–83 (1941). [CrossRef]  

16. M. Born and E. Wolf, “Basic properties of the electromagnetic field,” in Principles of Optics (Cambridge University Press, 1999), pp. 40–43.

17. W. W. Nichols, M. F. O’Rourke, and C. Vlachopoulos, “Properties of the arterial wall: practice,” in McDonald’s Blood Flow in Arteries: Theoretical, Experimental and Clinical Principles, W. W. Nichols, M. F. O’Rourke, and C. Vlachopoulos, eds. (Hodder Arnold Publishers, 2011), pp. 77–109.

18. J. Krejza, M. Arkuszewski, S. E. Kasner, J. Weigele, A. Ustymowicz, R. W. Hurst, B. L. Cucchiara, and S. R. Messe, “Carotid artery diameter in men and women and the relation to body and neck size,” Stroke , 37, 1103–1105 (2006). [CrossRef]   [PubMed]  

19. J. Mobley and T. Vo-Dinh, “Optical properties of tissue,” in Biomedical Photonics Handbook, T. Vo-Dinh, ed. (CRC Press, 2003), pp. 1–75.

20. W. -F. Cheong, S. A. Prahl, and A. J. Welch, “A review of the optical properties of biological tissues,” IEEE J. Quantum Electron. 26(12), 2166–2185 (1990). [CrossRef]  

21. P. V. Stroev, P. R. Hoskins, and W. J. Easson, “Distribution of wall shear rate throughout the arterial tree: A case study,” Atherosclerosis 191(2), 276–280 (2007). [CrossRef]  

22. S. K. Samijo, J. M. Willigers, R. Barkhuysen, P. J. E. H. M. Kitslaar, R. S. Reneman, P. J. Brands, and A. P. G. Hoeks, “Wall shear stress in the human common carotid artery as function of age and gender,” Cardiovascular Res. 39, 515–522 (1988). [CrossRef]  

23. K. Witzig, “Ueber erzwungene Wellenbewegungen zaeher, inkompressibler Fluessigkeiten in elastischen Roehren,” Inaugural Dissertation, University of Bern, Bern, (1914).

24. J. R. Womersley, “Section III: The equations of motion of the freely-moving elastic tube and the derivation of the pulse-velocity,” in An Elastic Tube Theory of Pulse Transmission and Oscillatory Flow in Mammalian Arteries, J. R. Womersley, ed. (Wright Air Development Center, 1957), pp. 19–29.

25. I. Mirsky, “Wave propagation in a viscous fluid contained in an orthotropic elastic tube,” Biophys. J. 7, 165–186 (1967). [CrossRef]   [PubMed]  

26. D. Ruh, S. Sherman, M. Theodor, J. Ruhhammer, K. Foerster, C. Heilmann, F. Beyersdorf, H. Zappe, and A. Seifert, “Determination of vessel wall dynamics by optical microsensors,” in Engineering in Medicine and Biology Society (EMBC), 2012 Annual International Conference of the IEEE, pp. 2359–2362 (2012).

27. G. B. Thurston, “Rheological parameters for the viscosity, viscoelasticity and thixotropy of blood,” Biorheology 16, 149–162 (1979).

28. D. Ruh, P. Reith, S. Sherman, M. Theodor, J. Ruhhammer, A. Seifert, and H. Zappe, “Stretchable optoelectronic circuits embedded in a polymer network, Adv. Mater. DOI: [CrossRef]  , (2013).

29. M. Keijzer, S. L. Jacques, S. A. Prahl, and A. J. Welch, “Light distributions in artery tissue: Monte Carlo simulations for finite-diameter laser beams,” Lasers Surgery Med. 9, 148–154 (1989). [CrossRef]  

30. R. R. Allison, H. C. Mota, V. S. Bagnato, and C. H. Sibata, “Bio-nanotechnology and photodynamic therapyState of the art review,” Photodiagnosis and Photodynamic Therapy 5, 19–28, (2008). [CrossRef]  

31. D. Ruh, S. Sherman, H. Zappe, and A. Seifert, “Optoelectronics on flexible substrates for biomedical applications,” in Optical MEMS and Nanophotonics (OMN), 2012 International Conference on, pp. 186–187 (2012).

32. J. Fiala, P. Bingger, D. Ruh, K. Foerster, C. Heilmann, F. Beyersdorf, H. Zappe, and A. Seifert, “An implantable optical blood pressure sensor based on pulse transit time,” Biomed. Microdev. 15(1), 73–81 (2013). [CrossRef]  

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (8)

Fig. 1
Fig. 1 Comparison of the simulation geometry between a) MCML and b) extMCML. a) In MCML, a stack of infinitely extended planes is illuminated by a Dirac beam in the z-direction. b) In extMCML, the Dirac beam strikes an infinitely long multilayered cylinder in the lateral direction.
Fig. 2
Fig. 2 Simulations by extMCML and MCML. The left plot shows the radiative flux obtained by MCML, the plot on the right presents the radiative flux determined by extMCML. The two different detectors are sketched in red, both detector areas are 1 mm2.
Fig. 3
Fig. 3 Simulation geometry used to calculate the photon energy distribution inside a static, large artery. The 3D plot of the simulation geometry introduces the thicknesses of the various layers forming the arterial wall, adventitia ta, media tm and intima ti; the length of the investigated arterial segment l; inner and outer radii of the arterial wall ri and ro; the Cartesian coordinate system x,y,z; and the Dirac beam striking the arterial wall in z direction.
Fig. 4
Fig. 4 Decay of the radiative flux inside the model artery at rest. left: normalized radiative flux at y = 0μm as a function of z calculated by integrating along the x direction. The blood layer is marked between the two arterial walls. The three different arterial layers forming the arterial wall are labeled as adventitia ta, media tm, and intima ti (dashed, black). right: Normalized radiative flux in x direction integrated over z at y = 0μm. Since there are no optical boundaries in direction of x, the plot shows no discontinuities. The insets in both graphs show the location of the plotted radiative flux inside the cylinder.
Fig. 5
Fig. 5 Schematic of the arrangement for calculating the photoplethysmograms in large, distending arteries. The upper part shows two photon energy distributions in a cross-sectional view for two different points in time t0 and t1. The cross-sectional view at t0 is to scale while the cross-sectional view at t1 is magnified to point out the change in diameter. The exact change in diameter is presented in the lower part of this figure. Here, the temporal behavior of the inner (red, dashed) and outer (red, dotted) arterial wall distensions ur are plotted. The dotted vertical lines indicate the two states above, where t0 and t1 represent the arterial state at the diastole and systole, respectively.
Fig. 6
Fig. 6 Measured pressure waveform p(t) and calculated inner (red, dashed) and outer (red, dotted) arterial displacements. The displacement waveforms are calculated according to Equation 10 with respect to the biomechanical parameters of Table 5.
Fig. 7
Fig. 7 Normalized transmission mode photoplethysmograms (blue, circle) versus arterial strain ε(%). left: data based on the intrinsic optical parameters of blood given by Friebel et al. [8]; right: data based on the intrinsic optical parameters reported by Enejder et al. [9]. The red error bars show the standard deviations of the normalized transmission photoplethysmograms with respect to variations of the intrinsic optical parameters in the wall shear rate range of 200 s−1WSR ≤ 1000 s−1.
Fig. 8
Fig. 8 Derived photoplethysmogram and the corresponding intra-arterial pressure. The PPG is inverted as per common clinical practice.

Tables (6)

Tables Icon

Table 1 Optical parameters to compare extMCML with the original MCML code.

Tables Icon

Table 2 Verification of the extMCML code by comparison of total transmittance and diffuse reflectance with the MCML code provided by [11]. Input power is 1 W, detector area is 1 mm2.

Tables Icon

Table 3 Relative error between MCML and extMCML for various anisotropy factors.

Tables Icon

Table 4 Optical properties [19,20] of the three different layers of the arterial wall and blood at a wavelength of 633 nm.

Tables Icon

Table 5 Biomechanical parameters of the arterial wall defining the linear viscoelastic behavior given in Equation 8 optimized to yield arterial strains of roughly 5 %.

Tables Icon

Table 6 Intrinsic optical parameters of flowing blood for shear rates ranging from 200 s−1 to 1000 s−1. Bold values are taken from literature [8, 9], the other values are linearly interpolated.

Equations (12)

Equations on this page are rendered with MathJax. Learn more.

( s ) I ν ( r , s ) + ( μ a + μ s ) μ t I ν ( r , s ) = μ t 4 π 4 π I ν ( r , s ) p s ( s , s ) ) d Ω ,
p s ( θ ) = 1 4 π 1 g 2 ( 1 + g 2 2 g cos θ ) 3 / 2
g = 0 π p s ( θ ) cos θ 2 π sin θ d θ .
p r = 2 ( p i r ) p i + p i ,
p t = n i n t p i + ( 1 n i n t 1 ( p i r ) 2 n i p i r ) r
A ( x , y , z ) = N p h Δ W i N p h d V
Φ ( x , y , z ) = A ( x , y , z ) μ a .
μ * = μ 0 1 + i ω λ 2 1 + i ω λ 1 ,
p ( r , x , t ) = n N P n ( r ) exp [ i ( n ω t γ n x ) ] ,
u r = N n γ n P n n ω ρ ( Δ 13 Δ 11 J 1 ( k n r ) + Δ 14 Δ 11 Y 1 ( k n r ) + Δ 15 Δ 11 J 1 ( i γ n r ) + Δ 16 Δ 11 Y 1 ( i γ n r ) ) exp [ i ( n ω t γ n x ) ] ,
k n 2 = n 2 ω 2 ρ w μ * γ n 2 ,
WSR = v x r | r = r i .
Select as filters


Select Topics Cancel
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.