Open Access
Add to BibSonomy
Abstract
PDT-SPACE is an open-source software tool that automates interstitial photodynamic therapy treatment planning by providing patient-specific placement of light sources to destroy a tumor while minimizing healthy tissue damage. This work extends PDT-SPACE in two ways. The first enhancement allows specification of clinical access constraints on light source insertion to avoid penetrating critical structures and to minimize surgical complexity. Constraining fiber access to a single burr hole of adequate size increases healthy tissue damage by 10%. The second enhancement generates an initial placement of light sources as a starting point for refinement, rather than requiring entry of a starting solution by the clinician. This feature improves productivity and also leads to solutions with 4.5% less healthy tissue damage. The two features are used in concert to perform simulations of various surgery options of virtual glioblastoma multiforme brain tumors.
© 2023 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. Introduction
Photodynamic therapy (PDT), or the use of light activated drugs, called photosensitizers (PS), is making continuous inroads for oncological indications which are proximal to tissue surfaces as the case of skin [1], esophageal [2] or bladder cancer [3,4]. In these situations, the photosensitizer is systemically administered to the patient and a large tissue surface area, comprising the target tumor and some surrounding normal is irradiated by the activation light resulting in a one-dimensional photon density gradient as function of depth into the tissue.
The photon gradient is governed by the PDT activation light’s effective attenuation coefficient, and the attenuation to $1/e$ ($\sim 37\%$) of the surface irradiance is comparable to the thickness of the targeted tissue. For skin, esophagus and bladder the effective attenuation coefficients are in the 0.3 to 2 $mm^{-1}$ ranges [5]. The shallow treatment depth and the generally homogeneous illumination of the target by a wide field source allow the use of empirically derived treatment protocols based on the administered photosensitizer concentration and radiant exposure [$J \cdot cm^{-2}$] or photon density [$hv \cdot cm^{-2}$].
While surgery remains the dominant primary cancer therapy not all tumors are eligible for surgery, in particular brain tumors and prostate cancers, due to the high probability for surgery related morbidity and mortality. Additionally, less invasive localized interstitial therapies are being thought to improve patient quality of life. The clinical development of interstitial PDT (iPDT) for solid tumors is ongoing, but remains lacking in clinical acceptance. Ongoing and completed phase II or III clinical trials are few, focusing on the prostate, with individual clinical trials ongoing for hepatocellular carcinoma, head and neck tumors, locally advanced lung cancers and glioma [6].
The slow progress of iPDT translation into the clinic is partially due to the difficulty of light delivery planning, with additional impediments being hypoxia in tumors and the limited tumor selectivity of classic photosensitizers [7–9]. This work focuses on light delivery planning, which is complex because for a given tumor geometry multiple PDT activation photon sources, here optical fibers with isotropic or cylindrical diffuser tips, need to be placed at predetermined positions. Each optical fiber generates a two-dimensional photon density gradient given its emission profile and the tumor and host tissues’ effective attenuation at the PDT wavelength. The photon gradients of multiple implanted fiber sources overlap to generate sufficient photon density within the tumor target. Hence, iPDT treatment planning and execution are significantly more challenging for the surgeon than trans-surface illumination.
Currently, guidelines for photon source fiber placement are often based on empirical treatment protocols with fixed inter-fiber distances [10–13]. Based on these guidelines, confining the treatment to the targeted tumor and avoiding over-treatment of the surrounding host tissue remains challenging. Previous work from this group presented interstitial photon source placement optimization protocols [14–17] based on Monte Carlo [18–20] simulations of light propagation in a 3D space combined with photosensitizer uptake and PDT responsivity in the tumor and host tissue values from the literature.
The previous work realized automatic PDT treatment planning with a software tool called PDT-SPACE that selects both the light source locations and their optical powers. PDT-SPACE demonstrated improved complete coverage of the target tumor with the required minimum PDT-dose, given by the photosensitizer concentration and photon density product, while limiting the volume of host tissue where the maximum permissible PDT-dose is exceeded resulting in collateral damage.
Clinical translation of this treatment planning tool is currently hindered by two significant shortcomings. First, the current approach does not consider surgically permissible access to the target volume, but instead assumes ubiquitous surgical access. Hence, the potential PDT-SPACE treatment planning solutions may not be surgically implementable due to a high risk of mechanical damage and associated morbidity to the patient. An example of this limitation includes producing unacceptable treatment plans which require crossing eloquent areas of the brain, such as the speech and motor cortex, with the needles or catheters holding the optical fibers. Hence, clinical iPDT treatment planning solutions need to consider these surgical limitations and be limited to fiber access directions which do not cross these areas nor risk damage to important blood vessels along their path. The second shortcoming is that currently the surgeon needs to provide an initial guess for the source fiber placement. This necessitates a skilled user with a prior empirical understanding of the appropriate number of source fibers, potentially limiting the attainable treatment planning optimization as well as increasing the manual effort to create a plan. In this work we enhance PDT-SPACE to remove both of these limitations: clinicians can specify surgical feasibility constraints that limit the solution space, and PDT-SPACE can generate high-quality initial solutions with appropriate numbers of sources.
This work is organized as follows. The next section lays out necessary background on iPDT planning in general and PDT-SPACE in particular. Section 3 outlines the test cases used to evaluate the algorithms, while Section 4 details the algorithmic enhancements. Section 5 evaluates the effectiveness of the new features and uses them to explore planning trade-offs for the test cases. The significance of these results are discussed in Section 6 and Section 7 summarizes conclusions.
2. Background
iPDT delivered by the placement of the light emitting optical fibers [21] in a predetermined grid pattern covering the tumor target volume can be optimized when the prescript photon energy is adjusted by the local photosensitizer concentration. The photosensitizer concentration can be quantified using the emitting fibers, their position and known tumor and host tissue response.
Current clinical interstitial PDT treatment planning is commonly based on placing optical fibers in a predetermined grid pattern independent of a PDT tissue response model, as in the case of TooKad mediated vascular PDT for prostate cancer [22], and using predetermined administered photosensitizer dose [$mg/kg$] and optical energy density [$J/cm$-diffuser-length]. The iPDT treatment planning approach based on administered photosensitizers and light energy can be largely agnostic about the tissue optical properties as fiber placement is dense enough compared to the effective activation light attenuation in the target tissue. Light emitters also acting as detectors for the PDT activation light and the photosensitizer’s fluorescence, as a quantifiable measure of the tissue response measure, enables monitoring of the treatment progression. In the latter approach, the photosensitizer concentration, are quantified over 1-2 $cm^{3}$ large tissue volumes. These approaches have been demonstrated in the prostate [23] and the brain when targeting glioma [24]. Adjusting the activation energy delivery at each light emitter based on the local quantification of the total photosensitizer concentration and its bleaching necessitates the knowledge of the tissue optical properties on similar length scales [25,26].
iPDT treatment plans based on predetermined grid pattern source placement can provide visualization of the iPDT treatment outcome based on the forward calculation of the photon- or iPDT- dose simulations by either finite element solution of the diffusion equation [27,28] or by Monte Carlo simulation [19,20,29–31].
The major advantage of grid patterns for particular indications such as the prostate or the brain is the use of alignment templates placed at positions avoiding critical structures during fiber source placement.
Treatment planning taking tumor geometry and PDT responsiveness of all tissues in the treatment volume into consideration must be based on a numerical or analytical dose description. Models of varying complexity have been proposed. Models are based either only on the photon propagation, [27,29–31], explicitly calculating the spatial distribution of the cytotoxic singlet oxygen [32–36] using photochemical and photobiological considerations, or models based on known tissue damage responses such as the Photodynamic Threshold Model [37–41].
The more complex treatment planning approaches require knowledge of the photosensitizer accumulation, the tissue’s responsiveness, and the anatomy of the target and host tissues as utilized in PDT-SPACE [14–17]. PDT-SPACE utilizes Monte Carlo simulations of photon propagation through the FullMonte software tool and optimizes for each fiber source’s position, length, and optical power to achieve a desired target PDT dose throughout the target while minimizing the risk and morbidity to the adjacent host tissues, or organs at risk (OARs). In this approach a 3D tetrahedral mesh is created from the clinical images which act as the virtual input mesh for PDT-SPACE. In this work we focus on treatment planning of GBM tumors, but PDT-SPACE is a general tool that can be used for PDT planning for different indications. However, to date, placement of the light sources and the transparent catheter guiding them does not consider clinical or surgical access limitations, restricting the utility of these planning tools.
Access to different locations in the human brain is challenging for PDT and other interventions such as Deep Brain Stimulation [42]. Clinical mapping of eloquent brain areas using MRI and electrocorticography [43] is driven by Epilepsy [44] and glioma [45] surgery. These studies have identified areas not to be traversed by the catheters with the needles or catheters carrying the optical fibers, restricting the clinically implementable solutions for PDT-SPACE. Hence, access points for burr holes in the skull need to be determined.
PDT-SPACE Optimization Capabilities and Approach. The PDT-SPACE software tool consists of two major modules: a power optimization module that optimizes for source power allocations with a given source position allocation, and a position optimization module that optimizes for source positions.
Power optimization makes use of convex optimization [14] by defining a linear program minimizing the weighted sum of the overdosages of OARs and under-coverage of tumor. The weight for tumor underdosage is iteratively adjusted to achieve 98% tumor dose coverage. 98% cytoreduction of a GBM is required to provide the patient with a meaningful survival advantage; the median survival was 13 months for a cytoreduction of 98% or more, but only 4.2 months for less than 98% tumour elimination [46]. PDT-SPACE uses FullMonte [18] to perform Monte Carlo simulations of the light propagation within the treatment planning volume; simulating once for each light source with unit power to obtain a fluence dosage matrix, with rows containing photon dose received at each volume element and columns representing each source. Each column of the matrix can then be multiplied by the assigned power to the respective source and summed by the law of superposition to arrive at the final photon dose resulted from the power allocation. This matrix is used in the linear program definition to optimize for the respective powers assigned to each light source. The linear program also includes a user-specified constraint on the maximum power emitted by each source to avoid overheating tissues.
Position optimization uses simulated annealing (SA) to iteratively refine the light source positions [16]. Starting from a given reasonable initial source placement, the program computes the placement cost by performing power optimization and obtaining the resulting over- and underdosages for the current placement. PDT-SPACE then performs iterative refinements of source positions via a reinforcement learning algorithm that repeatedly chooses from among the following perturbation moves:
- 1. Random Moves: randomly relocates a source.
- 2. Displace Moves: randomly displaces a source without changing its orientation.
- 3. Force-based Directed Moves: each OAR volume element and tumor volume element generates a repelling or attracting force, respectively, on a source. The source is then moved accordingly.
Not all proposed moves are accepted; the placement cost change produced by each is evaluated and its acceptance (update of the light source positions) is determined using the simulated annealing acceptance function. Please see [16] for further details.
3. Simulation methodology and evaluation metrics
The simulation setups are based on those used in the team’s prior works [15,16], with slight modifications to test the new features.
3.1 Geometry test set
The brain geometry used is a segmented 3D tetrahedral brain mesh taken from the Colin27 [47,48] brain atlas. The tumor shapes were approximated by elongated ovoid or several connected spheres to represent clinical glioblastoma multiforme (GBM) images obtained from the cancer imaging archive [49]. Note that these tumor shapes are virtual (constructed), and are intended to test the functionalities of PDT-SPACE across a range of tumor sizes and shapes.
The tetrahedral meshes (T1 $\sim$ T9) each consist of a total of 423K tetrahedra representing five regions: skull, cerebrospinal fluid (CSF), gray matter, white matter, and tumor. Figure 1 shows the front view of the tumor models visualized in Paraview. In addition, Table 1 lists geometric information about each of the 9 virtual tumors, along with the number of light sources used to treat it when an initial placement of light sources is given to PDT-SPACE. The light sources are modeled as infinitely thin line sources that emit light isotropically along a finite length that is chosen by PDT-SPACE to be less than or equal to a user-specified maximum length. This is a reasonable model of the thin (1 mm diameter) cylindrical diffusers generally used in PDT; these diffusers are commercially available in 5 mm length increments [50,51] and other techniques exist to realize finer gradations in the diffuser length [52,53].
Fig. 1. Paraview visualizations of all tumor models viewed from the front of the brain. All meshes are 175 $mm$ in height.
Table 1. Geometry information for each tumor model. The OBBs are oriented such that $l>w>h$. The number of sources listed were used to constrain the source numbers in initial source placement generation unless otherwise specified.
Table 2. Optical properties for brain tissues (tumor, gray matter, white matter, skull, and cerebrospinal fluid) used in this work. Table obtained from [16].
Table 3. Relative tissue dose threshold values with respect to the minimum tumor dose threshold. $d_{min}$ refers to the minimum dosage that must be achieved by a tissue, and is set to 0 for all non-tumor tissues; $d_{max}$ refers to the maximum dosage a tissue can receive without overdosing, set to $\infty$ for tumor. Note the gray and white matter $d_{min}$ thresholds here are $0.1\times$ the actual tissue death thresholds. Table obtained from [16].
Table 4. Geometric averages across 9 tumor models of the average values and standard deviations of runtime and OAR overdosage taken over 10 PDT-SPACE runs with and without injection constraints.
Table 5. Results for generated initial source placement. Geometric averages across 9 tumor models of the average values of OAR overdosage, runtime, and number of sources taken over 10 PDT-SPACE runs.
Table 6. Results of additional optimizations. This table reports the geometric average across 9 tumor models.
3.2 Photosensitizer, optical properties and dose thresholds
The optical properties (Table 2) of brain tissues used in the simulations are at a PS activation wavelength of 635 $nm$ [54–57]. 5-ALA (5-aminolevulinic acid) induced PpIX mediated PDT is modeled here [58].
The tissue PDT-threshold doses used in the tool, defined by the threshold light dose received leading to cell death for a given tissue type, are determined by the relative PS uptake ratio of a tissue to the tumor along with the number of photons absorbed by the photosensitizer per unit volume, assuming a uniform photosensitizer distribution per tissue. The resulting fluence threshold doses, i.e. energy per unit area passing through the surface of the tetrahedrons representing the problem geometry, are detailed in Table 3. A detailed description and computation of tissue death thresholds is in [14,41]. To evaluate OAR overdose at a given level of the tumor treatment the relative dose threshold values per tissue are the key metric; accordingly the threshold doses in Table 3 are normalized such that the minimum tumor dose threshold is 1 $J/mm^2$. The thresholds of gray and white matters have a 10$\times$ guard-band vs. the tissue death threshold; this leads to a forced increase of the over-treatment penalty during power allocation optimizations to maximize safety or minimizing morbidity. It also means that the overdosed gray and white matter volumes demonstrated in Section 5 are all evaluated with a dose threshold that is $0.1\times$ the tissue death thresholds, so it is a very conservative overestimate of the actual OAR damage.
3.3 Computational platform
The algorithms were implemented in C++ and tested on Ubuntu 16.04 machines with 128 $GB$ of RAM and 12-core 2.2GHz Xeon E5-2650v4 CPUs; each core can run 2 threads simultaneously for a total of 24 hyper-threads.
The source code of the PDT-SPACE tool is available on Gitlab at http://FullMonte.org/.
3.4 Evaluation metrics
The $v_\alpha$ value of a tissue refers to the volume of the tissue receiving $\alpha \%$ of its dose threshold. Results used for comparisons in Section 5 are the runtimes along with the gray matter and white matter $v_{100}$ values, while targeting the tumor $v_{100}$ value to be $98\%$ of the tumor volume. For each feature, we will compare 1) average runtime and overdosages over 10 independent runs for each tumor model, and 2) geometric averages of the mean and standard deviations over the 10 runs across all tumor models.
4. Features and algorithms
The new features this work introduces into PDT-SPACE fall into two categories. First, we extend PDT-SPACE to allow specification of constraints on the tissue insertion location of all light sources. This can address surgical complexity by, for example, limiting the number and size of burr holes through the skull, and can ensure that key structures, such as eloquent areas of the brain, are not mechanically penetrated by the inserted sources. Second, PDT-SPACE now optimizes light source placement with less manual effort and higher quality. Prior versions of PDT-SPACE required an initial solution of source locations which was then iteratively refined, this required the clinician to estimate the required number of sources and their rough locations to begin optimization, increasing manual effort and the expertise required. PDT-SPACE now includes automatic algorithms to generate an intelligent initial source placement and number of sources. The iterative improvement algorithms that refine this initial solution have also been enhanced to make better use of the light sources and to eliminate any light sources not providing significant benefit in order to simplify surgery and PDT delivery.
4.1 Injection constraint
The addition of an injection constraint provides an option to limit the valid locations of the insertion points of light sources. The clinician is able to provide the valid injection locations by specifying a circle on the surface of the input mesh; the program will then generate only those source placements that can be inserted through the circle. Figure 2 shows an example of a brain tumor model with a sample source placement and injection constraint on the surface of the skull. In this work we implement injection constraints that are circular, as these are a good model for the burr holes of cranial surgery. However, if other shapes (e.g. oval or square) are desired the formulation could be extended to support these additional constraint shapes.
Fig. 2. Example of injection constraint on a brain tumor model with a parallel placement of 19 sources satisfying the constraint.
4.1.1 Injection validation
The algorithm determining whether a valid injection direction consists of two steps: 1) compute the intersection of the extended line form the diffuser, representing the optical fiber shaft or catheter, with the plane where the injection constraint circle lies, and 2) compute the distance of the intersection from the center of the circle and compare with the circle radius. See Fig. 3.
The intersection of a line and a plane can be computed geometrically. Assume the line is defined by the set of $\textbf {p}_l \in \mathbb {R}^3$ satisfying Eq. (1), where $\textbf {p}_\textbf {A}\in \mathbb {R}^3$ is a point on the line, $k\in \mathbb {R}$ is any real number and $\vec{\mathbf {l}}\in \mathbb {R}^3$ is the direction vector, and the plane is defined by the set of $\textbf {p}_p \in \mathbb {R}^3$ satisfying Eq. (2), where $\textbf {p}_\textbf {B} \in \mathbb {R}^3$ is a point on the plane, and $\vec{\mathbf {n}}\in \mathbb {R}^3$ is the normal vector of the plane. The intersection point $\textbf {p}_i$ (setting $\textbf {p}_l = \textbf {p}_p$) can then be calculated from Eq. (3).
4.1.2 Generation of valid random moves
With the addition of an injection constraint, the PDT-SPACE iterative improvement algorithm will often generate invalid random source moves that do not satisfy the injection constraint, wasting runtime. The current PDT-SPACE Simulated Annealing engine contains two types of random moves, complete random move and random displacement move; both are modified to generate only moves that satisfy the injection constraint.
Random Moves To generate a complete random move, the program first randomly picks a valid source length satisfying any source length constraints specified, then chooses a point in space to place the first source endpoint (proximal end of the diffuser), and finally randomly rotates in space to place the second source endpoint (distal end of the diffuser).
The modification added is during the final rotation step. With a fixed proximal endpoint and a fixed diffuser length, the rotation is generated by choosing a random point in the valid circle for injections, then extending from this point to the first endpoint and placing the distal endpoint at the appropriate distance from the proximal endpoint.
Displace Moves. To randomly displace a source without changing its orientation, the program first randomly picks a displacement direction perpendicular to the source, then picks a random distance to displace the source. With an injection constraint, after determining the displacement direction $\hat {\textbf {dir}}$, the program finds and constrains the maximum distance to displace without moving out of the valid injection circle. See Fig. 4 for a visualization of the problem definition. The target is to find the maximum legal displacement distance $d$, given $\textbf {c}$ as the center of the circle, $\hat {\textbf {n}}$ as the normal unit vector, $r$ as the circle radius, $\hat {\textbf {dir}}$, and the current source position denoted by its endpoints.
The solution consists of four steps: 1) find the intersection $\textbf {p}_\textbf {i}$ from Eq. (3); 2) find the projection of $\hat {\textbf {dir}}$ onto the plane and normalize as $\hat {\textbf {d}_\textbf {1}}$; 3) find the distance $l_1$ from Eq. (8); and 4) calculate $d$, the perpendicular distance from $\textbf {p}_\textbf {edge}$ to the line, using trigonometry from Eq. (9).
4.2 Initial source placement generation
An initial source placement generation engine was added to provide a more automated usage option that removes the need for the manual entry; it can operate with or without an injection constraint specified.
The initial source placements used in prior PDT-SPACE testing [16], as well as those used in clinical studies and PDT plannings [11–13], were created following the same heuristic approach. Sources are to be placed in parallel, at least 8 $mm$ $\sim$ 10 $mm$ from the tumor boundaries, and at least 10 $mm$ apart from each other. The distances were obtained based on the effective light penetrations for the specific clinical indication used for testing [12], and will be used throughout this work as an approximation of the treatment planning heuristic.
The source placements used in prior works provide a baseline comparisons in which the sources are injected at the most proximal exterior surface to the target volume. As the exact heuristics are difficult to be reproduced by software the following implementations approximate this placement heuristic. The descriptions are separated by whether the source injections are constrained, since it affects the direction where the parallel diffusers can be inserted.
4.2.1 With injection constraint
A specific injection constraint implicitly limits the solid angle for the source placement direction. The program first draws a reference line $l_{ref}$ from the center of the injection circle to the center of the tumor (center of the axis-aligned bounding box of the tumor shape), as shown in Fig. 5(a). A reference coordinate system is created with the z-axis aligned with this reference direction, denoted by $z_{ref}$. It then displaces the reference line within the coordinate system and completes a parallel source placement, as shown in Fig. 5(b). For each cross-point on the grid, the program searches in the $z_{ref}$ direction to determine if a source can be placed there - a source will only be placed if it satisfies the injection constraint, passes through the tumor volume, and the two endpoints can be placed at a certain defined minimum distance from the tumor boundary. The pseudo-code for placing a single source along a provided line is given in Algorithm 1, where it takes as input a line with the distal and proximal endpoints oriented in the $z_{ref}$ direction. The distance between sources, $d_p$, is adjusted based on the available number of sources; the program will run a binary search algorithm until it reached the user defined sources number. The pseudo-code for the overall initial source placement procedure is shown in Algorithm 2
Fig. 5. Example visualization of initial source placement generation with an injection constraint on T1. a) Definition of $l_{ref}$ and $z_{ref}$. b) A 2D coordinate system in a plane perpendicular to $z_{ref}$. sources are placed at each grid cross point along the $z_{ref}$ direction.
Algorithm 1. Function to place a source along a specified line.
Algorithm 2. Generating initial source placement with injection constraint
4.2.2 No injection constraint
If no injection constraint is provided, it is assumed that the light sources and their optical fibers do not traverse critical tissue regions and injection is permitted at all surface points. In this case, minimizing the number of light sources becomes the most important factor to be considered; therefore ideally light sources will be placed along the tumor’s long axis. In this implementation, the tumor long axis is approximated by computing the Oriented Bounding Box (OBB) of the tumor shape.
The Oriented Bounding Box Problem The OBB of a finite set of 3-dimensional points is defined by the smallest rectangular volume circumscribing all of the points [59]. Figure 6 shows an OBB of an ellipsoid object. The authors of [60] surveyed methods of solving this problem and compared their computational complexity, implementation difficulty, and OBB quality. The Principle Component Analysis (PCA) [61] approach produces a high quality, but generally not optimal, solution with a low time complexity of $O$($N$), where $N$ is the number of 3D points in the volume to be circumscribed. Since the OBB in this case is only a starting point for initial source placement generation, a fast and reasonable-quality solution is more appropriate than a slow and optimal one, making PCA a good choice.
Algorithm Implementation The PCA approach computes the OBB from a set of points from the following steps: 1) compute the 3x3 covariance matrix of the positions for the data points to be considered, 2) diagonalize the matrix to find its Eigenvectors to form a basis to which the bounding box will be aligned, and lastly 3) find the boundaries of the OBB projected onto the Eigen basis. Diagonalization of the covariance matrix is done by iterative rotations of the matrix using Jacobi Transformations of symmetric matrices [62]. The source placement is similar to that as described in Section 4.2.1 with ($l_{ref})$ and $z_{ref}$ aligned with the longest OBB edge, rather than computed from an injection constraint, as shown in Fig. 7. Figure 8 shows the final OBB and a sample placement result from one of the tumor models used for testing (a detailed description of the testing methodology is given in Section 3.1).
Fig. 7. Example visualization of initial source placement generation with no injection constraint on T1.
Fig. 8. OBB and sample initial source placement result for T1. The longest and shortest dimensions of the OBB define the major and minor axes of the tumor, respectively.
4.3 Additional optimizations
4.3.1 Initial placement to final result correspondence
The initial source placement algorithm is fairly fast, taking from 30 seconds to 10 minutes, depending on the number of sources and tumor size. Iterative refinement of the source placement, on the other hand, typically takes several hours to complete. Consequently, if we could quickly determine that some initial placements would lead to a higher quality result after iterative refinement, it would be advantageous to generate multiple initial solutions and pick the best one for refinement.
Figure 9 shows the correlation between the OAR damage when we optimize the power allocation of the initial source placement to achieve 98% tumor coverage (x-axis) vs. the OAR damage when we optimize the power allocation on the final source placement after iterative refinement (y-axis). The figure combines the results for the 9 virtual tumors T1 to T9, described in Section 3.1. Observing the figure, some of the placements do show a strong correlation between initial placement and final placement, while some other placements are affected by randomization during optimization but still showing weaker correlation. Overall, performing 3 initial placements and choosing the placement with the lowest OAR damage after power allocation as input to the simulated annealing optimization process improved result quality without significantly impacting run time; we call this the Multi-Init approach and test its utility in Section 5.3.
Fig. 9. Overall (a) and tumor model specific (b) correspondence graphs between initial source placement overdosage ($mm^3$) (horizontal axis) and final placement overdosage ($mm^3$) (vertical axis). Note this graph consists of the combination of all points obtained for T1$\sim$T9.
4.3.2 Perturbation moves
The performance of the simulated annealing refinement algorithm is strongly impacted by the moves, or perturbations, that it proposes. In addition to the moves as described in [16], two moves were added to optimize for source utilization.
Source Removal Move The source removal move eliminates the source with the least power from further optimization steps as a 0 W emitting source may block the program from placing more useful sources near it. Additionally this leads to a simpler surgery when PDT dose refinement is not impacted.
Zero-Powered Source Relocation Move This move identifies all sources unused in the current placement (assigned zero powers), and then partitions the underdosed tumor volumes with a naive partition strategy and tries to place parallel sources into the partitions. This move, especially when combined with source removal moves, improves the source utilization by removing less useful sources and ensuring no largely under-treated volume remains.
The current partition strategy projects the centroids of the underdosed tetrahedra onto the injection constraint plane, or, if there is no such constraint, onto the plane perpendicular to the injection direction of the initial source placement. This produces a set of 2-dimensional points on a single plane. PDT-SPACE then loops through the set and merges the points that are close to each other by computing their joint centroid, the merging distance is heuristically set according to the optical properties governing light penetration [12]. The final partitions are defined by the points remaining after merging is completed; a weight factor is accumulated at each partition as the number of points merged into it.
After obtaining the partitions, the program sorts the weights and places the zero-powered sources starting from the partitions with the largest weights. Sources are placed along the line perpendicular to the projection plane and passing through the partition points with Algorithm 1.
5. Results
5.1 Injection constraint
The test-cases used to evaluate the quality of result obtained when the injection site is constrained vs. unconstrained were generated based on the manually generated initial source placements used in prior PDT-SPACE studies [16]. These manual initial placements are created aligned to the most proximal outer surface, with a clinical heuristic of 10 $mm$ center to center spacing of the optical diffusers in a triangular pattern [10]. This spacing was shown to increase the fraction of glioma patients to survive more than 30 months versus the standard 13-15 month survival for surgery and standard chemo and radiation therapy [63,64]. The number of light sources for each tumor is given by Table 1. Injection constraints were created as circles on the skull that allow these initial light source placements, with about 5 $mm$ extra radius to allow source perturbations. Figure 10 shows the manual initial source placement and the injection constraint used for T1; other tumor test cases were created in a similar manner.
Fig. 10. The T1 tumor (red volume), manual initial source placement (white lines inside the tumor) and injection constraint (white circle) used in testing.
Figure 11 shows a comparison of computational runtimes and gray or white matter overdosage, with and without an injection constraint. The values plotted for each tumor are the averages taken from 10 PDT-SPACE runs with different random seeds (starting points) in order to reflect optimization result variation. Table 4 lists the geometric averages across the 9 brain tumor models of the per-tumor average and standard deviation results. Constraining the injection point reduces the average result quality by a modest amount (10.4% increase in total overdosage) and the standard deviation is slightly reduced, which is expected given the constrained solution space. The per tumor results of Fig. 12 show that for most tumors specification of an injection constraint causes only small changes in white and gray matter damage; tumor T9 is an outlier where OAR damage significantly increases. In fact on average over tumors T1 - T8 an injection constraint produces almost no change in OAR overdose (only a 0.2% decrease), while for T9 OAR overdose more than doubles. Without an injection constraint, PDT-SPACE produces a light source placement for T9 where one of the four light sources is oriented in a direction almost perpendicular to the other three. This solution is surgically complex (it would require two burr holes) and is far from satisfying the injection constraint, which accounts for the large OAR damage change when an injection constraint is specified. An in-depth look of the T9 case is included in Supplementary Material section 2.
Fig. 11. Runtime ($h$) and OAR overdosage ($mm^3$) results with and without injection constraints for 9 tumor cases. Each result for a tumor is the average of 10 PDT-SPACE runs with different random seeds.
Fig. 12. OAR overdosage ($mm^3$), runtime ($h$) and number of sources of manually generated initial source placement (MANUAL), automatically generated initial placement using the same number of sources and source lengths as manual initial placement (AUTO-D), and automatically generated initial placement with unconstrained number of sources (AUTO-U). a) shows results with no injection constraint and b) shows results with an injection constraint. Each result for a tumor is the average of 10 PDT-SPACE runs with different random seeds.
5.2 Initial source placement generation
In this section we evaluate the quality of results PDT-SPACE obtains when it automatically generates an initial light source placement vs. when a manual initial placement is provided. The results from Section 5.1, denoted by Manual in Fig. 12 and Table 5, are compared with simulation results from running the same experimental settings – including the same injection constraints, if applicable – but with initial placements automatically generated by the algorithms described in Section 4.2. Two cases of automatic initial placement generation are evaluated. AUTO-D chooses initial source locations, but uses the same maximum number of sources and maximum source lengths as the Manual solution. AUTO-U is less constrained, and can change the maximum number of sources and their maximum lengths. The computing time to automatically create an initial placement is small, typically taking between 30 and 60 seconds, with the longest time over any test case being less than 10 minutes.
We first experimented by generating initial source placements with the same number of sources and source lengths used in manual source placements to see if the automatic process is able to preserve result quality, denoted by Auto-D. The average OAR damage was reduced for the cases with an injection constraint (4.5% decrease overall), and increased in the case of no injection constraint (17.7% increase overall). We believe this reduction in the treatment plan quality is due to over constraining PDT-SPACE; both the number of sources and their lengths were created by fine-tuning an initial manual placement, and rediscovering a high quality source placement that matches these lengths is challenging when there is no injection constraint to guide the search. To justify this explanation, a second experiment was performed by letting the program generate the initial source placements with a large number of sources and unconstrained lengths, denoted by Auto-U. Note, in order for this experiment to arrive at valid results without impacting runtime, we added a combination of the source removal and relocation moves described in Section 4.3.2 to optimize for source numbers, as well as added a hard constraint on the minimum source spacing $d_p$ to strictly replicate the original clinical heuristic described in Section 4.2. Algorithm 2 was modified such that, instead of adjusting $d_p$ to use all provided sources, the program decreases the target number of sources iteratively if minimum $d_p$ is not met and changes the reference $x-y$ coordinate system. The final targeted number of sources is then set to the maximum generated number of sources obtained in this iterative process.
The results from Auto-U show that PDT-SPACE can consistently reduce OAR damage (overall decrease of 17.8% without an injection constraint and 27.3% with an injection constraint) if it can personalize the treatment with a larger number of sources and can choose their lengths, at an average cost of 2 to 3 more sources and 20% to 30% more runtime. Overall the results show that a skilled user who can accurately estimate the number of sources and their lengths, may improve result quality and save runtime by providing an initial source placement. For use cases where the user is less experienced or where what-if analysis of quality vs. source count are desired, the automatic creation of the initial solution by PDT-SPACE is more appropriate.
5.3 Additional optimizations
The additional optimizations described in Section 4.3 were tested with an injection constraint and automatically generated initial source placements. The respective simulation results obtained in Section 5.2, Auto-D with injection constraint, are used as baseline comparisons in this section. As can be seen from Table 6, the resulting average $v_{100}$ values are only marginally affected by the optimization attempts. Generating 3 random source placements and taking the one with the least overdosage (denoted by Multi-Init) slightly reduced the average OAR damage and effectively decreased the standard deviation caused by randomness. Neither the addition of the removal move (SA-Removal) or the combination of removal and source relocation moves (SA-Relocate) significantly reduced overall OAR damage from the baseline, but they reduced the number of sources used to achieve the same results, simplifying the iPDT delivery.
As a result, both Multi-Init and SA-Relocate are enabled in the final PDT-SPACE implementation.
5.4 Investigations on injection constraint direction and diameter
The introduction of the new PDT-SPACE features allows in-silico studies to evaluate the impact of various treatment plan options on plan quality and surgical complexity. An example study on the effect of the position and size of injection constraints was performed on a large tumor model with an obvious longest axis T1, and a small tumor model with a more cubic OBB T4.
This experiment examines the results generated from injection constraints along the major and minor axes (along the longest and shortest sides of the OBB), for both the case of the smallest possible injection constraint and the case of a large injection constraint, for both tumor models. PDT-SPACE was first used to determine the smallest possible injection constraint required to achieve 98% tumor volume damage for each case. As can be seen from Table 7, T1 requires a 10 $mm$ larger injection constraint to obtain a 98% tumor coverage on the minor axis than the major axis due to the difference in area of the tumor volume projection onto the injection plane. Treating a tumor with probes inserted along its minor axis also significantly increases the number of sources required, by 30% to 60% in most cases. This increased source count does lead to lower OAR damage; for example T1 required 7 more sources on average and reduced damaged OAR volume by 35% when treated on the minor axis instead of the major axis. Table 7 also shows that larger injection constraints lead to PDT-SPACE using more sources (as it has more freedom to place them) and achieving a plan with significantly lower damaged OAR volume (over 50% reduction in all cases). An example visualization of final source placements obtained for T1 on the major axis and different injection diameters is included in Supplemental Material section 3. Note the total energy shown in Table 7 is the total energy that must be administered to all sources, normalized to the geometric mean of the total energy of all 8 test cases to facilitate comparison. The actual total energy will be scaled based on clinical assumptions and the choice of PS and its administered dose and selected light diffusers.
Table 7. Results for case study on injection constraint impact. The values of overdosage, runtime, and number of sources are taken as the average over 10 PDT-SPACE runs with different random seeds. The relative total energy is a unit-less value normalized to the geometric mean of all energy values for comparison purpose.
Figure 13 shows the correspondence between the total number of required sources and the total required energy (log scale). A clear relationship can be observed; when using more sources, the total energy that must be administered to treat the tumor is reduced, which will lead to significantly shorter treatment times.
Fig. 13. Relationship between the number of sources required (horizontal axis) vs. the total energy required (vertical axis). Vertical axes shown in log-scale. Showing unit-less relative total energy; actual total energy should be further scaled with practical PS and light diffuser assumptions.
To further examine the impact on treatment plan quality of the injection constraint size when the number of light sources allowed is held constant, another experiment was performed with the larger injection diameters (90 $mm$ diameter for T1 and 50 $mm$ diameter for T4) with the treatment plans constrained to use a maximum of 12 sources for T1 and 4 sources for T4. Observe the results in Table 8; with the same injection constraint size and a similar number of sources, on average the overdosage obtained by aligning the sources with the minor axis is lower than that obtained by aligning the sources with the major axis.
The results from Table 7 and Table 8 indicate that aligning the light sources along the tumor minor axis generally allows the program to utilize a larger cross-sectional area to generate initial source placements and results in more freedom for optimization that will eventually lead to reduction in OAR damage. The trade-off here is that more sources will be used by the optimizer, and the chance for the program not able to find a valid solution will increase. Treating T1 along its minor axis, for example, requires an injection constraint diameter of at least 50 $mm$ and the optimizer has a high chance of failing to find any plan that can treat 98% of the tumor using a small number of sources ($\leq 10$). If treating along the major axis of T1, on the other hand, valid solutions can be generated with a 40 $mm$ injection diameter and 8-9 sources.
Table 8. Results for case study on injection constraint impact with constrained number of sources (T1 constrained to maximum 12 sources, T4 constrained to maximum 4 sources). The values of overdosage, runtime, and number of sources are taken as the average over 10 PDT-SPACE runs with different random seeds.
6. Discussion
Surgeons and physicists have recently recognized the need to develop stricter procedures to improve iPDT efficacy and hence patient outcome, in particular for glioma therapy [10,63,64]. The procedures to plan iPDT treatment are now in the foreground of research efforts and clinical translation.
In this work, we have introduced new features to prior work in [16] for usages in clinical settings, as well as proposed several improved approaches to reduce damage to healthy tissue while keeping the number of sources at a minimum. This iPDT treatment planning tool can now be used in a more realistic setting where the point of injection of the sources can be suggested by the user, and it is able to produce reasonable or better starting source placements than manually-generated ones in its automated planning process to reduce human effort without affecting planning runtime. These new capabilities enable experimentation with many possible source counts, injection locations and burr hole sizes to find a plan that achieves the best result quality with a feasible surgery.
The results shown in this work are all obtained by running the same experiments 10 times and recording the averages and standard deviations for further analysis. Through the experiment, it was observed that the iPDT treatment outcome prediction for a single experiment varies by a large amount due to randomness, as shown by the geometric averages of the standard deviations shown in Tables 4 and 5. Further, in extreme cases a single run of the tool might give results with more than 2$\times$ the OAR damage of the average run, or with $\sim 0.5\times$ the average OAR damage. This is also true for the planning runtimes: the larger tumor models (T1 and T2) requiring more sources typically take longer to run compared to the smaller tumor models; if the annealer makes a move that escapes the current local minima, the runtime will significantly increase while the program explores the new placement. During the experiments there are cases where the tool took more than 24 hrs to run, but these usually lead to a reduction in the OAR volumes.
One way to make use of this randomness would be to run the tool a few times and take the best placement if time permits. However, considering the average $\sim 4.4 h$ runtime across all tumor models, the runtime differences between large and small tumors, and the runtime variations caused by randomness, the planning process would be significantly extended, especially for the larger tumors. One way to reduce runtime is to target the $\sim 57\%$ [16] of runtime spent on running FullMonte to model the forward light propagation in tissues, using accelerated Monte Carlo simulations with FPGA [65], GPU [66], or cloud computational platforms [67]. The majority of the remaining $\sim 43\%$ of runtime is spent solving the convex optimization power allocation problem and adjusting the tissue weights to target 98% tumor coverage; this process is sequential and only one thread is active. Future research could investigate methods to make this step parallel to reduce runtime. The total runtime is also affected by the tumor shape, the number of sources, and the number of moves made.
As more options are introduced, we decided to organize all inputs to the Simulated Annealing engine in a single, self-documenting, and easily modifiable XML file. An example of the file is included in the Supplemental Material and further descriptions can be found on the project Gitlab Wiki at https://gitlab.com/FullMonte/pdt-space/-/wikis/Running-PDT-SPACE. Future work on PDT-SPACE will extend this xml input file to allow fine-tuning of the plan to satisfy any detailed trajectory constraints.
7. Conclusion
This work extends the PDT-SPACE tool to further automate treatment plan optimization while taking into account additional clinical surgery constraints. The PDT-SPACE optimizer can now automatically create an initial source placement if none is provided, simplifying its use. Light sources can be constrained in variety of ways, including limits on the locations through which the sources are inserted, allowable source lengths, and maximum source count. Several improvements to the iterative source location optimization algorithms also reduce the damage to OAR. Taken together, these enhancements enable what-if analysis of various treatment plan options and surgery paths for both pre-clinical studies and personalized PDT treatment planning.
Funding
Ontario Ministry of Economic Development, Job Creation and Trade (ORF 08-023); Ontario Ministry of Health and Long-Term Care.
Acknowledgments
The authors would like to thank Abdul-Amir (Abed) Yassine for his contributions on the PDT-SPACE open-source software tool, as well as insights provided during the completion of this work. We would also like to thank Dr. Nardin Samuel for useful advice related to practical constraints in neurosurgery.
Disclosures
The authors declare no conflicts of interest.
Data Availability
Data underlying the results presented in this paper are available in Ref. [68].
Supplemental document
See Supplement 1 for supporting content.
References
1. C. Morton, R.-M. Szeimies, A. Sidoroff, and L. Braathen, “European guidelines for topical photodynamic therapy part 1: treatment delivery and current indications–actinic keratoses, bowen’s disease, basal cell carcinoma,” J. Eur. Acad. Dermatol. Venereol. 27(5), 536–544 (2013). [CrossRef]
2. C. Bennett, N. Vakil, J. Bergman, et al., “Consensus statements for management of barrett’s dysplasia and early-stage esophageal adenocarcinoma, based on a delphi process,” Gastroenterology 143(2), 336–346 (2012). [CrossRef]
3. K. Inoue, “5-aminolevulinic acid-mediated photodynamic therapy for bladder cancer,” Int. J. Urol. 24(2), 97–101 (2017). [CrossRef]
4. L. D. Lilge, J. Wu, Y. Xu, A. Manalac, D. Molenhuis, F. Schwiegelshohn, L. Vesselov, W. Embree, M. Nesbitt, V. Betz, A. Mandel, M. A. S. Jewett, and G. S. Kulkarni, “Minimal required PDT light dosimetry for nonmuscle invasive bladder cancer,” J. Biomed. Opt. 25(6), 068001 (2020). [CrossRef]
5. L. Lilge, M. Olivo, S. Schatz, J. MaGuire, M. Patterson, and B. Wilson, “The sensitivity of normal brain and intracranially implanted vx2 tumour to interstitial photodynamic therapy,” Br. J. Cancer 73(3), 332–343 (1996). [CrossRef]
6. G. Shafirstein, D. Bellnier, E. Oakley, S. Hamilton, M. Potasek, K. Beeson, and E. Parilov, “Interstitial photodynamic therapy–a focused review,” Cancers 9(12), 12 (2017). [CrossRef]
7. G. Gunaydin, M. E. Gedik, and S. Ayan, “Photodynamic therapy–current limitations and novel approaches,” Front. Chem. 9, 691697 (2021). [CrossRef]
8. Y. Wan, L.-H. Fu, C. Li, J. Lin, and P. Huang, “Conquering the hypoxia limitation for photodynamic therapy,” Adv. Mater. 33(48), 2103978 (2021). [CrossRef]
9. S. Kwiatkowski, B. Knap, D. Przystupski, J. Saczko, E. Kédzierska, K. Knap-Czop, J. Kotlińska, O. Michel, K. Kotowski, and J. Kulbacka, “Photodynamic therapy–mechanisms, photosensitizers and combinations,” Biomed. & pharmacotherapy 106, 1098–1107 (2018). [CrossRef]
10. M. Aumiller, C. Heckl, S. Quach, H. Stepp, B. Ertl-Wagner, R. Sroka, N. Thon, and A. Rühm, “Interrelation between spectral online monitoring and postoperative t1-weighted mri in interstitial photodynamic therapy of malignant gliomas,” Cancers 14(1), 120 (2021). [CrossRef]
11. A. Johansson, F. Faber, G. Kniebühler, H. Stepp, R. Sroka, R. Egensperger, W. Beyer, and F.-W. Kreth, “Protoporphyrin ix fluorescence and photobleaching during interstitial photodynamic therapy of malignant gliomas for early treatment prognosis,” Lasers Surg. Med. 45(4), 225–234 (2013). [CrossRef]
12. A. Rendon, J. Beck, and L. Lilge, “Linear feasibility algorithms for treatment planning in interstitial photodynamic therapy,” Proc. SPIE 6845, 68450O (2008). [CrossRef]
13. A. Johansson, J. Axelsson, J. Swartling, T. Johansson, S. Pålsson, J. Stensson, M. Einarsdóttír, K. Svanberg, N. Bendsoe, K. M. Kälkner, S. Nilsson, S. Svanberg, and S. Andersson-Engels, “Interstitial photodynamic therapy for primary prostate cancer incorporating real-time treatment dosimetry,” Proc. SPIE 6427, 64270O (2007). [CrossRef]
14. A.-A. Yassine, W. Kingsford, Y. Xu, J. Cassidy, L. Lilge, and V. Betz, “Automatic interstitial photodynamic therapy planning via convex optimization,” Biomed. Opt. Express 9(2), 898–920 (2018). [CrossRef]
15. A.-A. Yassine, L. Lilge, and V. Betz, “Optimizing interstitial photodynamic therapy with custom cylindrical diffusers,” J. Biophotonics 12(1), e201800153 (2019). [CrossRef]
16. A.-A. Yassine, L. Lilge, and V. Betz, “Optimizing interstitial photodynamic therapy planning with reinforcement learning-based diffuser placement,” IEEE Trans. Biomed. Eng. 68(5), 1668–1679 (2021). [CrossRef]
17. A.-A. Yassine, L. Lilge, and V. Betz, “Machine learning for real-time optical property recovery in interstitial photodynamic therapy: a stimulation-based study,” Biomed. Opt. Express 12(9), 5401–5422 (2021). [CrossRef]
18. J. Cassidy, A. Nouri, V. Betz, and L. Lilge, “High-performance, robustly verified Monte Carlo simulation with FullMonte,” J. Biomed. Opt. 23(08), 1–11 (2018). [CrossRef]
19. J. Jönsson and E. Berrocal, “Multi-scattering software: part i: online accelerated monte carlo simulation of light transport through scattering media,” Opt. Express 28(25), 37612–37638 (2020). [CrossRef]
20. D. Frantz, J. Jönsson, and E. Berrocal, “Multi-scattering software part ii: experimental validation for the light intensity distribution,” Opt. Express 30(2), 1261–1279 (2022). [CrossRef]
21. T. M. Baran and T. H. Foster, “Comparison of flat cleaved and cylindrical diffusing fibers as treatment sources for interstitial photodynamic therapy,” Med. Phys. 41(2), 022701 (2014). [CrossRef]
22. P. Dasgupta and B. Challacombe, “Technological Innovation in the BJUI,” (2013).
23. A. Johansson, J. Axelsson, S. Andersson-Engels, and J. Swartling, “Realtime light dosimetry software tools for interstitial photodynamic therapy of the human prostate,” Med. Phys. 34(11), 4309–4321 (2007). [CrossRef]
24. G. Hennig, H. Stepp, and A. Johansson, “Photobleaching-based method to individualize irradiation time during interstitial 5-aminolevulinic acid photodynamic therapy,” Photodiagn. Photodyn. Ther. 8(3), 275–281 (2011). [CrossRef]
25. A. Kim, M. Khurana, Y. Moriyama, and B. C. Wilson, “Quantification of in vivo fluorescence decoupled from the effects of tissue optical properties using fiber-optic spectroscopy measurements,” J. Biomed. Opt. 15(6), 067006 (2010). [CrossRef]
26. T. M. Baran, M. C. Fenn, and T. H. Foster, “Determination of optical properties by interstitial white light spectroscopy using a custom fiber optic probe,” J. Biomed. Opt. 18(10), 107007 (2013). [CrossRef]
27. S. R. H. Davidson, R. A. Weersink, M. A. Haider, M. R. Gertner, A. Bogaards, D. Giewercer, A. Scherz, M. D. Sherar, M. Elhilali, J. L. Chin, J. Trachtenberg, and B. C. Wilson, “Treatment planning and dose analysis for interstitial photodynamic therapy of prostate cancer,” Phys. Med. Biol. 54(8), 2293–2313 (2009). [CrossRef]
28. G. Shafirstein, D. A. Bellnier, E. Oakley, S. Hamilton, M. Habitzruther, L. Tworek, A. Hutson, J. A. Spernyak, S. Sexton, L. Curtin, S. G. Turowski, H. Arshad, and B. Henderson, “Irradiance controls photodynamic efficacy and tissue heating in experimental tumours: implication for interstitial pdt of locally advanced cancer,” Br. J. Cancer 119(10), 1191–1199 (2018). [CrossRef]
29. C. Dupont, G. Baert, S. Mordon, and M. Vermandel, “Parallelized monte-carlo dosimetry using graphics processing units to model cylindrical diffusers used in photodynamic therapy: from implementation to validation,” Photodiagn. Photodyn. Ther. 26, 351–360 (2019). [CrossRef]
30. C. Dupont, N. Betrouni, S. Mordon, N. Reyns, and M. Vermandel, “5-ALA photodynamic therapy in neurosurgery, towards the design of a treatment planning system: a proof of concept,” IRBM 38(1), 34–41 (2017). [CrossRef]
31. C. Dupont, N. Reyns, S. Mordon, and M. Vermandel, “Photodynamic therapy in neurosurgery: a proof of concept of treatment planning system,” Proc. SPIE 10047, 1004714 (2017). [CrossRef]
32. K. W. Beeson, E. Parilov, M. Potasek, M. M. Kim, and T. C. Zhu, “Validation of combined Monte Carlo and photokinetic simulations for the outcome correlation analysis of benzoporphyrin derivative-mediated photodynamic therapy on mice,” J. Biomed. Opt. 24(03), 1 (2019). [CrossRef]
33. M. M. Kim, R. Penjweini, X. Liang, and T. C. Zhu, “Explicit macroscopic singlet oxygen modeling for benzoporphyrin derivative monoacid ring A (BPD)-mediated photodynamic therapy,” J. Photochem. Photobiol. B: Biol. 164, 314–322 (2016). [CrossRef]
34. M. M. Kim, A. A. Ghogare, A. Greer, and T. C. Zhu, “On the in vivo photochemical rate parameters for pdt reactive oxygen species modeling,” Phys. Med. Biol. 62(5), R1–R48 (2017). [CrossRef]
35. A. Izumoto, T. Nishimura, H. Hazama, N. Ikeda, Y. Kajimoto, and K. Awazu, “Singlet oxygen model evaluation of interstitial photodynamic therapy with 5-aminolevulinic acid for malignant brain tumor,” J. Biomed. Opt. 25(06), 1 (2019). [CrossRef]
36. K. K.-H. Wang, J. C. Finlay, T. M. Busch, S. M. Hahn, and T. C. Zhu, “Explicit dosimetry for photodynamic therapy: macroscopic singlet oxygen modeling,” J. Biophotonics 3(5-6), 304–318 (2010). [CrossRef]
37. L. B. Rocha, H. T. Soares, M. I. P. Mendes, A. Cabrita, F. A. Schaberle, and L. G. Arnaut, “Necrosis depth and photodynamic threshold dose with redaporfin-PDT,” Photochem. Photobiol. 96(3), 692–698 (2020). [CrossRef]
38. K. T. Ramadan, C. McFadden, B. Gomes, F. Schwiegelshohn, R. V. Ribeiro, H. H. Chan, V. Betz, M. Cypel, and L. Lilge, “Determination of optical properties and photodynamic threshold of lung tissue for treatment planning of in vivo lung perfusion assisted photodynamic therapy,” Photodiagn. Photodyn. Ther. 35, 102353 (2021). [CrossRef]
39. R. Ferraz, J. Ferreira, P. Menezes, C. Sibata, O. C. e Silva Jr, and V. S. Bagnato, “Determination of threshold dose of photodynamic therapy to measure superficial necrosis,” Photomed. Laser Surg. 27(1), 93–99 (2009). [CrossRef]
40. T. J. Farrell, B. C. Wilson, M. S. Patterson, and M. C. Olivo, “Comparison of the in vivo photodynamic threshold dose for photofrin, mono-and tetrasulfonated aluminum phthalocyanine using a rat liver model,” Photochem. Photobiol. 68(3), 394–399 (1998). [CrossRef]
41. L. Lilge and B. C. Wilson, “Photodynamic therapy of intracranial tissues: a preclinical comparative study of four different photosensitizers,” J. clin. laser med. & surgery 16(2), 81–91 (1998). [CrossRef]
42. A. M. Lozano, N. Lipsman, H. Bergman, P. Brown, S. Chabardes, J. W. Chang, K. Matthews, C. C. McIntyre, T. E. Schlaepfer, M. Schulder, Y. Temel, J. Volkmann, and J. K. Krauss, “Deep brain stimulation: current challenges and future directions,” Nat. Rev. Neurol. 15(3), 148–160 (2019). [CrossRef]
43. G. Raffa, A. Scibilia, A. Conti, G. Ricciardo, V. Rizzo, A. Morelli, F. F. Angileri, S. M. Cardali, and A. Germano, “The role of navigated transcranial magnetic stimulation for surgery of motor-eloquent brain tumors: a systematic review and meta-analysis,” Clin. Neurol. Neurosurg. 180, 7–17 (2019). [CrossRef]
44. A. Trébuchon and P. Chauvel, “Electrical stimulation for seizure induction and functional mapping in stereoelectroencephalography,” J. Clin. Neurophysiol. 33(6), 511–521 (2016). [CrossRef]
45. L. Bello, A. Castellano, E. Fava, G. Casaceli, M. Riva, G. Scotti, S. M. Gaini, and A. Falini, “Intraoperative use of diffusion tensor imaging fiber tractography and subcortical mapping for resection of gliomas: technical considerations,” Neurosurg. focus 28(2), E6 (2010). [CrossRef]
46. M. Lacroix, D. Abi-Said, D. R. Fourney, Z. L. Gokaslan, W. Shi, F. DeMonte, F. F. Lang, I. E. McCutcheon, S. J. Hassenbusch, E. Holland, K. Hess, C. Michael, D. Miller, and R. Sawaya, “A multivariate analysis of 416 patients with glioblastoma multiforme: prognosis, extent of resection, and survival,” J. Neurosurg. 95(2), 190–198 (2001). [CrossRef]
47. D. L. Collins, A. P. Zijdenbos, V. Kollokian, J. G. Sled, N. J. Kabani, C. J. Holmes, and A. C. Evans, “Design and construction of a realistic digital brain phantom,” IEEE Trans. on Med. Imaging 17(3), 463–468 (1998). [CrossRef]
48. Q. Fang, “Mesh-based monte carlo method using fast ray-tracing in plücker coordinates,” Biomed. Opt. Express 1(1), 165–175 (2010). [CrossRef]
49. K. Clark, B. Vendt, K. Smith, J. Freymann, J. Kirby, P. Koppel, S. Moore, S. Phillips, D. Maffitt, M. Pringle, L. Tarbox, and F. Prior, “The Cancer Imaging Archive (TCIA): Maintaining and Operating a Public Information Repository,” J. Digit. Imaging 26(6), 1045–1057 (2013). [CrossRef]
50. T. J. Beck, F. W. Kreth, W. Beyer, J. H. Mehrkens, A. Obermeier, H. Stepp, W. Stummer, and R. Baumgartner, “Interstitial photodynamic therapy of nonresectable malignant glioma recurrences using 5-aminolevulinic acid induced protoporphyrin ix,” Lasers Surg. Medicine: The Off. J. Am. Soc. for Laser Medicine Surg. 39(5), 386–393 (2007). [CrossRef]
51. L. M. Vesselov, W. Whittington, and L. Lilge, “Performance evaluation of cylindrical fiber optic light diffusers for biomedical applications,” Lasers Surg. Medicine: The Off. J. Am. Soc. for Laser Medicine Surg. 34(4), 348–351 (2004). [CrossRef]
52. L. Vesselov, W. Whittington, and L. Lilge, “Design and performance of thin cylindrical diffusers created in ge-doped multimode optical fibers,” Appl. Opt. 44(14), 2754–2758 (2005). [CrossRef]
53. A. Rendon, J. C. Beck, and L. Lilge, “Treatment planning using tailored and standard cylindrical light diffusers for photodynamic therapy of the prostate,” Phys. Med. Biol. 53(4), 1131–1149 (2008). [CrossRef]
54. A. Yaroslavsky, P. Schulze, I. Yaroslavsky, R. Schober, F. Ulrich, and H. Schwarzmaier, “Optical properties of selected native and coagulated human brain tissues in vitro in the visible and near infrared spectral range,” Phys. Med. Biol. 47(12), 2059–2073 (2002). [CrossRef]
55. N. Honda, K. Ishii, Y. Kajimoto, T. Kuroiwa, and K. Awazu, “Determination of optical properties of human brain tumor tissues from 350 to 1000 nm to investigate the cause of false negatives in fluorescence-guided resection with 5-aminolevulinic acid,” J. Biomed. Opt. 23(07), 1 (2018). [CrossRef]
56. V. N. Du Le, J. Provias, N. Murty, M. S. Patterson, Z. Nie, J. E. Hayward, T. J. Farrell, W. McMillan, W. Zhang, and Q. Fang, “Dual-modality optical biopsy of glioblastomas multiforme with diffuse reflectance and fluorescence: ex vivo retrieval of optical properties,” J. Biomed. Opt. 22(2), 027002 (2017). [CrossRef]
57. H. Soleimanzad, H. Gurden, and F. Pain, “Optical properties of mice skull bone in the 455-to 705-nm range,” J. Biomed. Opt. 22(1), 010503 (2017). [CrossRef]
58. S. W. Cramer and C. C. Chen, “Photodynamic therapy for the treatment of glioblastoma,” Frontiers in Surgery 6, 81 (2020). [CrossRef]
59. J. O’Rourke, “Finding minimal enclosing boxes,” Int. J. Comput. Inform. Sci. 14(3), 183–199 (1985). [CrossRef]
60. C.-T. Chang, B. Gorissen, and S. Melchior, “Fast oriented bounding box optimization on the rotation group SO (3, $\mathbb {R}$),” ACM Trans. Graph. 30(5), 1–16 (2011). [CrossRef]
61. I. T. Jolliffe, Principal component analysis for special types of data (Springer, 2002).
62. W. Ames and C. Brezinski, Numerical recipes in Fortran (The art of scientific computing): WH Press, SA Teukolsky, WT Vetterling and BP Flannery, Cambridge Univ. Press, Cambridge, 1992. 963 pp., ISBN 0-521-43064-X. (North-Holland, 1993).
63. H.-A. Leroy, L. Guérin, F. Lecomte, G. Baert, A.-S. Vignion, S. Mordon, and N. Reyns, “Is interstitial photodynamic therapy for brain tumors ready for clinical practice? A systematic review,” Photodiagn. Photodyn. Ther. 36, 102492 (2021). [CrossRef]
64. H. Stepp and W. Stummer, “5-ALA in the management of malignant glioma,” Lasers Surg. Med. 50(5), 399–419 (2018). [CrossRef]
65. T. Young-Schultz, L. Lilge, S. Brown, and V. Betz, “Using OpenCL to enable software-like development of an FPGA-accelerated biophotonic Cancer treatment simulator,” in Proceedings of the 2020 ACM/SIGDA International Symposium on Field-Programmable Gate Arrays, (Association for Computing Machinery, New York, NY, USA, 2020), FPGA ’20, p. 86–96.
66. T. Young-Schultz, S. Brown, L. Lilge, and V. Betz, “FullMonteCUDA: a fast, flexible, and accurate GPU-accelerated Monte Carlo simulator for light propagation in turbid media,” Biomed. Opt. Express 10(9), 4711–4726 (2019). [CrossRef]
67. S. Wang, X. Y. Dai, S. Ji, T. Saeidi, F. Schwiegelshohn, A.-A. Yassine, L. Lilge, and V. Betz, “Scalable and accessible personalized photodynamic therapy optimization with fullmonte and pdt-space,” J. Biomed. Opt. 27(08), 083006 (2022). [CrossRef]
68. S. Wang, “PDT-SPACE Data,” Gitlab, 2022, https://gitlab.com/FullMonte/pdt-space/-/tree/master/data.
