Analytic function for predicting light fluence rate of circular fields on a semi-infinite turbid medium

: Accurate determination of in-vivo light fluence rate is critical for preclinical and clinical studies involving photodynamic therapy (PDT). The light fluence distribution in tissue depends on both the tissue optical properties and the incident field size. This study compares the longitudinal light fluence distribution inside biological tissue in the central axis of circular uniform light field with different radii for a range of in-vivo tissue optical properties (absorption coefficients (µ a ) between 0.01 and 1 cm − 1 and reduced scattering coefficients (µ s ’) between 2 and 40 cm − 1 ). This was done using Monte-Carlo simulations for a semi-infinite turbid medium in an air-tissue interface. The end goal is to develop simple analytical expressions that would fit the results from the Monte Carlo simulation for circular beams with different radii. A 6-parameter model ( can be used to fit MC simulation. Each of these parameters ( b , C 2 , C 3 , λ 1 , λ 2 , and λ 3 ) is expressed as a function of tissue optical properties and beam radius. These results can then be compared against the existing expressions in the literature for broad beam for analysis in both accuracy and applicable range. The analytical function can be used as rapid guide in PDT to calculate in vivo light fluence distribution for known tissue optical properties.

can be used to fit MC simulation. Each of these parameters (b, C 2 , C 3 , λ 1 , λ 2 , and λ 3 ) is expressed as a function of tissue optical properties and beam radius. These results can then be compared against the existing expressions in the literature for broad beam for analysis in both accuracy and applicable range. The analytical function can be used as rapid guide in PDT to calculate in vivo light fluence distribution for known tissue optical properties.

Introduction
Light propagation in tissue has been extensively studied for the past three decades to understand how light energy is distributed and absorbed within the target tissue. The in vivo light fluence distribution is influenced by numerous factors including the tissue optical properties, the tissue geometry and the irradiation geometry. Accurate determination of light fluence distribution in tissue is important to ensure successful photodynamic therapy (PDT) treatment as the rate of photochemical reactions depends on the light dose delivered to the tissue.
In vivo light fluence distribution can be directly measured by inserting an isotropic detector interstitially [1]. Due to the invasiveness, this technique is limited to measurement of light fluence rate at selected points or along a catheter rather than in a full three-dimensional volume. Model calculations such as Monte Carlo simulation [2] and diffusion approximation theory [3] are therefore more practical and have been widely employed to estimate the light distribution in tissue noninvasively, given that the 3D distribution of the tissue optical properties is known. Diffusion theory is attractive due to its simplicity in calculating light fluence rate quickly and accurately in scattering media. It has been found useful in various tissue optics applications as its validity condition of a μ << ' s μ is satisfied for most biological tissues at visible to near infrared wavelengths [4]. However, for regions that are close to sources or boundaries, diffusion theory is not appropriate as photons have a preferential direction of travel and cannot be approximated as a diffusion process. The light distribution described by diffusion theory is therefore less accurate in the first several hundred microns of tissue which can be an important region for PDT light-tissue interactions for surface irradiation [4]. Monte Carlo modeling, on the other hand, can be used to precisely simulate the light distribution in tissue with complex geometries and inhomogeneities and is often used as the gold standard method to model light transport in turbid media. However, the requirement of intensive computation load to achieve results with high accuracy has become the main limitation of MC method in application where speed is desirable [5]. Simple analytical or empirical expressions which allow for fast calculation of light distribution that mimic the accuracy of Monte Carlo simulations are therefore desirable in PDT treatment planning to optimize light fluence rate distribution during treatment.
Literature search reveals only analytical expression for broad beams, when the lateral dimension is much larger than the mean-free path [6,7]. The light distribution of broad beams in tissue is only dependent on the two dimensional tissue optical properties ( = + [6]. For small diameter circular light fields, the distribution of light fluence rate becomes dependent on the field size when the lateral dimensions are comparable to the effective mean-free-path of photons for optical properties relevant in PDT. Although this field size effect is not of great concern in most clinical PDT in which expanded or broad beams are usually used, it has to be taken into careful consideration in the treatment planning of experimental PDT using small animals where small diameter beams are often used to treat surface lesion as small as 1-cm in diameter [8,9]. For clinical cases such as PDT of skin and recurrent breast cancers on the chest wall, small circular light fields of ≤ 2 cm diameter are sometimes used. It is common that the size of the treatment beam is determined based on the tumor size and the identical light fluence rate is used for the PDT treatment irrespective of the beam size [10].
MC simulation of light distribution clearly demonstrated that in vivo light fluence distribution changes with incident field size with significant reduction of light fluence rate at the central axis for small diameter circular beam from that for a broad beam [11][12][13]. Light scattering causes the light fluence distribution in tissue to be dependent on the incident beam diameter. Photons are scattered and diffuse radially away from the source. Closer to the edge of the beam, there are reduced scattering photons thus creating a large radial gradient in fluence rate at the beam perimeter. The edge effects overlap at the center of the beam when the incident beam diameter is narrow causing the central light fluence rate to be significantly reduced [12]. As the beam diameter increases, the central axis is less affected by the edge effects and the radius dependence of the in vivo light fluence distribution reduces. When the incident beam diameter is broader than 3 cm, the fluence rate at the center of beam approaches that for an infinitely wide beam with uniform incident irradiance [12]. It is therefore important for PDT treatments using small irradiance beam to be calibrated carefully for radius dependence in order to deliver optimal light dose which would improve the PDT treatment outcome.
In this paper, we will determine the central longitudinal light fluence rate distribution of circular fields with radii of 0.5 cm to 8 cm for a range of tissue optical properties which are relevant in PDT and many medical applications (absorption coefficients ( a μ ) between 0.01 and 1 cm −1 and reduced scattering coefficients ( ' s μ ) between 2 and 40 cm −1 ) based on a review on the in-vivo tissue optical properties [14]. These results are expressed as simple analytical functions in term of beam radius and tissue optical properties ( ' , a s μ μ ), and will be compared with Monte-Carlo simulations in term of maximum relative deviations as summarized in Tables 1 and 2. We used Monte-Carlo simulation because the diffusion theory is not valid when the lateral dimension of beam geometry becomes comparable to the meanfree-path of the photons or when the region of interest is near the air-tissue interface. These results will then be compared against the existing expressions in the literature for both accuracy and applicable range.

MC Simulation
The setup geometry to be calculated by Monte-Carlo simulation for a semi-infinite medium with uniform optical properties, i.e. the absorption coefficient µ a , the scattering coefficient µ s , the scattering anisotropy g ( = 0.9), and the index of refraction (n 2 = 1.4 for tissue) is shown in Fig. 1(a). The outside medium is air (n 1 = 1). The light field is parallel and uniform inside a circle with radius R and incident toward the air-tissue interface. (R = 0.5, 0.75, 1, 1.5, 2, 3 and 8 cm.) Most tissues have a preferential scattering in the forward direction with a scattering anisotropy g = 0.90 [15]. In the diffusion approximation, the tissue scattering is converted to isotropic conditions (g = 0) and the tissue scattering property is described by a reduced scattering coefficient To improve the scoring efficiency, Monte-Carlo simulation is performed for a pencilbeam incident on a semi-infinite tissue phantom as shown in Fig. 1(b). The photon density ( We have also scored the fractional photons escaping the air-tissue interface into the air as a function of lateral radius (with a resolution of 0.05 cm) and emitting angle, θ (with a resolution of 1°). The reciprocity theorem is used to calculate the angular distribution of differential diffuse reflectance for circular fields as a function of radius and emitting angle. The light fluence rate above the tissue surface is calculated using the diffuse reflectance r d for circular field with radius r, as previously described [11,17,18], Equation (2) states that the light fluence at a tissue surface is composed of 1 part from the incident photons and 2r d part from the backscattered photons as described in the literature [19].
The Monte Carlo algorithm is similar to approaches used previously in the literature [2,11,20,21]. The effect of reflection at the air-tissue interface, resulting from the refractive index mismatch, was modeled according to the Fresnel reflection coefficient for unpolarized light [20]. The initial weight of each launched photon is assigned to be unity. Photons are propagated through the medium step by step with random step size and scattering angle sampled based on the optical properties of the tissue model and their respective probability distributions. Each step size, s Δ , is determined by ln where ξ is a uniform distributed random number and l is the mean free path. The scattering angle is selected from the phase function determined by the Henyey-Greenstein phase function [21] with anisotropy g = 0.90. At the end of each step, the weight of the photon is reduced according to the absorption probability. At least 2 millions incident photons are launched normal to the interface along the z direction for each Monte Carlo simulation. We followed each of them until it escapes from the tissue model or when it is completely absorbed. The MC result provides an approximation to light distribution in tissue with a typical statistical uncertainty of less than 0.1%. The small standard variation for 2 millions photons is possible because a convolution of pencil-beam result is used to calculate the broad beam parameters.

Analytical expression and fitting algorithm
For small diameter circular fields, we found that Eq. (3) can fit the resulting MC data very well, typically to a depth of at least 0.5 cm or a few mean-free-path, l, whichever is larger. This fitting equation is termed 6-parameter model since it includes 6 fitting parameters: λ 1 , λ 2 , is called in-air fluence rate that can be calculated using the total power, P (in mW), and the treatment area, A, of the collimated beam. The 6-parameter model ( ) describes the attenuation of the scattered photons. However, a single exponential decay is found to be insufficient to describe the attenuation of the scattered photons for small field sizes, thus we have added a second exponential decay term. λ 1 is the rate of buildup of scattered photons, 1-b is the ratio of fluence rate at surface to in-air fluence rate, air φ φ (d = 0), C 2 is the backscatter coefficients for scattered photons of broad beam and λ 2 is the effective attenuation coefficient of scattered photons, 3 3 d C e λ − term describes the reduction of the backscatter coefficients for scattered photons near the surface for small field sizes.
For a broad beam, existing literature shows that a 4-parameter model (λ 1 , λ 2 , b, C 2 ) is sufficient to fit the MC data [6,7]: We have modified the original form ( The objective function for minimization is defined as the least square difference between analytical fit of fluence rate, fit φ , and the MC simulation, MC φ , for a range of tissue optical properties μ a between 0.01 to 1 cm −1 and μ s ' between 2 and 40 cm −1 and beam radius between 0.5 and 8 cm. where fit φ is based on either Eq. Multivariable optimization algorithm uses a differential evolution algorithm described elsewhere [22]. The non-linear optimization routine was performed using Matlab.  Tables 3 to 7 in the appendix section. Tables 3, 4, and 5 show the summary results of 6-parameter fit to the MC simulation data for small circular fields with diameter of 1 cm, 1.5 cm, and 2 cm, respectively, for depth up to 3.5 cm deep using Eq. (3). For 1 cm and 1.5 cm diameter circular beams, C 3 approaches zero for a subset of data when μ eff > 2 cm −1 while for 2 cm diameter circular beam, C 3 approaches zero when μ eff > 1 cm −1 so that only a 4-parameter fit is necessary. In these cases, the value of λ 3 is unnecessary and its value is expressed as "-" in Tables 3, 4 and 5. Tables 6 and 7 show the summary results of 4-parameter fit using Eq. (4) to the broad beam MC simulation data (diameter of 8 cm) for an air-tissue interface (n 2 /n 1 = 1.4) and water-tissue interface (n 2 /n 1 = 1.053), respectively.

Broad beam
For light fluence rate of broad beams, we have found that all four fitting parameters in Eq. (4) are one-dimensional function of R d as used by Jacques [6], where R d is the diffuse reflectance of broad beam. The analytic expression of R d can be found in the discussion section of this paper. The two fitting parameters (λ 1 and λ 2 ) should be proportional to μ eff . We have found the following relationships: Inside tissue: On tissue surface

Small circular field
For small circular field with beam diameter smaller or equal to 2 cm, we found that all the 6 parameters are as a function of radius, r. The parameters b, C 2 and λ 1 are two dimensional functions of R d and radius, r and the parameters λ 2 is two dimensional function of μ eff and radius, r. C 3 Table 6 and Table 7, respectively. Red triangles represent parameters using formulas from Ref [6]. R 2 refers to the goodness of fit between the fitting results of each parameter to their respective equation. Figure 4 shows the fitting results for the 4-parameters Eq. (4) for broad beam, and Eqs. (6) - (9) and Eqs. (20)-(23) for their respective parameters at different refractive index mismatch. The fitting is performed using the curve fitting tool from Matlab global optimization toolbox (cftool.m). R 2 of the fitting results are all very good (R 2 > 0.96). For comparison, the fitting results for broad beam using the original form of Eq. (4) reported by Jacques [6] are included in Fig. 4 Tables 3 to 7 and Fig. 5.
Analytical expression based on the same form as a diffusion approximation [3,24] can be determined between R d and the albedo, a' as follows: where the albedo, a' is defined in Eq. (18). Figure 5 shows the goodness of the fit between R d and '/ (1 ') a a − using Eqs. (19) and (24), for air-tissue and water-tissue interface, respectively.  Figure 6 shows the fitting results for the broad beam using 4-parameters, Eqs. (6)-(9) and small circular beams with diameter of 1 cm, 1.5 cm and 2 cm using 6-parameters, Eqs. (11)-(16). Notice that for small circular field, the scaling theorem is no longer in existence that most of the parameters have to be expressed as at least two dimensional function of diffuse reflectance, R d or tissue optical properties ( ' , a s μ μ ) and beam radius, r. R 2 of the fitting results are also very good (R 2 > 0.95, not shown in Fig. 6). To minimize the potential of overfitting in the multi-variable optimization, we started the fitting using 4 parameters (Eq. (4)) to as deep a depth as possible that still fit the data, e.g., 0.5 cm for μ a = 0.1 cm −1 and μ s ' = 10 cm −1 , of the MC data to get the best fit for parameters b and λ 1 , which account for the buildup region of the light fluence rate, and C 2 and λ 2 . Then, by making parameters b and λ 1 constants with values obtained from the first fit and C 2 value to be the same as that for broad beam, parameters λ 2 , C 3 and λ 3 are refitted using Eq. (3) to MC data to a greater depth, at least 3.5 cm. By using this method, we are able to obtain consistent fitting results. A differential evolution algorithm [22] is used to perform the optimization and ideally find a global minimum for either 6 parameters or 4 parameters fitting results to the MC data. The differential evolution algorithm is a method that optimizes a problem by iteratively trying to improve a candidate solution with regard to a given measure of quality. This method is also known as metaheuristic as it makes few or no assumptions about the problem being optimized and can search very large spaces of candidate solutions. We found that 6 parameters are necessary to fit the data for small circular field data, particularly those with beam diameters smaller than 3 times light penetration depth, δ. Light penetration depth is the reciprocal of effective atten than zero and ' , a s μ μ ). C 3 te results as it a cylinder of th to the increas to or larger t tissue optical beam with di results as C 3 t    Tables 3 to 6. 2 , λ 1 and λ 2 to b auses a depleti n beam radius is smaller than b increases by 2 the increase in . λ 1 /µ eff increa dependence o s is larger than and (4), we c small and a fas dence on beam n be described adius and appro m radius ≤ 2.5 = 0.5 cm, λ 2 i C 3 is found to b sue optical pro circular field to uence rate at th e optical prope en beam diame of beam diam scattering decr d to fit the dat eff , and λ 3 ) and the 1 cm and 8 cm and e the red, blue and th their respective beam radius u ion of fluence is smaller than n 1 cm for albe 23% when bea n b is about 10% ases rapidly w f λ 1  µ eff is 0.5 cm −1 , 2.5 times when µ eff is 1 cm −1 , 1.5 times when µ eff is 2.5 cm −1 and the ratio approaches to 1 when µ eff is larger than 7.5 cm −1 .   Table 6, as well as Eqs. (6)- (9). For comparison, we have also included the fitting results in the literature [6]. We can see that the fitting results from Ref [6]. are good for μ a ≥ 0.1 cm −1 and it is unable to fit the buildup region accurately for μ a = 0.01 cm −1 . The fitting results using 4 parameters and Eqs. (6)-(9) are very good for all three optical properties presented in Fig. 8(a) for depth up to 3.5 cm. Equation (4), which is modified from the original form from Ref [6]. is able to improve the fitting crosstalks between parameters b and λ 1 , which correspond to the buildup region, as shown in Fig. 4. Figure 8(b) shows the comparison between MC simulation results of / air φ φ of small diameter circular beam with radius of 0.5 cm, 0.75 cm and 1 cm using the fitting parameters in Tables 3, 4, and 5, respectively, as well as Eqs. (11)-(16). We can see greater radius dependence of / air φ φ in small circular diameter beams when μ a and μ s ' decreases. In fact, Zhu et. al. [11] has reported that the radius dependence of / air φ φ is a function of μ a and μ s '. The relative deviations of both broad beam and small diameter circular beams for these three tissue optical properties over a depth of 3.5 cm are presented in Fig. 9. Relative deviation is defined as the difference between the fit and the MC result relative to the maximum fluence rate.  Dashed line are 4 parameters fits for (a) and 6 parameters fits for (b). Black lines in (b) correspond to beam with radius = 0.5 cm, blue lines correspond to beam with radius = 0.75 cm, and red lines correspond to beam with radius = 1 cm. Fig. 9. Relative deviation of the fitting results to the MC simulation data for (a) broad beam and (b) small diameter circular beam for three tissue optical properties (μ a , μ s ') = (i) (1, 40), (ii) (0.1, 10), and (iii) (0.01, 2) cm −1 . Data for (i) and (ii) are vertically shifted for clarity: + 25% for (ii) and + 30% for (i) for (a) and + 10% for (ii) and + 16% for (i) for (b), respectively. Black lines in (b) correspond to beam with radius = 0.5 cm, blue lines correspond to beam with radius = 0.75 cm, and red lines correspond to beam with radius = 1 cm. From Fig. 9(a), we can see that Eqs. (6)-(9) fits all broad beam MC data to a relative deviation of ± 3% for three optical properties up to 3.5 cm depths, while parameters from Table 6 fits the MC data to within 1% relative deviation. Parameters from Ref [6]. fit most of the data except for = 0.01 cm −1 , where 20% error is observed. Among the difference between our MC data and Jacque's fit, 3% can be attributed to difference in the value of index of refractions used for the MC simulation: 1.4 and 1.38 between us and Jacques, as well as other subtle differences in MC simulation. Figure 9(b) shows the relative deviation of small diameter circular beam fits to the MC simulation results. For optical properties = 1 cm −1 and = 40 cm −1 , both 6 parameters fitting and Eqs. (11)-(16) are able to fit the MC calculated / air φ φ accurately to within 1.5% relative deviation. It should be noted that C 3 and λ 3 term can be ignored for this optical properties as they approach to zero for all the beam radii studied. A 4parameters fitting using Eq. (4) should give similar accuracy to the results presented. For optical properties μ a = 0.1 cm −1 and μ s ' = 10 cm −1 , the maximum relative deviation of 6 parameters fitting and Eqs. (11)-(16) to the MC calculated / air φ φ is 4.5%. For optical properties μ a = 0.01 cm −1 and μ s ' = 2 cm −1 , 6 parameters fitting and Eqs. (11)-(16) are able to fit the MC calculated / air φ φ to within 4% and 9% relative deviation. The maximum relative deviations of the fittings using Eqs. (6)-(9) and Eqs. (11)-(16) to broad beam and small diameter circular beams for all optical properties studied, up to a depth of 3.5 cm, are summarized in Tables 1 and 2.   Table 1 shows the maximum relative deviation of the 4 and 6 parameters fitting to the MC simulation results for optical properties used in model training (µ a = 0.01, 0.1, 0.5, 1 cm −1 and µ s ' = 2, 10, 23.5, 34.5, 40 cm −1 ) and beam radii = 0.5, 0.75, 1, 3, 8 cm. Independent MC simulations using different combinations of optical properties (µ a = 0.05, 0.3 and 0.8 cm −1 ; µ s ' = 5, 15, 25 cm −1 ) and same beam radii were run to compare the performance of the 6parameter model. The maximum relative deviations of the fitting at these additional optical properties to the MC simulation data are shown in Table 2. It should be noted that 4 parameters fitting using Eqs. (11)- (14) are performed for small circular fields with beam diameter larger than 3 times light penetration depth, δ. The maximum relative deviation is defined as the maximum difference between the fit and the MC result relative to the maximum fluence rate. 4-parameters fitting using Eqs. (6)-(9) is sufficient to fit for broad beam for all optical properties to within 4.5% relative deviation of the MC results. For small diameter circular beam (radius = 0.5 cm, 0.75 cm, 1 cm and 3 cm), 6-parameters fitting using Eqs. (11)-(16) is able to fit for higher scattering coefficient (μ s ' ≥ 23.5 cm −1 ) to within 6% relative deviation and for lower scattering coefficient (μ s ' ≤ 10 cm −1 ) to within 10.6% relative deviation of MC results. For small diameter beam, when beam diameter is greater than 3 times of light penetration depth, C 3 approaches zero and both C 3 and λ 3 terms has minimal impact on the overall fitting to the MC calculated / air φ φ . When beam diameter is smaller than 3 times of light penetration depth (or 3/µ eff ), C 3 and λ 3 terms are significant and has greater impact on the overall fitting to the MC calculated / air φ φ as beam diameter decreases.
For instance, the light penetration depth at optical properties μ a = 0.01 cm −1 and μ s ' = 2 cm −1 is 4.07 cm. The diameter of a 1-cm beam at this optical properties is smaller than 3/µ eff and the maximum relative deviation for 6-parameters fit to MC results is 3.5%. The maximum relative deviation increases to 140% (result not shown in Table 1) when the C 3 and λ 3 terms are excluded from the fitting.  Another advantage of our expressions, Eqs. (3) and (4), is that the fitting based on diffusion approximation can be obtained by dividing them by ( )  Figure 11 shows the relative deviation between the fit without the buildup term and the MC simulation, normalized to the maximum MC calculated fluence rate. Beyond the buildup region, the maximum error is 1% for the broad beam and around 5% for the small diameter circular beam. Within the buildup region, the maximum error is 18% for the broad beam and 25% for the small diameter circular beam. The buildup term has minimal impact on the light fluence distribution on the central axis for both the broad beam and small diameter circular beam when the absorption and scattering coefficients are high (μ a ≥ 0.1 cm −1 , μ s ' ≥ 10 cm −1 ) but has a larger impact on the central light fluence rate distribution for low absorption and scattering coefficients (μ a = 0.01 cm −1 , μ s ' = 2 cm −1 ).

Effect of refraction index mismatch
Intracavitory clinical PDT cases such as intraperitoneal and pleural PDT often involve watertissue interface with a different refractive index mismatch than air-tissue interface [17,[25][26][27]. At water tissue interface the refractive index mismatch is about 1.053 with a critical angle of 71.8° compared to air-tissue interface with refractive index mismatch of 1.4 and critical angle of 48.8°. Fewer diffusely backscattered light that reaches the surface of the water-tissue boundary is reflected back into the medium. The reduction in total internal reflection and increment in escaping reflectance augments the behavior of light fluence rate near the surface. The expressions for the 4 parameters in Eqs.
The light fluence rate on the tissue surface can be calculated using Eq. (10) where it is composed of 1 part from the incident photons and 2 R d part from the backscattered photons. The value of R d depends on both the optical properties and also the refractive mismatch at the tissue surface boundary. At water-tissue interface, the analytic expression of R d can be determined as follows. The goodness of the fit between R d using Eq. (24) and MC data as a function of a'/(1-a') can be found in Fig. 5. It should be noted that the expression of diffuse reflectance R d for Eqs. (19) and (24) depend only upon the transport albedo a' and the ratio of the index of refraction, n 2 /n 1 , and one can calculated the total diffuse reflectance for any given optical properties and refractive indices using the analytical expression based on the same form as a diffusion approximation given by Farrell [3]. The radius dependence of R d can be calculated using Eq. (17).

Conclusion
We have performed MC simulation for air-tissue interface and provided a 4-parameters and 6-parameters expressions for the longitudinal light fluence rate distribution on the central axis with an accuracy of better than a relative standard deviation of 4.5% for broad beams and 10.6% for the small diameter circular beams for the full range of the in-vivo tissue optical properties studied. Simple analytical expressions are provided to determine the parameters as a one dimensional function of diffuse reflectance, R d , for broad beam and as at least two dimensional functions of tissue optical properties, μ a and μ s ', or diffuse reflectance and beam radius, r. We have also independently validated the analytical fitting result for broad beam and small diameter circular beams. We found that 4-parameters expression is able to fit accurately for circular beam with beam diameters larger than 3 times optical penetration depth for all optical properties studied. 6-parameters expression is necessary to fit for small circular diameter beams with beam diameters smaller than 3 times optical penetration depth. The additional two terms in 6-parameters expression, C 3 and λ 3 , account for the field size effect for small diameter circular beams which reduces the light fluence rate at the central axis of the beam near the surface. Both 4 and 6 parameters expressions are simple to use and can accurately mimic the Monte Carlo simulations of light transport including at small depths close to air-tissue interface where diffusion theory fails. The accuracy and simplicity of these radial dependent analytic expressions would be extremely helpful to PDT in optimizing treatment plan involving various irradiance beam sizes that would improve the treatment outcome. Table 3. Summary of the 6 fitting parameters (b, λ 1 , λ 2 , C 2 , λ 3 , C 3 ) using Eq. (3) to the MC calculated φ/φ air for 1-cm diameter circular beam for optical properties (µ a , µ s ') and index of refraction, n 2 /n 1 = 1.4. MC calculated diffuse reflectance for broad beam, R d , is also listed.