Investigation of non-uniformly emitting optical fiber diffusers on the light distribution in tissue

: The inﬂuence of inhomogeneities in the emission characteristics of optical ﬁber diﬀusers on the light distribution within biological tissue was evaluated by Monte Carlo (MC) simulations and by experiments on optical phantoms. Due to the strong scattering of light within biological tissue, inhomogeneities in the emission proﬁle become blurred within a short light propagation distance, so that the light distribution within the tissue approaches that of a homogeneous diﬀuser. The degree of feature vanishing in the light distribution is mainly determined by the width of the inhomogeneities. It was shown that the inﬂuence of local inhomogeneities on top of a homogeneous light distribution fades away very eﬀectively within 1 mm of tissue depth, which results in a light distribution very close to that for a homogeneously emitting diﬀuser. Emission proﬁles composed of multiple narrow peaks distributed over the full diﬀuser length with a peak-to-peak distance of less than 2 mm result in an almost homogeneous light distribution after approximately 1 mm of tissue depth. While this article is focused on the impact of diﬀuser inhomogeneities on the light distribution within the tissue, the importance of further investigations on the related thermal eﬀects is also discussed.


Introduction
Interstitial photodynamic therapy (PDT) is an evolving method for the treatment of malignant brain cancer [1][2][3][4]. Much research is currently focussing on the translation from the laboratories to clinical application [5]. The strength of this method is its minimally invasive surgical approach in combination with the selective targeting of individual tumor cells. The selectivity of the tumor treatment is achieved by the increased uptake of aminolevulinic acid (5-ALA) in high-grade brain tumor tissue and the subsequent transformation to protoporphyrin IX (PPIX) as part of the heme-biosynthesis pathway. By illuminating the affected region with light of a wavelength of λ = 635 nm, the present PPIX molecules become excited. The excited PPIX molecules transfer their energy to the present intracellular oxygen molecules, which triggers the conversion of oxygen into a reactive oxygen species. Reactive oxygen species may cause apoptosis and necrosis in the adjacent regions [6][7][8]. In order to treat all targeted malignant cells, a sufficiently intense and sufficiently homogenous light distribution within the targeted tumor volume is necessary.
The illumination of such tissue volumes within the human body could be achieved by means of optical fibers with attached cylindrical diffusely scattering fiber tips [9][10][11][12][13]. Depending on the volume to be treated, multiple fiber diffusers may have to be positioned simultaneously within the targeted treatment site. The inter-fiber distance is set to approximately 9 mm in order to avoid thermal side effects between neighbouring fibers, while still retaining a photodynamic response at the outer tumor boundary [14]. Regarding the light intensity limits, it had been determined that the minimal applied photon number per cm 3 to reach a PDT effect in cells is approx. 1.3 × 10 18 photons [15]. An upper limit, especially for the brain, is defined by the critical temperature limit of 42°C, where healthy tissue starts to degenerate [16,17].
Cylindrical diffuser tips can be distinguished technically into two groups: volume and surface scattering devices. Volume scatterers are commonly realized by attaching polymer tubes containing scattering particles (e.g. TiO 2 ) to the end of the light-delivering fiber [18], or by generating internal scattering centers inside the fiber core by laser processing [19][20][21][22][23]. Surface scatterers can be manufactured by roughening the fiber core, which had been proposed to be accomplished with different laser systems (CO 2 , Ultrafast, and Excimer) [24][25][26][27][28], mechanically or by HF-etching techniques [29].
The challenge, which all manufacturing techniques have in common, is to provide homogenous light emission along the complete diffuser length. This is difficult to achieve, and, thus, many diffusers with inhomogeneous emission patterns are on the market. For example, when examined with a laboratory setup, intensity peaks may be present at the distal and/or proximal diffuser end [27,30,31], or the spatial distribution of the scattering elements may lead to alternating intensity maxima and minima along the diffuser [19,25]. This leads to the general question, how intensity variations in the emission profile of fiber diffusers influence the light distribution within the tissue.
This study investigates the impact of intensity inhomogeneities of diffusers on the light distribution within the surrounding tissue. A model based on Monte Carlo ray tracing was developed to simulate the light emission from a cylindrical light diffuser. This model was calibrated by means of the intensity distribution of a diffuser measured on a real cylindrical light diffuser by means of an imaging and a non-imaging camera setup [31]. In addition, the real diffuser was positioned surface-parallel orientation at different depths within an artificial tissue phantom with known optical properties, and the light distribution transmitted to the surface of this phantom was measured and investigated for each diffuser depth. The measurements on the real diffuser were then compared to the respective simulated light distributions in the correspondingly simulated tissue phantom. By varying the modelling parameters of the inhomogeneous emission profiles, the impact of height and width of single intensity peaks, as well as the interplay between two adjacent intensity peaks were investigated regarding depth dependent intensity changes in the modelled surrounding environment.

Experimental Part
As real diffuser with inhomogeneous emission profile, a junk diffuser (based on a CD-403, LifePhotonic GmbH, Bonn, Germany) was selected. It will be termed reference diffuser in the following. The emission profile of the reference diffuser is characterized by two local intensity peaks at the distal end, while the radial intensity distribution of the remaining diffuser length is measured to be approximately homogeneous. The diffuser region consists of a 30 mm long polymer tube, filled with TiO 2 polymers acting as scattering centers, attached to the end of a light waveguide with a core diameter of 400 µm. For all subsequently presented measurements and simulations, a diffuser length of 30 mm was used. To reduce the forward transmission of light, a mirror is attached at the distal diffuser end. Thereby, the primary transmitted light is reflected back into the diffuser, which serves to improve the homogeneity of the emission profile.
As sketched in Fig. 1(a), the spatially resolved intensity distribution on the diffuser surface was determined by means of an imaging setup, while the angular resolved emission profile was measured using a non-imaging camera setup as shown in Fig. 1(b). In addition, the light distribution on the surface of an artificial tissue phantom in which the diffuser itself is placed at different depths, as depicted in Fig. 1(c), was measured using the imaging technique. Furthermore, the experiment in Fig. 1(c) was reconstructed and investigated using MC-simulations ( Fig. 1(d)).
For all experiments, light of a LED (M625D3, Thorlabs, Newton, NJ, USA) emitting at a wavelength of λ = 635 nm was coupled into the proximal end of the diffuser fiber. In this way, the whole NA of the fiber was used. The emission power measured at the diffuser end was set to 80 The surface emission profile measured with setup (a) was allocated to the modelled surface emitter. A cone-shaped emission was assumed from each spot of the diffuser surface, with full cone opening angle 2α and cone tilt angle β as adjustable model parameters, which were calibrated using the results obtained with setup (b). The resulting light distribution within virtual receiver planes at distances of 2 mm, 5 mm, 10 mm, and 15 mm from the diffuser surface were evaluated. mW measured by a calibrated integrating sphere (IS200-4, Thorlabs, Newton, NJ, USA). This power level was sufficient to obtain signals from the maximum investigated fiber depth of 15 mm within the tissue phantom.
Emission profiles from cylindrical diffusers and intensity distributions based on subsequent light penetration through tissue were determined experimentally using three setups: (a) The imaging setup ( Fig. 1(a)) for characterization of the emission profile of a cylindrical diffuser consisted of a CCD camera (DMK 41BU02, The Imaging Source Europe GmbH, Bremen, Germany) and a lens (f = 12 mm, TCL 1216 5MP, The Imaging Source Europe GmbH, Bremen, Germany). The distance between lens and the diffuser surface was 30 cm while the diffuser was aligned horizontally in the imaging plane. For the determination of the intensity profile along the longitudinal diffuser axis the exposure time of the camera had to be carefully adjusted to avoid oversaturated pixels. Then, an image was captured and the grey-scale values along the diffuser axis were extracted. These intensity values were derived by averaging over 10 adjacent pixels in the direction perpendicular to the diffuser axis. This procedure was used to determine the emission intensity profile in air ( Fig. 1(a)) and on the artificial tissue phantom tissue surface ( Fig. 1(c)).
(b) The non-imaging camera setup ( Fig. 1(b)) [32] for characterization of the emission profile of a cylindrical diffuser consisted of a CCD camera chip (ICX205AL, Sony, Tokyo, Japan, dimension: (4.5 × 6.0) mm 2 , 1280 × 960 pixel, pixel size (4.65 × 4.65) µm 2 ) that was positioned at a distance of 2 mm from the surface of the fiber diffuser. As the diffuser length exceeds the chip sensor's dimension, the sensor was continuously moved by means of a computer controlled stepping motor (LTS300/M, Thorlabs, Newton, NJ, USA) along the diffuser surface with a velocity of 0.1 mm/s, while recording a video-sequence at 15 frames per second. The principle of the raw data acquisition and evaluation is sketched in Fig. 2. Mimicking a line sensor perpendicular to the moving direction, only one pixel line per frame was evaluated. The analysed pixel line on the sensor is fixed (marked yellow in Fig. 2). Thereby, each frame, with frame number n, corresponds to a distinct position x n along the longitudinal diffuser axis. Only the intensity value of the brightest pixel within this pixel line (marked by an "X" in Fig. 2) was extracted, as intensity value I n to build up the emission profile I n (x n ) along the full diffuser length. By means of this non-imaging method, also forward or backward emitted light could be detected at adjacent pixel positions x n , which is not possible in the imaging setup.
(c) The investigation of the light distribution in tissue arising from the reference diffuser was performed with the imaging setup sketched in Fig. 1(c). The diffuser was submerged in an artificial phantom mimicking the optical properties of human brain tissue. The phantom was prepared using a standard recipe [33]. In brief, 434 ml of water were mixed with 0.6 Vol.-% Agar (3 ml). The solution was heated to 100°C, while continuously stirring. After cooling to 60°C, 12.5 Vol.-% of Lipovenös (Lipovenös MCT 20%, Fresenius Kabi Deutschland GmbH, Bad Homburg, Germany) (62.5 ml) as scatterer and 0.1 Vol.-% of black ink (0.5 ml) as absorber were added to the solution. The scattering and absorption coefficients are defined by the used concentrations of Lipovenös and ink, respectively. The optical properties of this specific artificial phantom were measured to have an absorption coefficient µ a = 0.09 mm −1 , a reduced scattering coefficient µ s = 2.2 mm −1 , and an anisotropy factor g = 0.95 [33]. After further cooling to 40°C, the liquid phantom was poured into a metal container as illustrated in Fig. 1(c). The container allows for positioning the test diffuser in defined phantom depth positions below the phantom surface by means of pre-drilled lateral horizontal boreholes. The boreholes were set at a depth of 2 mm, 5 mm, 10 mm, and 15 mm below the artificial phantom surface. Preliminary tests revealed that measuring the light distribution for tissue depth smaller than 2 mm was not feasible as rupturing of the artificial phantom surface by the diffuser was not avoidable in this case. After preparation, the artificial phantom was solidified by cooling to room temperature (21°C). The measurement of the light distribution on the surface of the artificial tissue phantom was performed five times for each depth position as follows: The diffuser was placed through one of the boreholes into the phantom. The phantom surface was imaged using the configuration in Fig. 1(a). The intensity distribution on the phantom surface in the direction parallel to the longitudinal diffuser axis was evaluated in analogy to the abovementioned imaging procedure applied to the diffuser surface in air. This procedure was repeated for diffuser placements in each borehole of the container (n = 5) independently and in random order.

MC-simulation part
Simulation runs were performed using Monte Carlo ray-tracing software (LightTools, Version 8.7.0, Synopsys, Mountain View, CA, USA). In the simulation model, the fiber diffuser is described as cylindrical geometry with a diameter of 1.1 mm and a length of 30 mm. Intensity values were associated to a pixel grid positioned on the cylindrical lateral surface of the modelled diffuser. The grid was composited from 3600 angular pixels (angular resolution: 10 pixel per degree) and 4000 longitudinal pixels (longitudinal resolution: 133 pixel per mm). A rotationally symmetric intensity distribution was assumed, as a worst-case scenario in terms of tissue heating for fixed overall maximum intensity value. For rotationally asymmetric distributions with the same maximum intensity, the overall impact on the surrounding tissue would be reduced. The longitudinal intensity values and the emission direction were adjusted according to the emission profiles measured by the imaging ( Fig. 1(a)) and non-imaging setup ( Fig. 1(b)), respectively. The emission angle distribution of the light rays was assumed to be Lambertian within limited angular ranges, revealing a good agreement between simulation and experiment in initial testings. As shown in Fig. 1(d), the angles 2α and β define the full opening angle of the respective emission cone and its tilt in relation to the surface normal, respectively. All simulation model runs were carried out with 10 8 rays.
The angles α and β were determined as follows. The emission profile of the reference diffuser measured with the imaging setup ( Fig. 1(a)) was imported into the simulation. In analogy to the non-imaging measurement illustrated in Fig. 1(b), in the simulation model, a virtual receiver plane was positioned in the simulation model parallel to the longitudinal diffuser axis at a distance of 2 mm from the diffuser surface, as sketched in Fig. 1(d). The spatial resolution on the virtual receiver plane was defined according to the resolution of the camera sensor with a pixel size of (4.65 × 4.65) µm 2 . Simulation runs were performed for α = 1°, 20°, 25°, 30°, and 90°and β = 0°. For α = 25°, additional runs were performed to define the fraction of forward and backward emitted rays with β = −35°, 0°, and 35°. The parameter set of α and β, which provided the best agreement between non-imaging experiment ( Fig. 1(b)) and simulation, was selected for all further simulations presented in this study.
Further MC-simulations mimicking the geometries of the phantom experiments ( Fig. 1(c)) were performed as illustrated in Fig. 1(d). For this purpose, the simulated cylinder diffuser was embedded in a medium to which the known optical properties of the artificial tissue phantom (brain tissue) were assigned. In analogy to the measurement of the light distributions for different depth positions of the diffuser below the phantom surface ( Fig. 1(c)), the light distribution in the simulation model was recorded at the corresponding distances from the diffuser surface (2 mm, 5 mm, 10 mm, 15 mm) using virtual receiver planes. For reasons of simplification, the refractive index mismatch between tissue phantom and air in the experiment was neglected in the simulation. However, initial tests have shown that this only minorly influences the results in the presented setups.
While the knowledge of the light distribution at a few selected distances from the cylindrical diffuser inside a surrounding tissue volume contains valuable information about the light distribution induced in the tissue by the diffuser, it only provides rather coarse information about the continuous evolution of the light distribution as a function of distance to the diffuser surface.
To investigate this evolution with higher resolution in the distance coordinate, an additional virtual receiver plane was introduced, in this case oriented perpendicular to the diffuser axis. In the longitudinal direction, it was placed at the position of the intensity maximum to be investigated. The resolution on the virtual receiver plane was defined by a pixel size of (100 × 100) µm 2 .
In addition, the impact of height and width of intensity peaks in the emission profile of cylindrical diffusers on the light distribution in tissue, as well as the interplay between two adjacent intensity peaks were studied by modelling the light distribution in tissue for artificially generated emission profiles. To model a single local intensity peak, a homogeneous intensity distribution over the entire diffuser length, f (x), was combined with a Gaussian distribution, f ∧ (x). The total emission profile f (x) was described according to Eq. (1).
I 0 is the background intensity, I max is the height, x 0 the position, and w FWHM the full width at half maximum of the added inhomogeneity on top. The position of the intensity maximum was set constant at x 0 = 15 mm. The intensity of the homogeneous emission was set to I 0 = 1, the initial maximum intensity I max was set to I max = 1.5, 2.0, 3.0, and 6.0, and the width parameter w FWHM to w FWHM = 0.2 mm, 1.0 mm, and 5.0 mm.
The impact of two adjacent intensity peaks on the light distribution in tissue was studied combining f with two Gaussian distributions f ∧ . The total emission profile was described by Eq. (2).
The distance between the two maxima was set to ∆x = 2 mm, 3 mm, 5 mm, 8 mm, and 10 mm. The Gaussian function parameter were set to I R (0) = 2.0 and w FWHM = (0.2 mm, 1.0 mm). Thereby, the influence of the width and the distance between the two intensity peaks on the light distribution in tissue was investigated.

Analysis part
The influence of inhomogeneities in the emission profile of a cylindrical diffuser on the light distribution in the surrounding tissue was analysed by determining the difference between light distributions from homogeneous emission profiles with and without added local intensity peaks at certain longitudinal positions. For this purpose, it is necessary to separate the additional contribution in the light distribution that is due to the local intensity maximum (A) from the one that is due to the perfectly or at least approximately homogeneous baseline emission (B). In this context, different evaluation procedures had to be applied for the light distributions measured in the experiment and those obtained from the simulations.
In case of intensity profiles measured experimentally after penetration through tissue phantom, the following approach was used to separate the fraction (A) attributed to local intensity maxima and the fraction (B) attributed to an approximately homogeneous baseline intensity. The experimentally measured intensity profile was modelled as a sum of five Gaussian curves, where four Gaussian curves were used to describe the slowly varying intensity distribution interpreted as fraction (B) and one Gaussian curve to describe the local intensity peak interpreted as fraction (A). A sum of five Gaussian curves was chosen, because five was the lowest number of Gaussian curves that provided a mean deviation of less than 1% between fit curve and measured curve. The sum of the five Gaussian functions that was finally used to fit the experimentally measured intensity profile can be written by Eq. (3) as where a denotes the peak height, b the position, and c the width of each Gaussian curve. In Fig. 3, an example of a fit curve according to Eq. (3) to an experimentally measured emission profile for a diffuser depth position of d = 5 mm below the phantom surface is shown. The blue curve describes the composition of four of the five Gaussian curves, fraction (B), attributed to the approximately homogeneous baseline intensity. The fifth Gaussian curve, fraction (A), models a local intensity maximum at the location x with a peak intensity I max . Dividing the maximal intensity I max,x (d) of the intensity peak at the location x by the corresponding baseline intensity I 0,x (d) at the same location x yields the intensity ratio I R (d) = where d denotes the diffuser depth below the phantom surface. In this way, the intensity maximum can be evaluated in comparison to the baseline intensity. Such evaluations of the measured intensity distributions at the phantom surface were performed for four different depth positions d of the reference diffuser below the phantom surface, d = 2 mm, 5 mm, 10 mm, and 15 mm below the phantom surface (n = 5 measurements each). Average value and standard deviation of the ratios I R (d) were calculated for each depth value d from the repeated experiments and evaluations. In this way, the fading of such intensity peaks in the experimental intensity distribution with increasing diffuser depth below the phantom surface was analysed. In case of intensity profiles, arising from simulated emission profiles after penetration through tissue, the fraction (A) attributable to local intensity peaks in the emission profile on the diffuser surface was separated from the background, fraction (B), by simulating the light distribution after tissue penetration for two scenarios: a perfectly homogeneous emission profile with and without artificially added intensity peaks of Gaussian shape.
For a single intensity peak at the position x, the fading of intensity at the position x with increasing depth d of the diffuser below the phantom surface was evaluated in comparison to the corresponding decrease in intensity at the same position x in the situation without added peak in the emission profile. Dividing the two intensity values at the position x yields the intensity ratio I R,single (d) = I max, x (d) Fig. 4, the evaluation process for determining I max, x and I 0, x is illustrated. A "50% tissue depth" d p,single at which the excess of the intensity ratio I R,single (d) over 1 drops to 50% of its initial value was defined by the criterion I R,single (d p,single )−1 I R,single (0)−1 = 0.5. By choosing the value 0.5 instead of the more common value 1/e to describe the decay, the corresponding depth value d p,single could be determined for all measured and simulated scenarios, even for parameter configurations resulting in a rather slow decay.  I max (d p,double ) becomes identical to I center (d p,double ) and, thus, the two intensity peaks are no longer discernible.

Results
In Fig. 6, 7, and 8 the radial emission profile along a diffuser of 30 mm length (0 mm: proximal diffuser part and entrance of light, 30 mm: distal end of diffuser with mirror) is shown for the different experiments on the reference diffuser and corresponding simulations. Fig. 6. Simulated emission profiles for a cylindrical light diffuser in air for a cone tilt angle β = 0°and five different cone opening angles α (α = 1°, 20°, 25°, 30°, 90°). The experimental emission profiles measured on the reference diffuser in air by the non-imaging and imaging setup ( Fig. 1(a) and (b)) are shown as black dashed and dotted line, respectively. The emission profiles at the surface of the reference diffuser measured by the imaging (dotted black line) and non-imaging (dashed black line) method according to Fig. 1(a) and (b) are depicted in Fig. 6. Although, the two intensity decays between the longitudinal positions 0 mm and 25 mm are very similar for the two measurement techniques, the amplitude and shape of the two intensity peaks observed between the longitudinal positions 28 mm and 32 mm differ dramatically. For comparison, five different emission profiles (five solid lines of different colours) were simulated according to Fig. 1(d), for an emission cone perpendicular to the optical axis (i.e.  Fig. 1(a)) and on the surface of the artificial tissue phantom (colored lines, according to Fig. 1(c)), with the reference diffuser submerged in the tissue phantom at depth positions of 2 mm, 5 mm, 10 mm, and 15 mm. β = 0 • ) and five different cone opening angles α. It can be perceived that the two intensity peaks at the distal diffuser end measured experimentally with the imaging method can be reproduced relatively well by simulation with α = 1 • (yellow line), while the measurement result of the non-imaging method can be reproduced better, at least qualitatively, with α = 25°(blue line). By comparing the five solid lines, it becomes obvious, that by increasing the cone opening angle α for fixed cone tilt angle β = 0 • , the intensity peaks broaden and decrease in amplitude.
Further, the influence of the cone tilt angle β on the emission profile was examined for a fixed cone opening angle α = 25°. In Fig. 7, the simulated profiles obtained for perpendicular emission (β = 0 • ), preferential forward (β = + 35 • ) and backward (β = − 35 • ) emission (three solid lines of different colours) are compared to the experimental result of the non-imaging method (dash-dotted black line). The scenario of emission perpendicular to the diffuser surface (β = 0 • , blue line) leads to rather broad intensity peaks. A preferential forward emission (β = +35 • , black line) results in intensity peaks with reduced FWHM and a tiny intensity redistribution into the forward direction in the longitudinal position range 0 -2 mm, as well as a more pronounced intensity redistribution for positions above 31.5 mm. By contrast, preferential backward emission (β = −35 • , red line) generates a shoulder in the emission profile near the longitudinal position 27 mm, as well as an intensity redistribution into the backward direction at the diffuser end near the longitudinal position 31 mm. It can be deduced that a variation of the cone tilt angle β in the simulation model results in minor influences on the shape of the simulated emission profile. Thus, further simulations were simplified by fixing the tilt angle to β = 0 • , thus mimicking an average emission direction perpendicular to the diffuser surface.
Based on this iterative process, the shape of the emission profile of the reference diffuser can be described by MC-simulations according to Fig. 1(d) with cone tilt angle β = 0 • and cone opening angle α = 25°. In all further investigations, the light distribution generated by the reference diffuser in the surrounding tissue phantom according to Fig. 1(c) was simulated with these fixed values for the simulation parameters α and β.
In Fig. 8 the intensity profile measured on the surface of the tissue phantom with the diffuser placed at a depth of 0 mm (i.e. in air), 2 mm, 5 mm, 10 mm and 15 mm, measured with the imaging method is shown. The results reveal that the intensity distributions at the phantom surface appear more homogeneous than the emission profile of the diffuser in air, and that they become broader with increasing depth of the diffuser below the phantom surface. Characteristic intensity peaks, like the two intensity peaks at the distal diffuser end, are strongly reduced already for a diffuser position 2 mm below the surface.
In Fig. 9, the simulated (black line) and calculated (blue dots) intensity ratios I R,single (d) are shown as a function of the diffuser depth d below the phantom surface. The comparison shows a close relationship between measurement and simulation within the error limits of the measurement. The simulation reveals that the excess of intensity ratio I R,single (d) over 1 drops to 50% (from I R,single (0) − 1 = 4.61 to I R,single (d) − 1 = 2.31) within the first 100 µm in tissue and reaches a value of 10% (I R,single (d) − 1 = 0.46) after about 2 mm. Thus, it can be concluded that local intensity peaks in the emission profile of the diffuser in air fade away quite effectively within the first 1 mm of tissue depth in the case the diffuser is inserted in an artificial brain tissue phantom. Simulations for virtual diffusers with artificial emission profiles, placed into an infinite tissue phantom according to Fig. 1(d), were performed to show the influence of height and width of an intensity peak on top of a homogeneous emission profile and the interplay of two adjacent intensity peaks on the light distribution in tissue. For this purpose, artificial emission profiles were generated and implemented into the MC-simulation scenario mimicking the situation of a cylindrical diffuser submerged in an artificial tissue phantom as shown in Fig. 1(c).
First, the width w FWHM and height I R,single (0) of a single intensity peak in the center position of the diffuser, x 0 = 15 mm, were varied. Figure 10 shows the relative peak intensity I R,single (d) as a function of the tissue depth d from the diffuser for intensity peaks with two different widths of w FWHM = 0.2 mm (Fig. 10(a)), w FWHM = 1.0 mm (Fig. 10(b)) and w FWHM = 5.0 mm (Fig. 10(c)). For each width w FWHM , four scenarios with four different initial values of the relative peak intensity were examined, for I R,single (0) = 1.5, 2.0, 3.0, and 6.0. All curves are characterised by a steep drop in intensity within the first millimetres. The corresponding "50% tissue depths" d p,single , defined as Table 1. This comparison shows that the width of the intensity peak, in contrast to the height, has a strong influence on d p,single . d p,single was below 100 µm for a width of w FWHM = 0.2 mm, while it was determined to be 0.79 mm -0.91 mm and 1.91 -2.11 mm at a width of w FWHM = 1.0 mm and w FWHM = 5.0 mm, respectively. This finding suggests that narrow intensity peaks lead to more superficial tissue damage close to the diffuser surface, while broader ones penetrate deeper into tissue.
In Fig. 11(a) simulated intensity profiles are shown for an artificial diffuser with two local intensity peaks of width w FWHM = 1.0 mm separated by a peak-to-peak distance of ∆x = 5 mm. Results are shown for several tissue depths from the diffuser. The diagram shows that the intensity Fig. 10. Evolution of the intensity ratio I R,single (d) with increasing tissue depth for artificially generated emission profiles with one local maximum. The intensity peak was varied in width w FWHM and initial height I R,single (0). Results for a width w FWHM of 0.2 mm, 1.0 mm, and 5.0 mm are illustrated in panels (a), (b), and (c), respectively. For each width, four scenarios with four different initial intensity ratios were simulated, for I R,single (0) = 1.5, 2.0, 3.0, and 6.0. Table 1. "50% tissue depth" d p,single at which the intensity ratio I R,single (d) drops to 50% of its initial value at the diffuser surface. at the positions of both maxima decreases with increasing tissue depth, but the intensity in the center between the two maxima increases. For a tissue depth of d = 5 mm, the two intensity peaks have already merged into a single, broad intensity peak. To provide a deeper insight into the interplay of two adjacent intensity peaks, their distance and width and w FWHM was varied. The intensity ratio I R,double (d) = I max (d) I center (d) was calculated and plotted as function of the tissue depth d for w FWHM = 0.2 mm and 1.0 mm (Fig. 11(b) and (c)). The intensity peaks become indistinguishable when the intensity ratio I R,double reaches 1.
The corresponding "peak coalescence tissue depths" d p,double are listed in Table 2 for different separations ∆x of the two peak maxima. It becomes obvious, that for one parameter setting (w FWHM = 0.2 mm and ∆x = 2 mm), the two peaks become indistinguishable within 1 mm of tissue depth, while for an initial peak distance of ∆x = 10 mm, the two intensity peaks remain distinguishable over the full range of investigated tissue depths, independently of the initial peak width w FWHM . Comparing the differences between narrower (Fig. 11(b)) and wider intensity peaks (Fig. 11(c)), it can be noted that narrower intensity peaks lead to a faster attenuation of I R,double . The tissue depth d p,double at which both peaks merge does not differ for the two different peak width values investigated. Overall, it can be suggested that emission profiles consisting of several intensity maxima with a peak-to-peak distance below 2 mm lead to a homogeneous light distribution within the first 1 to 2 mm of tissue depth.

Discussion
In this investigation, the influence of local inhomogeneities in the emission profile of fiber diffusers on the light distribution within tissue after interstitial insertion was studied. This was assessed by dedicated MC-simulations based on reproducing experimental measurements.
The experimental measurements were performed using an imaging and a non-imaging camera method [32]. This made it possible to describe the emission behaviour of the diffuser more precisely. While the imaging method could describe the intensity distribution on the surface, the direction of the emission was measured with the non-imaging method. Both measurements were performed in air and the light distribution in tissue was deduced by means of a MC-simulation. Although the simulation was validated by measuring the light distribution in tissue, the question of a measuring method continuously and easily accessing the light distribution in tissue still arises.
Although the MC simulation presented was calibrated and validated by measurements, some assumptions were made. The emission angle was set to β = 0°. However, this provides a good approximation only for volume-scattering diffusers with a reflecting mirror at the fiber end. The mirror causes light to propagate in both fiber directions, away from and towards the light source. The emission should therefore be almost symmetrical about the angle β = 0°. For diffusers without mirrors or surface scattering diffusers, this assumption loses its validity and forward directed emission (β > 0°) must be expected.
Furthermore, the refractive index mismatch between tissue phantom and the surrounding medium air was neglected in the simulation. However, it can be estimated that the angle α = 25°u sed here would be increased to 33°by a refractive index mismatch. As a consequence, the emission profile, including potential inhomogeneities, would appear broadened. However, the calibration performed in Fig. 6 shows that the intensity distribution for an emission angle of α = 30°deviates only slightly from that for the selected angle α = 25°.
Other approaches to air-based measurement include integrating sphere measurement setups [20,32], goniometric approaches [27,34] or photodiode-based measurement methods [35]. A description of the radiation profile by an integrating sphere and photodiode-based method would theoretically lead to similar results as the imaging and non-imaging camera method. An experimental validation will be left to further investigations. A complete description of the light distribution in air can only be obtained by goniometric imaging of the radiation pattern. In addition to the intensity distribution on the fiber surface, the goniometric representation can also directly measure the directional radiation [34]. The disadvantage of this method is the difficulty in aligning the diffuser in the measurement setup coupled with sophisticated mathematical evaluation procedures [34].
The vanishing of single and multiple intensity maxima during light penetration through a scattering tissue phantom were examined. It was shown, that the light distribution is mainly influenced by the width w FWHM . Intensity peaks with a width in the range of one millimetre have an approximate "50%" and "peak coalescence" tissue depth d p of one millimetre in tissue. An explanation might lie in the light propagation originating from point and line sources. Since the light intensity flux of point and line sources decreases geometrically by I ∼ 1 r 2 and I ∼ 1 r [36], respectively, narrow intensity peaks, representing point sources, decay faster with increasing distance to the diffuser surface compared to wide peaks, representing line sources [33,37]. An analytical solution proves to be difficult, as the transition between point source to elongated line source is continuous. The deeper penetration depth of extended light sources is also used in Diffuse Reflectance Spectroscopy [38][39][40]. Here, light emitted by one fiber is detected by another. If the distance between the fibers, or the fiber radius itself, is increased, deeper tissue regions can be monitored [38][39][40].
Comparing the simulated results to emission profiles in literature, an evaluation of the profiles with regard to the width w FWHM of the intensity peaks has to be performed. The emission profiles presented in [27,30,31] have pronounced, however, narrow intensity peaks at the distal and/or proximal diffuser end. While the height I R (0) can reach up to three times the average emitted intensity, the width w FWHM of the intensity peaks does not exceed a width of w FWHM = 1 mm. Therefore, it can be concluded, that the influence towards the light distribution is limited to the first millimetre of tissue depth from the diffuser surface.
The light distribution of a diffuser in tissue phantom shown in [41] slightly differs from the results presented here. An emission of a width of 5 mm showed a reduction of the intensity to 50% of the original intensity after only 0.9 mm. The faster decrease might be due to the fact that sinusoidal intensity peaks were used and thus the geometry differs from the here presented Gaussian inhomogeneities.
The investigation on the influence of two narrow peaks in the emission profile on the light distribution within the tissue showed that two peaks with an initial peak-to-peak distance of ∆x = 2 mm become indistinguishable after approximately 1 mm tissue depth. For larger peak-to-peak distances ∆x, the two peaks remain separate, and each of them acts as a single intensity peak propagating into the tissue. Applying the gained knowledge to the field of diffuser production, it can be concluded that after only a few millimetres of tissue, an emission profile composed from multiple close intensity peaks over the full diffuser length exhibits an equivalent behavior comparable to homogeneous emission profiles. Other emission profiles, composed from multiple intensity peaks over the full diffuser length, showed a peak-to-peak distance between the intensity maxima of ∆x = 0.5 mm [19,25]. Therefore, the light distribution has a homogeneous intensity distribution after less than 1 mm of tissue depth.
The light distribution in tissue and the attenuation of single and multiple intensity peaks on emission profiles was extensively discussed. While it is difficult to make a statement about the biological influences, a conclusion about spatial allocation of potential tissue damage can be done. The local excess intensity within the tissue induced by local maxima in the emission profile is largest close to the diffuser surface. Thereby, the first cell layers close to the diffuser are affected by strong intensity peaks. This can affect about 2 to 20 cell layers since each cell layer has an approximate size of 5 -50 µm. However, the simulation model lacks a description of actual damage to tissue or the diffuser itself due to thermal overtreatment. Therefore, the thermal properties of the targeted tissue has to be examined, the calculated light intensities transferred to temperature, and the thermal diffusion considered. Exceeding 42°C in brain tissue can lead to cellular damage causing side effects to the treatment [17].
While in the present work the influence of inhomogeneities on the light distribution in brain tissue was presented, their influence on other tissue types must be critically discussed, since the optical as well as the thermal properties of different tissue types strongly differ from each other. A decrease in absorption and/or scattering may lead to an enhanced penetration of inhomogeneities to deeper tissue regions, which is the case for skin and prostate tissue at a wavelength of λ = 635 nm [42]. Furthermore, a variation in thermal properties can influence the temperature distribution in the targeted tissue region. For increased heat capacity, a lower temperature rise is expected, while for increased thermal conductivity, a locally induced temperature rise fades away effectively into the surrounding tissue.
Further investigations will evaluate the effect of certain light doses on the temperature increase in various types of tissue. Mathematical approaches could be built on the solution of Penne's bio-heat equation [16,27,43,44]. Thereby, an applied diffuser power could be defined at which the temperature increase is limited to a total temperature of 42°C. Thus, quality criteria for the homogeneity of diffusers could be defined.

Conclusion
Inhomogeneities in the measured emission profile of cylindrical light diffusers in air affect the light distribution in the surrounding tissue. The observed effects are mainly dependent on the width of the local intensity peaks.
Narrow intensity peaks with a width smaller than one millimeter in the emission profile have only a minor effect on the light distribution somewhat deeper within the tissue, while broader or a multitude of intensity maxima can maintain a strong effect also at larger distances from the diffuser. Potential overheating near narrow intensity maxima should be avoided, in particular in case of tissue with strong absorption and scattering and/or low heat capacity and/or low heat conductivity. Clinically, strong local intensity maxima constitute a potential risk factor for sensitive tissue structures in this area.