2.1.

Integrated Computational Device

The Dosie™ hardware and proprietary software incorporates, in one integrated computational device, Monte Carlo (MC) simulations of light transport, light fluence rate, and light dose (fluence) in target regions as well as PK simulations of singlet oxygen dose (${[{{}^1\mathrm{O}}_2]}_{rx}$) and PDT dose. The MC approach is ideal for simulating light transport in biological tissue.^41,42 This method has been greatly expanded over the years and has been successfully ported⁴³ onto graphical processor units (GPU) to improve its performance. The subsequent PK simulations use the MC results and PS properties as inputs to calculate the resulting PDT dose and ${[{{}^1\mathrm{O}}_2]}_{rx}$.

A schematic diagram of Dosie™ is shown in Fig. 1. The host process is an independent component of the software that is responsible for scheduling all PDT-related numerical calculations. It provides access to and controls two modules: light transport (MC-module) and PK-module. A graphical user interface (GUI) enables users to define PDT session simulation parameters, to launch and control simulations, and to visualize the output. The parameters needed for the numerical calculations are added from a database or from direct optical measurements before treatment. They include the 3-D target shape (prepared by third-party software, generated by measurements, such as CT or MRI), the optical parameters of the target (light scattering and light absorption coefficients), PK parameters for the PS (see Table 1 for BPD parameters), and the initial oxygen concentration. The variable parameters include laser power and location(s), PS (drug) concentration, and treatment time. The first calculation step is to use the MC-module to determine the light intensity (fluence rate) and light dose (fluence) at every microscopic point in the tumor. CUDA calculations are performed on a massively parallel GPU processor to simulate the paths of tens of millions of light rays. The calculated light intensity field is transferred from GPU memory to the host memory to launch the PK-module on the CPU processers to calculate the PDT PK that includes the light-PS-excitation (and photobleaching), the PS-oxygen excitation to generate singlet oxygen, and the singlet oxygen reaction with the target (including singlet-oxygen induced photobleaching). Comprehensive proprietary graphics, developed by Simphotek in the Dosie™ device, are used to visualize 2-D and 3-D outputs of the light fluence (light dose) and fluence rate, PDT dose, and singlet oxygen dose as well as concentrations of PS, ground-state oxygen, and cancer killing the ${[{{}^1\mathrm{O}}_2]}_{rx}$ at every microscopic point in the 3-D target for the duration of the treatment. This enables the researcher or physician to localize areas of possible undertreatment or overtreatment in the operating room while the patient is undergoing treatment and make corrections.

Fig. 1

Schematic diagram for Dosie™ integrated computer and software computational device. The computer includes both a GPU and a multicore CPU. The software first performs a MC simulation of light transport that is followed by a PK simulation.

Table 1

BPD simulation parameters.10,44

2.2.

MC-Module (Light Transport)

The MC-module is based on MC photon transport that allows modeling light propagating in a heterogeneous translucent medium. Such a medium is characterized by high likelihood of (mostly forward) light scattering events occurring within short distances (submillimeter scale) and by a high albedo (absorption is typically two orders of magnitude smaller than the scattering). The MC-module numerically estimates the resulting 3-D maps of light fluence rates and light doses (fluences) within the PDT target translucent region (target volume) during a planned PDT session for a given arrangement of the incident or interstitially applied light sources.

For BPD simulation, a collimated light source is located in air and directed perpendicular to the surface of the target volume. Our modified proprietary version of MC models’ light propagation within voxel-based geometry, where the target volume is subdivided into voxels (cubes of identical sizes) so that the graphical output 3-D map of Dosie™ assigns one fluence rate value for each voxel. The target volume consists of heterogeneous media with known optical properties (i.e., absorption coefficient, scattering coefficient, scattering anisotropy factor, and refractive index) so that each voxel can be assigned a unique set of properties.

During MC simulations, millions of photon packets are launched toward the target volume, according to the light excitation distribution. Each packet—we will be using word “photon” instead for the rest of this section—has an initial weight of 1. Specular reflection/refraction is taken into account at the air-target surface interface by applying the Fresnel model for unpolarized light with the refractive index 1.4. Each photon can undergo multiple scattering events in the target volume until the photon either escapes the volume or is absorbed within the volume when its weight drops below a specified threshold. Within the medium, a photon propagates free, voxel by voxel, till the next scattering event. The length of such free propagation is randomly determined according to the scattering coefficient ${\mu }_{\mathrm{s}}$. If the optical properties change while the photon travels from one voxel to another, the free path is properly adjusted. At each scattering event, a trajectory direction is determined by applying the Henyey–Greenstein phase function. The photon weight is gradually reduced while traveling from one voxel to another by a factor $e^{-{\mu }_{a}l}$, where $l$ is the portion of its trajectory within each such a voxel. All fractions of the weight lost that occurred during propagation of each photon are deposited at every voxel it passes so that by the end of the simulation, we can estimate the fluence from the accumulated weights.

A voxel-based implementation of MC is relatively straightforward, as it does not require adding potentially time-consuming ray-boundary intersection tests. However, for geometries with multiple curved surfaces and/or for less scattering media, a boundary-based MC implementation may result in more accurate solutions due to more precise reflection/refraction calculations along the boundaries (e.g., as in this tetrahedral-based⁴⁰ MC implementation).

MC simulations can be very time consuming (i.e., hours) for translucent media. The main challenge is to make the simulations run in near real-time (i.e., a minute or less). At the current state of the computer hardware, only parallel calculations are capable of reaching near real-time performance. MC-proprietary module software in Dosie™ is a GPU CUDA code (an extension to C++ programming language that allows scheduling identical compute kernels to run as threads on GPU Core processors in parallel; CUDA compute capability 5.0 has been used in Dosie™) with a performance target of completing target light transport simulations within minutes. Each CUDA thread performs a simulation of a photon propagation accessing the voxels’ optical properties of each encountered voxel from the device constant memory and depositing the weight’s losses to the device shared memory. The final fluence calculations are done in the host process after the accumulated weights are transferred from the device memory to the host memory using the proprietary Dosie™ device, including experimental hardware configuration of Core i7-6820HQ (16G memory, 64-bit) CPU and NVIDIA Quadro M2000M (640 cores, 4G memory, CUDA 5.0) GPU, the near-real time target performance has been achieved.

2.3.

PK-Module

The PK-module of Dosie™ is a complicated and essential part of modeling PDT. Unlike the light transport, this is a quasi-quantum mechanical process in which the laser light energy (photon) is transformed into reactive chemical agents that kill cancer cells. The PK calculations are done for each voxel independent from other voxels. This allows the calculations to run in parallel by taking advantage of a multicore CPU architecture. PK Eqs. (1)–(3) are used to calculate the time evolution of the (ground state) PS concentration, $[S_0]$, the ground state oxygen concentration, $[{{}^3\mathrm{O}}_2]$, and the reacted singlet oxygen concentration, ${[{{}^1\mathrm{O}}_2]}_{rx}$, where the value $[S_0]$ is the BPD concentration and $\varphi$ is the fluence rate. This is an initial value problem with the given initial fields $[S_0]$ $(x,y,z,t=0)\ge 0$ and $[{{}^3\mathrm{O}}_2]$ $(x,y,z,t=0)>0$ while ${[{{}^1\mathrm{O}}_2]}_{rx}$ $(x,y,z,t=0)=0$, and with satisfying a non-negativity requirement at all times: $[S_0]$$(t>0)\ge 0$ and $[{{}^3\mathrm{O}}_2]$ $(t>0)\ge 0$. The initial $[S_0]$ for each mice group is listed in Table 2. A variant of an embedded Runge–Kutta formulation RK5(4), adjusted to satisfy a non-negativity condition on a solution, is implemented in PK-module as a numerical proprietary C++ code to solve the system Eqs. (1)–(3). Previously, Penn solved^10,13 a similar system using MATLAB internal solver. Simphotek’s PK-module has been validated against MATLAB calculations⁴⁵ resulting in a perfect match for the published parameters. The equations are as follows:^10,13

Table 2

Summary of the results for all mice groups. Penn data are taken from Ref. 10. The optical power (column 6) used for the Dosie™ simulations equals the area of the 10-mm diameter emitting disc (0.7854  cm2) times the Penn in-air fluence rate (column 2). The [O12]rx simulated results for 1.00-, 0.50- and 0.25-mm voxels are shown in columns 8, 9, and 10, respectively. The values in parentheses for columns 8, 9, and 10 are the mismatch to the Penn values (column 7).

In these equations, $g$ is the oxygen supply rate to the tissue, $\delta$ is the low-concentration correction, $\beta$ is the oxygen quenching threshold concentration, $\sigma$ is the specific photobleaching ratio, and $\xi$ is the macroscopic maximum oxygen supply rate. The parameter values specific for the PS, BPD, are given in Table 1. The initial oxygen concentration, [${}^3{\mathrm{O}}_2$] $(x,y,z,t=0)$, is assumed to be $40\mu \mathrm{M}$ for all simulations.

2.4.

Graphical User Interface

The calculations are broken into individual sessions. Using the GUI created in Dosie™, users can create sessions, archive current sessions, browse and search for previous sessions, define PDT simulation parameters, and launch simulations. The visualization component in Dosie™ includes tools to render 2-D and 3-D simulation and geometry data. It enables users to visualize and analyze dose maps resulting from MC and PK calculations (see Fig. 2) to visually inspect under- and overdosed regions. Visualization view has a control panel that lists all dose maps currently available for visualization for each completed session or allows uploading maps from the file system and creating groups of maps for intersession analysis.

Fig. 2

Dosie™ GUI screenshots: MC fluence rate results using 0.50-mm voxels for an incident fluence rate of $50{\mathrm{mW}/\mathrm{cm}}^2$. Volume rendering is done to visualize 3-D fluence rate map in the upper left view, whereas other views show three 2-D cross-sections (slices) of the map—sagittal, coronal, and transverse views.

The GUI shows four conventional views for visualization of the target data: one model view with a 3-D map and three 2-D slice views of the individual 2-D slices resulting in a 2-D sagittal (side) view, a 2-D coronal (front) view, and a 2-D transverse (horizontal) view. The slicing planes can be moved to visualize different cross-sections of the 3-D calculations shown in model view. This enables the user to investigate the fine, near microscopic details of the treatment in the target.

In this study, we will focus on the MC fluence rate and ${[{{}^1\mathrm{O}}_2]}_{rx}$ singlet oxygen dose—the quantity that showed a strong correlation with PDT outcome in recent animal studies.¹⁰ The voxel-based MC simulations of light transport run over $100\times 100\times 100$ grids of either ${(1\mathrm{mm})}^3$, ${(0.5\mathrm{mm})}^3$, or ${(0.25\mathrm{mm})}^3$ voxels. From the resulting fluence rate maps, the subsequent PK calculations determine the singlet oxygen dose or ${[{{}^1\mathrm{O}}_2]}_{rx}$ for each voxel. We show that Dosie™ and the Penn preliminary calculations give nearly equivalent results.

3.1.

Experimental Procedures Summary

Superficial PDT was done at Penn using mice having radiation-induced fibrosarcoma (RIF) tumors. The cancerous mice were injected with the PS benzoporphyrin derivative (BPD, Visudyne^®) followed by PDT treatments.¹⁰ The authors define a cure index for the mice and analyze the correlation of the cure index to the three types of doses: light fluence dose, PDT-dose, and singlet oxygen dose. Penn denotes the singlet oxygen dose as the reacted singlet oxygen concentration, or ${[{{}^1\mathrm{O}}_2]}_{rx}$.

Penn researchers divided the mice into 15 groups, each containing 3 to 5 mice. PDT was delivered to 14 groups. The 15th group was the control group that did not receive either BPD or PDT. RIF cells were cultured and $30\mu \mathrm{l}$ (at $1\times 10{}^7\text{cells}/\mathrm{ml}$) were injected in the right shoulders of 6- to 8-week old female C3H mice. The resulting tumors were treated with PDT when the tumors were $\sim 3$ to 5 mm in diameter and $\sim 3\text{-}\mathrm{mm}$ thick. The mice were injected with BPD 3 h before treatment. The mice were then treated with a 10-mm diameter, 690 nm, collimated light beam from an 8 W diode laser (B&W Tek Inc., Newark, Delaware). The incident in-air fluence was 30 to $350{\mathrm{J}/\mathrm{cm}}^2$, which corresponded from 50 to $150{\mathrm{mW}/\mathrm{cm}}^2$ fluence rate.

Simulations of singlet oxygen dose require knowledge of the BPD concentration in the tumor, the optical properties of the tumor, and the initial ground-state oxygen concentration, $[{{}^3\mathrm{O}}_2]$, in the tumor. Penn did fluorescence measurements on each mouse using a custom multifiber spectroscopic contact probe before treatment to determine the initial BPD concentration.⁴⁶ The average initial BPD concentrations, ${[\mathrm{BPD}]}_{\mathrm{pre}}$, for each group of mice, are listed in Table 2. Tissue optical properties for the tumors, ${\mu }_{\mathrm{a}}$ and ${\mu }_{\mathrm{s}}^{\prime }$, were also measured using a multifiber contact probe.⁴⁴ For the simulations, it is assumed that the tissue has a single set of optical properties, ${\mu }_{\mathrm{a}}=0.69{\mathrm{cm}}^{-1}$ and ${\mu }_{\mathrm{s}}^{\prime }=11{\mathrm{cm}}^{-1}$, which are typical values for tumors in all mice groups.

To determine a cure index, CI, Penn measured tumor dimensions daily for 14 days after PDT and calculated tumor volumes using the formula:⁴⁷ $V=\pi a^{2}b/6$, where $a$ and $b$ are diameters of the width and length axes. A tumor regrowth rate for each tumor was determined by the best fit to an exponential of the form, $e^{kt}$, where $t$ is the time (in units of days) after PDT. Then CI is obtained as follows:¹⁰

3.2.

Calculations of Fluence Rate and ${[{{}^1\mathrm{O}}_2]}_{rx}$

The bulk of Penn simulations have been performed in MATLAB. To determine fluence rate versus tissue depth at the center of the light beam, Penn used the previously determined, one-dimensional (1-D) analytical expression to approximate MC results at the center of a 10-mm diameter treatment light beam:⁴⁸

Full 3-D MC simulations of fluence rate versus depth were done for the same mice data using MC-module of Dosie™. The 3-D simulations of Dosie™ allow the visualization of the results in all three dimensions rather than just the 1-D results obtained by Penn only at the center of the treatment beam. Each MC simulation launched 20 million collimated light rays directed perpendicular to the target surface from a 10-mm-diameter disc (the area is $0.7854{\mathrm{cm}}^2$).

4.1.

Fluence Rate Calculations

In Figs. 2 and 3, we show Dosie™’s partial screenshots of MC light transport results for the fluence rates inside the target volume using 0.50 mm voxels, ${(0.50\mathrm{mm})}^3$ cubes, an incident fluence rate of $50\mathrm{mW}/{\mathrm{cm}}^2$ assumes that the tissue has a single set of optical properties, ${\mu }_{\mathrm{a}}=0.69{\mathrm{cm}}^{-1}$ and ${\mu }_{\mathrm{s}}^{\prime }=11{\mathrm{cm}}^{-1}$. The Dosie™ MC runtime was 21 s for 20 million photons.

Fig. 3

YZ slice of Dosie™’s MC fluence rate results using 0.50-mm voxels through the center of the target volume for an incident fluence rate of $50{\mathrm{mW}/\mathrm{cm}}^2$.

There are four views in the MC output graphics in Fig. 2. In the upper left, a portion of the target volume with significant fluence values is rendered. The view is interactive allowing the user to rotate and magnify the target volume. In this example, the volume is tilted to view, where the light is incident on the target surface. The other three views are an XY slice (upper right), an YZ slice (lower left), and an XZ slice (lower right). The positions of each slice can be moved by the user within the target volume. Each slice can be expanded to fill the graphical image space. The graphics showing the expanded YZ slice (with added labels) of fluence rate values is shown in Fig. 3.

Note that although the incident fluence rate is $50{\mathrm{mW}/\mathrm{cm}}^2$, the peak fluence rate induced just below the surface is $168{\mathrm{mW}/\mathrm{cm}}^2$ (as shown on the scale at the right), which is 3.36 times larger than the incident fluence rate. This enhanced fluence rate occurring just under the surface of the target volume is typical for scattering media, where light can be scattered in all directions and can undergo total internal reflection under the target surface, which adds significantly to the magnitude of the primary incident light.

A comparison of the analytical formula, Eq. (5), to MC calculations of Dosie™ for the ratio (fluence rate in tissue)/(incident fluence rate in air) versus depth in tissue is shown in Fig. 4. For MC, the results at the center of the beam using 1.0, 0.5, and 0.25 mm grid steps are done for an incident fluence rate of $50{\mathrm{mW}/\mathrm{cm}}^2$. For the 1.0-mm grid, the first voxel (centered at 0.5 mm depth) misses the initial higher ratio of fluences at $\sim 0.25\mathrm{mm}$ depth in the target. The 0.5-mm and 0.25-mm grids more accurately resolve the fluence ratios for the smaller depths. Consequently, the MC results for such grids are approximately identical to the analytical approximation.

Fig. 4

The ratio (fluence rate in tissue)/(fluence rate in air) versus depth in tissue that compares an analytical fit⁴⁸ to Dosie™’s MC results.

4.2.

PK Calculations

After the MC calculation of the fluence rate for each voxel of the target, the PK-module performs PK numerical calculations to find ${[{{}^1\mathrm{O}}_2]}_{rx}$ for each such voxel of the target volume for a treatment time of 600 s. The YZ slice at the center of the full 3-D simulation is shown in Fig. 5. The maximum ${[{{}^1\mathrm{O}}_2]}_{rx}$ just below the surface at the center of the beam is $711\mu \mathrm{M}$ or 0.711 mM. The Dosie™ PK runtime was 12 s using 4 CPU processor cores.

Fig. 5

YZ slice of Dosie™’s PK ${[{{}^1\mathrm{O}}_2]}_{rx}$ results through the center of the target volume shown for an incident fluence rate of $50\mathrm{mW}/{\mathrm{cm}}^2$ for 600 s and for 0.50-mm voxels.

To insure each tumor gets a sufficient dose through the full 3-mm tumor thickness, we determined the singlet oxygen dose at a distance of 3 mm below the surface of the target at the center of the illumination. For PK simulations with a 1-mm grid, one needs to look at the third and fourth layers of 1-mm-thick voxels at the center of illumination and to average the two results to get the singlet oxygen dose 3 mm below the target surface. For a 0.50-mm grid, as shown in Fig. 5, this was done by looking at the sixth and seventh layers and by averaging the two results to get the singlet oxygen dose at the interface 3 mm below the target surface. For a 0.25-mm grid, this was done by averaging the values from the 12th and 13th layers.

The final Penn and Dosie™ simulation results for all the mice groups are shown in Table 2. Columns 8, 9, and 10 include Dosie™ PK calculated values of ${[{{}^1\mathrm{O}}_2]}_{rx}$ at 3 mm below for different grid resolutions: 1.0, 0.5, and 0.25 mm voxels, respectively. The numbers in parentheses show the mismatch with the corresponding Penn results (see column 7). Except for mice groups 1 and 2, the differences between the PK results and the Penn results are less than about 2%.

Plots showing the differences between the PK results and the Penn results are shown in Fig. 6. In going from 1.00-mm voxels to 0.50-mm voxels, the differences with the Penn results change on the order of 1%. However, in going from 0.50-mm voxels to 0.25-mm voxels, the differences only change $\sim 0.2\%$, indicating that reducing the grid step below 0.50 mm is unnecessary.

Fig. 6

Comparison of Dosie™’s PK results to Penn results for 1-mm, 0.50-mm, and 0.25-mm voxels.