PARAMETRIC STUDY OF THE DERMAL TEMPERATURE PROFILES DURING CRYOLIPOLYSIS By Dillon H. McClintock A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Mechanical Engineering – Master of Science 2021 ABSTRACT PARAMETRIC STUDY OF THE DERMAL TEMPERATURE PROFILES DURING CRYOLIPOLYSIS By Dillon H. McClintock Cryolipolysis is a noninvasive clinical procedure to locally reduce adipose tissue. During the procedure, paddles that are maintained at a prescribed temperature (commonly -14 ‰ and 3.1 ‰) are placed in good thermal contact with the skin. Cryolipolysis was inspired by panniculitis observed in a pediatric patient. The goal is to cool the adipose tissue to 10.4 ‰ in order to induce apoptosis in adipocytes, presumably by the crystallization of the intracellular triglycerides. The dermal cells between the cooled paddle and the adipocytes are colder than the adipocytes, but are less susceptible to death due to cooling. Thus, when properly administered, cryolipolysis leaves the dermis and epidermis unharmed. There are some adverse outcomes, e.g. transient neuropathy and rare adipose hyperplasia. The clinical cooling protocol has been based largely on animal experiments and subsequent clinical experience. A mathematical model could aid clinicians by providing insight into the temperature history of the tissues, potentially allowing optimization of the procedure. A model is here presented based on the Pennes bioheat transfer equation. Scaling the equation reveals two parameters in the Pennes equation: the Fourier number and a blood perfusion parameter. For optimizing treatments, the relationship between temperature and time required to trigger apoptosis in adipocytes needs to be quantified and coupled with the temperature model. ACKNOWLEDGEMENTS I would like to thank my thesis advisor Dr. Neil T. Wright for all of his support and genuine investment in my academic career. His drive and passion for learning makes him a great teacher and mentor. He made completing a whole thesis virtually not just possible, but also enjoyable. I would additionally like to thank the rest of my committee, Sara Roccabianca and Christina Chan, for their interest in my work. iii TABLE OF CONTENTS LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi KEY TO SYMBOLS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi CHAPTER 1 INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 CHAPTER 2 PROBLEM FORMULATION . . . . . . . . . . . . . . . . . . . . . 10 2.1 Anatomy of Dermis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.2 Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.2.1 Pennes Bioheat Equation . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.2.2 Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.2.3 Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.3 Solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.3.1 Early times, no perfusion – Error Function solution . . . . . . . . . . 14 2.3.2 Early Temperature Response of Perfused Tissue . . . . . . . . . . . . 15 2.3.3 Longer term solution for perfused tissue – Finite Depth, X11 solution 17 2.3.4 Green’s Function Solution . . . . . . . . . . . . . . . . . . . . . . . . 18 CHAPTER 3 RESULTS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 CHAPTER 4 DISCUSSION AND CONCLUSIONS . . . . . . . . . . . . . . . . . 39 BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 iv LIST OF TABLES Table 2.1: Model skin layer thicknesses. Nominal value of 13 mm used for adipose tissue for 15 mm over all model length [56]. . . . . . . . . . . . . . . . . . 10 Table 2.2: Property values used during study besides during sensitivity studies were one parameter was varied. Properties assumed uniform throughout all three layers [53, 46]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 Table 3.1: Minimum and maximum heat flux during treatment. Values are negative for this model but can be assumed to be positive for direct comparison to the manufacturer heat flux. Maximum value is the largest transient value with the minimum heat flux is the smallest flux during the steady state portion of both the temperature profile and the heat flux. The maximum heat flux is not measured at time zero, but one second into the time interval due to the heat flux going to infinity at zero. . . . . . . 36 v LIST OF FIGURES Figure 1.1: The anatomy of the skin that the problem model is based upon. The layers of skin are not distinct as assumed in the model. It can also be seen that there are other tissues and mechanisms present besides the skin tissue in the layers. Created with BioRender.com [10]. . . . . . . . . 2 Figure 1.2: Schemas of the two applicators used during CLL (a) paddle and (b) vacuum assisted. The vacuum assisted applicator pulls the skin into the applicator in order to increase contact conductance. The vacuum assisted may also restrict the local blood perfusion by increasing the pressure to the skin could increase efficacy and better possible cooling than the paddle applicator. . . . . . . . . . . . . . . . . . . . . . . . . . 4 Figure 1.3: The proposed CIF for each heat flux described by Zeltiq [41]. The relation between the CIF and heat flux is linear. CIF values: 21.5, 24.4, and 42 corresponding to heat flux values :-36.8 mWcm-2 , -43.8 mWcm-2 , -72.9 mWcm-2 [17, 42]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 Figure 2.1: Diagram model given 3 distinct layers of skin and a defined arterial tem- perature of Ta and a paddle temperature of T0 . L is the characteristic length with can be divided into Li where i=1,2,3 for the epidermis, der- mis, and subcutaneous, respectively. La denotes the length of muscle tissue included in the model until arterial temperature is reached. The dimensional coordinate is x̂. . . . . . . . . . . . . . . . . . . . . . . . . . 11 Figure 3.1: Comparison of emperature profiles as a function of treatment time for the error function. The Zeltiq paddle temperature of -14 ◦ C is used [41]. The length of the model is set to be 3 cm with all other properties uniform. A thermal conductivity of k=0.52 Wm-1 K-1 is used and blood perfusion is not included in the model. The time varies from 1 minute to 90 minutes. Due to the semi-infinite distant internal boundary con- dition, the profile does not reach steady state. . . . . . . . . . . . . . . . 21 Figure 3.2: Comparison of temperature profiles as a function of treatment time for the error function. The Lipocryo paddle temperature of 3.1 ◦ C is used [45]. The length of the model is set to be 3 cm with all other properties uniform. A thermal conductivity of k=0.52 Wm-1 K-1 is used and blood perfusion is not included in the model. The time varies from 1 minute to 90 minutes. Due to the semi-infinite distant internal boundary condition, the profile does not reach steady state. . . . . . . . 22 vi Figure 3.3: Comparison of the temperature profiles as a function of treatment time for the X10 Green’s function solution. The Zeltiq paddle temperature of -14 ◦ C is used [41]. The length of the model is set to be 3 cm with all other properties uniform. A thermal conductivity of k=0.52 Wm-1 K-1 and a blood perfusion value of wb = 0.0005 s-1 are used. The time varies from 1 minute to 90 minutes. The model gets closer to steady state approsximately 15 minutes with the 60- and 90-minute treatments being similar, but true steady state is not reached due to the semi-infinite internal distant boundary condition. . . . . . . . . . . . 23 Figure 3.4: Comparison of the temperature profiles as a function of treatment time for the X10 Green’s function solution. The Zeltiq paddle temperature of -14 ◦ C is used [41]. The length of the model is set to be 3 cm with all other properties uniform. A thermal conductivity of k=0.52 Wm-1 K-1 and a blood perfusion value of wb = 0.0005 s-1 are used. The time varies from 1 minute to 90 minutes. The model gets closer to steady state around 15 minutes with the 60- and 90-minute treatments being similar, but true steady state is not reached due to the semi-infinite internal distant boundary condition. . . . . . . . . . . . . . . . . . . . . 24 Figure 3.5: Comparison of the treatment time for the X11 Green’s function solution. The Zeltiq paddle temperature of -14 ◦ C is used [48]. The length for the distant boundary condition is set to be 3 cm with all other properties uniform. A thermal conductivity of k=0.52 Wm-1 K-1 and a blood per- fusion value of wb = 0.0005 s-1 are used. The time varies from 1 minute to 90 minutes. The model nearly reaches steady state after 15 minutes of treatment and does not change post 60 minutes. . . . . . . . . . . . . 25 Figure 3.6: Comparison of the treatment time for the X11 Green’s function solution. Lipocryo paddle temperature of 3.1 ◦ C is used [45]. The length for the distant boundary condition is set to be 3 cm with all other properties uniform. A thermal conductivity of k=0.52 Wm-1 K-1 and a blood per- fusion value of wb = 0.0005 s-1 are used. The time varies from 1 minute to 90 minutes. The model nearlyreaches steady state by 15 minutes of treatment and does not change after 60 minutes. . . . . . . . . . . . . . . 26 vii Figure 3.7: Comparison of the temperature solutions for the error function and Green’s functions of X10 and X11 boundary conditions at 60 minute cooling and a characteristic length of 3 cm. The Zeltiq paddle tem- perature of -14 ◦ C is used [41]. All parameter values are uniform and perfusion was set to the nominal value of wb = 0.0005 s-1 for the Green’s function solutions. The thermal conductivity was set to be k=0.52 Wm-1 K-1 for all solutions. The error function and the X10 Green’s function have a semi-infinite internal boundary condition resulting in a model with a distant boundary condition colder than the internal body temperature. The Green’s function of the X11 case keeps the distant boundary condition at the set internal arterial temperature resulting in a more realistic model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 Figure 3.8: Comparison of the temperature solutions for the error function and Green’s functions of X10 and X11 boundary conditions at 60 minute cooling and a characteristic length of 3 cm. The Lipocryo paddle tem- perature of 3.1 ◦ C is used [45]. All parameter values are uniform and perfusion was set to the nominal value of wb = 0.0005 s-1 for the Green’s function solutions. The thermal conductivity was set to be k=0.52 Wm-1 K-1 for all solutions. The error function and the X10 Green’s function have a semi-infinite internal boundary condition resulting in a model with a distant boundary condition colder than the internal body temperature. The Green’s function of the X11 case keeps the distant boundary condition at the set internal arterial temperature resulting in a more realistic model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 Figure 3.9: Study of sensitivity of the Green’s function of the X11 case to the com- putational length of the model given nominal perfusion of wb = 0.0005 s-1 and a treatment time of 60 minutes. The thermal conductivity is set to be k=0.52 Wm-1 K-1 . The Zeltiq paddle temperature of -14 ◦ C is used [41]. The grey portion of the graph highlights the adipose tissue layer ending at 1.5 cm into the skin. The lengths of 3 cm, 4.5 cm, and 6 cm show a similar temperature profiles up until the end of the adipose tissue layer. The temperature profile for the Green’s function of the X11 case is insensitive to model length after 3cm. . . . . . . . . . . . . . 29 viii Figure 3.10: Study of sensitivity of the Green’s function of the X11 case to the com- putational length of the model given nominal perfusion of wb = 0.0005 s-1 and a treatment time of 60 minutes. The thermal conductivity is set to be k=0.52 Wm-1 K-1 . The Lipocryo paddle temperature of 3.1 ◦ C is used [45]. The grey portion of the graph highlights the adipose tissue layer ending at 1.5 cm into the skin. The lengths of 3 cm, 4.5 cm, and 6 cm show a similar temperature profiles up until the end of the adipose tissue layer. The temperature profile for the Green’s function of the X11 case is insensitive to model length after 3cm. . . . . . . . . . 30 Figure 3.11: Sensitivity of the Green’s function with X11 boundaries to perfusion. Three values of perfusion were used of either half or double the nominal wb = 0.0005 s-1 . The Zeltiq paddle temperature of -14 ◦ C is used [41]. The treatment time varies from 1 minute to 60 minutes. The length is 3 cm and all other properties are uniform with the thermal conductivity being set to k=0.52 Wm-1 K-1 . The transient temperature profile is in- sensitive to the blood perfusion, but the steady state temperature profile has a larger difference between the upper and lower blood perfusion limits. 31 Figure 3.12: Sensitivity of the Green’s function with X11 boundaries to perfusion. Three values of perfusion were used of either half or double the nom- inal wb = 0.0005 s-1 . The Lipocryo paddle temperature of 3.1 ◦ C is used [45]. The treatment time varies from 1 minute to 60 minutes. The length is 3 cm and all other properties are uniform with the thermal conductivity being set to k=0.52 Wm-1 K-1 . The transient temperature profile is insensitive to the blood perfusion, but the steady state temper- ature profile has a larger difference between the upper and lower blood perfusion limits. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 Figure 3.13: The sensitivity of the Green’s function of the X11 case to the thermal conductivity. Values ±0.1 Wm-1 K-1 to a nominal value of 0.5 Wm-1 K-1 were used to assess the sensitivity at treatment times of 1 minute, 5 minutes, and 60 minutes. The Zeltiq paddle temperature of -14 ◦ C is used [41]. The computational length is 3 cm and the nominal blood perfusion value is wb = 0.0005 s-1 . All other properties are uniform. The temperature profile shows little sensitivity to thermal conductivity with the variation between values being small and constant throughout the transient and steady-state profiles. . . . . . . . . . . . . . . . . . . . 33 ix Figure 3.14: The sensitivity of the Green’s function of the X11 case to the thermal conductivity. Values ±0.1 Wm-1 K-1 to a nominal value of 0.5 Wm-1 K-1 were used to assess the sensitivity at treatment times of 1 minute, 5 minutes, and 60 minutes. The Lipocryo paddle temperature of 3.1 ◦ C is used [45]. The computational length is 3 cm and the nominal blood perfusion value is wb = 0.0005 s-1 . All other properties are uniform. The temperature profile shows insensitivity to thermal conductivity with the variation between values being small and constant throughout the transient and steady-state profiles. . . . . . . . . . . . . . . . . . . . . . 34 Figure 3.15: Sensitivity study of the Green’s function solution for X11 boundary condition to the paddle temperature. The model uses the nominal per- fusion value of wb = 0.0005 s-1 and a treatment time of 60 minutes. All other properties are constant with a characteristic length of 3 cm and a thermal conductivity of k=0.52 Wm-1 K-1 . T0 =-20 ◦ C is when tissue starts to develop necrotic behavior [49]. Bernstein explored pad- dle temperatures of -10 ◦ C and -5 ◦ C [9] while the two manufacturers, Zeltiq and Lypocryo, used paddle temperatures of -14 ◦ C [41] and 3.1 ◦ C [45], respectively. The temperature solution is sensitive to paddle temperature causing variance in the model curve up until the specified distant boundary condition. . . . . . . . . . . . . . . . . . . . . . . . . . 35 Figure 3.16: Comparison of the surface heat flux between the three model solutions over a 60 minute treatment time interval. The Zeltiq paddle tempera- ture of -14 ◦ C is used [41]. The perfusion is not included in the error function. All other parameters are kept constant. . . . . . . . . . . . . . 37 Figure 3.17: Isolation of 2500 seconds to 3600 seconds of -14 ◦ C paddle temperature heat flux comparison of error function and the Green’s function of the X10 case and X11 case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 Figure 3.18: Comparison of the surface heat flux between the three model solutions over a 60 minute treatment time interval. The Lipocryo paddle temper- ature of 3.1 ◦ C is used [45]. The perfusion is not included in the error function. All other parameters are kept constant. . . . . . . . . . . . . . 38 Figure 3.19: Isolation of 2500 seconds to 3600 seconds of 3.1 ◦ C paddle temperature heat flux comparison of error function and the Green’s function of the X10 case and X11 case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 x KEY TO SYMBOLS Symbols c tissue specific heat [J kg-1 ◦ C -1 ] cb blood specific heat [J kg-1 ◦ C -1 ] k tissue thermal conductivity [W m-1 K-1 ] t non-dimensional time [-] t̂ dimensional time [s] T temperature [◦ C] Ta arterial temperature [◦ C] To paddle temperature [◦ C] x non-dimensional depth [-] x̂ dimensional depth [m] Greek Letters α thermal diffusivity [m2 s-1 ] γ perfusion parameter [m-2 ] η dimensionless error function variable [-] ρ tissue mass density [kg m-2 ] ρb blood mass density [kg m-2 ] φ scaled temperature [-] wb blood perfusion [s-1 ] Subscripts i tissue layer xi CHAPTER 1 INTRODUCTION Cryolipolysis (CLL) is a cosmetic procedure that uses local, superficial cooling of the epidermis to induce apoptosis in the subcutaneous adipose tissue. Manstein [41] proposed pressing a paddle maintained at -14 ◦ C against the skin to cool the subcutaneous adipose tissue to a temperature of 9 ◦ C for duration of 30 to 60 minutes [48]. A candidate patient for CLL is seeking local fat reduction and would reduce their weight by other means such as dieting [38, 11]. CLL was inspired by the clinical observation of “popsicle panniculitis” [20] and the apparent susceptibility of adipocytes to cold temperatures, as compared with dermal cells. Epstein and Oren [20] discussed a child who presented with a 2-cm red nodule in her cheek after having had a popsicle rest against the cheek there for 5 minutes. Biopsies showed healthy epidermis and dermis but inflamed subcutaneous tissue. Similar responses of adipose tissue to cooling had been observed by Haxthausen [24]. CLL is a non-invasive alternative to liposuction. Liposuction removes subcutaneous adi- pose tissue via suction through a cannula that has been inserted in the area to be treated. The surgical nature of liposuction suggests that infections, permanent numbness, or other complications are more likely to occur from liposuction as compared with the non-invasive fat removal treatments [2]. CLL tends to be less painful than other noninvasive procedures, such as ultrasound or photodynamic treatments that rely on heating of adipose tissue [2, 44]. The skin is comprised of layers (see Fig. 1.1) of the epidermis, the dermis, and the hy- podermis, also called the subcutaneous adipose tissue. The thickness of these layers varies between patients and between different parts of the same patient. If the epidermis is cooled sufficiently, the subcutaneous tissue will also be cooled, although not to the same tempera- 1 Figure 1.1: The anatomy of the skin that the problem model is based upon. The layers of skin are not distinct as assumed in the model. It can also be seen that there are other tissues and mechanisms present besides the skin tissue in the layers. Created with BioRender.com [10]. ture. CLL relies on the epidermis and dermis remaining undamaged even though they are colder than the adipose tissue, which is thought to undergo apoptosis following CLL. The difference in death rate following cooling results from the triglyceride-rich adipocytes being more susceptible to cooling than are the water filled cells in the outer layers of the skin. It has been proposed that crystallization of the intracellular triglycerides of the adipocytes induce apoptosis at 10.4 ◦ C [48], while the outer cells with an aqueous cytosol can tolerate sub-zero temperatures [45, 48]. Apoptosis is a programmed cell death that may be triggered by naturally occurring processes, such as mitosis, or through external stressors such as, in this case, cold temperatures. The minimum time required for the adipose tissue to be cooled to trigger apoptosis is unknown. Loap and Lathe suggest that rather than apoptosis, a thermogenic mechanism is responsible for the loss in adipose tissue following CLL [40], but most researchers consider apoptosis to be the mechanism of cell death [48, 2, 4, 11, 26]. 2 There is not a single skin cell death temperature because multiple variables affect cell survival [49]. The environmental temperature and the duration for which the skin experiences this temperature are the two most significant variables that determine if necrosis will occur. The dermis can withstand temperatures of -10 ◦ C or warmer for extended times without necrosis or serious skin damage. If the environment is at -50 ◦ C or colder [49], the skin may quickly go necrotic. Nonfreezing injuries can occur in wet and cold conditions before actually freezing. Between these two temperature extremes, skin biopsies have shown necrosis for various temperatures, especially as the time intervals increased [49]. Commercial Devices In addition to temperature and duration of treatment, design of the cooling paddle may play a role in outcome. Two common designs are a flat paddle and a vacuum assisted applicator. The vacuum assisted applicator may restrict local blood perfusion by applying a greater perceived pressure to the skin. Quantification of blood perfusion is challenging although some typical values are available. There is a correlation between the amount of pressure applied and the effectiveness of the treatment post procedure [41]. It was suggested that by increasing the perceived pressure, the vacuum paddle has better thermal contact with the skin or with enough pressure, perfusion is reduced in the skin. Either might account for better cooling of the adipose tissue. The relationship between perfusion and pressure increases the efficacy of the treatment [48, 4]. It was shown that the greater the pressure using during treatment the better results were and results were seen for a longer period post-treatment. In both cases, gel is typically applied to the applicator to increase thermal contact conductance between the skin and paddle. Two commercially available devices for CLL are the Lipocryo [45] and the Zeltiq [41]. Each has a cooled paddle that is pressed against the skin as shown in Figure 1.2. The Lipocryo is said to cool the dermal surface to 3.1 ◦ C within four minutes at a maximum cooling heat flux of 1400 Wm-2 which is comparable to 140 mW cm-2 . Zeltiq cools to -14 3 Figure 1.2: Schemas of the two applicators used during CLL (a) paddle and (b) vacuum assisted. The vacuum assisted applicator pulls the skin into the applicator in order to increase contact conductance. The vacuum assisted may also restrict the local blood perfusion by increasing the pressure to the skin could increase efficacy and better possible cooling than the paddle applicator. ◦ C. Zeltiq also defines a proprietary “Cooling Intensity Factor” (CIF), which is related to the heat flux extracted from the surface. Published values suggest that the CIF has units of 21.5 for -36.8 mWcm-2 , 24.5, for -43.8 mWcm-2 , and 42 for -72.9 mWcm-2 [17, 42]. This appears to be a linear relationship (see Fig.1.3). Effectiveness of CLL Mostafa and Elshafey [44] found that CLL was more effective than laser lipolysis in reducing abdominal fat in adolescents. Neither treatment reduced the body mass index (BMI1 ) of the subjects, but locally reduced adipose tissue. Avram and Harry [4] reviewed the safety and efficacy of CLL based on studies that used paddle temperatures of from -8 ◦ C to -5 ◦ C and treatment times of up to 10 minutes. This is much warmer and for shorter treatment times than protocol used by other researchers. They found it to be 1 BMI is defined as the ratio of the mass in kg of the subject to the square of the height of the subject in m [22]. 4 Figure 1.3: The proposed CIF for each heat flux described by Zeltiq [41]. The relation between the CIF and heat flux is linear. CIF values: 21.5, 24.4, and 42 corresponding to heat flux values :-36.8 mWcm-2 , -43.8 mWcm-2 , -72.9 mWcm-2 [17, 42]. safe and effective, even though the mechanism was unknown. Ingargiola et al. [26] examined the available studies on the effectiveness of CLL and concluded that it is generally safe and resulted in substantial local reduction in adipose tissue as measured by calipers. At 16-weeks post treatment, a 2.8 mm reduction of the mean fat layer had been measured [64]. Zelickson, Burns, and Kilmer note that a standard treatment protocol is yet to be established. Bernstein et al. found that 94.4% of the time physicians were able to correctly identify blinded clinical pre- and post- treatment pictures of patients [8]. Paddle temperatures of -5 ◦ C and -10 ◦ C for 10 minutes was used. Bernstein et al. found that the localized fat reduction was maintained 6 to 9 years post-treatment in patients whose weight remained essentially unchanged [9]. Measurements taken over a 12-week span post procedure have shown no uptake of fat in liver. This holds even though the apoptosis of adipocytes releases triglycerides, presumably to the blood stream [36]. 5 Post-treatment erythema, bruising, and transient numbness are common side effects, of CLL. These usually resolve within two weeks, as does the less common (c.a., 0.1%) complaint of late-onset pain [38]. Coleman et al. [14] found that CLL produced short-term, reversible reduction in peripheral sensory nerve function. This loss lasted less than 4 weeks on average. Rare cases of paradoxical adipose hyperplasia (PAH) occur following CLL [58, 34]. PAH is the local growth of adipose tissue at the treatment site. PAH is not subject to spontaneous resolution and may require surgical intervention (i.e., liposuction). PAH is reported to occur in 0.0051% of cases [27], but others have suggested greater incidence [52, 33]. Keaney and Naga [33] observed that PAH occurs more commonly in men. Although generally done safely, without damage in the epidermis or dermis, there are some examples of the treatment going wrong even when done by a certified technician. A beautician performed CLL on a woman who then presented at a hospital with burns that required surgical repair [7]. Benoit and Modarressi [7] describe treatment to correct full-thickness frostbite following cryolipolysis treatment by non-medically trained operators. Leonard et al. [39] discusses the case of a patient who suffered full thickness burns after self-administering CLL using dry ice. These injuries also required surgical intervention. The paddle temperature recommended for CLL is at the warmer end of the necrotic spectrum, as mentioned above [49]. Derrick et al. reviewed the available literature on the safety and efficacy of CLL [17]. They suggested that future studies should focus on quantifying treatment parameters as reflected in treatment outputs. Many of the reports on the efficacy of CLL lack rigor in methodology or reporting of results [11]. There is uncertainty in what constitutes success from a clinical perspective and a lack of rigorous study regarding the mechanism of CLL [2, 11]. Garibyan et al. [23] measured the change in tissue volume using calipers as well as an optical system, but the clinical standard of calipers is most commonly used [26]. 6 Modeling Tissue Temperature Many studies have modeled the temperature of the skin, but none have considered CLL. Those focused on hypodermis have focused on cryosurgery, for which the water-rich dermis undergoes freezing [28, 47]. Baldry et al. [5] modeled the cooling of the scalp in a procedure designed to prevent hair loss in patients undergoing chemotherapy using a numerical model. The target temperature for stimulating follicles was 22 ◦ C. The cooling power of thermoelectric coolers were kept to less than 3175 W m-2 . Modeling of perfused tissue is usually done starting with the Pennes bioheat equation [57, 60]. The Pennes bioheat transfer equation may be written for 1D rectangular coordinates and assuming uniform properties, as ∂T ∂ 2 (T ) ρc =k + (ρb cb ) wb (Ta − T ) (1.1) ∂ t̂ ∂ x̂2 | {z } blood perfusion term where x̂ is depth from surface, t̂ time, T and Ta are tissue and arterial temperatures, respec- tively, as seen in Figure 2.1, k, c, and ρ are the tissue conductivity, specific heat, and density, respectively, cb is the blood specific heat, ρb the blood density, and wb is the blood perfusion rate with net units of s-1 . The Pennes equation is the diffusion equation with a term included to model blood perfusion. This makes the equation a transient fin equation. Some modelers have suggested using the Catteneo constitutive model for conduction in tissue rather than the Fourier equation [55, 54]. This is based on the experiments in wet sand [55] and in processed meat [54]. This would change the Pennes equation to ∂ 2T ∂T ∂ 2 (T ) τ + ρc = k + (ρb cb ) wb (Ta − T ) (1.2) ∂ t̂2 ∂ t̂ ∂ x̂2 with τ being a relaxation time. Flaws in the experimental design of [55] and [54] have been demonstrated by Scott et al. [51] and Herwig and Beckert [50]. As yet, no reliable experiment 7 using tissue has demonstrated the appropriateness of a non-zero τ , so the Pennes equation will be used here. Project scope While multiple clinical and laboratory studies have examined the efficacy and safety of CLL, no study has quantified the sensitivity of the temperature field to the several clinical parameters. The clinically controllable parameters include paddle tempera- ture and treatment duration. Parameters that are not controlled include tissue properties of perfusion (which may be affected by the pressure of the applicator bit is a challenge to quan- tify), thermal conductivity, and thermal capacity. The sensitivity of the controllable and uncontrolled parameters will be examined here through calculations of the 1D temperature field calculated using Green’s functions. The Greens function solution will be validated using the usual error function solution for the temperature in a semi-infinite domain. In addition, published studies have assumed that the heat flux at the skin surface is constant, although it varies with time as the temperature field is established. Here, the heat flux removed from the skin surface will be examined as a function of time. The results of this study should guide clinicians in the application of CLL. It will provide a starting point for optimizing treatment based on patient specific parameters of such as skin thickness and hydration. There is no treatment protocol for CLL. This is dangerous for a few reasons. Prolonged exposure to low temperatures can cause permanent numbness, frostbite “burning” of the skin, hard lumps, and extreme redness and swelling [2, 64]. Without a set treatment temperature or procedure to follow it is possible for the treatment to not be as effective as it should be or damaging to the skin. A model of the thermal profile of the skin can be used to accurately design a treatment procedure to best fit the condition of the client and to ensure their safety. 8 Besides the implications of optimizing CLL, the creation of an accurate thermal profile of the skin allows for exploration into other medical procedures that are based upon cooling and non-freezing injuries. Frostbite, trench foot, and chilblain are a couple of common non- freezing injuries and frostbite is also a possible negative side effect of CLL. Frostbite does not require the skin to actually reach freezing, but to get to at least a cold enough temperature to send cells into apoptosis or if cold enough, necrosis. Local hypothermia is also used in many medical procedures or therapies so being able to derive and understand the temperature profile of the skin as well as how it would react is important. Improved understanding of the response of adipose tissue to environmental changes may allow for better understanding of fat metabolism which, is important for shedding light on obesity in terms of fat storage and reaction to environmental changes. Creating a thermal profile for the skin extends beyond the reach of CLL. An accurate thermal profile of the skin allows for exploration into both freezing and non-freezing injuries as well as medical procedures that use non-freezing tissue reactions. Sensitivity studies were done for the characteristic length, blood perfusion, paddle temperature, treatment time, and thermal conductivity. 9 CHAPTER 2 PROBLEM FORMULATION 2.1 Anatomy of Dermis The skin is composed of three primary layers: the epidermis, dermis, and subcuta- neous(also called the hypodermis). The adipose tissue that CLL aims to remove is in the subcutaneous layer, which is the deepest in the tissue (see Fig. 1.1). In CLL, a paddle is placed at the surface of the skin and is maintained at -14 ◦ C for Zeltiq and 3.1 ◦ C for Lipocryo. Because a thermally conducting gel is applied to the skin, it is assumed that there is perfect contact between the paddle and the outer most point of skin. Adipose tissue thickness ranges from 1.6 to 30.6 mm [56] for normally weighted to obese subjects. This is taken from randomly selected sites on the body giving an average for the whole body. 2.2 Model The farthest point of the model can either be assumed to be the same as the arterial of internal temperature of the patient since the model has a fixed length and is not semi- infinite. Inside each layer of skin, the same equation is used throughout assuming that each Range Value Used Layer (mm) (mm) Epidermis 0.05 - 1.5 1 Dermis 0.6 - 3 2 Hypodermis 1.6 - 30.6 13 Table 2.1: Model skin layer thicknesses. Nominal value of 13 mm used for adipose tissue for 15 mm over all model length [56]. 10 Figure 2.1: Diagram model given 3 distinct layers of skin and a defined arterial temperature of Ta and a paddle temperature of T0 . L is the characteristic length with can be divided into Li where i=1,2,3 for the epidermis, dermis, and subcutaneous, respectively. La denotes the length of muscle tissue included in the model until arterial temperature is reached. The dimensional coordinate is x̂. layer has constant perfusion, and the properties are not changing. This is based upon the Pennes Bio Heat Transfer equation. Between each layer of skin an energy balance equation was done assuming that there is also perfect contact between layers and the skin is the same temperature at the start and end of either tissue layers that is in contact. 2.2.1 Pennes Bioheat Equation The Pennes bioheat transfer equation is commonly used to model the temperature in blood perfused tissue [57, 60]. The equation was developed using the 1D diffusion equation with 11 the inclusion of a blood perfusion term and terms to account volumetric heat generation due to metabolic activity or external sources, such as therapeutic ultrasound. Pennes initially determined values of perfusion by fitting temperature measurements in the forearms of prone, resting subjects [57]. It was later determined that the perfusion values were in error due to an error in the assumed value of the thermal conductivity of the tissue [60]. While that sensitivity to the difficult to measure blood perfusion, along with some theoretical questions [12], the Pennes equation remains the equation used to model the temperature profile in perfused tissue [60]. The Pennes bioheat transfer equation may be written for 1D rectangular coordinates and assuming uniform properties, as ∂T ∂ 2 (T ) ρc =k + (ρb cb ) wb (Ta − T ) (2.1) ∂ t̂ ∂ x̂2 | {z } blood perfusion term with x̂ depth from surface, t̂ time, T and Ta are tissue and arterial temperature respectively, k is the tissue conductivity, c is the specific heat, ρ is the mass density, cb is the specific heat of the blood, ρb the mass density of blood, and wb is the blood perfusion rate, with net units of s-1 .Typical values of wb are between 0.0001 and 0.0005 s-1 [43]. Since the epidermis is a very thin layer and the distinction between the dermis and hypodermis layers is unclear, it is reasonable to model the skin as one layer. 2.2.2 Boundary Conditions The boundary and initial conditions may be written as x̂ = 0 : T = T0 x̂ → ∞ : T →Ta or x̂ = L : T =Ta t̂ = 0 : T =Ta 12 Parameter Range Value Used k (W m-1 K-1 ) 0.16 - 0.52 0.52 ρ (kg m-2 ) 850 - 1270 1200 c (J kg-1 ◦ C -1 ) 2288 - 3770 3400 ρb (kg m-2 ) 850 - 1270 1200 cb (J kg-1 ◦ C -1 ) 2288 - 3770 3400 wb (s-1 ) 0 - 0.00125 0.0005 Ta (◦ C) - 37 T0 (◦ C) - -14, 3.1 Table 2.2: Property values used during study besides during sensitivity studies were one parameter was varied. Properties assumed uniform throughout all three layers [53, 46]. The distant boundary can be treated as an infinite boundary during the initial response. As treatment times increase, specifying the core body temperature of 37 ◦ C at a finite length is the appropriate distant boundary condition [28]. The appropriate value of L is examined in the results section. 2.2.3 Properties All three layers of tissue would have slightly different parameters, but due to differences between patients and inconsistency in measuring parameters, uniform properties are were used. Values used throughout the model are seen in Table 2.2 besides when property values were varied for sensitivity studies. 2.3 Solutions The error function, and Green’s function with two different sets of boundary conditions were used. The error function does not include perfusion and, with a semi-infinite distant internal boundary condition, only functions for early times; however, this works well as a check for the other solutions. The Green’s function of the X10 boundary conditions includes the two parameters of α and γ once non-dimensionalized meaning that it includes perfusion. This solution technique still does not include a finite length and instead a semi-infinite distant 13 boundary condition meaning that it is still only reasonable for early temperature responses. The Green’s function solution with X11 boundary conditions includes perfusion and has a finite length with a specified temperature at the depth L̂. 2.3.1 Early times, no perfusion – Error Function solution For the case of no perfusion, wb is zero and the Pennes equation simplifies to the diffusion equation ∂T ∂ 2 (T ) ρc =k (2.2) ∂ t̂ ∂ x̂2 Because this solution considers the initial temperature response of the tissue, the range is considered semi-infinite with the distant boundary condition being the core body tempera- ture. The boundary and initial conditions may be written as x̂ = 0 : T = T0 x̂ → ∞ : T →Ta t̂ = 0 : T =Ta After a sudden change in the boundary temperature, the familiar error function solution is   T − T0 x̂ φ= = erf √ (2.3) Ta − T0 4αt̂ This has been left in terms of dimensional x̂ and t̂ because the analysis naturally forms   the dimensionless variable η = √ x̂ . The heat flux during this initial response may be 4αt̂ determined using Fourier’s Law. Thus,  2 ∂T k(T0 − Ta ) −x k(T − Ta ) q000 = −k = √ exp = √0 (2.4) ∂x x=0 παt 4αt x=0 παt This solution provides a check of the solutions that include perfusion. 14 2.3.2 Early Temperature Response of Perfused Tissue Including blood perfusion, that is a wb 6= 0, Eq. 2.1 has the form of the equation for the transient temperature response of a fin. The initial temperature response of perfused tissue retains the semi-infinite domain of the error function solution. For the perfused case, it is useful to scale Eq. 2.1 by defining a perfusion parameter as ρ c  b b γ= wb (2.5) k with the thermal diffusivity being α = k/ρcp ; γ has units of m-2 . Equation 2.1 becomes 1 ∂θ ∂ 2θ = − γθ (2.6) α ∂ t̂ ∂ x̂2 with the two parameters of α and γ. The boundary conditions are x̂ = 0 : T = T0 x̂ → ∞ : T →Ta t̂ = 0 : T =Ta Comparing Eq. 2.6, γ plays the role m2 in the usual notation for the fin equation (e.f. √ Bergman & Lavine). Scaling the independent variables of t = αγ t̂ and x = γ x̂ yields ∂φ ∂ 2φ = −φ (2.7) ∂t ∂x2 The boundary conditions are x = 0: φ=1 x → ∞: φ→0 t = 0: φ=0 Transformation of Fin Equation Using the transformation for a transient fin equation to a diffusion equation [15] φ = W exp [−t] 15 Eq. 2.7 can be rewritten as ∂W ∂ 2W = (2.8) ∂t ∂x2 The boundary conditions transform to x=0: W = exp [t] (2.9) x→∞: W =0 t=0: W =0 where the boundary condition at x = 0 becomes time dependent instead of constant. Equa- tion 2.8 may be solved using Green’s functions, where the GF equation [15] reduces to Z t  h i ∂G  W (x, t) = −α dτ − exp m ατ 2 (2.10) τ =0 ∂x for the case of no initial temperature profile, no heat generation, and only one inhomogeneous type I boundary condition at the origin. The appropriate Green’s function for the X101 case is [15] x2   ∂GX10 x − = exp − (2.11) ∂n0 x0 =0 n 3 1/2 o 4α (t − τ ) 4π [α (t − τ )] Substituting Eq. 2.11 into 2.10 and integrating, the solution for the non-dimensional tem- perature is √     1 −x x φ = e −1 − erf t− √ 2 2 t √    x − ex 1 − erf t+ √ (2.12) 2 t which can be written in terms of dimensional temperature as T (x, t) = T0 + (Ta − T0 ) √ √        1 −x x x x ∗ e −1 − erf t− √ − e 1 − erf t+ √ (2.13) 2 2 t 2 t 1 X10 is the numbering system used by Cole et al. citecolebeck2011 to categorize boundary value problems with the diffusion equation. The X refers to rectangular coordinates, the 0 refers to an infinite boundary, 1 refers to a boundary condition of the 1st kind (specified temperature), 2 to a specified heat flux, and 3 for a specified convective surface condition. 16 √ Substituting for the dimensional independent variables of x = x̂ γ and t = γαt̂ leads to T (x, t) = T0 + (Ta − T0 ) 1 −√γ x̂   q  x̂ ∗ e −1 − erf γαt̂ − √ 2 2 αt̂ √  q  x̂ − e γ x̂ 1 − erf γαt̂ + √ (2.14) 2 αt̂ The heat flux at the boundary may be calculated again by using Fourier’s law. The gradient of temperature is √ 2 √ 2 √ x̂ √ x̂ " − γ x̂− αγ t̂− √ γ x̂− αγ t̂+ √ ∂T 1 e 2 αt̂ e 2 αt̂ = (Ta − T0 ) √ + √ (2.15) ∂ x̂ 2 παt̂ παt̂  !# √ −√γ x̂  x̂ i √ √γ x̂  q hq x̂ + γe 1 + erf αγ t̂ − √ − γe erfc αγ t̂ + √ 2 αt̂ 2 αt̂ Then, the heat flux evaluated at the cooled surface, x̂ = 0, is ! e−(αγt) √  √  q000 = k(T0 − Ta ) √ + γ erf αγt (2.16) παt Equation 2.16 reduces to Eq. 2.4 when there is no perfusion. 2.3.3 Longer term solution for perfused tissue – Finite Depth, X11 solution As time progresses, the temperature disturbance from the surface of the X10 solution reaches deeper into the tissue. The human body has an active mechanism of thermal regula- tion and the tissue temperature in the body is kept nearly constant. Here, this is simulated by assuming that the temperature is fixed to the core-body temperature at some depth L̂ [28]. The analysis for the longer-term temperature response is similar to the X10 case but uses a different Green’s function. Because the depth of the core body temperature is not fixed, the sensitivity of the solution to the choice of L̂ will be examined. The governing 17 equation remains Eq. 2.1 and the initial and boundary conditions for the X11 case become x̂ = 0 : T = T0 x̂ = L̂ : T = Ta t̂ = 0 : T = Ta In terms of φ, the initial and boundary conditions become x̂ = 0 : φ = 0 (2.17) x̂ = L̂ : φ = 1 t̂ = 0 : φ = 0 2.3.4 Green’s Function Solution Equation 2.8 subject to conditions 2.17 may be solved using Green’s functions. All of the independent variables are still dimensional, but the carats are here omitted. The end result will show The GF equation (p. 66, CBHL) [15] becomes 2 Z t " # ∂GX11 θ0 eγατ X W (x, t) = −α dτ (2.18) τ =0 ∂x0i x0 =0 i=1 with the boundary conditions of x=0: W =0 (2.19) x=L: W = exp [t] t=0: W =0 For the X11 case, Cole et al. [15] gave ∞ ∂GX11 2π X −m2 π2 α(t−τ )/L2  x − = e m sin mπ (2.20) ∂n0 x0 =0 L2 L m=1 18 Since the boundary condition is homogeneous at x = 0, that is f (τ ) = 0 at x = 0, the integral vanishes there an is not written. Then, the derivative of the Green’s function at the x = 0 boundary can be written in dimensional coordinates ∞  x ∂GX11 2π X −n2 π2 α(t−τ )/L2 − = e n sin nπ (2.21) ∂n0 x0 =0 L̂2 L n=1 So, substituting into Eq. 2.18 gives Z t ∞ γατ 2π X −n2 π2 α(t−τ )/L2  x W (x, t) = α e e n sin nπ dτ (2.22) τ =0 L2 L n=1 Integrating gives ∞    nπx  X γαt −n 2 π 2 αt/L2 n W (x, t) = 2π e −e sin (2.23) L γL + n2 π 2 2 n=1 The inverse transformation gets to the temperature excess of φ(x, t) = e−γαt W (x, t) ∞   −γαt X γαt −n 2 π 2 αt/L2  nπx  n = 2πe e −e sin (2.24) L γL + n2 π 2 2 n=1 which may be simplified to ∞   X −γαt −n 2 π 2 αt/L2  nπx  n φ(x, t) = 2π 1−e e sin (2.25) L γL + n2 π 2 2 n=1 This may be written in terms of temperature as T (x, t) = T0 + (Ta − T0 ) ∞   ! 2 2 2  nπx  n 1 − e−γαt e−n π αt/L sin X 2π (2.26) L γL2 + n2 π 2 n=1 In order to calculate the heat flux at the surface of the cooled skin, the temperature gradient is taken to be ∞   ! ∂T 2 2 2 nπ  nπx  n 1 − e−γαt e−n π αt/L X = (Ta − T0 ) 2π cos (2.27) ∂x L L γL2 + n2 π 2 n=1 Then at x = 0, the heat flux becomes ∞  ! n2 π  00 X −γαt −n 2 π 2 αt/L2 q0 = k(T0 − Ta ) 2π 1−e e (2.28) γL3 + n2 π 2 L n=1 19 CHAPTER 3 RESULTS To aid in clinical outcomes, the temperature and duration of a CLL treatment are exam- ined in terms of sensitivity to various parameters such as blood perfusion, tissue thermal conductivity, and characteristic length. Due to the relationship between tissue density and specific heat, the product of these two parameters is kept constant. The error function does not include perfusion, but is well known and is used to verify the other models. The X10 case includes perfusion, but at longer times, has a non-physiological response because the infinite boundary condition neglects the core body temperature. The X11 case also includes perfusion and has a more realistic boundary condition deep in the tissue. While the skin has three layers, it is modeled as one in this case due to the thinness of the epidermis and the uncertainty in the distinction between the dermis and hypodermis layers in terms of thermal properties. This assumption will be checked in calculations. All of the analysis was done non-dimensionally and re-dimensionalized for paddle temper- atures of -14 ◦ C (Zeltiq [41]) or 3.1 ◦ C (Lipocryo [45]) for comparison of clinical protocols. Typical clinical protocols use treatment times of 30 to 60 minutes. For each plot, aside from direct comparison of paddle temperature, there is a -14 ◦ C and 3.1 ◦ C paddle temperature version to aid in examination of the depth of the 10.4 ◦ C isotherm. Time Evolution The solutions can be compared to one another over the clinically treat- ment period of 30 to 60 minutes total. Early treatment lengths and a 90-minute treatment length were included to watch the evolution of the temperature profile over and past the treatment time. The treatment time greatly affects the temperature profile for each of the three solutions examined in this study. It affects the error function and Green’s function 20 more since the distant boundary condition is taken to be semi-infinite and then results in an unreasonably cold temperatures in the deep tissue. Figure 3.1: Comparison of emperature profiles as a function of treatment time for the error function. The Zeltiq paddle temperature of -14 ◦ C is used [41]. The length of the model is set to be 3 cm with all other properties uniform. A thermal conductivity of k=0.52 Wm-1 K-1 is used and blood perfusion is not included in the model. The time varies from 1 minute to 90 minutes. Due to the semi-infinite distant internal boundary condition, the profile does not reach steady state. 21 Figure 3.2: Comparison of temperature profiles as a function of treatment time for the error function. The Lipocryo paddle temperature of 3.1 ◦ C is used [45]. The length of the model is set to be 3 cm with all other properties uniform. A thermal conductivity of k=0.52 Wm-1 K-1 is used and blood perfusion is not included in the model. The time varies from 1 minute to 90 minutes. Due to the semi-infinite distant internal boundary condition, the profile does not reach steady state. The error function predicts initial temperatures with the solution not reaching steady state due to the semi-infinite distant boundary condition (see Fig. 3.1 for -14 ◦ C and 3.2 for 3.1 ◦ C). The Green’s function of the X10 boundary condition case is similar to the error function with the semi-infinite internal boundary condition, but perfusion is included. Due to the distant internal boundary condition omitting the core body temperature, the overall temperature profile is not anatomically correct. The Green’s function of the X11 case also includes perfusion, but the core body temperature is assumed at the boundary distant from the skin surface. 22 Figure 3.3: Comparison of the temperature profiles as a function of treatment time for the X10 Green’s function solution. The Zeltiq paddle temperature of -14 ◦ C is used [41]. The length of the model is set to be 3 cm with all other properties uniform. A thermal conductivity of k=0.52 Wm-1 K-1 and a blood perfusion value of wb = 0.0005 s-1 are used. The time varies from 1 minute to 90 minutes. The model gets closer to steady state approsximately 15 minutes with the 60- and 90-minute treatments being similar, but true steady state is not reached due to the semi-infinite internal distant boundary condition. 23 Figure 3.4: Comparison of the temperature profiles as a function of treatment time for the X10 Green’s function solution. The Zeltiq paddle temperature of -14 ◦ C is used [41]. The length of the model is set to be 3 cm with all other properties uniform. A thermal conductivity of k=0.52 Wm-1 K-1 and a blood perfusion value of wb = 0.0005 s-1 are used. The time varies from 1 minute to 90 minutes. The model gets closer to steady state around 15 minutes with the 60- and 90-minute treatments being similar, but true steady state is not reached due to the semi-infinite internal distant boundary condition. 24 Figure 3.5: Comparison of the treatment time for the X11 Green’s function solution. The Zeltiq paddle temperature of -14 ◦ C is used [48]. The length for the distant boundary condition is set to be 3 cm with all other properties uniform. A thermal conductivity of k=0.52 Wm-1 K-1 and a blood perfusion value of wb = 0.0005 s-1 are used. The time varies from 1 minute to 90 minutes. The model nearly reaches steady state after 15 minutes of treatment and does not change post 60 minutes. 25 Figure 3.6: Comparison of the treatment time for the X11 Green’s function solution. Lipocryo paddle temperature of 3.1 ◦ C is used [45]. The length for the distant boundary condition is set to be 3 cm with all other properties uniform. A thermal conductivity of k=0.52 Wm-1 K-1 and a blood perfusion value of wb = 0.0005 s-1 are used. The time varies from 1 minute to 90 minutes. The model nearlyreaches steady state by 15 minutes of treatment and does not change after 60 minutes. For the Green’s function of the X11 case (see Fig. 3.4 and 3.5), the temperature profile fairly quickly cools to reach almost steady state by 15 minutes. After 30 minutes the tem- perature profile is nearly stationary with the 60- and 90-minute treatment times essentially matching. For Zeltiq, the 15-minute treatment time reaches the target temperature of 10.4 ◦C has been reached for the first 0.8 cm of the adipose tissue. This suggests that the treat- ment needs to be longer than 15 minutes, but the amount of time for crystallization to occur is uncertain. Model Comparison The longest treatment time is typically 60 minutes and is used here to compare the three model solutions. The Zeltiq paddle temperature results in a colder temperature profile of all solutions (see Fig. 3.7) as compare to the Lipocryo paddle tem- 26 perature (see Fig. 3.8), which is expected. The error function does not include perfusion in the solution and thus cools faster than either of the Green’s function solutions. The error function and the Green’s function of X10 case both have a semi-infinite boundary condition at the distant boundary resulting an internal temperature colder than the arterial temper- ature attained deep the body. Perfusion in the X10 case causes the profile to be warmer than the error function solution. The Green’s function of the X11 case allows for the distant boundary condition to be set at the arterial temperature and the profile maintaining a higher and more realistic temperature. Figure 3.7: Comparison of the temperature solutions for the error function and Green’s functions of X10 and X11 boundary conditions at 60 minute cooling and a characteristic length of 3 cm. The Zeltiq paddle temperature of -14 ◦ C is used [41]. All parameter values are uniform and perfusion was set to the nominal value of wb = 0.0005 s-1 for the Green’s function solutions. The thermal conductivity was set to be k=0.52 Wm-1 K-1 for all solutions. The error function and the X10 Green’s function have a semi-infinite internal boundary condition resulting in a model with a distant boundary condition colder than the internal body temperature. The Green’s function of the X11 case keeps the distant boundary condition at the set internal arterial temperature resulting in a more realistic model. 27 Figure 3.8: Comparison of the temperature solutions for the error function and Green’s functions of X10 and X11 boundary conditions at 60 minute cooling and a characteristic length of 3 cm. The Lipocryo paddle temperature of 3.1 ◦ C is used [45]. All parameter values are uniform and perfusion was set to the nominal value of wb = 0.0005 s-1 for the Green’s function solutions. The thermal conductivity was set to be k=0.52 Wm-1 K-1 for all solutions. The error function and the X10 Green’s function have a semi-infinite internal boundary condition resulting in a model with a distant boundary condition colder than the internal body temperature. The Green’s function of the X11 case keeps the distant boundary condition at the set internal arterial temperature resulting in a more realistic model. Characteristic Length Sensitivity The temperature profile of the Green’s function of the X11 case is sensitive to the computational length, up to a point. These are compared in Figs. 3.9 and 3.10. The 3 cm, 4.5 cm, and the 6 cm lengths share similar temperature profiles until about 1 cm into the skin, which is approximately the average end of the adipose tissue layer. A computational domain of 6 cm into the skin is much deeper than the 3 layers of skin, meaning that a lot of muscle is included in the domain. The internal tissue beyond the adipose tissue would be cooled in this case, but due to much greater perfusion, it can be assumed to be kept at the core arterial temperature. The 6 cm length is not realistic 28 but could be used since it affects the temperature profile minimally within the skin region. The 1.5 cm characteristic length is too short to adequately model the temperature profile the hypodermis in the 60-minute treatment time. Figure 3.9: Study of sensitivity of the Green’s function of the X11 case to the computational length of the model given nominal perfusion of wb = 0.0005 s-1 and a treatment time of 60 minutes. The thermal conductivity is set to be k=0.52 Wm-1 K-1 . The Zeltiq paddle temperature of -14 ◦ C is used [41]. The grey portion of the graph highlights the adipose tissue layer ending at 1.5 cm into the skin. The lengths of 3 cm, 4.5 cm, and 6 cm show a similar temperature profiles up until the end of the adipose tissue layer. The temperature profile for the Green’s function of the X11 case is insensitive to model length after 3cm. 29 Figure 3.10: Study of sensitivity of the Green’s function of the X11 case to the computational length of the model given nominal perfusion of wb = 0.0005 s-1 and a treatment time of 60 minutes. The thermal conductivity is set to be k=0.52 Wm-1 K-1 . The Lipocryo paddle temperature of 3.1 ◦ C is used [45]. The grey portion of the graph highlights the adipose tissue layer ending at 1.5 cm into the skin. The lengths of 3 cm, 4.5 cm, and 6 cm show a similar temperature profiles up until the end of the adipose tissue layer. The temperature profile for the Green’s function of the X11 case is insensitive to model length after 3cm. Blood Perfusion Sensitivity The temperature solution is highly sensitive to variation in blood perfusion regardless of model solution. In Figs. 3.11 and 3.12, the X11 Green’s function is used. During the semi-transient portion of the temperature profile for either paddle temperature, the blood perfusion difference has a small effect on the profile. Once the near steady state is reached, starting at 15 minutes until the 60 minutes in the plot, the difference in the temperature profile for each solution is significant. One possible source of perfusion reduction could be the difference in paddle temperatures. A colder paddle may decrease perfusion due to the vasoconstriction. There is significant uncertainty in reported perfusion values [43]. The two different paddle types (as seen in Fig. 30 1.2) could be another reason for varying perfusion. The vacuum assisted paddle is thought to reduce perfusion in the skin by increasing the perceived pressure. The less perfusion the tissue has, the quicker it cools in this model with an almost 10 ◦ C difference at the maximum point of variance. Figure 3.11: Sensitivity of the Green’s function with X11 boundaries to perfusion. Three values of perfusion were used of either half or double the nominal wb = 0.0005 s-1 . The Zeltiq paddle temperature of -14 ◦ C is used [41]. The treatment time varies from 1 minute to 60 minutes. The length is 3 cm and all other properties are uniform with the thermal conductivity being set to k=0.52 Wm-1 K-1 . The transient temperature profile is insensitive to the blood perfusion, but the steady state temperature profile has a larger difference between the upper and lower blood perfusion limits. 31 Figure 3.12: Sensitivity of the Green’s function with X11 boundaries to perfusion. Three values of perfusion were used of either half or double the nominal wb = 0.0005 s-1 . The Lipocryo paddle temperature of 3.1 ◦ C is used [45]. The treatment time varies from 1 minute to 60 minutes. The length is 3 cm and all other properties are uniform with the thermal conductivity being set to k=0.52 Wm-1 K-1 . The transient temperature profile is insensitive to the blood perfusion, but the steady state temperature profile has a larger difference between the upper and lower blood perfusion limits. If the target adipose tissue temperature is 10.4 ◦ C [48], a -14 ◦ C paddle temperature with reduced perfusion, possibly from the vacuum assisted paddle, results in almost all of the adipose tissue reaching this target temperature. Even with reduced perfusion, the 3.1 ◦ C paddle temperature results in an adipose tissue temperature of about 20 ◦ C at the farthest point of the hypodermis. Tissue Thermal Conductivity Sensitivity The variation in comparing the tissue ther- mal conductivity is over a small range near 0.52 Wm-1 K-1 . Thermal conductivity is largely a function of tissue hydration [13]. For either paddle temperature (Fig. 3.13 and Fig. 3.14), the thermal conductivity has a minimal effect on the temperature profiles. The difference in the temperature profile between thermal conductivity values also does not drastically change 32 regardless of treatment time. The solution is insensitive in changes in conductivity during the transient or steady state portion of the treatment. Figure 3.13: The sensitivity of the Green’s function of the X11 case to the thermal conductivity. Values ±0.1 Wm-1 K-1 to a nominal value of 0.5 Wm-1 K-1 were used to assess the sensitivity at treatment times of 1 minute, 5 minutes, and 60 minutes. The Zeltiq paddle temperature of -14 ◦ C is used [41]. The computational length is 3 cm and the nominal blood perfusion value is wb = 0.0005 s-1 . All other properties are uniform. The temperature profile shows little sensitivity to thermal conductivity with the variation between values being small and constant throughout the transient and steady-state profiles. 33 Figure 3.14: The sensitivity of the Green’s function of the X11 case to the thermal conductivity. Values ±0.1 Wm-1 K-1 to a nominal value of 0.5 Wm-1 K-1 were used to assess the sensitivity at treatment times of 1 minute, 5 minutes, and 60 minutes. The Lipocryo paddle temperature of 3.1 ◦ C is used [45]. The computational length is 3 cm and the nominal blood perfusion value is wb = 0.0005 s-1 . All other properties are uniform. The temperature profile shows insensitivity to thermal conductivity with the variation between values being small and constant throughout the transient and steady-state profiles. Paddle Temperature Sensitivity The temperature profile is sensitive to the paddle temperature at the surface of the skin where contact with the paddle is made (see Fig. 3.15). Setting the paddle temperature hotter than -14 ◦ C results in a curve that does not reach the target temperature of 10.4 ◦ C for the adipose tissue, even at a 60-minute treatment time. Colder than -14 ◦ C results in a cool enough temperature profile to have all of the adipose tissue at or less than 10.4 ◦ C in order to try and meet the conditions for possible crystallization of the adipose tissue in the subcutaneous layers. Necrosis in the dermis can occur for surface temperatures of -20 ◦ C, regardless of treatment time [49]. The paddle temperatures of -10 ◦ C and -5 ◦ C [9] may be reasonable for CLL. The Zeltiq paddle of -14 ◦ C is more likely to cause adverse side effects and tissue damage since it is 34 much colder than the Lipocryo paddles at 3.1 ◦ C. The 3.1 ◦ C is relatively warm and does not cause much adipose tissue to reach the target temperature for crystallization. If the paddles used by Bernstein to explore cooling of the skin were used in a clinical CLL treatment, either might work better than Lipocryo and be safer than Zeltiq. Figure 3.15: Sensitivity study of the Green’s function solution for X11 boundary condition to the paddle temperature. The model uses the nominal perfusion value of wb = 0.0005 s-1 and a treatment time of 60 minutes. All other properties are constant with a characteristic length of 3 cm and a thermal conductivity of k=0.52 Wm-1 K-1 . T0 =-20 ◦ C is when tissue starts to develop necrotic behavior [49]. Bernstein explored paddle temperatures of -10 ◦ C and -5 ◦ C [9] while the two manufacturers, Zeltiq and Lypocryo, used paddle temperatures of -14 ◦ C [41] and 3.1 ◦ C [45], respectively. The temperature solution is sensitive to paddle temperature causing variance in the model curve up until the specified distant boundary condition. Heat Flux The heat fluxed provided by the manufacturers state values were presumably for the steady-state portion of the temperature profile during treatment. While the heat flux is initially transient, after about 5 minutes the steady state values are reached. The final portion of the steady state heat flux can be seen in Fig. 3.17 and 3.19. The temperature profile is not sensitive to the tissue thermal conductivity meaning that the sensitivity studies 35 Max Heat Flux Min Heat Flux Solution (mWcm-2 ) (mWcm-2 ) erf -4191 -70 -14 ◦C GFX10 -4193 -168 GFX11 -3053 -173 erf -2786 -46 3.1 ◦ C GFX10 -2787 -112 GFX11 -2029 -115 Table 3.1: Minimum and maximum heat flux during treatment. Values are negative for this model but can be assumed to be positive for direct comparison to the manufacturer heat flux. Maximum value is the largest transient value with the minimum heat flux is the smallest flux during the steady state portion of both the temperature profile and the heat flux. The maximum heat flux is not measured at time zero, but one second into the time interval due to the heat flux going to infinity at zero. done for the heat flux with varying k made the heat flux a function of only the difference in k. For this reason, all heat flux plots were only done for the nominal thermal conductivity value. For either paddle temperature the X11 and X10 cases for the Green’s function are almost indistinguishable from each other (see Fig. 3.16 and 3.18), except for the maximum heat flux during the transient portion of the profile. It can be seen from the time evolution plots of the GF models for either paddle temperature that the X10 case cools much quicker in Figures 3.3 and 3.4. In the direct comparison of the models in Fig. 3.7 for the Zeltiq paddle temperature that the GF profiles are very similar near the surface of the skin. The error function has a smaller minimum heat flux value, which corresponds to the model comparison plot. 36 Figure 3.16: Comparison of the surface heat flux between the three model solutions over a 60 minute treatment time interval. The Zeltiq paddle temperature of -14 ◦ C is used [41]. The perfusion is not included in the error function. All other parameters are kept constant. Figure 3.17: Isolation of 2500 seconds to 3600 seconds of -14 ◦ C paddle temperature heat flux comparison of error function and the Green’s function of the X10 case and X11 case. 37 Figure 3.18: Comparison of the surface heat flux between the three model solutions over a 60 minute treatment time interval. The Lipocryo paddle temperature of 3.1 ◦ C is used [45]. The perfusion is not included in the error function. All other parameters are kept constant. Figure 3.19: Isolation of 2500 seconds to 3600 seconds of 3.1 ◦ C paddle temperature heat flux comparison of error function and the Green’s function of the X10 case and X11 case. 38 CHAPTER 4 DISCUSSION AND CONCLUSIONS CLL is a popular cosmetic procedure used to reduce subcutaneous fat in localized areas (e.g., under the chin) through superficial cooling of the epidermis, which is thought to trigger apoptosis in subcutaneous adipocytes. While there is uncertainty as to the mechanism by which the adipocytes become apoptotic, Manstein suggests that the intracellular triglycerides of the adipocytes crystallize at 10.4 ◦ C, thus triggering apoptosis [48]. CLL is usually a safe and effective procedure with few side effects. A common side effect is transient neuropathy that typically resolves shortly after the procedure. Paradoxical adipose hyperplasia is a rare outcome that may require surgical intervention. Improved outcomes after CLL might be possible if the mechanism of adipocyte loss and the temperature distribution during treatment were known. One step to reducing the uncertainty in the treatment is to develop a mathematical model of the tissue temperature during treatment. This study uses a target temperature of 10.4 ◦ C to estimate the extent of effective treat- ment. Additionally, the adipose tissue in the subcutaneous layer is assumed to end at 1.5 cm into the skin (see table 2.1). The temperature of the skin is modeled assuming a single layer with uniform properties, even though the skin has three different layers with varying properties. This assumption appears reasonable because the unperfused epidermis is thin and has little effect on the temperature profile. The dermis and hypodermis likely have different blood perfusion rates and thermal conductivities. There is, however, considerable uncertainty in blood perfusion rate of tissues. Blood rates are often published in qualitative terms and not quantified. In addition, blood perfusion can be affected by apparent pressure on a tissue and by the temperature of the tissue. Thus, a value often quoted for skin is 39 used. Various values are also quoted in the literature for the thermal conductivity of tissue. Tissue thermal conductivity varies considerable with the mass fraction of water in the tissue [13]. Some of the largest values quoted are likely to include the effects of blood perfusion. Sensitivity calculations for perfusion and conductivity justify these assumptions, although quantifying the perfusion rate would aid in calculation of the steady temperature. The transient temperature profiles are based on the Pennes bioheat equation 2.1 solved using the error function and the Green’s functions for the X10 and X11 boundary condition cases. The error function solution omits perfusion and has a semi-infinite boundary condi- tion for the boundary distant from the skin (i.e., the internal temperature). This solution is best for the initial response of unperfused tissue and is used to verify the Green’s function solutions. The Green’s function solutions include perfusion, making them more representa- tive of the tissue temperature. The X10 case includes blood perfusion and, similar to the error function, the distant internal boundary condition is the semi-infinite condition. This results in a colder internal temperature as time progresses than the core body temperature. The X11 case is also perfused, but has a set internal boundary condition of the core body temperature (i.e., 37 ◦ C), making it the most anatomically correct. Realistically, the internal body temperature fluctuates and would result in a temperature profile between the X11 and X10 solutions, but the X11 case is still closer to an anatomical reaction than the X10 case. Each Green’s function solution shows good agreement with the error function solution after setting the blood perfusion term to zero. With no blood perfusion, the X10 Green’s function solution equals the error function solution identically. The X11 Green’s function solution, while not identical to the error function solution, agrees with the error function solution during early temperature responses. These results verify the Green’s function solutions. 40 All temperatures were calculated in non-dimensional terms, and are presented in dimen- sional terms to aid comparison in terms of the 10.4 ◦ C target temperature for the adipose tissue. Except where noted, such as for studying paddle temperature sensitivity, results are presented for the paddle temperature of two prominent manufacturers Zeltiq (-14 ◦ C) [48] and Lipocryo (3.1 ◦ C) [45]. Aside from the paddle temperature sensitivity, all of the sensi- tivities of the properties are independent of paddle temperature. For discussion purposes, unless the paddle temperature is specified, the assumption can be made that it applies to the paddles of both manufacturers. A model domain of 3 cm was used for the X11 Green’s function solution. Figure 3.9 justifies this length. For referenced, the subcutaneous layer is assumed to end at 1.5 cm into the skin (see Table 2.1). A model length of 1.5 cm does not represent the susceptibility of the adipocytes to cooling. Temperatures for model depths of 4.5 cm and 6 cm to set the boundary temperature of the core body temperature show good agreement in the adipose layer. These depths are unrealistic for anatomically. The 3 cm depth is more realistic anatomically and compares well with the deeper boundaries. Both manufacturers suggest treatment times of from 30 to 60 minutes. Due to its formu- lation, the error function solution (see Fig. 3.1 and 3.2) does not reach steady state. For a 60-minute treatment duration, the temperature at the deepest point in the solution do- main is less than 20 ◦ C for the Zeltiq paddle temperature and approximately 25 ◦ C for the Lipocryo. This is much colder than the core body temperature, so even though the adipose tissue is colder than the target of 10.4 ◦ C, the error function solution does not representative of an accurate physiological response. As seen in Figure 3.3 and 3.4, inclusion of blood perfusion in the X10 Green’s function solution leads to a warmer profile than for the error function solution, albeit with the same 41 distant boundary condition. The temperature profile varies slowly after 30 minutes and the 60-minute treatment time results in a temperature at a depth of 3 cm of just cooler than 30 ◦C for the Zeltiq paddle and approximately 30 ◦ C for the Lipocryo paddle. Still colder than the core body temperature. The X11 Green’s function solution represents a more clinically accurate response as the body seeks to maintain homeostasis. With the internal boundary (dimensional) temperature set to 37 ◦ C at the distance computation boundary, the temperature profile is the warmest of the solutions. This means that the 10.4 ◦ C target temperature is not met for much of the hypodermis. For the -14 ◦ C paddle temperature and 60 minutes of treatment, the deepest point of the adipose tissue reaches approximately 20 ◦ C. For the 3.1 ◦ C paddle, the deepest adipose tissue reaches approximately 25 ◦ C. Figures 3.11 and 3.12 show the sensitivity of the X11 Green’s function solution to variation in blood perfusion. Perfusion was varied about the nominal value of 0.0005 s-1 between 0.00025 s-1 and 0.001 s-1 . The transient temperature profiles are less affected by the blood perfusion variation than is the steady state portion. The difference in temperature between the profiles for the minimum and maximum perfusion values is near 10 ◦ C for the -14 ◦ C paddle temperature and just below 10 ◦ C for the 3.1 ◦ C paddle temperature. This is a significant difference in temperatures suggesting that the model is sensitive to the values of blood perfusion. This suggests that quantifying the rate of blood perfusion is important for improved clinical predictions. Smaller blood perfusion rates lead to colder temperature profiles. At the minimum perfusion value, the temperature at the end of the assumed adipose tissue layer is just about the target 10.4 ◦ C if the -14 ◦ C paddle is used. This also suggests that the vacuum assisted paddle, which may reduce blood perfusion, would produce a colder adipose tissue temperature. 42 The temperature profile varied little with varying thermal conductivity over the range of 0.5 ± 0.1 W m-1 K-1 . Thermal conductivity of tissues varies about a small range changing most significantly with the mass fraction of water [13]. The maximum differences in tem- perature for this range is small (see Figures 3.13 and 3.14 for the two paddle temperatures). Unlike for blood perfusion, this difference in temperature profiles was similar for the transient and steady state. Each of the two manufacturers referenced in this study lists heat flux for the removal of energy during CLL. Lipocryo suggests a heat flux value of 140 mWcm-2 [45], while Zeltiq offers several values, although in terms of a proprietary cooling index factor or CIF. Based on available values the heat flux ranges from -36.8 mWcm-2 to -72.9 mWcm-2 [17, 42] and the relationship between heat flux and CIF can be seen in Figure 1.3. These values are presumably stated for the steady state, as initial transient values are much greater in mag- nitude. The initial maximum and steady state minimum values for each solution can be seen in Table 3.1. Comparing the X11 Green’s function steady state heat flux values with those stated by the manufacturers shows reasonable agreement. For Zeltiq, the calculated heat flux during steady state was -173 mWcm-2 , while the manufacture defined values are lower. The calculated heat flux during steady state for the Lipocryo is calculated to be -115 mWcm-2 , which is closer to the 140 mWcm-2 defined by the manufacturer. Of course, the magnitude of the heat flux should be compared in this instance, rather than the sign, be- cause the manufacturers may define it in a positive sense for the heat gain in the paddle. The manufacturers may also be measuring heat flux as the cooling power of the device, rather than the heat loss of the skin, accounting for some of the differences. In order to provide the safest and most effective CLL treatment, the temperature and du- ration of cooling necessary to trigger apoptosis in the adipose tissue, the assumed mechanism of CLL, should be quantified. The skin does not need to be exposed to freezing tempera- 43 tures in order to be injured, so exploration of nonfreezing dermal injuries is also important for a safe CLL treatment. The paddle temperature significantly changes the outcomes of the temperature profile due to changing the starting point of the curve at the outer surface boundary condition. In addition, the manufacturers should provide more accurate measure- ments of heat flux and the CIF in order to allow for the clinician to provide a safer treatment to their patient. The CIF and heat flux measurements used by the common producers of machines for CLL [45, 41, 17] do provide enough insight into CLL for thermodynamics and heat transfer analyses of the mechanisms of CLL. Additionally, more study is required for the mechanisms of apoptosis in adipocytes. The time required for exposure of tissue to 10.4 ◦ C [48], assuming that it is the crystallization of triglycerides that causes apoptosis, needs to be considered. Since the temperature profiles demonstrate that all of the adipose tissue is unlikely to be at 10.4 ◦ C, the role of intercellular communication in triggering cell death should also be examined. Moreover, CLL could benefit if the patient-specific tissue thicknesses and tissue properties were included in a pre-procedure analysis of CLL. This suggests the continuing need for a defined relationship between time and temperature of the adipocytes in order to crystallize the cells and when a dangerous environment is reached. 44 BIBLIOGRAPHY 45 BIBLIOGRAPHY [1] Oktay Akkus, Aytekin Oguz, Mehmet Uzunlulu, and Muhammed Kizilgul, “Evaluation of Skin and Subcutaneous Adipose Tissue Thickness for Optimal Insulin Injection”, Journal of Diabetes and Metabolism, 3:8, 2012. [2] Zahra Alizadeh, Farzin Halabchi, Reza Mazaheri, Maryam Abolhasani, and Mastaneh Tabesh, “Review of the Mechanisms and Effects of Noninvasive Body Contouring De- vices on Cellulite and Subcutaneous Fat”, International Journal of Endocrinology and Metabolism, 14, 2016. [3] Naim Alkhouri, Agnieszka Gornicka, Michael P. Berk, Samjhana Thapaliya, Laura J. Dixon, Sangeeta Kashyap, Philip R. Schauero, and Ariel E. Feldstein, “Adipocyte Apop- tosis, a Link between Obesity, Insulin Resistance, and Hepatic Steatosis”, The Journal of Biological Chemistry, 285, 3428-3438, 2010. [4] Mathew M. Avram and Rosemary S. Harry, “Cryolipolysis for Subcutaneous Fat Layer Reduction”, Lasers in Surgery and Medicine, 41, 703-708, 2009. [5] Mark Baldry, Victoria Timchenko, Chris Menictas, “Thermal modelling of controlled scalp hypothermia using a thermoelectric T cooling cap”, Journal of Thermal Biology, textbf76:8-20, 2018, doi: 10.1016/j.jtherbio.2018.06.008. [6] S. M. Becker, “One-Dimensional Transient Heat Conduction in Composite Living Per- fuse Tissue”, textitJournal of Heat Transfer, 135, 071002, 2013, doi: 10.1115/1.4024063. [7] Cédric Benoit and Ali Modarressi, “Severe frostbite complication after cryolipolysis: A case report”, JPRAS Open, 25, 46-51, 2020. [8] Eric F. Bernstein, Jason D. Bloom, Lisa D. Basilavecchio, and Jessica M. Plugis, “Non- Invasive Fat Reduction of the Flanks Using a New Cryolipolysis Applicator and Over- lapping, Two-Cycle Treatments”, Laser Surg. Med. 46:731-735, 2014. [9] Eric F. Bernstein, “Long-term efficacy follow-up on two cryolipolysis case studies: 6 and 9 years post-treatment”, Journal of Cosmetic Dermatology 15, 561-564, 2016. [10] BioRender.com, Reprinted from “Anatomy of the Skin”, 2021, https://app.biorender.com/biorender-templates. [11] Atiyeh, Bishara S., Romeu Fadul, Jr. and Fadl Chahine, “Cryolipolysis (CLL) for Re- duction of Localized Subcutaneous Fat: Review of the Literature and an Evidence-Based Analysis”, Aesth Plast Surg, 44:2163-2172, 2020. [12] Charney, C. K., “Mathematical Models of Bioheat Transfer”, Advances in Heat Transfer, Elsivier, 22:19-155, 1992. 46 [13] John Chato, “Selected Thermophysical Properties of Biological Materials”, Heat Trans- fer in Medicine and Biology:Analysis and Applications, 1, Editors: Eberhart, R.C., Shitzer, A. (Eds.),Appendix 2, 413-423, 1985. [14] Sydney R. Coleman, Kulveen Sachdeva, Barbara M. Egbert, Jessica Preciado, and John Allison, “Clinical Efficacy of Noninvasive Cryolipolysis and Its Effects on Peripheral Nerves”, Aesth Plast Surg, 33:482–488, 2009. [15] Cole, K. D., Beck, J. V., Haji-Sheikh, A., Litkouhi, B., Heat Conduction Using Green’s Functions, CRC Press, 2011. [16] Zhong-Shan Deng and Jing Liu, “Analytical Study on Bioheat Transfer Problems with Spatial or Transient Heating on Skin Surface or Inside Biological Bodies”, Journal of Biomechanical Engineering, 124, 638-649, 2002. [17] Derrick, Chase D., Shridharani, Sachin M., and Broyles, Justin M., “The Safety and Ef- ficacy of Cryolipolysis: A Systematic Review of Available Literature”,Aesthetic Surgery Journal, 35: 830-836, 2015, doi:10.1093/asj/sjv039. [18] J. W. Durkee Jr and P P Antich and C E Lee, “Exact solutions to the multiregion time- dependent bioheat equation. I: Solution development”, Phys. Med. Biol., 35, 847-867, 1990. [19] J. W. Durkee Jr and P P Antich and C E Lee, “Exact solutions to the multiregion time-dependent bioheat equation. II: Numerical evaluation of the solutions”, Phys. Med. Biol., 35, 869-889, 1990. [20] Epstein, E.H. and Oren, M.E., “Popsicle Panniculitis”, The New England Journal of Medicine, 282: 966-967, 1970. [21] fellrnr.com“CoolSculpting and DIY CoolSculpting (Cryolipolysis)”, 2019, https://fellrnr.com/wiki/CoolSculpting. [22] Keys A, Fidanza F, Karvonen MJ, Kimuru N, Taylor HL., “Indices of relative weight and obesity”, J Chron Dis, 25, 329–343, 1972. [23] Lilit Garibyan, William H. Sipprell III, H. Ray Jalian, Fernanda H. Sakamoto, Mathew Avram, and R. Rox Anderson, “Three-Dimensional Volumetric Quantification of Fat Loss Following Cryolipolysis”, Lasers in Surgery and Medicine, 46:75–80, 2014. [24] Haxthausen, H., “Adiponecrosis E Frigore”, British Journal of Dermatology, 53: 83-89, 1941. [25] D. A. Hodson, J. C. Barbenel, and G. Eason, “Modelling transient heat transfer through the skin and a contact material”, Physics in Medicine & Biology, 34:10, 1493-1507, 1989. 47 [26] Michael J. Ingargiola, Saba Motakef, Michael T. Chung, Henry C. Vasconez, and Gordon H. Sasaki, “Cryolipolysis for Fat Reduction and Body Contouring: Safety and Efficacy of Current Treatment Paradigms”, American Society of Plastic Surgeons, 135:6, 1581- 1590, 2015. [27] H. Ray Jalian, Mathew M. Avram, Lilit Garibyan, Martin C. Mihm, and R. Rox An- derson, “Paradoxical Adipose Hyperplasia After Cryolipolysis”, JAMA Dermatol., 150: 317-319, 2014 doi:10.1001/jamadermatol.2013.8071. [28] Latif M. Jiji and Peter Ganatos, “Freezing of blood perfused tissue: An improved quasi- steady solution”, International Journal of Heat and Mass Transfer, 51, 5679-5686, 2008. [29] Purva Joshi and Yoed Rabin, “Thermal analysis of marginal conditions to facilitate cryopreservation by vitrification using a semi-empirical approach”, Cryobiology, 91, 128-136, 2019. [30] Junghyo Jo, Oksana Gavrilova, Stephanie Pack, William Jou, Shawn Mullen, Anne E. Sumner, Samuel W. Cushman, and Vipul Periwal, “Hypertrophy and/or Hyperplasia: Dynamics of Adipose Tissue Growth”, PLos Computational Biology, 5:3, 2009. [31] Junghyo Jo, Juen Guo, Teresa Liu, Shawn Mullen, Kevin D. Hall, Samuel W. Cushman, and Vipul Periwal, “Hypertrophy-Driven Adipocyte Death Overwhelms Recruitment under Prolonged Weight Gain”, Biophysical Journal, 99, 3535-3544, 2010. [32] Junghyo Jo, Zeina Shreif, Jonathan R. Gaillard, Matilde Arroyo, Samuel W. Cushman, and Vipul Periwal, “Mathematical Models of Adipose Tissue Dynamics”, Studies in Mechanobiology, Tissue Engineering and Biomaterials, 16, 11-34, 2015. [33] Terrence C Keaney and Lina I Naga, “Men at risk for paradoxical adipose hyperplasia after cryolipolysis”, Journal of Cosmetic Dermatology, 15: 575-577, 2016. [34] Kelly, M. E., Rodrı́guez-Feliz, J., Torres, C., Kelly, E., “Treatment of Paradoxical Adi- pose Hyperplasia following Cryolipolysis: A Single-Center Experience”, Plastic and Reconstructive Surgery, 142: 17e-22e, 2018, doi: 10.1097/PRS. 0000000000004523. [35] Sepideh Khoshnevis, Natalie K. Craik, R. Matthew Brothers, and Kenneth R. Diller, “Cryotherapy-Induced Persistent Vasoconstriction After Cutaneous Cooling: Hysteresis Between Skin Temperature and Blood Perfusion”, Journal of Biomechanical Engineer- ing, 138, 2016. [36] Kenneth B. Klein, MD, Brian Zelickson, MD, Jeffrey G. Riopelle, MD, Eric Okamoto, MD, Eric P. Bachelor, MD, Rosemary S. Harry, MSBME, and Jessica A. Preciado, PhD, “Non-Invasive Cryolipolysis for Subcutaneous Fat Reduction Does Not Affect Serum Lipid Levels or Liver Function Tests”, Lasers in Surgery and Medicine, 41:785-790, 2009. 48 [37] Paulus A. Kroon and Monty Kriegerg, “The Mobility of Cholesteryl Esters in Native and Reconstitute Low Density Lipoprotein as Monitored by Nuclear Magnetic Resonance Spectroscopy”, The Journal of Biological Chemistry, 256:11, 5340-5344, 1981. [38] Nils Krueger, Sophia V. Mai, Stefanie Luebberding, and Neil S. Sadick, “Cryolipolysis for noninvasive body contouring: clinical efficacy and patient satisfaction”, Clinical, Cosmetic and Investigational Dermatology, 7, 201-205, 2014. [39] Leonard, C. D., Kahn, S. A., and Summitt, J.B., “Full-thickness wounds resulting from ‘do-it-yourself’ cryolipolysis: a case study”, Journal of Wound Care, 25: S30-S33, 2016, doi: 10.12968/jowc.2016.25.Sup4.S30. PMID: 27068348. [40] Suvaddhana Loap and Richard Lathe, “Mechanism Underlying Tissue Cryotherapy to Combat Obesity/Overweight: Triggering Thermogenesis”, Journal of Obesity, 2018, https://doi.org/10.1155/2018/5789647. [41] Dieter Manstein, Hans Laubach, Kanna Watanabe, William Farinelli, David Zu- rakowski, and R. Rox Anderson, “Selective Cryolysis: A Novel Method of Non-Invasive Fat Removal”, Lasers in Surgery and Medicine, 40, 595-604, 2008. [42] Mark Melkerson, ”510(K) Sumary of Safety and Effectiveness”, Letter to ZELTIQ Aes- thetics Inc, FDA, May 2, 2012. [43] Mudaliar, A.V., Brent E. Ellis, Patricia L. Ricketts, Otto I. Lanz, Charles Y. Lee, Thomas E. Diller, and Elaine P. Scott, “Noninvasive Blood Perfusion Measurements of an Isolated Rat Liver and an Anesthetized Rat Kidney”, Journal of Biomedical Engineering, 130:6, 2008. [44] Mostafa, M.S.E.M. and Elshafey, M.A., “Cryolipolysis Versus Laser Lipolysis on Ado- lescent Abdominal Adiposity”, Lasers in Surgery and Medicine, 48:365–370, 2016. [45] Hernán R. Pinto, Eduardo Garcia-Cruz, and Graciela E. Melamed, “A Study to evaluate the action of lipocryolysis”, CryoLetters, 33:176-180, 2012. [46] Ricardo Romero-Méndez, Francisco G. Pérez-Gutiérrez, Joseph J. Musacchia, and Wal- fre Franco, “Passive cooling of cutaneous and subcutaneous tissues using phase changing materials: feasibility study using a numerical model”, International Journal of Hyper- thermia, textbf34:4, 363-372, 2018. [47] Y. Rabin and A. Shitzer, “Combined Solution of the Inverse Stefan Problem for Suc- cessive Freezing/Thawing in Nonideal Biological Tissues”, Journal of Biomechanical Engineering, 119, 146-152, 1997. [48] Amir Y. Sajjadi, Dieter Manstein, and Stefan A. Carp, “Measuring Temperature In- duced Phase Change Kinetics in Subcutaneous Adipose Tissues Using Near Infrared Spectroscopy, MR Imaging and Spectroscopy and OCT”, Scientific Reports, 7, 177-186, 2017. 49 [49] Andrew A. Gage, Joseph A. Caruana, Jr., and Mario Montes, “Critical Temperature for Skin Necrosis in Experimental Cryosurgery”, Criobiology, 19, 273-282, 1982. [50] H. Herwig and K. Beckert, “Fourier Versus Non-Fourier Heat Conduction in Materials With a Nonhomogeneous Inner Structure”, J Heat Transfer,122:363–365, 2000. [51] Elaine P. Scott, Muluken Tilahun, and Brian Vick, “The Question of Thermal Waves in Heterogeneous and Biological Materials”, J. Biomechanical Engineering, 131: 074518– 12009, 2009. [52] Singh, S.M., Geddes, E. R.C., Boutrous, S. G., Galiano, R.D., and Friedman, P.M., “Paradoxical Adipose Hyperplasia Secondary to Cryolipolysis: An Underreported En- tity?”, Lasers in Surgery and Medicine 47:476-478, 2015. [53] Malgorzata A. Jankowska and Grazyna Sypniewska-Kaminska, “An Interval Finite Dif- ference Method for the Bioheat Transfer Problem Described by the Pennes Equation with Uncertain Parameters”, Mechanics and Control, 31:2, 77-84, 2012. [54] Mitra, K., Kumar, S., Vedavarz, A., and Moallemi, M. K., “Experimental Evidence of Hyperbolic Heat Conduction in Processed Meat”, J. Heat Transfer, 117: 568–573, 1995. [55] Kaminski, W., “Hyperbolic Heat Conduction Equation for Materials With a Nonhomo- geneous Inner Structure”, J. Heat Transfer, 112: 555–560, 1990. [56] Paul Störchle, Wolfram Müller, Marietta Sengeis, Sonja Lackner, Sandra Holasek, and Alfred Fürhapter-Rieger, “Measurement of mean subcutaneous fat thickness: eight stan- dardised ultrasound sites compared to 216 randomly selected sites”, Scientific Reports, 8, 2018. [57] Pennes, H. H., “Analysis of tissue and arterial blood temperature in the resting human forearm”, J. Applied Physiology, 85: 93–122, 1948. [58] Seaman, S. A., Tannan, S. C., Cao, Y., Peirce, S. M. Peirce, and Gampper, T. J., “Paradoxical Adipose Hyperplasia and Cellular Effects After Cryolipolysis: A Case Report”, Aesthetic Surgery Journal, Vol 36(1) NP6-NP13, 2016. [59] Reeta Vyas and M. L. Rustgi, “Green’s function solution to the tissue bioheat equation”, Medical Physics, 19, 1319-1324, 1992, doi: 10.1115/1.4024063. [60] Wissler, E.H., “Pennes 1948 paper revisited”, J. Applied Physiology, 85: 35–41, 1998. [61] Hui Xin, Bin Chen, Zhifu Zhou, Dong Li, and Jiameng Tian, “Numerical investigation of multi-pulsed cryogen spray cooling for skin cold protection in laser lipolysis”, Numerical Heat Transfer, Part A: Applications, 77:7, 730-742, 2020. [62] Lisa X. Xu, “New Developments in Bioheat and Mass Transfer Modeling”, Annual Review of Heat Transfer, 10, Begell House, Chang-Lin Tien, 1-23, 1999. 50 [63] Christopher J. Yeung and Ergin Atalar, “A Green’s function approach to local rf heating in interventional MRI”, Med. Phys., 28, 826-832, 2001, doi: 10.1118/1.1367860. [64] Brian D. Zelickson, A. Jay Burns, and Suzanne L. Kilmer, “Cryolipolysis for Safe and Effective Inner Thigh Fat Reduction”, Lasers in Surgery and Medicine, 47, 120-127, 2015. [65] Le Zhang, Weizhong Dai, and Raja Nassar, “A Numerical Method for Optimizing Laser Power in the Irradiation of a 3-D Triple-Layered Cylindrical Skin Structurel”, Numerical Heat Transfer, Part A, 48, 21-41, 2005. [66] Le Zhang, Weizhong Dai, and Raja Nassar, “A Numerical Method for Obtaining an Optimal Temperature Distribution in a 3-D Triple-Layered Cylindrical Skin Structure Embedded with a Blood Vessel”, Numerical Heat Transfer, Part A, 49, 765-784, 2006. 51