Depth Penetration of Light into Skin as a Function of Wavelength from 200 to 1000 nm

An increase in the use of light‐based technology and medical devices has created a demand for informative and accessible data showing the depth that light penetrates into skin and how this varies with wavelength. These data would be particularly beneficial in many areas of medical research and would support the use and development of disease‐targeted light‐based therapies for specific skin diseases, based on increased understanding of wavelength‐dependency of cutaneous penetration effects. We have used Monte Carlo radiative transport (MCRT) to simulate light propagation through a multi‐layered skin model for the wavelength range of 200–1000 nm. We further adapted the simulation to compare the effect of direct and diffuse light sources, varying incident angles and stratum corneum thickness. The lateral spread of light in skin was also investigated. As anticipated, we found that the penetration depth of light into skin varies with wavelength in accordance with the optical properties of skin. Penetration depth of ultraviolet radiation was also increased when the stratum corneum was thinner. These observations enhance understanding of the wavelength‐dependency and characteristics of light penetration of skin, which has potential for clinical impact regarding optimizing light‐based diagnostic and therapeutic approaches for skin disease.


INTRODUCTION
With the markedly increased incidence of skin cancer (1,2), and the increasing demand for noninvasive, light-based therapies and technologies (3)(4)(5)(6), it is more important than ever to understand the intricate interactions that occur between light and skin. A detailed knowledge of the penetration depth that light can reach within skin is fundamental to this understanding.
For example, one of the challenges faced by the clinicians who prescribe UV-phototherapies, psoralen-UVA (PUVA) photochemotherapy, photodynamic therapy (PDT) or laser therapy is that of gaining an accurate quantification of the most effective dose. This challenge stems from the difficulty in gaining a reliable picture of the fluence rate distribution that the phototherapeutic source will produce within an individual's skin as well as an accurate idea of where the light will be absorbed (7)(8)(9). A simple and accessible method to gain this information in order to produce this picture would therefore be of great benefit. On the other hand, measuring the level of protection skin has against ultraviolet (UV) radiation is also vital in determining the risk from devices such as sunbeds and germicidal UV lamps (10,11).
We use Monte Carlo radiative transfer (MCRT) to probabilistically predict the path of light through a 6-layer skin model of Fitzpatrick skin type I. Such methods have been applied previously for many different applications (6,12,13). By compiling a full set of optical properties for six separate skin layers and running MCRT simulations for the wavelength range 200-1000 nm, a detailed data set of fluence rate distribution within skin has been produced.
As light travels through the skin, it is absorbed or scattered, resulting in a gradual reduction of energy density with depth. Due to the varying optical properties of skin layers, this reduction is not linear. It is the case however, that much of the data published in the literature represents the penetration depth into skin portrayed as solid lines ending suddenly at some depth, often with no reference to the intensity of the decided end point. Many diagrammatic illustrations represent the penetration of light by plotting only a single line for full spectral categories such as UV, visible and infra-red light (5,14). In reality, there are large differences in the depths that wavelengths within these spectral groups can reach, rendering these illustrations uninformative and potentially misleading.
Using the data set we have produced, plots detailing the depth at which the fluence rate into skin reaches a certain percentage of the incident value could be produced on a wavelength-bywavelength basis. This allows more detailed and informative images of light penetration into skin to be constructed.
To test the effect of different setups and demonstrate the adaptability of the simulations, the model was initially run for both direct and diffuse incident light to examine the difference this might make to penetration depth. The effect of changing the incident angle of the light was also investigated. Simulations were then run for several different thicknesses of stratum corneum, the protective outermost layer of the skin, to highlight the difference that the small changes within this layer can make to penetration depth of radiation within the UV spectrum. Finally, the lateral spread of light within the model was investigated by setting the incident light to be a circular beam of 1 mm in diameter.

MATERIALS AND METHODS
Monte Carlo radiative transport. Monte Carlo Radiative Transport (MCRT) is used to produce simulations of light transport within human skin. The simulations are written in FORTRAN adapted from code developed by McMillan et al. (15) which was adapted from code developed for astronomy (16). The simulation starts by defining a 3D grid of voxels with total grid size is 1 mm × 1 mm × 5 mm and 100 × 100 × 800 voxels in each direction. The grid contains 5 different voxel volumes in order to provide higher resolution data in the thinner upper layers of the skin.
Simulations are undertaken with a spectral resolution of 1 nm over the wavelength range 200-1000 nm. For each simulation, a total of 1 billion Monte Carlo packets, or power packets, are sent in one at a time, evenly distributed through the top of the grid. Fresnel reflections make the proportion of light entering the grid vary depending on the chosen angle of incidence and the refractive index difference at the grid-air boundary. The energy deposited by each power packet is recorded by storing an energy value within every voxel that it enters, proportional to the distance the power packet has traveled through the voxel. The power packet undergoes scattering as it propagates through the grid before either escaping or becoming absorbed the packet is terminated, and the process starts again. The optical properties of each voxel determine the overall proportion of scattering to absorption events that occur by being used to weight probability distribution functions which are then randomly sampled to make decisions. This also helps to decide the average distance the power packet travels between events as well as the new direction it takes after being scattered. Repeating boundary conditions are implemented such that when a power packet exits the grid in either the x or y direction, it is returned at the opposite side with other properties retained. However, if it exits though either z plane (top or bottom of the grid) the power packet escapes. This helps to simulate a semi-infinite slab of skin, making the 1 mm × 1 mm × 5 mm grid effectively infinite in the x and y direction. The simulation output is a 3D array of fluence rate values which can then be used to plot the penetration depth via post processing.
Skin model. The 6-layer skin model is based on the one developed by Barnard et al. with the addition of a subcutaneous fat layer (10).
The layers are formed within the grid by assigning each voxel a set of wavelength-dependent optical properties based on the layer it is located in. Figure 1 shows the model layers and their corresponding thicknesses. The dimensions of the layers were chosen with reference to a combination of Iglesias-Guitian et al. (17) and Barnard et al. (10).
The Stratum Corneum (layer 1) is thin, flat and makes up the outermost layer of the skin.
Due to its protective function, penetration of light through the stratum corneum is a topic of notable interest and so it is given a z-resolution of 200 voxels with a thickness of 0.1 μm to allow higher resolution data to be obtained within it. The epidermis (layer 2) is thicker than the stratum corneum and most of the layer is contained within a z-depth of 100 voxels of thickness 0.53 μm. The bottom of the epidermis has an undulating structure due to the presence of dermal papillae and this is simulated within the grid using Eq. 1, taken from Barnard et al. (10). Here, z is the depth in mm of a voxel at each x and y position, and depth is the average depth of the layer from the top of the grid, given in Table 1.
The melanin and basal layers (layers 3 and 4) are thin layers at the bottom of the epidermis that follow the shape of the dermal papillae. Another 200 voxels with thickness 0.32 μm are assigned to these two layers, along with the parts of the epidermis and dermis that the basal and melanin layers protrude into. It should be noted that melanin itself is not a skin layer, but a pigment produced by melanocytes and distributed within the epidermis. For Fitzpatrick skin phototypes I and II, most of the melanin present in the skin is concentrated in the basal layer of the epidermis (10,18). This helps to protect the DNA present within the basal layer from damage caused by dangerous optical radiation. For simplicity, the melanin layer described in this model is chosen to be a thin layer of melanized epidermis located just above the basal layer. The dermis and subcutaneous fat layers (layers 5 and 6) are significantly thicker than the top four layers. However, they are only assigned a z-depth of 150 voxels each, with 13 μm voxel thickness for layer 5 and 20 μm voxel thickness for layer 6, as there is less need for high resolution data here due to the larger change in penetration depth per wavelength.
Optical properties. The depth and direction that light travels through skin is determined by optical properties that are wavelength-dependent and vary depending on the biological content at each position. For simplicity, the skin model optical properties are chosen to vary between skin layers but remain constant within each layer. These properties consist of absorption and scattering coefficients (μ a and μ s ), anisotropy factor (g) and the refractive index (n). The refractive index does not vary significantly between layers (17,19) and so, for simplicity, was chosen to be a constant value of 1.38 (20,21) for all wavelengths and all skin layers. As a result, the layers are distinguished only by their scattering and absorption coefficients. The anisotropy factor is described as the average of cos(θ), where theta is the angle between the direction of initial travel of a photon and the direction after scattering. The Henyey Greenstein phase function (Eq. 2) is used together with g to describe the probability distribution of the scattering direction within biological tissue where forward scattering is dominant (22,23).  Table 1. A melanin volume fraction of 2% is used to simulate Fitzpatrick phototype I. The top of the dermis is covered with dermal papillae which cause the shape of the melanin layer, basal layer and bottom of the epidermis to follow the "egg-box" undulating shape described by Eq. 1.
For every layer, g is wavelength-dependent and modeled using Eq. 3, presented by van Gemert et al. where λ is the wavelength of light in nm (22).
Stratum corneum. The absorption and scattering coefficients selected for the stratum corneum (layer 1) were taken from those presented by van Gemert et al. (22). These were obtained using reflectance and transmission measurements through a slab of skin, measured by Everett et al. (24).van Gemert reports the optical properties for the stratum corneum in the wavelength range 200-400 nm (22). In order to extend this range to 1000 nm, reflection and transmission data were extracted up to 700 nm from relevant plots from Everett et al. (24). At 700 nm, both the transmission and reflection data were flattening off and so it was assumed that they remained roughly constant from 700 to 1000 nm where only the absorption coefficient of water, of which there is very little present in the stratum corneum, changes significantly (19). Following the method used by van Gemert, μ a and μ s were then calculated for these wavelengths.
Epidermis. Scattering coefficients for the epidermis were taken from the data presented by van Gemert (22) from 200 to 800 nm and extending it to 1000 nm using the equations from van Gemert, the data given by Wan et al. (25) and the assumption that the reflection and transmission values beyond this point are constant. The absorption coefficients for the epidermis were also taken from van Gemert for the range 200-450 nm while the coefficients presented by Meglinski et al. (19) were used for the range 451-1000 nm. As the 6-layer skin model used assumes that melanin is only present in the melanin layer, the absorption coefficient of melanin calculated using 4 was subtracted from the epidermis absorption for layer 2. μ a;melanin λ ð Þ ¼ 6:6 Â 10 11 λ À3:33 F melanin (4) where λ is the wavelength in nm and F melanin is the fraction of melanin present in the skin which varies with skin type. Melanin layer. As mentioned above, the 6-layer model assumes that the melanin in the epidermis is concentrated into a single layer near the bottom of the epidermis. The model used was chosen to be unexposed Fitzpatrick skin phototype I, which has the lowest melanin content and so will help determine the worst-case scenario when looking at UV penetration. This is estimated to be equivalent to a melanin volume fraction of 2% (10,17,26). However, as the melanin from the entire epidermis is concentrated into a thin layer, the fraction used in the calculation for the absorption coefficient for this layer had to be adapted to represent this. The adapted fraction was calculated using Eq. 5 This adapted fraction was then used to find the melanin absorption using Eq. 4, replacing F melanin . This melanin absorption was then combined with the absorption coefficient for the epidermis in layer 2 to compute the final absorption coefficient for the melanin layer. The scattering coefficients for the melanin layer were assumed to be equal to that of the epidermis in layer 2.
Basal layer. The basal layer forms the bottom of the epidermal layers, and so its optical properties are identical to that of the melanin-less epidermis in layer 2.
Dermis. From 200 to 449 nm, the dermis absorption coefficients were taken from the data presented by van Gemert (22) and from 450 to 1000 nm, the coefficients presented by Meglinski et al. (19) were chosen.
The scattering coefficients for the dermis were calculated using Eq. 6 from Jacques (27).
Subcutaneous fat layer. The absorption coefficient for the subcutaneous fat layer was taken from the data presented by Meglinkski et al. (19).
Comparison tests. In order to demonstrate the power and adaptability of the simulations, a range of different setups were run, and the results compared.
Direct and diffuse incident light. The simulation was run using the setup described above for the full range of 200-1000 nm. First, direct radiation, evenly distributed across the top of the grid, was simulated by choosing the initial power packet direction of entry to be normal to the surface of the voxel grid. This was then repeated using a diffuse light source, which was simulated by randomly choosing the incident angle of each power packet individually from a uniform distribution of all possible angles incident onto the top of the grid.
Stratum corneum thickness. The simulation was run for a wavelength range of 200-400 nm at various thicknesses of the stratum corneum, the upper most layer of the skin, in order to investigate its effect on the penetration depth of light within the UV range. It can be seen from Figs. 2 and 4 that for wavelengths beyond the UV range, the stratum corneum becomes optically thin. It therefore made sense to only investigate wavelengths within the UV range. Simulations with various stratum corneum thicknesses were run using evenly distributed, direct radiation. Thicknesses of 30, 20, 10 and 14.8 μm were chosen with reference to the range of stratum corneum thicknesses commonly in the face, forearm, shoulder and buttock, respectively (29,30). A thickness of 0 μm was also chosen to investigate the effect of a compromised or nonexistent stratum corneum found in a range of skin conditions (31). Areas of typically higher sun exposure such as the cheeks, forehead and nose were shown by Chirikhina et al. to have larger average stratum corneum thicknesses at 27.6, 24.6 and 23.1 μm, respectively, which would result in lower penetration depth of UV radiation in these areas (30).
Incident angle. The effect of the angle of incidence on radiation penetration depth was investigated by running the simulations first with an angle of incidence of 45°from the normal and again with an angle of 65 degrees from the normal. These simulations were run at wavelengths of 222, 300, 350, 400, 630 and 800 nm which are specific wavelengths of interest in either UV radiation safety studies or light-based therapies.
Lateral spread. The lateral spread of incident light through the grid was investigated by illuminating the center of the top of the grid with a 0.1 mm diameter beam of uniform intensity. As for the incident angle investigation, the simulations were run for wavelengths of 222, 300, 350, 400, 630 and 800 nm.

RESULTS
The absorption and scattering coefficients selected for each skin layer are displayed in Fig. 2.

Comparison tests
Direct and diffuse incident light. Direct and diffuse incident light for both direct and diffuse incident light, the depth that light can reach at each wavelength is displayed in Fig. 3. A list of fluence rates within a central column of voxels spanning the full depth of the grid was extracted and used to plot the data. The plots display the depths at which the fluence rate reaches 90%, 50%, 10%, 1% and 0.1% of the incident light value for each wavelength. As expected, the penetration depth increases with increasing wavelength until water absorption starts to take effect in the infrared (IR) region where it then decreases. It can be seen that for all wavelengths, a larger percentage of the direct incident light can reach further into the skin than the diffuse incident light.
Stratum corneum thickness. Stratum corneum thickness for each stratum corneum thickness simulated, the results are compared in Fig. 4 by plotting the depth at which the fluence rate reaches 90%, 10% and 1% of the incident value at each wavelength. As expected, simulations with a thinner stratum corneum show deeper penetration of light at all wavelengths from 200 to 400 nm. There is reduced depth of radiation travel through the highly absorbing stratum corneum, particularly at shorter wavelengths.
Incident angle. Figure 5 compares the fluence rate down a central column of the grid as a function of depth for all three incident angles at each wavelength. The fluence rate is normalized to the incident irradiance. It can be seen that incident angles further away from the normal, result in a larger fluence rate at the surface but lower fluence rate further down. Fluence rate values higher than 1 are seen near the top of the grid due to scattering effects. A higher fluence rate at the top of the skin for incident angles further away from the normal is possibly due to the forward bias of scattering in biological tissue (22,23). As light enters the skin at more oblique angle, it has a higher chance of being scattered in the x and y directions and exiting through the  Fig. 1, as a function of wavelength. Data for the stratum corneum are taken either directly from van Gemert et al. or derived using the same methods using data from Everett (22,24). Epidermis and melanin layer absorption is taken part from van Gemert et al. (22) and part from Meglinksi et al. (19) with melanin absorption calculated using Eq. 4 subtracted from the epidermis. Scattering data for the epidermal layers is taken either directly from van Gemert et al. or derived using the same methods using data from Everett (22,24). Dermis absorption is taken part from van Gemert et al. (22) and part from Meglinksi et al. (19) and dermis scattering is calculated using Eq. 6 from Jacques (27). Subcutaneous fat absorption is taken directly from Meglinski et al. (19) and scattering is calculated using eqs. 7 and 8 from Bashkatov et al. (28). top of the grid rather than traveling straight down in the z direction as light entering along the normal would. This causes an increase in the number of photons traveling through voxels in the upper layers of the skin, resulting in the higher fluence rate. For the same reason, an increased fluence rate is also observed in the upper layers for the diffuse incident light results compared with the direct incident light.
For each plot in Fig. 5, each dataset intersects once near the top of the grid, within the stratum corneum. The depths of these intersection points are displayed in Table 2.
Lateral spread. Contour plots of a central slice of the fluence rate grid are produced for each wavelength. Contours are plotted for the points where the fluence rate reaches 90%, 50% 10%, 1% and 0.1% of the incident value (Fig. 6). The spread of light in the lateral direction from a direct beam occurs due to scattering within the skin. It can be seen from Fig. 6 that longer wavelengths have different lateral spread and the different scattering coefficients for each layer result in different rates of lateral spreading, which dictate the unique contour shape. We can also see from Fig. 6e and f that the lateral spread of the 630 and 800 nm radiation goes beyond the 1 mm grid width.

DISCUSSION
With so much variability in the choices of optical properties to use within the Monte Carlo simulations, an emphasis was made on finding properties that were realistic when considering the biological contents of each layer. Figure 2 displays the final chosen optical properties for each layer. Several distinct peaks can be seen in the absorption plots. These peaks correspond to particular absorption features of the chemical constituents within the layer. The small peak in the stratum corneum and epidermis at around 280 nm is due to the keratin present in each of these layers (32). Blood present in the dermis and subcutaneous fat produces several peaks. The first is at around 425 nm and is named the Soret band or Soret peak. There are then two more smaller peaks between 500 and 600 nm which form the Q-band (23,32,33). There is also a peak at around 970 nm in the epidermis, dermis and subcutaneous fat, which is due to the increased absorption of water at this wavelength range (19,33).
Discontinuities appear at the transition points between the plotextracted and calculated data in the stratum corneum scattering and absorption at 400 nm, as well as the epidermal scattering at 800 nm. These jumps are likely mostly due to small differences in the datasets or anisotropy (g) values used in calculation, as well as uncertainties introduced when extracting data from plots produced by van Gemert et al, Everett and Wan et al. (22,24,25).
The results show large differences in light penetration depth with wavelength which correspond closely with the differences in optical properties. This highlights a drawback of this method; its reliability is fully dependent on the availability of accurate optical properties. These can be hard to measure for skin, especially in the shorter-UV wavelength region where the absorption is high (34). There is also a large variability in the geometry and contents of skin, and therefore optical properties both between and within individuals (29,33,35,36). This may limit comparisons between studies and while it is unclear how applicable the results are at an individual level, the data are more easily interpreted collectively, and broad comparisons can be made. For example,    The different optical properties in each layer change the rate that the fluence rate depletes as the light travels through, dictating the shape of the contour. Panels (e) and (f) show the contours being cut off due to the lateral spread of these wavelengths being wider than the grid width. It should be noted then that repeating boundary conditions were not in place in these simulations to allow the lateral spread of a single beam to be investigated.
in vivo measurements of the penetration depth of UV radiation into forearm skin by Meinhardt et al. show around 37% of the light at 300 nm reaches 15-30 μm and 40-60 μm at 340 nm (34). This is consistent with the data used to make Fig. 3 (37). This is consistent with the results in Fig. 3(b) at 450 nm and 650 nm which show 1% of the light reaching around 1.6 mm and a little over 5 mm, respectively. This shows that the results obtained by the simulation are well within the expected range and provide useful information across a broad wavelength region.
One advantage of this model is its adaptability. Provided there are optical properties available, which, it should be noted, there is a current lack of, the simulation can be easily adapted to represent many different skin types and layer thicknesses, as well as different incident light setups. As the accuracy and diversity of available optical properties increases, the reliability and applicability of this model will also increase proportionally. Specifically, the models adaptability will be useful when it comes to disease targeting. It has been shown that changing the standard treatment bandwidth for some phototherapy treatable skin diseases based on the diagnosis and disease stage and severity has led to more effective treatment (38,39). This is, in part, due to selecting phototherapy sources with greater penetration depths when needed. The model data will therefore provide useful information when choosing the most appropriate treatment for an individual circumstance or even when developing new light-based treatments.

CONCLUSION
To the authors' knowledge, this is the first presentation of the penetration depth of optical radiation in skin over the wavelength range 200-1000 nm. The impact on penetration depth of varying parameters, including optical radiation incident angle and stratum corneum thickness was also investigated. Researchers and healthcare professionals can use this information when evaluating new and existing techniques for diagnostic and therapeutic applications as it provides a useful guide relating to light penetration. Targeted approaches can be produced by determining the optical properties specific to an individual and to a disease process and incorporating them within the model. The penetration depth of light into skin increases with wavelength from the UV to the visible light range before decreasing again in the IR range in accordance with the selected optical properties. This depth is further increased as the thickness of the stratum corneum is decreased. Additionally, the depth reached by light from a diffuse source is slightly less than that of a direct source, although the fluence rate at the top of the grid from light with more oblique incident angles is higher.
Acknowledgements-Louise Finlayson acknowledges financial support from EPSRC Industrial Doctorate Centre Scheme (2262922) and the Laser Research and Therapy Fund (registered charity SC030850).

SUPPORTING INFORMATION
Additional supporting information may be found online in the Supporting Information section at the end of the article: Dataset S1. Fluence rate data from this paper can be accessed and combined with a spectrum of the users' choice via a web application at: https://loufin.shinyapps.io/Depth_Pen_App/. Table S1. Input parameters to use within the web application (at https://loufin.shinyapps.io/Depth_Pen_App/) to reproduce the data used to produce Fig. 3(a) and (b).