A two-layer model of laser interaction with skin: A photothermal effect analysis

Optics & Laser Technology 43 (2011) 425–429

Contents lists available at ScienceDirect

Optics & Laser Technology journal homepage: www.elsevier.com/locate/optlastec

A two-layer model of laser interaction with skin: A photothermal effect analysis Kui-Wen Guan, Yan-Qi Jiang, Chang-Sen Sun, Hong Yu  Lab of Biomedical Optics, College of Physics and Optoelectronic Engineering, Dalian University of Technology, Dalian 116023, People’s Republic of China

a r t i c l e in fo

abstract

Available online 12 January 2010

In order to understand the photothermal effect mechanism of laser interaction with skin, we employed a two-layer model to describe the heat generation, transportation, and dispersion in the procedure of laser interaction with skin. A skin temperature distribution corresponding to the laser interaction direction is calculated to describe the time of skin gasification and the possible thermal injury. The magnitude of time is ms. This basic process provides a possible quantitative recognition of the applications of laser in clinical skin care. & 2009 Elsevier Ltd. All rights reserved.

Keywords: Thermal interaction Thermal properties Gasification

1. Introduction Lasers have been widely used in clinical skin care. However, the mechanism underlying these applications is still a problem for the researchers. Usually, laser–skin interaction can be categorized into five aspects [1], namely:

(1) Actinism: The irradiance could be continuous-wave or pulsed, from 10.3 to 1 W=cm2 or from 10.2 to 105 J=cm2 ; its temporal duration is very variable. (2) Thermal interaction [2]: The radiation is absorbed by the tissue transforming in internal energy producing a temperature increment, in the following manner: thermophysical, chemical, and biological. The irradiance could be continuouswave or pulsed, from 0.1 to 105 W=cm2 or from 0.1 to 106 J=cm2 ; its temporal duration is from 1 ms to hours. (3) Photoablation (ablative photodecomposition): UV radiation is used because the high energy of the photons can break molecular links and to ionize atoms. By this way only focalized atoms and molecules can be affected. The exposition times are very short, from 1 to 500 ns, and irradiance from 105 to 109 W=cm2 . (4) Plasma-induced ablatin: In this process the time focus is in ns. In this way irradiances in the magnitude of 1011 W=cm2 are obtained. (5) Photodisruption: In this process a very short exposition from ps to ns and high-power laser focused by a lens into the treated tissue are used. In this way irradiances in the magnitude of 1014 W=cm2 are obtained.  Corresponding author.

E-mail address: [email protected] (H. Yu). 0030-3992/$ - see front matter & 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.optlastec.2009.12.007

In order to understand each case a different model is tentatively proposed to formulate it. The most critical fact among these aspects is the photothermal effect [3,21,25,28], since these are particularly very complex, resulting from three distinct phenomena, namely, conversion of light to heat, transfer of heat and the tissue reactions, which are ultimately related to the temperature and the heating time. The models need to deal with how the heat generated, instantly transport, and possibly thermal damage during laser interaction. In the 1970s, Mainster et al. used the finite-difference method to solve the heat conduction equation for cylindrical symmetrical and thermal homogeneous media [8]. In 1971, Vassiliadis, Christian and Dedrick developed a Green’s function solution for symmetrical noncoherent sources [9]. Later, Welch combined a rate process model with the finite-difference model [10], and Takata expanded the model to relate to multiple layers, blood perfusion, blister formation, and damage [11]. Torres and Motamedi modified the Takata skin model to predict the temperature response of an in vitro aorta irradiated with an argon laser. In that study, the Beers law heat source of the Takata model was replaced by a Monte Carlo model of photon propagation to account for light scattering. Even with that modification, however, Torres and Motamedi found that the computed steady state and relaxation temperatures were much higher than the measured results. A second set of experiments demonstrated the existence of temperature-dependent surface cooling associated with the evaporation of water at the surface of the skin [12]. In this paper, we will use a two-layer model [17,19,23,27] to analyze the skin surface’s photothermal effect. The near-infrared or visible light laser radiation is delivered to the target area via an optical fiber with a cylindrical diffusively scattering surface. The targeted pathological tissue is destroyed by the immediate or

426

K.-W. Guan et al. / Optics & Laser Technology 43 (2011) 425–429

delayed hyperthermic and coagulative effects due to photon absorption and heat transfer in the skin.

2. Theory The near-infrared and visible light lasers are usually used as light sources for laser–interaction treatments [6,7]. The transport of the near-infrared range light in the tissue can be mathematically described by the transport equation [13,14]: Z dJðr; sÞ m ¼ ðma þ ms ÞJðr; sÞ þ s pðs; s0 ÞJðr; s0 Þ do0 ; ð1Þ ds 4p 4p where Jðr; sÞ is the radiant intensity at position r in the direction s. ma and ms are the absorption coefficients and the scattering coefficients of tissue for incident light as a function of wavelength (such as the absorption coefficient ma at wavelength li which is given by

ma ðli Þ ¼

N X

ei ðli ÞCi ;

ð2Þ

i¼1

where ei is the extinction coefficient and Ci is the concentration). ds and do0 are the differential lengths along the s direction and the differential solid angle in the s0 direction. pðs; s0 Þ is the phase function. The local fluence rate cðrÞ is defined as Z cðrÞ ¼ Jðr; sÞ do: ð3Þ 4p

The integro-differential equation in Eq. (1) cannot be solved analytically for complex geometries and boundary condition, so many analytical approximations such as the diffusion approximation and Kubella–Muck theory have been developed for conditions which are not matched in most in vivo tissue. Moreover, when the optical properties vary with position and time, analytical approaches become almost impossible. The laser light absorption results in internal energy increases in the local tissue. The various thermal models for blood perfused tissue have been reviewed perviously. The standard Pennes Bioheat equation is still the most commonly used model to describe the thermal behavior of living tissue without large, thermally significant blood vessels:

rc

@T ¼ r  ðkrTÞ þ rb cb ob ðTb TÞ þ Qm þQr ; @t

ð4Þ

where r, c and k are the tissue density, specific heat and heat conductivity. rb , cb and ob are the blood density, specific heat and perfusion rate per unit volume. T is the tissue temperature. Tb is the arterial blood temperature. Qr is the internal heat source term due to photo absorption and Qm is the metabolic heat generation per unit volume which is usually negligible. The theoretical model is shown in Fig. 1. In addition, we only discussed the normal incidence onedimensional case; Eq. (4) is simplified to the following form: @Ti ðz; tÞ ¼ r  ðki rTi ðz; tÞ; Ti ðz; tÞ ¼ Tðz; tÞT0 ; i ¼ 1; 2: ð5Þ @t Ti ðz; tÞ is defined as the time taken for the excess temperature, and T0 is the initial temperature. Though the heat conductivity is uniform, therefore, we can further simplify the form:

ri ci

@Ti ðz; tÞ @2 Ti ðz; tÞ ¼ ki ; i ¼ 1; 2: ð6Þ @t @z2 z-Axes normal to the free surface of the considered system, along the direction of the incident irradiance, are used to describe the problem, where the boundary z ¼ 0 represents the front surface, t is the transit time, and T0 is the initial temperature, defined as the time taken for the excess temperature to change from zero.

ri ci

Fig. 1. A schematic of the Cartesian coordinate system set up on two-layered tissues.

The corresponding initial conditions and boundary conditions are as follows: 1. The initial conditions: Ti ðz; 0Þ ¼ 0;

i ¼ 1; 2:

ð7Þ

2. The boundary conditions: k1

@T1 ð0; tÞ ¼ aQ0 ; @z

ð8Þ

a ¼ ma z0 ;

ð9Þ

T1 ðz0 ; tÞ ¼ T2 ðz0 ; tÞ; k1

ð10Þ

@T1 ðz0 ; tÞ @T2 ðz0 ; tÞ ¼ k2 ; @z @z

ð11Þ

T2 ð1; tÞ ¼ 0;

ð12Þ

where Q0 is the incident laser irradiance on the surface, a is the absorptance assumed to be temperature independent (we focus on the absorption of light), and z0 is the first layer’s thickness. We merger-related the constant, and further organize the form of the formula: @Ti ðz; tÞ @2 Ti ðz; tÞ ¼ ai ; @t @z2

ai ¼

ki

ri ci

;

i ¼ 1; 2:

ð13Þ

We assume that the incident light Q0 and the absorptance a is a constant, use the Laplace transform [4,5] to solve Eq. (13) and get the analytical solutions of two-layer system: ! rffiffiffiffiffiffiffiffiffiffi 1 aQ0 X 4a1 t ½2z0 ð1þ nÞz2 nþ1 T1 ðz; tÞ ¼ B exp  k1 n ¼ 0 p 4a1 t ! 2z0 ð1þ nÞz pffiffiffiffiffiffiffiffiffi ½2z0 ð1 þ nÞz erfc 4a1t ! rffiffiffiffiffiffiffiffiffiffi 1 aQ0 X 4a t ð2nz0 þzÞ2 1 n þ B exp  k1 n ¼ 0 p 4a1 t ! 2nz0 þ z ð2nz0 þ zÞ erfc pffiffiffiffiffiffiffiffiffiffi ; ð14Þ 4a1 t

K.-W. Guan et al. / Optics & Laser Technology 43 (2011) 425–429

T2 ðz; tÞ ¼

2aQ0 k1 ð1 þ eÞ n ¼ 0

B Bn B @

rffiffiffiffiffiffiffiffiffiffi 4a1 t

0

½ðzz0 Þ B expB  @ p

rffiffiffiffiffi 1 a1 þ ð2n þ 1Þz0 2 C a2 C A 4a1 t

rffiffiffiffiffi 31 a1 rffiffiffiffiffi   ðzz0 Þ þ ð2n þ 1Þz0 6 7C a1 a 7C; p2ffiffiffiffiffiffiffiffiffiffi  ðzz0 Þ þð2n þ1Þz0 erfc6 4 5A a2 4a t 2

1

ð15Þ where rffiffiffiffiffi t pffiffiffiffiffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi k2 r2 k2 c2 a2 k2 a1 ; e ¼ rffiffiffiffiffi ¼ pffiffiffiffiffi ¼ r1 k1 c1 k1 a2 t k1 a1 B¼

1e ; 1þ e

B o 1;

1. t = 0.05s 2. t = 0.1s 3. t = 0.2s

60 Temperature increase /K

0 1 X

50 40 30 20 1 0

ð16Þ

3

2

10

ð17Þ

427

0.2

0.4

According to Eqs. (14) and (15), we can calculate the temperature [18,20,22,26] in different organizations’ depth. For the palm of your hand, we use the epidermal thickness of 1.2 mm. When t ¼ 0:05, 0.1, 0.2 s, the effect is shown in Fig. 2, which is simulated by the software of Mathematica. For the cubital fossa, the thickness is 0.3 mm. When t ¼ 0:05, 0.1, 0.2 s, the effect is shown in Fig. 3 is showed the follow: For epidermal thickness of eyelid less than 0.1 mm, we choose a thickness of 0.07 mm; the temperature increase is shown in Fig. 4. For the normal epidermal thickness, z0 ¼ 0:1 mm. When t ¼ 0:05, 0.1, 0.2 s, the effect is shown in Fig. 5. Table 1 The relevant constants of all layers. i layer

ki ðW=m KÞ

ai ðm2 =sÞ

1

0.4

1:2  10-7

2

0.9

1:8  10-7

1.2

1. t = 0.05s 2. t = 0.1s 3. t = 0.2s

50 40 30 3

20 10

2 1 0.1

0.2 0.3 Distance /mm

0.4

0.5

Fig. 3. When t ¼ 0:05, 0.1, 0.2 s, the temperature increase of the cubital fossa.

45 Temperature increase /K

3.1. Internal temperature increase with time and depth changes

Temperature increase /K

60

Earlier, we calculated the heat transfer law in the skin radiated by laser. In this process, we approximate the epidermis as the first layer and the others under the epidermis as the second layer. Because of less blood flow in the epidermis, blood perfusion rate can be ignored, and the blood perfusion rate of the internal organization should be given consideration. The relevant constants [24] are shown in Table 1. The laser vaporization heat flux and absorption rate of epidermis is aQ0 ¼ 1:5  105 W=m2 . In different tissue, the thickness of skin is different. The thickness of skin is 0.5–4 mm (not including the subcutaneous adipose tissue). Epidermal thickness can vary to a great extent, the range being 0.07– 1.2 mm; the epidermis of the palm of your hand and foot are up to 0.8–1.4 mm, and the cubital fossa is 0.3 mm. Epidermal thickness of eyelid is less than 0.1 mm.

1

Fig. 2. When t ¼ 0:05, 0.1, 0.2 s, the skin temperature increase of the palm of your hand.

and erfc is the complementary error function.

3. Result and discussion

0.6 0.8 Distance /mm

1. t = 0.05s 2. t = 0.1s 3. t = 0.15s

40 35 30

3

25

2

20 15

1

0.05

0.10 Distance/mm

0.15

0.2

Fig. 4. When t ¼ 0:05, 0.1, 0.15 s, the temperature increase of the cubital fossa.

Analyzing the figures, we can summarize the differences as follows: (1) The thermal properties of biological tissue have a large effect on laser thermal interaction. The first layer’s thermal conductivity and thermal diffusivity are lower. The absorbed radiation cannot be quickly passed into an internal tissue, and resulted in the rapid increase in the temperature of the epidermal tissue; the second layer’s thermal diffusivity is higher than the first layer’s, and heat can be quickly delivered to the deep tissue, therefore, the epidermal tissue’s temperature increases slowly.

428

K.-W. Guan et al. / Optics & Laser Technology 43 (2011) 425–429

70 1. t = 0.05s 2. t = 0.1s 3. t = 0.2s

40 3

Temperature increase /K

Temperature increase /K

50

30 20

2 1

10 0

0.05

0.1

0.15 0.2 Distance /mm

0.25

0.3

Fig. 5. When t ¼ 0:02, 0.1, 0.2 s, the temperature increase of the normal epidermal thickness.

68 66 64 62

0.10

0.12

0.14

0.16 0.18 Time /s

0.20

0.22

0.24

Fig. 6. With changing time, different epidermal surface temperature changes: z0 ¼ 0:07 mm (red), z0 ¼ 0:1 mm (blue), z0 ¼ 0:3 mm (yellow), z0 ¼ 1:2 mm (gray).

(2) The epidermal thicknesses of different tissues are slightly different, using the same irradiation time; they are slightly different in terms of temperature increase, the thinner the epidermis, the higher the temperature increase. (3) The same epidermal thickness, and different irradiation times; the temperature increase on each layer is also quite different. With the increase of irradiation time, the temperature increase at all layers also increases significantly.

3.2. Determine the time of the beginning of the gasification of epidermis As the temperature increases, denaturation of proteins occurs, which leads to necrosis of tissue and cells. Thermal coagulation, including collagen hyalinization, collagen and muscle birefringence changes, tissue whitening and cell shrinkage, occurs at a higher temperature. When z ¼ 0, we can use Eq. (20) to get the expression of the epidermal surface excess temperature: ! rffiffiffiffiffiffiffiffiffiffi 1 aQ0 X 4a1 t ½2z0 ð1 þ nÞ2 T1 ð0; tÞ ¼ Bn þ 1 exp  k1 n ¼ 0 p 4a1 t ! 2z0 ð1 þ nÞ ½2z0 ð1þ nÞ erfc pffiffiffiffiffiffiffiffiffiffi 4a1 t " # ! rffiffiffiffiffiffiffiffiffiffi 1 X aQ0 4a1 t ð2nz0 Þ2 2nz0 2nz0 erfc pffiffiffiffiffiffiffiffiffiffi ; Bn þ exp  k1 n ¼ 0 p 4a1 t 4a1 t ð18Þ We assume that the surface temperature of biological tissue is 36 3 C and gasification begins at 100 3 C. The excess temperature of the palm of your hand, cubital fossa and the eyelids at the beginning of gasification is given as: T1 ð0; tÞ ¼ 100 3 C36 3 C ¼ 64 3 C. Fig. 6 shows the temperature rise with the changing time on the surface. From Fig. 6, we can get the time of different tissues at the beginning of gasification: z0 ¼ 0:07 mm (red), t ¼ 193 ms; z0 ¼ 0:1 mm (blue), t ¼ 174 ms; z0 ¼ 0:3 mm (yellow), t ¼ 188 ms; z0 ¼ 1:2 mm (gray), t ¼ 190 ms. We found: the epidermis of different tissues to achieve gasification need different times, and with the increase in epidermal thickness, the increase in time, when epidermal surface is gasification, is not monotonic. There is a critical time, which is the shortest time.

Fig. 7. 532 nm internal photon probability at different depths.

3.3. Distribution of heat injury The standard rate process model of tissue damage was introduced by Henriques and Moritz [15,16]. The damage parameter, which indicates the level of damage, was computed using the Arrhenius equation:

OðzÞ ¼ A

Z 0

1

 exp

 E0 dt; RTðz; tÞ

ð19Þ

where A ¼ 3:1  1098 is the molecular collision frequency factor, E0 ¼ 6:28  105 J=mol is denaturation activation energy, and R ¼ 8:3145 J=mol K is the universal gas constant. Henriques and Moritz assigned O ¼ 0:53 corresponding to a threshold of firstdegree burn (persistent but reversible erythema), O ¼ 1 to the threshold of second-degree burn (irreversible partial-thickness injury), and O ¼ 10; 000 corresponded to a threshold of thirddegree burn (irreversible full-thickness injury). For some proteins, when their temperatures reach 60 3 C, degeneration and irreversible thermal damage occur. According to the skin temperature change with time, we can determine the thermal damage of regional tissue at different times. When t ¼ 0:2 s, different tissues achieve different skin depths of thermal injury; we use the h instead of the skin depth: z0 ¼ 1:2, 0.3, 0.1 mm, hC 150, 120, 90 mm; t ¼ 0:15 s, z0 ¼ 0:07 mm, h C 90 mm.

K.-W. Guan et al. / Optics & Laser Technology 43 (2011) 425–429

3.4. Energy distribution In the derma, there are many independent vessels, we have investigated the change of energy aggradation density along with the depth of skin and the energy distribution of the two-layer skin model, which has cylindrical vessels in the derma. The thickness from epidermis to the dermis is 150 mm. This is shown in Fig. 7. The energy distribution of the epidermal tissue is similar to drawing a line in the figure. It starts gradually increasing from the organization and the surrounding environment of the contact interface and the energy absorption of the first order approximation curve of exponential function in the way of dramatic increase; there is a sharp decline in dermal tissue. The energy absorption of the epidermal tissue is far greater than the energy absorption dermis. Because of the difference of absorptance of the dermis and the epidermis, the energy distribution changes drastically.

4. Conclusion In this paper, learning from the multi-layer model of laser and biological tissue established by the predecessors, using the twolayer model to simulate the thermal interaction between laser and skin heat transfer, and studying the temperature changes of the skin, we get the following conclusions: (1) Through the Laplace transform, solution with time and depth change of a two-layer system’s biological tissue temperature distribution were obtained, and applied on the skin. (2) Through the analysis of skin surface excess temperature changes, we can deduce the temperature changes that have relation with the optical parameters and the thickness of skin, have a direct bearing on the incident laser intensity and duration and have a natural contact. (3) Different organizations have different epidermal thicknesses; we estimated the time when gasification is beginning and the situation of heat injury by the analytical solution, and found the time when gasification begins does not increase monotonically with the increase of epidermal thickness; there is a critical thickness. (4) The skins of different organizations use the same irradiation time; the depths of thermal injury is obviously different. However, they also do increase monotonically.

Acknowledgements This work was supported by funding from the Chinese National Science Foundations (No. 30470416; No. 30870582) and the National ‘985’ Project to C. Sun. References [1] Niemz MH. Laser-tissue interactions: fundamentals and applications, 2003. [2] Chen B, Thomsen SL, Thomas RJ, Welch AJ. Modeling thermal damage in skin from 2000-nm laser irradiation. J Biomed Opt 2006;11(6). 064028(1-15).

429

[3] Eggleston TA, Roach WP, Mitchell MA, Smith K, Oler D, Johnson TE. Comparison of two porcine (Sus scrofa domestica) skin models for in vivo near-infrared laser exposure. Comput Med 2000;50(4):391–7. [4] El-Adawi MK, Abdel-Naby MA, Shalabyt SA. Laser heating of a two-layer system with constant surface absorption: an exact solution. Int J Heat Mass Transfer 1995;38(5):947–52. [5] El-Adawi MK, Shalaby SA, Abdel-Naby MA. Laser heating of a two-layer system with temperature dependent front surface absorptance. Vacuum 1995;46(1):37–42. [6] De Felici M, Tartari A, et al. Effects induced by the fine structure of biological tissues on microbeam dose deposition: a numerical study. Nucl Instrum Meth Phys Res A 2005;548:59–64. [7] Handapangoda CC, Premaratne M. An approximate numerical technique for characterizing optical pulse propagation in inhomogeneous biological tissue. J Biomed Biotechnol 2008;2008:1–5. [8] Mainster MA, White TJ, Allen RG. Spectral dependence of retinal damage produced by intense light sources. J Opt Soc Am 1970;60(6):848–55. [9] Vassiliadis A, Christian HC, Dedrick KG. Ocular laser threshold investigations. USAF School of Aerospace Medicine, Stanford Research Institute, Brooks Air Force Base, TX, 1971. [10] Takata AN. Thermal model of laser-induced eye damage. Engineering Mechanics Division, IIT Research Institute, 10 West 35th, Chicago, IL, 1975. [11] Torres JH, Motamedi M. Experimental evaluation of mathematical-models for predicting the thermal response of tissue to laser irradiation. Appl Opt 1993;32(4):597–606. [12] Kim AD. Light propagation in biological tissue containing an absorbing plate. Appl Opt 2004;43(3):555–63. [13] Ma N, Gao X, Zhang XX. Two-layer simulation model of laser-induced interstitial thermo-therapy. Lasers Med Sci 2004;18:184–9. [14] Jiang SC, Zhang XX. Dynamic modeling of photothermal interactions for laser-induced interstitial thermotherapy: parameter sensitivity analysis. Lasers Med Sci 2005;20:122–31. [15] Henriques FC, Moritz AR. Studies of thermal injuries I: The conduction of heat to and through skin and the temperatures attained there in a theoretical and an experimental investigation. A J Pathol 1947;23:531C549. [16] Henriques FC, Moritz AR. Study of thermal injuries V., The predictability and the signifcance of thermally induced rate processes leading to irreversible epidermal injury. Arch Pathol 1947;43:489C502. [17] Tasakia A, Yamanakaa N, Kubo M, et al. Three-dimensional two-layer collagen matrix gel culture model for evaluating complex biological functions of monocyte-derived dendritic cells. J Immunol Meth 2004;287:79–90. [18] Loze MK, David Wright C. Temperature distributions in laser-heated biological tissue with application to birthmark removal. J Biomed Opt 2001;6(1):74–85. [19] Das M, Xu C, Zhu Q. Analytical solution for light propagation in a two-layer tissue structure with a tilted interface for breast imaging. Appl Opt 2006;45(20):5027–36. [20] Angell-Petersen E, Hirschberg H, Madsen SJ. Determination of fluence rate and temperature distributions in the rat brain; implications for photodynamic therapy. J Biomed Opt 2007;12(1). 014003(1-9). [21] Ao H, Xing D, Wei H, Gu H. In vitro study of thermal changes in optical properties of myocardium tissue with diffuse reflectance imaging. IEEE/ICME international conference on complex medical engineering, 2007, pp. 1042–7. [22] Stec B, Dobrowolski A. Paper Estimation of internal distribution of temperature inside biological tissues by means of multifrequency microwave thermograph. J Telecommun Inf Technol 2002;1:39–42. [23] Kim AD, Moscoso M. Light transport in two-layer tissues. J Biomed Opt 2005;10(3). 034015(1-10). [24] Cheong W-F, Prahl SA, Welch AJ. A review of the optical properties of biological tissues. IEEE J Quantum Electron 1990;26(12):2166–85. [25] Deng Z-S, Liu J. Analytical study on bioheat transfer problems with spatial or transient heating on skin surface or inside biological bodies. J Biomech Eng 2002;124:638–49. [26] Zhou J, Liu J, Yu A. Numerical study on the thawing process of biological tissue induced by laser irradiation. J Biomech Eng 2005;127:416–31. [27] Alexandrakis G, Busch DR, Faris GW, Patterson MS. Determination of the optical properties of two-layer turbid media by use of a frequency-domain hybrid Monte Carlo diffusion model. Appl Opt 2001;40(22):3810–21. [28] Diaz-Valdes SH, Aguilar G, Basu R, Lavernia EJ, Wong BJ. Modeling the thermal response of porcine cartilage to laser irradiation. In: Laser-tissue interaction XIII: photochemical, photothermal, and photomechanical, 2002. p. 47–56.