Numerical Modeling of CompressionControlled Lowlevel Laser Probe for Increasing Photon Density in Soft Tissue
 Author: Kwon Kiwoon, Son Taeyoon, Yeo Changmin, Jung Byungjo
 Organization: Kwon Kiwoon; Son Taeyoon; Yeo Changmin; Jung Byungjo
 Publish: Current Optics and Photonics Volume 15, Issue4, p321~328, 25 Dec 2011

ABSTRACT
Various methods have been investigated to increase photon density in soft tissue, an important factor in lowlevel laser therapy. Previously we developed a compressioncontrolled lowlevel laser probe (CCLLP) utilizing mechanical negative compression, and experimentally verified its efficacy. In this study, we used Bezier curves to numerically simulate the skin deformation and photon density variation generated by the CCLLP. In addition, we numerically modeled changes in optical coefficients due to skin deformation using a linearization technique with appropriate parameterization. The simulated results were consistent with both human in vivo and porcine ex vivo experimental results, confirming the efficacy of the CCLLP.

KEYWORD
Lowlevel laser therapy , Bezier curve , Photon density , Numerical modeling

I. INTRODUCTION
Ultraviolet (UV), visible, and near infrared (NIR) light have been extensively used in therapeutic and diagnostic modalities [13]. However, the depth of light penetration is still limited to a few centimeters because of the light scattering property of tissue, resulting in low photon densities.
Various methods such as increasing the power and diameter of the incident beam [35], application of optical clearing agents [6,7], and skin compression [810] have been suggested to increase the photon density in tissues. In a previous study, we developed a compressioncontrolled lowlevel laser probe (CCLLP) to increase the photon density in tissues by applying negative compression [11]. Figure 1 represents the CCLLP, which consists of a vacuum hole to supply negative compression, a hollow brace for guiding and supporting the optical fiber, and an acrylic probe body to maintain the negative compression applied to the tissue. We experimentally validated the efficacy of the CCLLP by measuring the photon density in ex vivo porcine skin samples [11].
In another previous study, we mathematically modeled and simulated the effects of tissue compression [10]. We simulated tissue compression in two ways: (1) nonuniform skin deformation with constant optical coefficients and (2) uniform skin deformation with variable optical coefficients depending on changes in skin thickness. However, neither of these approaches accurately simulated the mechanism of the CCLLP.
Although several numerical simulations for nonflat boundaries or layers have been reported [12,13], the geometry
of the outer boundaries and inner layers similar to the skin affected by the CCLLP has not been described. The skin region affected by the CCLLP swells outward except for the area under the optical fiber probe. In the previous uniform compression study [10], the photon density did not increase in the swollen areas but increased in the region of skin under the optical fiber probe. In this study, which is a fundamental generalization of the simulation reported in [10], we numerically simulated such competing effects between the swollen and compressed regions in the CCLLP.
In three subjects, we experimentally measured the change in skin thickness due to negative compression and used the values for the numerical simulations. We represented the deformation of the skin surface associated with the CCLLP using three quadratic Bezier curves. We also numerically modeled the changes in optical coefficients with changes in skin thickness using a linearization method. We obtained simulation results as functions of the modeling parameters, radius of the CCLLP, and negative compression, and compared these with results from an
ex vivo experiment using a scaling factor.II. METHODS
The photon density in tissues was modeled using the radiative transfer equation or its diffusion approximation [14]. The diffusion approximation in the frequency domain was represented by the following elliptic partial differential equation in the region of interest Ω :
where
r_{0} is some point on the boundary of Ω . For the notations used in the equations, see the Nomenclature section before the Reference section at the end of this paper.To solve (1) using the finite element method (FEM), the geometrical domain Ω was assumed to be bounded, and appropriate boundary conditions were assigned to the boundary
∂ Ω . Let us assume that Ω = [w ,w ] × [d, 0]. That is to say, the region of interest is a rectangular region in the twodimensionalxz plane modeling a few skin layers with depthd , deposited in the direction of the zaxis. Then∂ Ω is composed of four straight lines:∂top Ω =[w, w]×{0}: the top of the skin contacting the air
∂bottom Ω =[w, w]×{d}: the bottom of the skin region of interest
∂left Ω ={w}×{d,0}: the left of the skin region of interest
∂right Ω ={w}×{d,0}: the right of the skin region of interest
Note that only
∂^{top} Ω contacts air, while the other boundaries contact skin. Thus, we imposed the Neumann boundary condition on∂^{left} Ω ,∂^{right} Ω , and∂^{bottom} Ω and the Robin boundary condition on∂^{top} Ω to reflect differences in refractive index mismatching [14]. Therefore, we obtained the following boundary value problem for modeling light propagation in skin:Let us divide Ω into small triangles or rectangles, which are called finite elements, and {
u }_{i}=1,…_{i} Nb be the corresponding finite element bases spanning piecewiselinear or bilinear function spaces. The finite element formulation for (1) with the boundary condition (2) using the finite element bases is then given as follows:The finite element solution for (3) was found using COMSOL 3.4 (COMSOL AB, Stockholm, Sweden). In the simulation, the photon density, or light fluence rate, was measured at the bottom ∂^{bottom}Ω of the domain of interest and at a distance 7.5 cm away. However, we could neglect the 7.5 cm distance in the air for the photon density measurement since light is mostly absorbed and scattered in the tissue rather than in the air.
The optical properties of human tissues have been reported in several studies [3,15]. Based on these studies, a fivelayer skin model was used for a 633 nm light [3].
2.1. Skin Deformation by Negative Compression
For simplicity, we modeled the domain of interest Ω in two dimensions by neglecting the yaxis. When there is no compression, Ω is given as follows:
where Ω
_{i} anddi are the domain of interest and the thickness of thei th skin layer, respectively,dx is the dimension of the skin layer along thex axis, and t_{5} =d . The thicknesses of layersd_{i} (i = 1,2,3,4,5) are presented in Table 1 with their optical coefficients. The overall fivelayer structure is depicted in Figure 2.We then formulated the layers Ω
_{i} (i = 1,2,3,4,5) with negative compression. Let 2K and 2O be the length of the region contacting the entire acrylic probe body of the CCLLP andthe length of the region contacting the guide for the optical fiber, respectively. Let
C be the midpoint betweenK andO . If negative compression is applied throughout the regions [K ,O ] and [O ,K ] on the skin surface, these regions become swollen by the maximum height difference shown in Figure 3 for the applied negative compression. On the other hand, in the region [O ,O ], the skin surface should be retracted between the points (O , 0) and (O , 0) because of the swollen surfaces surrounding it. For more detailed evaluations, elastic or viscoelastic analysis could be used [16,17]; however, in this study, we used three quadratic Bezier curves for the regions [K , O ], [O ,O ], and [O ,K ]. The data shown in Figure 4 were used as the maxima of the two Bezier curves for [K , O ] and [O ,K ], and the Bezier curve for [O ,O ] was constructed to join the three curves without any discontinuity of tangential lines. Letl_{i} be the upper surface ofi th layer (i = 1,2,3,4,5), which is a straight line before the application of negative compression. For this purpose, control points were chosen as follows:where C was chosen as
K +O/2 to cause the two Bezier curves for [K , O ] and [O ,K ] to be symmetric;S = 2O/ (K O ) was chosen so that the three Bezier curves joined without any discontinuity of tangent lines, especially at O andO ; andsi(i = 1,2,3,4,5) is the maximum skin surface change fori th layer. We chose w = 3.5 cm,K = 2 cm, andO = 0.2 cm, based on the sizes of the skin sample and CCLLP.Using the definition of a quadratic Bezier curve, we obtained the following detailed representations for
l_{i}, i = 1,2,3,4,5:Therefore, the thicknesess of the i layers are as follows:
Assuming the thickness changes for the third, fourth, and fifth layers are relatively small and the change for second layer is half of the change in the first layer, we let
s_{3} = s_{4} = s_{5} = 0 and
2.2. Change in Optical Coefficients with Negative Compression
In the CCLLP, negative compression is applied to the skin surrounding the central area
[？O,？K ]×{0}？[K,O ]×{0} in contact with the probe, and positive compression is applied at the center [?K,K ]×{0} due to skin tension. The change in optical coefficients with uniform and partial compression of the skin was previously studied and reported [10]. Letand
μ _{a},μ _{s}'be the optical coefficients for the ith layer (
i = 1,2) when the thickness of the layer isd_{i} andrespectively. Using an analysis similar to that in [10], we assumed that the relationship between the thickness and the absorption coefficient is:
The coefficient
α is the optical modeling coefficient between 0 and 1, which related the reduction of skin penetration depth to the increase of optical coefficients. The more fundamental meaning of the coefficient is discussed in the “Discussion” section.Equation (8) is linearly approximated as follows:
Similarly, we obtained the following relationship between the thickness and the reduced scattering coefficient:
Given that the diffusion coefficient is defined as a constant multiplied by the reciprocal of the sum of the absorption and reduced scattering coefficients, using (9) and (10) we obtained:
In (9), (10), and (11), we used firstorder Taylor series approximations.
Using (7) and (11), we obtained an explicit representation of the diffusion coefficient as follows:
We chose the value for s1 based on the data shown in Figure 4, and
si(i = 1,2,3,4,5) were chosen such thatdecreased with the ratio of the depth from the surface to the entire depth of the skin. Thus, (12) becomes:
and (9) becomes:
III. RESULTS
Figure 1 shows CCLLPs with radii of 1, 1.5, and 2 cm. In Figure 3, Skin deformation by negative compression was investigated in the left arms of three subjects. Figure 4 shows the change in skin thickness averaged across the three subjects during negative compression as a function of the radius of the CCLLP. As described in the section II.1, the change in skin thickness was modeled using three Bezier curves with continuous tangential lines.
For the FEM analysis of the CCLLP, a finite element mesh was applied to the modeled geometry of the skin. We used meshes shown in Figure 5 for the skin deformation geometry to obtain the photon density distribution by FEM. The mesh is finer in the upper layers and coarser in the lower layers because skin deformation is expected to be less pronounced in the lower layers.
The optical coefficients for the five layers used in modeling the skin structure are shown in Table 1 [3]. The resulting photon density distributions at a depth of 0.2 cm below the skin surface are presented in Figure 6 as functions of negative compression, the thicknesscoefficient
relation parameter
α which is represented in (8), and the radius of the CCLLP.In our previous study [11], we investigated the efficacy of CCLLP by measuring the photon density distribution in
ex vivo porcine samples. In Figure 7, the photon density ratio, which is defined as the ratio of the photon density distribution between a given negative compression and no compression, was compared to the experimental data for a CCLLP with a radius of 1 cm. To compare the results of this study with theex vivo data, the maximum change in skin thicknessin vivo for human subjects was halved in the numerical simulation. The experimental and simulation results for various values ofα (0.5,0.6,0.7, and 0.8) was shown in Figure 7.IV. DISCUSSION
The photon density in the tissue due to negative compression was influenced by two successive and competing physical phenomena: (1) the change in skin thickness due to compression and (2) the change in optical coefficients due to the change in skin thickness. In previous studies, the correlation between the change in skin thickness and negative compression was mathematically modeled using elastic equations [18,19]. In this study, the correlation was experimentally obtained from three subjects, and the data (Figure 4) was used to mathematically simulate the mechanism of the CCLLP. The changes in skin thickness for all three subjects were approximately proportional to the logarithm of the negative compression, which corresponds well with the results of previous studies [18,19].
Skin surface deformation by the CCLLP was experi
mentally measured as a function of negative compression (Figure 3). The region of skin contacting the optical fiber probe was compressed for two reasons: (1) this region of the skin was pressed down and thus could not swell up, and (2) the skin in the contacted region was pushed aside to the surrounding region, which swells outward. In other words, the skin inside the CCLLP, except for that in contact with the optical fiber probe, swells outward, whereas the contacted region is compressed inward. Thus, the skin surface can be modeled by: (1) two straight lines representing the two flat areas of skin outside the CCLLP, (2) two concave Bezier curves representing the raised skin region inside the CCLLP not contacting the optical fiber probe, and (3) one convex Bezier curve representing the compressed skin region contacting the optical fiber probe. The two concave Bezier curves and the convex Bezier curve were constructed to satisfy the C1 assumption at the two points where the curves meet, representing the boundary of the optical fiber probe. The C1 assumption states that the tangential lines representing the skin surface, as well as the skin surface itself, are continuously connected across the optical fiber probe inside the CCLLP. In other words, the skin surface deformation due to negative compression should be smooth and not have any discontinuities along the skin surface. The resulting skin surface deformation was shown in Figure 5.
The optical coefficients are known to be inversely proportional to skin penetration depth [8,9]. However, the relationship has been experimentally investigated with only a limited number of samples and is not well established for human or animal skin samples. Therefore, the relationship was mathematically modeled [10] when uniform compression was applied to the skin surface, corresponding to
α = 0.5 in (8). Analysis of the changes in the optical coefficients with skin thickness was generalized from a uniform geometry [18] to a more complex geometry (Figure 5). Detailed changes in the optical coefficients with the generalized parameter α were computed using (13) and (14).The skin thickness at the center decreased because of negative compression, resulting in increased photon density. At the same time, by the result of [8,9,10], we assumed that the optical coefficients (including the absorption coefficient (8),(9), the reduced scattering coefficient (10), and the diffusion coefficient (11)) per unit skin penetration depth increased through the optical modeling coefficient
α as the skin thickness decreased, and thus interfered with the increase in photon density. Previous experimental [8,9] and simulation [10] results have shown that the increase in photon density due to the reduction in skin thickness is greater than the decrease due to changes in optical coefficients. However, this study modeled a more complicated phenomenon than the uniform compression case. In the skin deformation model, the skin regions inside the CCLLP and those contacting the optical fiber probe are raised and compressed, respectively. Therefore, the effects of skin deformation on photon density in these two regions are in opposite directions. A previous study usingex vivo porcine samples [11] presented the efficacy of the CCLLP in enhancing photon density. These results were verified by numerical simulation results (Figure 6) which used data onin vivo changes in human skin thickness. From these results, we conclude that the decrease in skin thickness directly under the optical fiber probe is more important than the increase in skin thickness in the surrounding areas. In addition, the negative compression, value forα , and radius of the CCLLP in the present results were all larger than in the previous study.The peak intensities of photons in the
ex vivo porcine samples [11] were compared with the simulated results. Changes in skin thickness were not obtained in the previous study because of the instability of the physical state of theex vivo porcine skin samples. Therefore, skin deformation was measured in the arms of three human subjects, and changes in the skin thickness were halved to enable comparisons with theex vivo porcine skin samples.V. CONCLUSION
The CCLLP was developed to increase the photon density in tissues by applying negative compression. Skin deformation due to the CCLLP was numerically modeled using three Bezier curves. The finite element simulation used experimental data on changes in skin thickness from three human subjects, skin deformation modeling, and changes in optical coefficients as a function of the optical modeling coefficient
α . The photon density increased when the negative compression and radius of the CCLLP increased. The simulated results at various values ofα were similar to theex vivo experimental results. Although we obtained favorable result using the optical modeling coefficient, we have not experimental verification about the optimal choice for the coefficient, which will be our next study subject.> NOMENCLATURE
Φ : The photon density profile [W/cm^{2}]
q : The light source function [W/cm^{3}]
R : The Robin function, which is the photon density profile with respect to the Dirac delta source [W/cm^{2}]
μ _{a} : The absorption coefficient [cm^{1}]
μ _{s} : The scattering coefficient [cm^{1}]
g : The anisotropy factor
μ _{s}' : The reduced scattering coefficient [cm^{1}] [μ _{s}' = (1g) μ s]
κ : The diffusion coefficient [cm] [κ = (3(μ _{a} μ _{s}'))1]
n : The refractive index
a : The reflection coefficient on the surface
c : The speed of light in a vacuum [cm/sec] (c = 2.99×10^{10} cm/sec)
ω : The angular modulation frequency of the source intensity [sec^{1}]
Ω : The region of interest (skin structure)
∂ Ω : The boundary of Ω (∂ Ω = ∂
^{top} Ω ∂^{left} Ω ∂ right Ω + ∂^{bottom} Ω )∂
^{top} Ω : The top boundary of Ω∂
^{bottom} Ω : The bottom boundary of Ω∂
^{left} Ω : The left boundary of Ω∂
^{right} Ω : The right boundary of ΩN
_{b} : The number of finite element basis functions

[FIG. 1.] The compressioncontrolled lowlevel laser probe (CCLLP) constructed with radii of (a) 1, (b) 1.5, and (c) 2 cm.

[TABLE 1.] Fivelayer skin model for 633 nm light used to numerically model the photon density profile in tissue

[FIG. 2.] The comparison of the fivelayer skin structure before (black line) and after deformation (red line on the surfaces of the first and second layers).

[FIG. 3.] In vivo skin deformation due to negative compression at (a) 0, (b) 5, (c) 10, (d) 15, and (e) 20 kPa applied to the skin surface of a subject using the CCLLP with a radius of 1 cm.

[FIG. 4.] The changes in skin thickness as a function of negative compression and CCLLP radius (1, 1.5, and 2 cm).

[FIG. 5.] Mesh grids used for finite element analysis of (a) the entire skin sample and (b) the region around the optical fiber probe with a negative compression of 5 kPa and a CCLLP radius of 2 cm.

[FIG. 6.] Photon density distribution at the base of the skin sample as a function of the (a) negative compression (α = 0.1, radius = 1 cm), (b) parameter α (radius = 1 cm, negative compression = 5 kPa), and (c) CCLLP radius (α = 0.1, negative compression = 5 kPa).

[FIG. 7.] The comparison of experimental photon density ratios [11] and modeling results with respect to the various values of the optical coefficient modeling parameter α . The photon density ratio is defined as the ratio of the photon density distribution with negative compression (0, 5, 10, 15, or 20 kPa) to that without compression (0 kPa)