Edited by: Ulas Sunar, University at Buffalo, USA
Reviewed by: Adam Thomas Eggebrecht, Washington University School of Medicine, USA; Rongxiao Zhang, Dartmouth College, USA
*Correspondence: Jeffrey Cassidy, Edward S. Rogers Sr., Department of Electrical and Computer Engineering, University of Toronto, 10 King's College Road, Toronto, ON M5S 3G4, Canada e-mail:
This article was submitted to Biomedical Physics, a section of the journal Frontiers in Physics.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
Monte Carlo (MC) simulation is recognized as the “gold standard” for biophotonic simulation, capturing all relevant physics and material properties at the perceived cost of high computing demands. Tetrahedral-mesh-based MC simulations particularly are attractive due to the ability to refine the mesh at will to conform to complicated geometries or user-defined resolution requirements. Since the chosen MC method does not inherently require approximations on the materials, it is applicable to the broadest set of biophotonic simulation problems. MC methods also have other implementation features including inherent parallelism, and permit a continuously-variable quality-runtime tradeoff. We demonstrate here a complete MC-based prospective fluence dose evaluation system for interstitial PDT to generate dose-volume histograms on a tetrahedral mesh geometry description. To our knowledge, this is the first such system for general interstitial photodynamic therapy employing MC methods and is therefore applicable to a very broad cross-section of anatomy and material properties. We demonstrate that evaluation of dose-volume histograms is an effective variance-reduction scheme in its own right which greatly reduces the number of packets required and hence runtime required to achieve acceptable result confidence. We conclude that MC methods are feasible for general PDT treatment evaluation and planning, and considerably less costly than widely believed.
Photodynamic Therapy (PDT) is an emerging treatment modality for cancerous lesions in a variety of indications. It is widely used for indications which are superficial in terms of both the method of light application and the lesion depth, so currently mostly actinic keratosis [
Photodynamic cell killing mechanisms depend on cell physiology (response to radical production), and require the overlap of photosensitizer, tissue oxygen, and photon fluence. Selectivity is achieved through influencing these factors by methods including: designing photosensitizer biochemistry for preferential uptake by the tumor; exploiting pharmacokinetic differences between diseased and healthy tissue; controlling oxygen delivery dynamics; artificially inducing or using naturally-occurring differential sensitivity between tumor and nearby tissue; and by controlling the shape of the light fluence field to achieve spatially selective photosensitizer activation. We focus here on control of the light fluence field, as it is largely independent of the other factors and all else being equal an improvement in light targeting always improves the outcome. The ability to simulate the light fluence field is also critical to understanding optical measurements made (including for implicit or explicit dosimetry parameters), and for optimal placement of monitoring probes to detect optical properties
MC is also the standard of care for ionizing radiation therapy (IRT), where the coefficients of absorption and scattering are lower resulting in significantly fewer interaction events per photon and hence much lower compute requirements. Ionizing radiation therapy (IRT) has become a main-line treatment for cancers, with well-developed workflows, numerical optimization based on Monte Carlo simulation, and rigorous quality-assurance procedures. Such practices are still emerging for PDT treatment planning and monitoring, with some challenges due to the unique characteristics of PDT. Compared to IRT, the red or near-infrared photons used for PDT have significantly higher coefficients of interaction with the target tissue. The greater propensity to scatter leads to light confinement within much smaller volumes around the implanted light source which facilitates treatment targeting. Such selectivity comes at a cost, though; since each photon undergoes a very large number of scattering events prior to absorption (and hence causing cell damage), the cost of computational modeling is correspondingly increased. A complete PDT treatment plan covers a large range of parameters, some under control of the treating physician and some outside as summarized in Table
Patient anatomy | ✓ | ||
Patient optical properties | ✓ | ||
Activation wavelengths | ✓ | ✓* | ✓ |
Number of sources | ✓ | ✓ | |
Source emission profile | ✓ | ✓ | |
Photosensitizer choice | ✓ | ||
Photosensitizer dose | ✓ | ||
Drug-light interval | ✓ | ||
Light time and power | ✓ | ✓ | |
Additional protective or sensitizing techniques | ✓ |
Based on simulations of the fluence field, we derive dose-volume histograms (DVH) and related metrics which are extremely widely used in radiation therapy planning. Operating within a defined planning volume, the dose-volume histogram displays the percentage of each organ within the planning volume receiving a dose (vertical axis) which exceeds a given dose level (independent, horizontal axis). The DVH permits the ability to query what fraction of an organ receives a dose exceeding ΦJ/cm2, or conversely to find the minimum dose received by the
One advantage of IRT over PDT is that the imaging (CT) and treatment use similar photon energies and hence the imaging modality used to build the geometric model directly provides measurements of the necessary propagation properties. On the other hand, the patient properties relevant to PDT (coefficients of absorption μ
Figure
Source placement | Derive plan which delivers the desired fluence dose | ≈ 102 − 103 | Low | High | Low |
Signoff | Verify plan remains acceptable under expected variation of optical properties and source placements | ≈ 101 − 102 | Moderate | Moderate | Low |
Monitoring placement | Find optimal sensor probe placements to monitor optical properties and/or fluence | ≈ 101 | High | Low | Low |
Online adjustment | (1) Verify configuration as-placed is acceptable, adjust power/time, (2) Infer and adjust to optical property changes | ≈ 101 − 102 | Moderate | Moderate | Critical |
From this workflow, we can derive a set of requirements for a simulation environment which could be applied to any indication for which PDT is suitable:
Gives valid results in low-scattering areas
Able to accommodate curved surfaces
Refinable geometry model to provide high detail where needed
Able to model varied light emitters and probes, including customized extended sources
Predictable, quantifiable result uncertainty
Modest run-time requirements to support iterative optimization
Broad range of quality-effort tradeoffs, capable of high precision when needed
We demonstrate below a software implementation “FullMonte” which meets all these criteria using a tetrahedral mesh description of the planning volume, and a Monte Carlo simulation to solve for the light distribution. We also show that the use of dose-volume histograms for assessing a proposed treatment geometry acts as a variance reduction scheme which greatly reduces the computational cost of Monte Carlo simulations. Taken together with our work on hardware acceleration, we conclude that MC-based prospective PDT treatment evaluation, planning, and online monitoring are both feasible and desirable.
Of the workflow illustrated in Figure
Simulation of light propagation through living tissue is complicated by the very high scattering coefficients within the optical window generally used for photosensitizer activation. Physically, a photon will scatter many times before absorption, with the expected number of interactions depending on the material's albedo. Tissue has a higher albedo for interactions with optical photons than x-rays, meaning that considerably more scattering must be modeled for PDT compared to IRT. The Radiative Transfer Equation (RTE, or sometimes ERT) expresses the conservation relationship which must hold for light transport in a scattering medium. Solution techniques can be classified into three groups: analytic and finite-element numerical solutions using the diffusion approximation (the simplest and least accurate method); other numerical methods to solve the RTE directly; and Monte Carlo methods.
The Diffusion Approximation (DA) [
However, there are several conditions which would invalidate the diffusion approximation in a large fraction of the treatment volume:
Use of extended sources where an appreciable volume of tissue is located near (within a few mean free paths of) a light source (Jacques and Pogue [
Heterogeneous structures with a large number of material interfaces (see previous point)
Presence of air-tissue interfaces with reflection and refraction (e.g., sinuses, esophagus in head and neck; see Figures 1, 7 in Jacques [
Presence of low-scattering voids (e.g., cerebral-spinal fluid in brain/spine; bladder) or high-absorbing regions (discussed with references in Boas et al. [
For many indications, one or more of these conditions will be violated, yielding incorrect simulation results which would jeopardize the safety and efficacy of the treatment.
Other solution methods with less onerous approximations have been applied to the RTE. Klose and Larsen, for instance, have shown an SP3 solution [
Monte Carlo methods are generally acknowledged as the “gold standard” for biophotonic simulations, covering the widest set of problem descriptions with accuracy which can be arbitrarily refined, at the perceived cost of long run times required to produce sufficient result confidence. Particularly for air-tissue interfaces (e.g., head-and-neck cancers) and volumes with low-scattering voids (e.g., bladder cancers), the approximations required for other solution techniques are likely to invalidate the results, while MC methods remain able to deliver a quality solution.
A diversity of prior work exists using Monte Carlo methods to solve biophotonic problems with a variety of geometry descriptions. The basic concept started in the early 1980's with Wilson and Adam [
Diffuse imaging applications [
To address refractive interfaces correctly, several 3D tetrahedral-mesh-based simulators have been created [
Offsetting the many benefits of Monte Carlo simulation are the perceived computational costs. Despite a high amount of effort expended on software optimization, required run times for high-quality simulations of complex geometries remain long. Shen and Wang's TIM-OS made efficient use of multi-core processors and Intel's automatic vectorization capabilities for double-precision floating point arithmetic. We improved further on those results when writing FullMonte by moving to single precision floating-point (including validating that there was no loss of quality), and using hand-tuned SIMD instructions for further performance gains. We have also demonstrated a pre-existing but previously-unknown quality-runtime tradeoff (the
For superficial applications, PDT treatment planning generally consists of illuminating the surface with a set energy per unit area, after giving a standard photosensitizer dose in mg/kg. For deeper-seated lesions optical diffusers may be used to increase the penetration depth of the light [
Some of the first algorithmic PDT planning was demonstrated by Altschuler et al. [
Researchers at Lund University and SpectraCure AB [
A series of clinical trials vascular-targeted photodynamic therapy of prostate cancer using Tookad at the Ontario Cancer Institute [
While a few works have used MC methods, for instance calculating the dose received at the surface of a cylindrical area around a light source in the prostate [
We divide the treatment evaluation process into three stages discussed individually below. First, the problem description consisting of patient geometry, optical properties, target volume, and light sources must be defined. Next, the description is given to the Monte Carlo simulator to run a user-specified number of photon packets. Finally, the results are interpreted in terms of dose-volume histograms and fluence field visualizations.
Future work would naturally use a numerical optimization method coupled with the plan evaluation methods shown here to iteratively refine the source descriptions toward user-prescribed goals. Given sufficient computing resources, different wavelengths may also be contemplated in the optimization by assigning wavelength-appropriate optical properties to each tissue type.
Variability of optical properties is one of the most important factors in safe and effective delivery of PDT. Computationally, our method can address this by evaluating multiple plans spanning the anticipated range of optical properties, an approach made more feasible by our short run times. A thorough study of the effect of such variability on PDT planning and delivery would be an important future contribution outside the scope of the present discussion, but enabled by the methods presented here. Given such an understanding, we would use these plan evaluation tools to develop treatment planning methods which are robust to the major sources of PDT variability.
The planning system accepts as input a tetrahedral mesh in which each mesh element is assigned a tissue type. Generally this would be done by a radiologist using CT or MRI information, followed by software generation of a volumetric mesh from the provided curves. Pending creation of such infrastructure, we demonstrate our capability using Digimouse [
Entire mouse | 306774 | 21354.2 |
Tumor | 436 | 4.2 |
Lung | 23148 | 179.7 |
Heart | 759 | 31.5 |
Liver | 1557 | 116.5 |
Muscle | 31873 | 579.5 |
Total DVH bounds | 57817 | 912.2 |
Monte Carlo simulation was carried out using our FullMonte software, which uses a tetrahedral mesh geometry description and sums the total energy
FullMonte is able to handle arbitrary combinations of a number of clinically relevant light sources: isotropic points, cut-end fibers (cone beam emitted from a disc of finite extent), volume sources (the interior of a tetrahedral element), face sources (a triangular face of a tetrahedral element), line sources, and pencil beams. For this instance, we used a single cylindrical diffusing fiber which emits isotropically in the plane normal to the fiber axis. Extension to sources with a tailored emission profile [
For each tetrahedral element which lies entirely within the bounding box, the fluence is calculated from the energy absorbed within the volume
Finally, the fluence field is visualized as both a set of isofluence contours to evaluate the spatial distribution. In this case, the source is close enough to the surface that the isofluence contours are truncated at the skin since fluence is not scored beyond the outer limit of the geometry. When using MC simulations for planning, we aim to use many low-quality runs to quickly eliminate poor candidates using the DVH variance reduction scheme, and produce a short list of several good plans for consideration. Having narrowed the list, we can spend more time on those few plans to produce high-quality results for visualization and sign-off.
Two treatment plans were tested to illustrate the method: a subjectively “good” plan, in which the line source was positioned along the long axis of the tumor [(11.4, 42.4, 10.4) − (11.7, 41.6, 11.0)], and a “bad” plan in which the source extended beyond the edges of the tumor into both lung and liver tissue [(11.5, 40.0, 10.4) − (12.5, 45.1, 12.1)]. The organs at risk which had significant volumes within the bounding box and received meaningful doses were the tumor, lung, liver, and surrounding muscle (the default material in Digimouse where no organ is assigned). The optical properties used in these simulations are listed in Table
Muscle, tumor | 10.0 | 0.23 | 0.9 | 1.0 | 0.9775 |
Heart | 11.0 | 0.11 | 0.9 | 1.0 | 0.9775 |
Stomach | 17.0 | 0.21 | 0.9 | 1.7 | 0.9878 |
Lung | 23.0 | 0.35 | 0.9 | 2.3 | 0.9850 |
Liver | 20.0 | 0.45 | 0.9 | 2.0 | 0.9780 |
Skeleton | 10.0 | 0.23 | 0.9 | 1.0 | 0.9775 |
The dose-volume histogram (DVH, as described in the introduction) results of one simulation run with one million packets for each of the plans (good/bad) is shown in Figure
As a statistical numerical method, Monte Carlo simulation shows variability in its output results, with the standard deviation scaling proportional to
The absorbed energy integrated over volume is not directly useful for calculating photodynamic dose because it does not separate photosensitizer absorption from other chromophores; that requires knowledge of photosensitizer concentration (outside the scope of the optical simulator), absorption cross-section (a material property), and fluence (derived from optical simulation). Conversion of absorbed energy into fluence requires simply multiplication by a constant (
The variance of an individual element in isolation is of greatest interest when looking at fluence maps, or calculating the value that would be measured by a very small isotropic point probe. For evaluating the merits of a specific treatment plan, though, examining and extracting values from the dose-volume histogram (DVH) is the customary method. The DVH is built from a set of simulated fluence values through sorting. We can expect, and indeed find, significant reduction in the variance of the order statistics represented by the DVH compared to the variance of its individual constituents. This is shown graphically in Figure
The same effect is noticeable when trying to determine
One interesting point is that the organs at risk (OAR), which have much
We posit that increasing the spatial resolution of the mesh may sometimes paradoxically decrease the variance in the DVH. The increased mesh resolution may modestly increase the compute time per packet, but a reduction in DVH variance will require the simulation of fewer packets to achieve an acceptable result. It may therefore be possible to reduce the required compute time by working with a finer mesh, increasing spatial resolution and the confidence in the DVH despite increasing the variance of the individual constituent elements. Figure
One advantage of the Monte Carlo method for biophotonic simulations is that it allows a continuously-variable tradeoff of run time and result variance. Due to statistical averaging in the output, the uncertainty (standard deviation) for the fluence in any individual output element is inversely proportional to
Simulations were conducted on a high-performance laptop (MacBook Pro 15″, mid-2014 model) with a latest-generation Intel
Desired result quality and problem geometry are the primary factors affecting required run time. As discussed in the DVH section above, for evaluation of treatment plans, several factors come into play to determine the required precision for individual elements to yield an acceptable measure of plan merit. Given the American Association of Physicists in Medicine (AAPM) absolute calibration standards of optical power meters for PDT to within 5% [
When conducting PDT treatment planning, a large number of candidate plans will be examined to determine which is the best under a given measure of plan merit. Comparing two plans will thus have some probability of yielding an incorrect result based on the confidence bounds of each of the solutions, and thus mis-directing the optimizer. One can and should try to use optimization techniques that are robust to noisy inputs. Ideally, though, we would like to run simulations which are “good enough” with respect to the degree of confidence in the output without incurring excessive computational cost (run time or financial cost of hardware).
A different degree of confidence may be required depending on the end goal of a simulation run. As described below, the generation of dose-volume histograms gives a variance reduction effect. If, on the other hand, the user is looking for the sensitivity of fluence within a single element to a small change of input parameters, for instance in looking to place a small fluence-rate monitoring probe, a much greater number of packets may be required to reduce the variance to an acceptable level. For such smaller regions, the options for variance reduction are more limited but spatial averaging remains a viable option. When designing a complete workflow one can achieve significant savings by customizing the level of effort to the level of quality needed, for instance running a large number of plans quickly to rapidly explore the space before running more refined calculations. It should also be feasible to separate source planning and monitor plan execution as discussed in the introduction such that many candidate placements can be examined at low packet counts, prior to running fewer high-resolution simulations for deriving a monitoring plan.
We have shown that Monte Carlo methods are computationally tractable for PDT treatment evaluation, permit continuously-variable quality-runtime tradeoffs, and are amenable to hardware acceleration. Additionally, the computation of dose-volume histograms (DVHs) from MC simulations functions as a significant variance-reduction scheme which greatly reduces the runtime required to achieve a given level of output confidence. Current run times can be sufficiently low, even on a laptop computer, that tetrahedral-mesh MC is a viable option for evaluating PDT treatment plans. Higher-precision simulations with a larger number of packets which will permit accurate sensitivities (e.g., to probe placements or to personalized optical properties) can also be calculated. For iterative treatment planning and algorithmic optimization, accelerated hardware enables a very large volume of simulations to be run quickly. In short, given an appropriate understanding of the requirements for acceptable output quality, and of the effects of DVH generation, Monte Carlo simulations for PDT treatment planning may be considerably (integer factors) less expensive than is generally perceived.
This work demonstrates the capability and low computational cost of using Monte Carlo methods to evaluate the dose delivered by a given light source configuration into a specified tissue geometry with known optical properties. In the practice of interstitial PDT delivery, however, patient optical properties vary widely [
Using the infrastructure developed, we aim to create an automated numerical optimization system which works toward user-defined treatment plan goals (minimum dose for tumor, maximum dose to organs at risk, and a penalty function for deviations). It will need to be robust to misdirection due to the output variance, though to what degree remains to be determined. One may contemplate a highly robust algorithm (e.g., simulated annealing) which works on very short, low-quality/high-variance simulation runs; or, one may choose to employ an algorithm which is less demanding in the number of sample plans to evaluate but requires higher-quality input.
In previous work, we have also demonstrated the capability to perform identical Monte Carlo simulations using Field-Programmable Gate Arrays (FPGA) custom digital logic [
While the present work focuses entirely on delivering a evaluating the light fluence delivered, definitions of photodynamic dose encompassing more factors exist [
The connection between mesh size and DVH variance needs further investigation. As discussed in the results section, there are contrary pressures when mesh element size decreases: individual element variance would tend to increase because of a smaller number of arrivals to be summed; DVH variance would however tend downwards due to the variance-reduction properties of sorting the fluence elements. Mesh refinement also adds a mild computational burden, which may or may not be offset by the decreased number of packets which will be required to achieve an acceptable level of confidence in the DVH.
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
We thank Dr. Robert Weersink of the Princess Margaret Cancer Centre for helpful discussions regarding treatment planning practices in radiotherapy. The authors are also grateful to KitWare Inc for creating and maintaining the Visualization Toolkit (VTK) which was used to create the interactive planning and visualization software.