Photodynamic therapy outcome modelling for patients with spinal metastases: a simulation-based study

Spinal metastases often occur in the advanced stages of breast, lung or prostate cancer, resulting in a significant impact on the patient’s quality of life. Current treatment modalities for spinal metastases include both systemic and localized treatments that aim to decrease pain, improve mobility and structural stability, and control tumour growth. With the development of non-toxic photosensitizer drugs, photodynamic therapy (PDT) has shown promise as a minimally invasive non-thermal alternative in oncology, including for spinal metastases. To apply PDT to spinal metastases, predictive algorithms that optimize tumour treatment and minimize the risk of spinal cord damage are needed to assess the feasibility of the treatment and encourage a broad acceptance of PDT in clinical trials. This work presents a framework for PDT modelling and planning, and simulates the feasibility of using a BPD-MA mediated PDT to treat bone metastases at two different wavelengths (690 nm and 565 nm). An open-source software for PDT planning, PDT-SPACE, is used to evaluate different configurations of light diffusers (cut-end and cylindrical) fibres with optimized power allocation in order to minimize the damage to spinal cord or maximize tumour destruction. The work is simulated on three CT images of metastatically involved vertebrae acquired from three patients with spinal metastases secondary to colorectal or lung cancer. Simulation results show that PDT at a 565 nm wavelength has the ability to treat 90% of the metastatic lesion with less than 17% damage to the spinal cord. However, the energy required, and hence treatment time, to achieve this outcome with the 565 nm is infeasible. The energy required and treatment time for the longer wavelength of 690 nm is feasible (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\sim }\,40$$\end{document}∼40 min), but treatment aimed at 90% of the metastatic lesion would severely damage the proximal spinal cord. PDT-SPACE provides a simulation platform that can be used to optimize PDT delivery in the metastatic spine. While this work serves as a prospective methodology to analyze the feasibility of PDT for tumour ablation in the spine, preclinical studies in an animal model are ongoing to elucidate the spinal cord damage extent as a function of PDT dose, and the resulting short and long term functional impairments. These will be required before there can be any consideration of clinical trials.


Scientific Reports
| (2021) 11:17871 | https://doi.org/10.1038/s41598-021-97407-z www.nature.com/scientificreports/ or nerve root compression which may result in motor or sensory deficits and bowel or urinary incontinence in severe cases, significantly impacting quality of life. Patients with advanced metastatic disease are typically poor surgical candidates, and receive treatment aimed primarily at palliative pain control. Current treatment options for spinal metastases include localized, targeted approaches such as surgical stabilization (for cases with vertebral instability requiring cord decompression) 2 and radiation and thermal therapies 3,5 as well as systemic treatments such as chemotherapy, immunotherapy and bisphosphonates 4 . While systemic treatments like bisphosphonates have shown benefit in relieving metastatic bone pain and delaying complications, they usually cause adverse effects on the gastrointestinal and haematopoietic system 6 . Conventional radiation therapy has been shown to achieve an overall response rate (palliation of pain symptoms) of approximately 60%, and complete abrogation of pain (complete response rate) in only ∼25% of patients 7 . In contrast, focal radiation therapy (i.e. stereotactic body radiation therapy, SBRT) is a well-established non-invasive approach that precisely targets metastatic lesions in bone and provides pain relief in the majority of patients (> 90% ). Thibault et al. 8 reported that patients undergoing SBRT had sustained pain relief of 86% and local control of 88% (using CT criteria) at a median follow-up of 21 months. Yet, SBRT's repeated use is limited by toxicity to the spinal cord (radiation-induced myelopathy) and an incidence of fracture post treatment [9][10][11][12] . Additionally, SBRT requires sophisticated hardware and software for treatment planning, patient setup, and careful patient selection 13 , which currently limit its widespread use.
Radiofrequency ablation (RFA) is also utilized clinically for localized treatment of spinal metastases 5 . RFA is a thermal modality that utilizes high-frequency alternating current (by placing needle electrodes into the surrounding tissues) to heat and eventually ablate tumours. RFA may also be coupled with vertebral cement augmentation (VCA) to provide vertebral stabilization and extended pain relief. Mayer et al. 14 recently reported that 80% of patients who underwent bipolar RFA with VCA achieved favourable pain relief (3 points reduction on the visual analogue scale, VAS) at a mean follow-up of 3.4 months. Neurologic injuries during RFA of spinal metastasis may occur if performed too close to critical structures, as such RFA generally is limited for tumours in the posterior vertebral body an in cases with a breach of the posterior vertebral body wall 5 .
Photodynamic therapy (PDT) is an emerging non-thermal minimally invasive modality that offers the potential to precisely target spinal metastasis. In preclinical rodent models, PDT has been shown to successfully ablate spinal metastases and improve vertebral mechanical stability with increased osteoid formation, particularly when combined with systemic bisphosphonates [15][16][17][18][19] . A recent phase I trial further demonstrated the safety and feasibility of BPD-MA mediated (Benzoporphyrin derivative mono-acid photosensitizer, trade name: Visudyne, Novartis, QC, Canada) PDT as a tumour-ablative adjunct modality prior to VCA in patients with spinal metastases 20 . The study evaluated various treatment parameters including the energy delivered and drug-light interval. The results suggested that vertebral PDT as an adjunct to VCA is safe from a pharmaceutical and neurological perspective. The 50 J cm −1 and 100 J cm −1 treatment groups showed a clinically significant reduction in pain 20 .
A potential advantage of PDT over radiation or thermal therapies is the ability to repeatedly treat the same site 21 , without risking toxicity to the spinal cord or other critical structures. High light scattering in intact vertebral bone confines the excitation photons and limits the light reaching the spinal cord. In cases of posterior vertebral body metastatic involvement, particularly if the cortical shell is compromised, establishing light scattering through the tumour and remaining bone is critical to safety. For this, personalized PDT treatment planning could be used to reduce any potential injury risk in targeting the malignancy and avoiding impact to adjacent spinal cord or nerve roots. Yet, PDT treatment planning is challenging due to the lack of established 3D modelling tools and framework for optimizing the source configuration. Such planning requires visualizing the final fluence distribution and evaluating the quality of the resulting treatment plan in a highly heterogeneous geometry.
Here, we investigate the feasibility of BPD-MA mediated PDT 22 in patients with spinal metastases by presenting a framework for PDT treatment modelling and simulation leading towards systematic PDT planning. We incorporate the complexity and heterogeneity of the 3D tumour geometry using original pre-treatment CT datasets and contours for patients treated with SBRT. Using an interstitial PDT (iPDT) planning optimization tool called PDT-SPACE 23,24 , the optimal source power allocation for two types of light sources embedded within the metastatic lesion is determined, demonstrating the ability to tailor the light dose (3D fluence distribution) to the tumour geometry while minimizing damage to the spinal cord using PDT dose threshold values (considering a PDT threshold model 25 ) derived from preclinical models. We simulate the attainable efficacy of PDT at two different activation wavelengths, 690 nm and 565 nm. Both of these wavelengths are peak absorbance bands for BPD-MA 26 . We report the predicted tumour coverage and damage to the spinal cord in both cases.

Methods
Treatment planning framework. Computer tomography (CT) images of patients with spinal metastases were used to generate the virtual PDT treatment plans used in this work (Fig. 1). Patient identifiers were removed prior to using these imaging datasets from three patients with spinal metastases treated with stereotactic body radiation therapy (Department of Radiation Oncology, Sunnybrook Health Sciences Centre, Toronto, ON, Canada). This is a retrospective study utilizing the imaging datasets. Informed consent from the patients was obtained. The study followed all applicable guidelines and regulations. The study was approved by the Sunnybrook Health Sciences Ethics Review Board (Toronto, ON, Canada). Table 1 describes the three cases along with their locations. Guided by the contours from the SBRT plan, target structures and organs at risks were segmented using ITK-SNAP 27 , including the metastasis, spinal cord, normal surrounding bone, and muscle. Note that the latter two structures were manually segmented as they were not explicitly contoured. Other surrounding structures were omitted as they are not expected to be significantly impacted by PDT, because they were distal to the target and were exposed only to a low photon density, or they have a limited vasculature and hence do PDT-SPACE 23,24 , an open-source optimization software for iPDT treatment planning, was used to optimize the power allocation across the source geometry, including a scenario of 4 cut-end fibres with a 400 µ m core diameter and NA = 0.22 (placed into the spinal metastasis, pointing diagonally at 45 • away from the spinal cord) as well as a scenario of two cylindrical diffusers with a 500 µ m radius inserted diagonally into the spinal metastases. Different weights were applied (called tumour weight) to drive the optimization algorithm towards favouring tumour targeting at the expense of critical structures such as the spinal cord or favouring no damage to the spinal cord at the expense of less damage to the tumour. Dose-volume histograms were generated for each tissue type by PDT-SPACE to assess the given PDT plan's quality. The final PDT treatment plan at two different activation wavelengths ( 1 = 690 nm and 2 = 565 nm) with the corresponding iso-fluence contours are visualized using ParaView 5.6.0 29 . Table 2 summarizes the tissue optical properties (absorption coefficient µ a , scattering coefficient µ s , anisotropy g, and refractive index n) at the two activation wavelengths. The optical properties of the osteolytic tumours, spinal cord, bone, and muscle are based on literature values [31][32][33][34][35][36] (interpolated for 565 nm), while those of the sclerotic spinal metastasis are assumed to be similar to the bone. Table 3 summarizes the PDT dose constraints considered in the current model. The spinal Figure 1. Overview of PDT treatment planning framework for a female patient with a metastatic T8 sclerotic lesion: (a) extraction of contours from stereotactic body radiation therapy plan using the Dicompyler tool, including the clinical target volume or metastasis (green) and spinal cord (blue), (b) segmentation of the original CT dataset using ITK-SNAP to delineate the metastasis (purple), normal bone (green), spinal cord (red), and muscle (yellow) at the T8 level, (c) generation of the 3D mesh geometry using MeshTool for light dose simulation with FullMonte Methods.  Table 2. Tissue optical properties for each region in the 3D model at two activation wavelengths ( 1 = 690 nm and 2 = 565 nm). a Assumed to have similar tissue optical properties as the bone. www.nature.com/scientificreports/ cord is assigned 1/10 the PDT dose threshold of normal bone, while the metastasis is assumed to be 10× more resistant than normal bone. Based on the photosensitizer uptake data (BPD-MA at 15 min with a dose of 2 mg/ kg) from an earlier preclinical study 37 , the threshold fluence, threshold (J cm −2 ) , at the boundary of necrosis is computed using the PDT threshold formula 38  Ethics approval and consent to participate. This is a retrospective study utilizing imaging datasets obtained from patients for whom informed consent was obtained. The study followed all applicable guidelines and regulations. The study was approved by the Sunnybrook Health Sciences Ethics Review Board (Toronto, ON, Canada).

Results
Using the final 3D mesh geometries and the dose constraints specified in Table 3, two light source configurations were evaluated to demonstrate the use of the 3D PDT treatment planning framework for the first T8 metastatic lesion in Table 1. For this model, the volumes of the different tissue types in the 3D model are as follows: spinal metastasis is 23.21 cm 3 , the spinal cord is 4.72 cm 3 , normal bone is 31.14 cm 3 and the muscle is 171.84 cm 3 . The corresponding iso-fluence contours based on the optimized power allocation are shown to compare the quality of the predicted treatment outcomes for all scenarios. Additionally, we present dose-volume histograms (DVHs) to see effect of both wavelengths on the spinal cord damage. In the first set of optimization simulations the power allocation was optimized to attain at least 90% necrosis in the metastatic lesion. To achieve this, the metastasis weighting parameters were automatically adjusted in the optimization framework. In the second set of optimization simulations, the weighting parameters were varied to prevent any necrotic damage to the spinal cord while maximizing the impact on the metastatic lesion.
Simulations for cut-end fibre light delivery. Figure 4 shows the iso-fluence contours at the centre slice of the metastatic lesion (see Fig. 2) for both the 690 nm and 565 nm treatments. The yellow contour shows the necrotic threshold dose of the spinal cord (0.02 J mm −2 ), while the orange contour indicates the necrotic threshold for the metastasis (0.4 J mm −2 ). For the 690 nm treatment, the overall energy needed is estimated to be 1247.6 J. However, the potential overall damage to the spinal cord was around 4.7 cm 3 , which is almost the entire segmented spinal cord tissue. The 565 nm treatment shows significantly reduced damage to the spinal cord (2.4 cm 3 , representing ∼ 50% damage reduction). However, the estimated required energy for the 565 nm treatment is infeasible at 167.3 kJ. Figure 3 compares the spinal cord dose-volume histograms of the two treatments. While the 565 nm wavelength can significantly reduce the damage to the spinal cord, it would require a much longer treatment duration.
Simulations for cylindrical diffusers light delivery. Similar to the cut-end fibres configurations, Fig. 5 shows the iso-fluence contours for a 2-cylindrical diffusers scenario at the two wavelengths. A similar trend is seen in that the 565 nm treatment significantly decreases the damage to the spinal cord. The damage reduction is more substantial than in the cut-end fibres scenario; the 565 nm wavelength yields a damage volume of only 1.4 cm 3 (70% damage reduction). However, the required energy (9.3 ×10 3 kJ) would necessitate a much longer treatment duration.

Summary of simulations.
The above simulations were executed for the three metastatically involved vertebrae shown in Table 1. The results of the simulations for the metastatic lesion (volume and damage fraction (%)), the spinal cord (volume at risk and the damaged volume) and the required energy at 690 nm and 565 nm for a treatment optimized for 90% tumour destruction are shown in Table 4. Notice that for the small metastatic lesion    To simulate a palliative treatment goal, the optimization was run focusing on spinal cord preservation. In this, the source placement remained fixed but the tissues' weighting parameters in PDT-SPACE were adjusted to attain near-zero damage on the spinal cord. Table 5 reports the limited (0.04% to 15%) attainable bone metastasis destruction fraction under this scenario.

Discussion
PDT is increasingly being utilized as a focal oncology therapy, with recent approvals for prostate and brain tumours in various jurisdictions including Europe, the USA and Japan. This is due to the non-toxicity of the photosensitizer drugs in the absence of light. Additionally, the short penetration of light in biological tissues makes PDT less invasive to the surrounding critical structures, making repeatability of the treatment in case of recurrent pain possible. In considering PDT for the treatment of spinal metastasis, the advantages of a localized non-thermal therapeutic approach must be weighed against any potential risk to the spinal cord. Predictive algorithms that minimize and quantify the risk of spinal cord damage based on each patient's anatomy and tumour volume will be essential to achieve a broad acceptance for PDT in these critical situations and guide appropriate patient selection.
Modelling of light propagation in complex geometries with heterogeneous optical property distributions requires high-resolution clinical imaging and precise delineation of these structures. In this study, the tissue regions of interest were manually delineated (as is the convention in the clinical workflow of radiation therapy); however, several automated and semi-automated approaches for the segmentation of spinal images have been proposed and can be implemented to streamline the workflow for clinicians 40 .
The simulations provided here are based on assumed photodynamic threshold data from prior publications 38 considering a least favourable responsivity difference between the spinal cord and the target metastasis, given the current lack of tissue specific PDT responsivity threshold values from the literature. The spinal cord mainly consists of a peripheral region that contains neuronal white-matter tracts. Internal to this region, there are neuronal grey-matter-like structures. White matter is rather resistant to PDT, while grey-matter is less resistant 38 . To be conservative, we assume in this work that the PDT threshold dose for the spinal cord is similar to that of grey-matter. In general, the threshold values of normal tissues are lower than those of malignancies, due to the latter's ability to neutralize higher reactive oxygen species (ROS) concentration 41 . As such, the potential damage to the spinal cord reported in this study is likely an overestimate when optimizing to a level of 90% destruction of the metastatic lesion. Similarly, the maximum attainable malignancy destruction is likely underestimated in this model when optimizing based on complete spinal cord preservation. Experimental determination of PDT threshold values is underway in ongoing preclinical rat models 42 , in which the extent of the treatment effect is assessed through histology 38 . Preliminary results on n = 3 T10 vertebrae cases show that the photodynamic threshold ranges between 0.61-2.36 ×10 18 photons/cm 3 , which is higher than the threshold assumed in this study. Using established tissue optical properties 31 , the PDT dose gradients per source can be calculated independent of the uncertainty in the PDT threshold values. However, the absolute power requirements, the total delivered energy and the energy distribution between the multiple sources would vary.
The following conclusions are not affected by the unknown PDT threshold values. Exploiting the strong light attenuation of the 565 nm excitation light aids in limiting damage to the spinal cord; however this required 134 and 3196 times more total energy to achieve the same fraction of tumour destruction for the cut-end and the cylindrical emitters, respectively (Figs. 4,5). The higher optical energy delivery is partially due to the 5-times lower molar extinction coefficient at 565 nm versus 690 nm. Assuming the same power delivery, the added energy translates directly into a proportional increase in exposure time. Using cylindrical fibres instead of the cut-end fibres resulted in a 130% increase in total energy requirements (2.91 kJ versus 1.25 kJ). However, the power per optical fibre can be higher when using cylindrical diffusers versus cut-end fibres. Assuming a power delivery of 200 mW cm −1 across the diffuser length, treatment time of ∼ 41 min would be required. It is worth mentioning, however, that depending on the blood volume in the tissues, the absorption coefficient, µ a , can vary significantly 35 . Our simulation results showed that for high blood content in the tumour, a solution may not be  43 .
In any case, the PDT-SPACE automated plan can be tailored to clinical plans focused on high levels of tumour destruction (Table 4) or accommodate less comprehensive, palliative treatment, which prioritizes preservation of the spinal cord volume at risk (Table 5). For example, Table 5 shows that treatment at 565 nm for case number 1 with around the same total energy (0.36 kJ) given to the 690 nm treatment-that kills 90% of the tumour (case 1 in Table 4)-reduces the metastatic lesion volume by at least 15% with no damage to the spinal cord.
Ultimately, treatment planning using the PDT-SPACE platform may allow a combination of different treatment wavelengths to be assessed. Future investigations may consider a 565 nm source positioned in the posterior aspect of the vertebral body (closer to the spinal cord), and a 690 nm source in a more anterior location. This would reduce the overall optical energy requirements and treatment time, enabling a safe, feasible and effective approach to PDT tumour ablation in the spine, specifically when there is a risk of cortical breach of the vertebral body towards the spinal cord. Another possible future direction is the use of a combination of BPD-MA and lipid-anchored BPD to achieve photo-damage at a lower light dose by targeting mitochondria, ER and lysosomes simultaneously as shown by Rizvi et. al. 44 in an in-vitro ovarian cancer 3D model. In principle, this would allow achieving metastatic damage at lower energy levels even for the shorter wavelength treatment. However, more pre-clinical and clinical studies are needed to determine the feasibility of this approach and to evaluate the lower PDT dose thresholds required by PDT-SPACE.  Optimized treatment plans for wavelengths of (a) 690 nm and (b) 565 nm using a 2-source configuration (2 cylindrical diffusers with a radius of 500 µm). The resulting iso-fluence contours are overlaid on the 3D model, highlighting the the minimum necrotic threshold within the metastatic region (0.4 J mm −2 ) as per the colour map shown. The total energy at the 690 nm wavelength ( E total = 2914 J) was distributed across the 2 light-emitting cylinders as: E 1 = 1981.8 J and E 2 = 932.4 J, for the left and right sources respectively. At 565 nm, the total energy was E total = 9.3 × 10 6 J, distributed as: E 1 = 8.35 × 10 6 J and E 2 = 0.95 × 10 6 J.