Mathematical modeling of the optimum pulse structure for safe and effective photo epilation using broadband pulsed light

The objective of this work is the investigation of intense pulsed light (IPL) photoepilation using Monte Carlo simulation to model the effect of the output dosimetry with millisecond exposure used by typical commercial IPL systems. The temporal pulse shape is an important parameter, which may affect the biological tissue response in terms of efficacy and adverse reactions. This study investigates the effect that IPL pulse structures, namely free discharge, square pulse, close, and spaced pulse stacking, has on hair removal. The relationship between radiant exposure distribution during the IPL pulse and chromophore heating is explored and modeled for hair follicles and the epidermis using a custom Monte Carlo computer simulation. Consistent square pulse and close pulse stacking delivery of radiant exposure across the IPL pulse is shown to generate the most efficient specific heating of the target chromophore, whilst sparing the epidermis, compared to free discharge and pulse stacking pulse delivery. Free discharge systems produced the highest epidermal temperature in the model. This study presents modeled thermal data of a hair follicle in situ, indicating that square pulse IPL technology may be the most efficient and the safest method for photoepilation. The investigation also suggests that the square pulse system design is the most efficient, as energy is not wasted during pulse exposure or lost through interpulse delay times of stacked pulses. PACS number: 87.10.Rt

types, namely free discharge, square pulse, close pulse stacking, and spaced pulse stacking, as illustrated in Fig. 1.
Using a computer simulation, this study investigates the effect these categories of pulse structure have on hair removal for critical evaluation of system performance. Our challenge in mathematical modeling the broad wavelength distribution in IPL systems is the requirement to incorporate the numerous optical parameters listed in Table 1, such as absorption coefficients, scattering coefficients, and anisotropy factors on light propagation within tissues. The absorption coefficients are tissue-and wavelength-dependant and thus cannot be averaged; therefore, we modeled the propagation in tissue of an IPL source simultaneously for each wavelength comprising the IPL spectrum. Mathematical models used to optimize selective photothermolysis in laser treatments have been widely published. (8)(9)(10)(11)(12) Monte Carlo simulations of light propagation in tissue are seen as the 'gold standard' to describe the transport of photons in biological turbid tissue (13) and have worked satisfactorily for some years, where each photon is followed through a propagation medium and its journey is directed by probability distributions of its interaction of scattering angles and absorption within the various tissue layers. In particular, Bjerring et al. (14) showed the importance of the pulse shape of a pulsed dye laser (PDL) using mathematical and clinical evidence. Steiner et al. (15) showed the effect of altering pulse shape and wavelengths of a laser system for epilation, and demonstrated that pulse shape bore significance in relation to efficacy and safety. Shafirstein et al. (16) presented a study in 2006 using mathematical modeling to aid vascular treatments using a dye laser, and recognized the temporal pulse shape as an important parameter affecting treatment effectiveness and safety outcomes/issues. Two recent studies have shown IPL technology to be a competitor to the laser in dermatology applications using modeling of blood vessels and hair follicles, respectively. (17)(18) Fig. 1. Time-resolved spectral images of the four categories of IPL systems currently available (intensity is arbitrary). (7)

Selective photothermolysis
Anderson and Parrish (19) first introduced the principle of selective photothermolysis applicable to the laser treatment of cutaneous lesions in 1983. Selective photothermolysis describes how a target can be selectively destroyed by using the appropriate wavelength, pulse duration, and energy, and that the optimum level of thermal absorption in the target should be achieved with the least interference to the surrounding non-target tissue. (19) In tandem with the principle of selective photothermolysis is the concept of thermal relaxation time. (20) Thermal relaxation time describes the time period for an object to lose 50% of its induced thermal energy to surrounding tissues. With the right combination of wavelength, energy, and pulse duration, it is therefore possible to target a hair follicle precisely without causing injury to the surrounding structures. One way to achieve greater injury to the hair follicle is by extending the pulse duration of the light exposure. The thermal relaxation time for hair follicles that are 200-300 μm in diameter is approximately 40 to 100 milliseconds. If pulse duration was the only factor, then the ideal pulse duration should lie between the thermal relaxation time for melanosomes found in the epidermis, which is approximately 3 to 10 milliseconds, and the thermal relaxation time for hair follicles. (19) Applying the theory of selective photothermolysis using a laser for cutaneous disorders is the same when using IPL technology. With respect to lasers, IPL systems emit multiple wavelengths of light within the range of 500 nm to 1100 nm at different intensities. (21) These could be modeled as 600 wavelengths of 1 nm bandwidth and varying radiant exposure simultaneously.

II. MAtErIALS And MEtHodS
The two-dimensional mathematical simulation described here was previously utilized for an earlier study modeling blood vessels, (14) but was modified to simulate the two-dimensional selective photothermolysis of a buried hair follicle using a typical IPL spectrum with a uniform spatial distribution. The authors of this study modeled the tissue matrix using a 5 μm cell size as a Cartesian coordinate, an 80 μm melanin rich epidermis with 5% melanin concentration, a melanin free dermis, and a buried hair follicle 200 μm below the surface of 200 μm diameter with a 30% melanin concentration extending to 2000 μm in length. The skin temperature was measured 20 μm below the skin surface; the hair follicle temperature was measured 10 μm into the side of the hair follicle at a depth of 2000 μm, as depicted in Fig. 2. The simulation is comprised of two interacting subroutines -a Monte Carlo photon distribution and a heat diffusion calculation -to compare four different pulse structures as produced by various IPL systems (shown in Fig. 3 and listed in Table 2). The model was calibrated against thermography measurements of human skin exposed to IPL exposure. These measurements used various skin tones (Caucasian, Indian, and Afro Caribbean) with ranging IPL exposure (2-20 J/cm 2 ) at a known distance from skin. The calculation of the energy distribution was achieved by using a Monte Carlo process. This computer simulation allows observation of comparative data that cannot be performed in skin. The Monte Carlo simulation of photons is a numerical calculation to compute the radiative transport in scattering and absorbing media. In the computations of photon pathways, the absorption coefficients (μ a ), the scattering coefficient (μ s ), and the anisotropy factor (g) of each material at the specific wavelength is considered. These respective values were taken from the literature (22,23) and the respective absorption coefficients at these wavelengths are shown in Fig. 4. The wavelength dependency of the optical parameters of the IPL spectrum complicates mathematical   modeling of light tissue interaction for IPL, as compared to monochromatic lasers. (17) For comparison of system design and efficiency, the total radiant exposure of the 1×10 9 photons that constitute the spectrum shown in Fig. 4 equals 10 J/cm 2 . This radiant exposure is above the published 5J/cm 2 threshold for follicular damage. (24) When filtered broadband IPL is delivered to the subject's skin, a fraction of the optical energy is reflected from the skin's surface while the rest is scattered beneath the skin's surface and is then absorbed in biological cells. Such absorption causes the hair shaft and hair bulb to heat up, resulting in damage to the follicle, shedding of the hair, and prevention of future regrowth. A custom Monte Carlo simulation package was created by the Faculty of Design and Engineering at Swansea Metropolitan University, UK, and implemented on a personal computer. Within this Monte Carlo simulation, the absorption coefficient of melanin, μ a , is computed and graphically shown in Fig. 4. The scattering coefficient (cm -1 ) is measured from experimental techniques originating from Mie scattering due to collagen fibres and from Rayleigh scattering due to small tissue structures, respectively. (25) Scattering in tissue by photons is characterized by the Henyey-Greenstein scattering phase function, (26,27,28,29) which is mathematically expressed in the form (1) P( ) 1 g 2 (1 g 2 2 g Cos( )) 3/2 The longitudinal angle of scattering, (θ), is characterized by the Henyey-Greenstein phase function and, since this phase function is forwardly biased, θ will be randomly influenced to reflect this characteristic. Applying the probability density function for the random number R to the above equation and solving for Cos(θ) leads the following expression: Cos( 2 1 g 2 ) 1 2g 1 g 2 1 g 2gR for g ≠ 0 (2) Fig. 4. Absorption coefficients of melanin, oxyhaemoglobin, and water. The IPL emission spectrum used for this evaluation is overlaid for Reference (6).
The azimuthal angle, ϕ, (Fig. 5) is uniformly distributed within the interval [0, 2π]. Its probability density function is constant and equals . Hence ϕ takes the form: A predefined number of photons are directed onto a skin model. These photons are given a weight (W), which is equal to 1. A photon starts at a boundary position with the tissue interface with a step size, anisotropy angle of rotation, and a deflection angle. At each interaction the photon delivers a portion of its energy, preset to values determined of tissue optical properties. The loss in energy at each of the ongoing interaction sites is determined by the weight that the photon deposits. The change in weight is defined by the following equation: (30,31,32) The photons are terminated when the value of ΔW is below a threshold. In this simulation the photons are terminated when one hundredth of the original weight remained after multiple interactions. Then after this photon is deposited in its final position, a new photon is released and the process is repeated. A sufficient number of photons are required to generate an absorption energy density matrix for the defined tissue configuration for absorption of the IPL power output. (30,32,33) This energy density matrix is used with a factor of total energy output within this defined area. (30,31) The absorbed energy at each Cartesian cell can be used as a heat source for thermal diffusion approximation. The time dependant heat flow equation is reduced to its 2-D form and subsequently evaluated over the discretized domain using the alternating direction implicit (ADI) method: where H is the volumetric distribution of energy calculated in the Monte Carlo simulation and stored in numabs array, and y and z are the radial and axial coordinates, respectively.

III. rESuLtS
The effect of delivered IPL pulse structure was studied on a buried hair follicle using the Cartesian model depicted in Fig. 2. The skin temperature was measured 20 μm below the skin surface; the hair follicle temperature was computed 10 μm into the side of the hair follicle at a depth of 2000 μm. This location is where pluripotential stem cells responsible for hair growth are thought to reside. It has been postulated that follicular damage would occur at temperatures higher than 70°C. (19) The data presented in Fig. 6(a) showing skin temperatures for all four systems indicates that the peak temperature for the free discharge system produced the highest absorbed temperatures and may cause greater patient discomfort as a result (22% greater than square pulse and pulse stacking techniques).Typically, free discharge systems tend to implement active or parallel cooling of the skin to prevent adverse reactions, such as prolonged erythema or hyperpigmentation. The higher peak power dissipated into the tissue matrix of the free discharge systems may explain the greater number of adverse reactions seen with these devices in practice. (34) Figure 6(b) shows free discharge, square pulse, and close pulse stacking systems to attain temperatures above 70ºC within the hair follicle, which is generally considered sufficient to denature pluripotential stem cells surrounding the hair follicle and supporting regrowth. The modeling of the spaced pulse stacking system suggests that the device would need to deliver a greater radiant exposure to attain similar results than the other systems, as much energy is lost during the prolonged off time between pulses. As a result of the increased energy necessary, spaced pulse stacking systems may require more internal water cooling and active or parallel skin surface cooling to generate and deliver energy safely.

IV. dIScuSSIon
Both the square pulse and close pulse stacking temporal profiles provide sustained heating to the follicle in a controlled repeatable and consistent delivery of optical energy. This controlled delivery should also provide repeatable clinical results shot-to-shot. Whereas the free discharge systems, although simple in their technology, inherently vary shot-to-shot, and may cause disparity during treatment and during the lifetime of the device. (35) Fig. 6. Computer modeling of the epidermis (a) using four types of commercially available temporal profiles with pulse durations of those measured in Reference (7), as depicted in Fig. 3. The temperature profile (b) of absorbed light of the four pulse profiles using computer modeling of an in situ hair follicle. The Monte Carlo model used for this study is two-dimensional in its representation of the light-tissue interaction. However, in reality the interaction of light with hair follicles is in three dimensions, where radial heat conduction is much faster than axial heat conduction. This may explain the retained temperatures within the temperature modeling and during exposure. The model is representative for a broad field as approximately the same amount of photons jump into as out of the azimuthal plane.
It has long been assumed that the optical properties of the various tissue layers do not change during exposure to light since the temperature increase is not sufficiently large, so constant thermal properties may be assumed. This may not be the case as during exposure, absorption of optical energy by chromophores may cause heating and mechanically modify some biological targets (for example, coagulation of blood vessels close to the surface). (36,37) The value of being able to compare the physics of light-tissue interaction with the ability to change modeled parameters individually is of significance to clinicians and light-based system designers. Repeatable simulations of light-tissue interaction may be the true value of Monte Carlo modeling; however, allowance should be made for some inaccuracies compared to reality. For example, the values of absorption (μ a ) and scattering (μ s ) of tissue layers are obtained experimentally and contain an uncertainty. (6) The mathematical model is a tool that can assist researchers and clinicians to identify positive theoretical models and avoid unnecessary experimental repetition. It facilitates the examination of models with hundreds of different parameters that would be practically impossible to conduct in a clinical setting. Mathematical modeling also allows investigations to be performed that cannot be conducted practically, such as viewing the deposition of photons through the skin and their final position depths in tissue.

V. concLuSIonS
Monte Carlo modeling is a technique that solves the radiative transfer equation using the probabilistic nature of photon interactions and has been used to simulate many such interactions. It may be efficiently implemented for complicated tissue geometries and without restrictions in optical properties. Furthermore, this technique provides accurate results for highly absorbing and scattering media at positions close to the surface, and can effectively handle anisotropic scattering.
The computer modeling of the various pulse profiles presented here indicates that square pulse and close pulse stacking systems are efficient in delivering an optimum dose of optical radiation to induce a sufficient thermal transient in the follicular structure for effective hair reduction, in accordance with the theory of selective photothermolysis. Pulses of light that are longer than the thermal relaxation time of melanosomes, approximately 3 to 10 ms, and shorter than the thermal relaxation time for hair follicles, approximately 40 to 100 ms, will efficiently heat the hair follicle whilst the epidermis is conducting heat to surrounding tissue. Therefore, as free discharge systems have short pulse durations of circa 4-6 ms, heat has insufficient time to transfer from the epidermis during the pulse, resulting in higher peak temperature induced in the skin compared with longer pulse duration pulse types.