2.1.

Simulation Overview

The simulation procedure to estimate ALA-iPDT outcomes consists of three major calculations: insertion conditions for cylindrical light diffusers and calculation of spatial distribution of light fluence rate at positions $(x,y,z)$, $\phi (x,y,z)$, calculation of spatial distribution of generated ${{}^1\mathrm{O}}_2$ concentration, $D_{\mathrm{SO}}(x,y,z)$, and estimation of cell death region, as shown in Fig. 1. A numerical three-dimensional (3-D) brain model with optical properties and PpIX concentrations was prepared. Cylindrical light diffusers were virtually inserted in the numerical 3-D brain model and the light propagation was calculated to obtain $\phi (x,y,z)$. Using $\phi (x,y,z)$, the $D_{\mathrm{SO}}(x,y,z)$ was calculated. By setting the threshold of the generated ${{}^1\mathrm{O}}_2$ concentration, $D_{\mathrm{SO}(\mathrm{th})}$, above the level at which cell death is induced, the treatment area in the numerical 3-D brain model was calculated. In this study, for quantitative comparison of several treatment conditions, treated volume (TV) and damaged volume (DV) values are defined as the volume of cell death region by ${{}^1\mathrm{O}}_2$ in tumor tissues and normal tissues. Also, tumor coverage (TC) is defined as the ratio of TV to total tumor volume.

Fig. 1

Schematic diagram of the calculation of ALA-iPDT outcomes. Here, ${\mu }_{\mathrm{a}}$, ${\mu }_{\mathrm{s}}$, and $g$ are the absorption coefficient, scattering coefficient, and anisotropy of scattering, respectively; $\phi (x,y,z)$ is the spatial distribution of fluence rate at position $(x,y,z)$; ${\phi }_{\mathrm{sur}}$, $T$, and $\lambda$ are the light fluences on the surfaces of inserted diffusers, irradiation time, and wavelength of light, respectively; $D_{\mathrm{SO}}(x,y,z)$ is the spatial distribution of generated ${{}^1\mathrm{O}}_2$ at position $(x,y,z)$; $\epsilon$, $C(x,y,z)$, $\beta$, and $\mathrm{\Phi }$ are the molar extinction coefficient, initial PS concentration at position $(x,y,z)$, photobleaching coefficient, and quantum yield for ${{}^1\mathrm{O}}_2$ generation, respectively.

2.2.

Numerical Tissue Model

To simulate ALA-iPDT outcomes, an anonymized MRI data set consisting of 156 two-dimensional slices (with resolution $512\times 512\text{pixels}$, pixel size $468\mu \mathrm{m}$, and slice thickness 1.25 mm) from a patient with a malignant brain tumor was used. The MRI slices were manually segmented into normal and tumor regions by a clinician using Synapse Vincent software (Fujifilm, Tokyo, Japan). The normal tissue regions around the tumor regions were assumed to be the white matter regions because most adult malignant brain tumors arise in white matter.³³ From the segmented MRI slices, a numerical 3-D voxel model having the segmented information was constructed³⁴ and optical property parameters were assigned to each voxel according to its tissue type. The absorption coefficient, ${\mu }_{\mathrm{a}}$ (${\mathrm{cm}}^{-1}$), scattering coefficient, ${\mu }_{\mathrm{s}}$ (${\mathrm{cm}}^{-1}$), and anisotropy factor of scattering, $g$, were set as $1.7{\mathrm{cm}}^{-1}$, $365{\mathrm{cm}}^{-1}$, and 0.9, respectively, for the tumor regions and $0.7{\mathrm{cm}}^{-1}$, $951{\mathrm{cm}}^{-1}$, and 0.9, respectively, for the white matter regions at the wavelength of 635 nm.³⁵ For iPDT outcome estimations, each voxel was divided into thirds in a direction perpendicular to the slice plane to give an isotropic voxel size of $468\times 468\times 417{\mu \mathrm{m}}^3$.

For detailed analysis of the influence of photobleaching and the initial PpIX concentration, a uniform brain tumor tissue model ($512\times 512\times 512\text{pixels}$, $100\times 100\times 100{\mu \mathrm{m}}^3$), which consisted of tumor only, was used. Although PpIX is induced heterogeneously in tissues, the PpIX distributions were assumed to be uniform in the tumor region of the model for simplicity. The ${\mu }_{\mathrm{a}}$, ${\mu }_{\mathrm{s}}$, and $g$ were set as $1.7{\mathrm{cm}}^{-1}$, $365{\mathrm{cm}}^{-1}$, and 0.9, respectively, for the wavelength of 635 nm.³⁵ The photobleaching coefficient and initial PpIX concentration values were varied in evaluation.

2.3.

Light Propagation Calculation

To obtain $\phi (x,y,z)$ for the constructed numerical 3-D models, the Monte Carlo-based method, which is regarded as the gold standard to calculate light propagation in nonhomogeneous scattering media, such as biotissues, was used.³⁴ Here, the cylindrical fiber diffuser was assumed to have a radiation section that uniformly emits photons. The outer diameter was 1.1 mm and the radiation length was 40 mm. To calculate light propagation from the inserted cylindrical light diffuser, a 3-D Monte Carlo program (mcxyz.c)³⁴ was customized. The launch points of photon packets were located inside the radiation section of the inserted cylindrical light diffuser according to Saccomandi et al.³⁶ As a laser source, a multiport laser system (ML 7710i, Modulight, Finland) with a wavelength of 635 nm and eight fiber ports was assumed.³⁷ Calculations of light propagation were performed for $1\times {10}^7\text{photons}$. Irradiation time and irradiation power density were changed according to the simulation conditions.

2.4.

Singlet Oxygen Generation Calculation

To calculate $D_{\mathrm{SO}}(x,y,z)$, we adopted the SOQY model. PDT works well at oxygenation levels over 1% dissolved oxygen saturation rate (DO), but the effect decreases rapidly at lower levels of oxygenation.³⁸ The oxygen concentrations of tissues during PDT did not fall below 1% DO.³⁹ The SOQYs of ALA-induced PpIX in aqueous solutions have been reported and are very similar in 1% and 20% DOs, where the 1% and 20% DOs are calculated as 13 and $260\mu \mathrm{M}$, respectively.³⁰ In the calculation, the same SOQY value were used for the normal and tumor regions.

Using the SOQY of PpIX, $\mathrm{\Phi }$, $D_{\mathrm{SO}}(x,y,z)$ can be expressed as

2.5.

Cell Death Region Calculation

In the ${{}^1\mathrm{O}}_2$ threshold models,^26,41 cell death is induced at the region where $D_{\mathrm{SO}}(x,y,z)$ exceeds $D_{\mathrm{SO}(\mathrm{th})}$. Normal tissues were assumed to be more resistant to ${{}^1\mathrm{O}}_2$ toxicity than tumor tissues because of superoxide dismutase, which exerts an important protective function against oxygen toxicity.⁴² However, to avoid underestimation of ${{}^1\mathrm{O}}_2$ toxicity in normal tissues, the $D_{\mathrm{SO}(\mathrm{th})}$ values for normal tissue damage were set to be equal to tumor tissue values. For estimation of the treatment region, the $D_{\mathrm{SO}(\mathrm{th})}$ value of 0.56 mM was used.²⁶

3.1.

Estimation of ALA-iPDT Outcomes Using Clinical Magnetic Resonance Imaging Data

To estimate the iPDT outcomes using clinical data, the TV and DV values were calculated with a cylindrical light diffuser inserted as shown in Fig. 2. Figure 2(a) shows the segmented tumor region and the insertion position. The total volume of the tumor region was $30.8{\mathrm{cm}}^3$. According to clinical studies,^12,13 $T$ and the light fluence at the surface of the diffuser, ${\phi }_{\mathrm{sur}}$, were set to 3600 s and $580{\mathrm{mW}/\mathrm{cm}}^2$ (linear density: $200\mathrm{mW}/\mathrm{cm}$), respectively. The $C_0(x,y,z)$ and $D_{\mathrm{SO}(\mathrm{th})}$ were set to $5.8\mu \mathrm{M}$ and 0.56 mM, respectively.^26,43 Figures 2(b)–2(d) show the $\phi (x,y,z)$ calculated by the Monte Carlo-based method, the $D_{\mathrm{SO}}(x,y,z)$ calculated by inputting the calculated $\phi (x,y,z)$, and the treatment region where $D_{\mathrm{SO}}(x,y,z)$ exceeds $D_{\mathrm{SO}(\mathrm{th})}$. TV and DV were calculated as 2.48 and $0{\mathrm{cm}}^3$, respectively. No damaged tissue in the normal tissue region was found under these irradiation conditions. When the PpIX uptake ratio of the tumor tissue region to white matter tissue (TN ratio) was 95,⁷ the DV was $0{\mathrm{cm}}^3$ even when the initial PpIX concentration in the tumor region was the reported maximum value⁴³ ($28.2\mu \mathrm{M}$). This result is consistent with a previous clinical study that showed almost no damage to white matter in ALA-PDT.⁷ However, when the TN ratio was set to 12, the value of the brain tissue adjacent to the tumor,⁷ damage in the white matter region occurred at higher PpIX initial concentrations ($>10.6\mu \mathrm{M}$) in the tumor region. This result shows that damage to normal tissue adjacent to the tumor is possible in clinical application.

Fig. 2

Calculation of treatment region under a set of irradiation treatment conditions. (a) Tumor regions segmented by a clinician and insertion position of a cylindrical light diffuser. Tumor regions are demarcated with white lines and the surface of the cylindrical light diffuser is shown in red. (b) Spatial distribution of light fluence, $\phi (x,y,z)$, calculated by a Monte Carlo-based algorithm, displayed as ${\log }_{10}[\phi (x,y,z)]$ ($\mathrm{mW}/{\mathrm{cm}}^2$). (c) Spatial distribution of generated ${{}^1\mathrm{O}}_2$, $D_{\mathrm{SO}}(x,y,z)$ (mM), calculated from $\phi (x,y,z)$. (d) Estimated treatment region (indicated in blue) setting an ${{}^1\mathrm{O}}_2$ concentration threshold to induce cell death at 0.56 mM.

3.2.

Evaluation of Irradiation Conditions

In ALA-iPDT, multiple cylindrical light diffusers are inserted depending on the tumor geometry.¹² To evaluate the dependence of the ALA-iPDT outcome on the insertion positions and ${\phi }_{\mathrm{sur}}$, TV and DV values were compared under several insertion conditions. The minimum distance between the diffusers (interfiber distance, $L$) has been determined as 0.9 cm to avoid a temperature increase above 42°C.¹² The distance $L$ was, therefore, chosen as 0.9 cm to concentrate photons between the diffusers as much as possible. Figure 3(a) shows the treated regions when $L$, $T$, and ${\mathrm{\Phi }}_{\mathrm{sur}}$ were 0.9 cm, 1 h, and $580{\mathrm{mW}/\mathrm{cm}}^2$ (linear density: $200\mathrm{mW}/\mathrm{cm}$), respectively. TV and TC were calculated as $17.4{\mathrm{cm}}^3$ and 0.56, respectively. The treated regions surrounding each diffuser overlapped. To improve the TC by avoiding region overlap, $L$ was increased. However, TC was calculated as 0.55 when $L=1.25\mathrm{cm}$ [Fig. 3(b)]. Setting ${\mathrm{\Phi }}_{\mathrm{sur}}$ as $2320{\mathrm{mW}/\mathrm{cm}}^2$ (linear density: $800\mathrm{mW}/\mathrm{cm}$), the TCs were 0.66 when $L=0.9\mathrm{cm}$ [Fig. 3(c)] and 0.71 when $L=1.25\mathrm{cm}$ [Fig. 3(d)]. The iPDT outcomes are strongly dependent on the relationship between insertion positions and ${\mathrm{\Phi }}_{\mathrm{sur}}$. These results demonstrate that both the insertion position and the ${\mathrm{\Phi }}_{\mathrm{sur}}$ should be designed by estimating the ALA-iPDT outcomes before treatment.

Fig. 3

Estimated treatment regions with eight cylindrical light diffusers with $D_{\mathrm{SO}(\mathrm{th})}$ set to 0.56 mM. The tumor region is inside the white line. The treatment region is shown in red and the centers of the inserted fibers are indicated by black points. [Interdiffuser distance (cm), fluence rate at the diffuser surface ($\mathrm{mW}/{\mathrm{cm}}^2$)] = (a) (0.9, 580), (b) (1.25, 580), (c) (0.9, 2320), and (d) (1.25, 2320).

The relationships between generated ${{}^1\mathrm{O}}_2$ concentration, $D_{\mathrm{SO}}$, and light fluence were investigated. As shown in Fig. 4, $D_{\mathrm{SO}}$ reaches a maximal value with increasing light fluence because of photobleaching. This indicates that an increase in light fluence does not always improve the ALA-iPDT outcomes. Moreover, in the ${{}^1\mathrm{O}}_2$-based models, $D_{\mathrm{SO}}$ is required to merely exceed the $D_{\mathrm{SO}(\mathrm{th})}$ value to provide treatment effects. When $D_{\mathrm{SO}(\mathrm{th})}$ was 0.4, 0.56, and 0.72 mM,²⁶ the required light fluences were 1.6, 2.3, and $3{\mathrm{J}/\mathrm{cm}}^2$, respectively. To avoid excessive irradiation, $D_{\mathrm{SO}(\mathrm{th})}$ should be considered when determining light fluence.

Fig. 4

Relationship between generated ${{}^1\mathrm{O}}_2$ concentration, $D_{\mathrm{SO}}$, and the light fluence for 0 to $100{\mathrm{J}/\mathrm{cm}}^2$.

3.3.

Analysis of Photobleaching Influence

To evaluate the photobleaching effect, $D_{\mathrm{SO}}$ was compared between the proposed model with and without the photobleaching process. The parameters $\beta$, $C_0$, which is PpIX initial concentration, $T$, and ${\mathrm{\Phi }}_{\mathrm{sur}}$ were $13.5{\mathrm{J}/\mathrm{cm}}^2$, $5.8\mu \mathrm{M}$, 3600 s, and $580{\mathrm{mW}/\mathrm{cm}}^2$, respectively. Figure 5 shows $D_{\mathrm{SO}}$ with and without photobleaching. When $D_{\mathrm{SO}(\mathrm{th})}$ was assumed as 7.9 mM,²⁶ the $D_{\mathrm{SO}}$ exceeded the threshold only in the model without photobleaching. When $D_{\mathrm{SO}(\mathrm{th})}$ was 0.4, 0.56, and 0.72 mM,²⁶ the $D_{\mathrm{SO}}$ exceeded the threshold to induce cell death in both models. Therefore, the parameter combinations of photobleaching and $D_{\mathrm{SO}(\mathrm{th})}$ effect the iPDT outcomes. Next, using $C_0=5.8\mu \mathrm{M}$, $T=\mathrm{3,600}\mathrm{s}$, ${\mathrm{\Phi }}_{\mathrm{sur}}=580{\mathrm{mW}/\mathrm{cm}}^2$, and $D_{\mathrm{SO}(\mathrm{th})}=0.56\mathrm{mM}$, TV was obtained to analyze the dependency of $C_0$, as shown in Fig. 6. When $\beta$ was 4.5, 13.5, and $33{\mathrm{J}/\mathrm{cm}}^2$,^27–29 the ALA-iPDT outcomes were induced when each $C_0$ was more than 0.4, 0.9, and $2.7\mu \mathrm{M}$. In the model without photobleaching, cell death was induced regardless of $C_0$. As $\beta$ becomes smaller, more $C_0$ is necessary for treatment. The optical properties of malignant brain tumor vary from patient to patient.³⁵ To estimate the possible ranges of the iPDT outcomes, TVs were calculated with the tissue optical properties, which maximized and minimized the light penetration depths. The maximum and minimum light penetration depths were defined using the standard deviation of the tissue optical properties. TV at the maximum light penetration depth was $20.0{\mathrm{cm}}^3$ when the ${\mu }_{\mathrm{a}}$, ${\mu }_{\mathrm{s}}$, and $g$ were set as $0.92{\mathrm{cm}}^{-1}$, $113{\mathrm{cm}}^{-1}$, and 0.9, respectively. TV at the minimum light penetration depth was $2.1{\mathrm{cm}}^3$ when the ${\mu }_{\mathrm{a}}$, ${\mu }_{\mathrm{s}}$, and $g$ were set as $2.5{\mathrm{cm}}^{-1}$, $617{\mathrm{cm}}^{-1}$, and 0.9, respectively. The TV difference was found since the reported standard deviations of optical properties are large.³⁵ A detailed analysis of the optical properties of the brain normal and tumor tissues is expected to increase the precision of the iPDT outcome estimations.

Fig. 5

Profiles of generated ${{}^1\mathrm{O}}_2$ concentration, $D_{\mathrm{SO}}$, from the center of the inserted diffuser for light fluence at the surface of the diffuser of 50, 75, 100, and $150{\mathrm{J}/\mathrm{cm}}^2$ (a) without and (b) with photobleaching when the coefficient was $13.5{\mathrm{J}/\mathrm{cm}}^2$.

Fig. 6

Relationship between treatment volume and PpIX initial concentration, $C_0$, of 0 to $10\mu \mathrm{M}$ in models with and without photobleaching for various values of photobleaching coefficients (4.5, 13.5, and $33{\mathrm{J}/\mathrm{cm}}^2$).

3.4.

Analysis of PS Parameter Sensitivities

The sensitivities of $C_0$ and $\beta$ were evaluated using the uniform living tissue model with optical properties of malignant brain tumors. Figure 7(a) shows the $D_{\mathrm{SO}}$ when $C_0$ was varied, assuming $\beta$ and $D_{\mathrm{SO}(\mathrm{th})}$ were $13.5{\mathrm{J}/\mathrm{cm}}^2$ and 0.56 mM, respectively. When $C_0$ was above $13.8\mu \mathrm{M}$, $1{\mathrm{J}/\mathrm{cm}}^2$ light irradiation was enough to induce cell death. When $C_0$ was 6.32 to $13.8\mu \mathrm{M}$, more light irradiation was necessary for treatment. Figure 7(b) shows the relationship between $D_{\mathrm{SO}}$ and $\beta$ when $C_0=5.8\mu \mathrm{M}$. The larger $\beta$ led to higher $D_{\mathrm{SO}}$ because of a smaller decrease of $C$. When $\beta$ was above $3.52{\mathrm{J}/\mathrm{cm}}^2$, $3{\mathrm{J}/\mathrm{cm}}^2$ light irradiation was enough for cell death. The $D_{\mathrm{SO}}$ varied depending on the light fluence and also on $C_0$ and $\beta$, which indicates that estimation of ALA-iPDT outcomes based only on light fluence, such as photon dose,^21,22 is less accurate and that ${{}^1\mathrm{O}}_2$ dose should also be considered.

Fig. 7

Generated ${{}^1\mathrm{O}}_2$ concentration, $D_{\mathrm{SO}}$, versus (a) PpIX initial concentration, $C_0$, from 0 to $30\mu \mathrm{M}$ ($\beta =13.5{\mathrm{J}/\mathrm{cm}}^2$) and (b) photobleaching coefficient, $\beta$, of 0 to $50{\mathrm{J}/\mathrm{cm}}^2$ ($C_0=5.8\mu \mathrm{M}$). ${{}^1\mathrm{O}}_2$ concentration threshold, $D_{\mathrm{SO}(\mathrm{th})}$, was 0.56 mM and light fluences were (i) $3{\mathrm{J}/\mathrm{cm}}^2$, (ii) $\text{2\hspace{0.17em}\hspace{0.17em}}{\mathrm{J}/\mathrm{cm}}^2$, and (iii) $\text{1\hspace{0.17em}\hspace{0.17em}}{\mathrm{J}/\mathrm{cm}}^2$.