The clinical radiobiological parameters (i.e. \(T_}}\), \(T_}}\), \(\alpha\), \(\alpha \beta\)) for twenty-five patients were derived using the diagnostic (FDG PET/CT) and therapeutic (90Y-SIRT PET/CT) images to compute the BED map, TCP and FDG PET total liver lesion glycolysis (TLG) for each patient. Furthermore, the relationship between the radiobiological parameters and the calculated dosimetric quantities was investigated.

The data for twenty-five patients with liver malignancy including pancreatic neuroendocrine tumours (PNET), colorectal cancer (CRC), pancreatic ductal adenocarcinoma (PDAC), small bowel neuroendocrine tumours (SBNET), hepatocellular carcinoma (HCC), neuroendocrine Carcinoma (NEC) and other metastatic tumours who were treated with [90Y]-SIR-Spheres between July 2019 and April 2021 were used in this study. All patients gave informed consent at the time of the procedure for their clinical and image data to be used for further research, education, training, and audit. For each patient a complete imaging set suitable for lesional analysis was available. Similar to our previous study [22], an individual imaging set consisted of baseline FDG PET/CT (acquired ≤ 28 days prior to radioembolisation), 90Y-SIRT PET/CT (acquired within 24 h of radioembolisation), and follow-up FDG PET/CT (acquired ≤ 80 days post-radioembolisation). Treatment time is the time interval between the baseline and follow-up FDG, and each patient’s treatment time is listed in Table 1. All patients underwent pre-treatment interventional arterial mapping and abdominal Technetium‐99 m macroaggregated albumin ([99mTc]MAA) single-photon emission computed tomography SPECT/CT imaging prior to treatment to assess the lung shunt fraction and possible extrahepatic uptake to determine the feasibility, safety, and number of injections required for selective treatments.

All imaging data were acquired using similar scanners and protocols that were used in our previous study [22]. Images were acquired on a Siemens Biograph mCT-S (64) PET/CT system (Knoxville, TN, USA) with 550 picosecond timing resolution time-of-flight (ToF) capabilities, an axial field of view of 21.8 cm and 78 cm crystal ring diameter. Images (with voxel size = 4.072 × 4.072 × 2 mm) were reconstructed using the standard OSEM (with 3 iterations and 21 subsets) reconstruction method in conjunction with ToF modelling and point spread function recovery. Our ‘low-dose protocol’ (i.e. as two 10 min frames over the liver and reconstructed with 3i21s and a 5 mm Gaussian filter) was used to acquire baseline and follow-up PET/CT data. The quantitative liver 90Y PET/CT data were reconstructed with 1i21s with no filtering. Siemens’ Intevo-6 or Symbia.T16 were used to acquire the MAA planning SPECT/CT data with low energy parallel hole collimators and standard CT-based attenuation correction were used for reconstruction. To avoid breakdown of the 99mTc-MAA in vivo, in all cases acquisition was performed within 1 h following implantation. The [90Y]-microsphere dosimetry navigator software (RapidSphere®) within a commercial platform (Velocity, Varian Medical Systems, Palo Alto, USA) was used to analyse the patient images including the absorbed dose. The absorbed dose was calculated using the local deposition method [30,31,32].

The 90Y-SIRT PET/CT images were first registered to the pre‐ and post‐implantation FDG PET/CT images using the deformable image registration package provided. The Velocity deformable registration algorithm applies multiresolution free‐form deformations and an intensity-based B-spline multipass algorithm to provide a high-level deformable image registration accuracy [32, 33]. The deformable registration provided a one‐to‐one correlation between voxels on different images and time points, allowing for mapping of anatomical data and structure sets from [90Y]-microsphere PET/CT, pre-treatment FDG PET/CT to follow‐up FDG PET/CT.

The 90Y dosimetry navigator was used to calculate the lesion absorbed dose distribution from the 90Y PET/CT images. Next, the dose and the normalised standard uptake value (SUV) Volume Histogram of the 90Y PET/CT and registered pre- and post-treatment FDG images were computed and imported into MATLAB (R2020a) software. In our previous in vitro study [23], it was demonstrated that the measured radiation survival by metabolic cell assay is comparable to that of clonogenic assays. Additionally, the metabolic assay measures all viable cells thus representing cells from a true tumour population rather than just clonogenic cells [34, 35]. The FDG scan uses a glucose analogue and is the most commonly used PET tracer to assess tumour metabolism. Due to increased glucose metabolism in most types of tumours, the FDG PET is widely used clinically for tumour imaging [36]. The FDG uptake in PET imaging is a measure of the tissue glucose metabolism and is usually high in high-grade tumours (e.g. maximum SUV = 11) and relatively low in low-grade (e.g. maximum SUV = 7) tumours [37]. Additionally, in our previous study [23] we demonstrated that in vivo metabolic imaging such as FDG PET can be used to assess the metabolic dose response as well as prognostic factors for radioembolisation of liver metastases from colorectal cancer. Furthermore, the change in tumour voxel SUV ratio at a specific dose level from a serial FDG PET/CT imaging (i.e. pre- and post-treatment FDG images) can be used to model the tumour voxel dose response [38]. In this study, similarly we have proposed that the ratio of voxel SUV from pre- and post-treatment FDG images should represent the metabolic radiation survival fraction (SF) due to [90Y]-microspheres irradiation. Therefore, the ratio of voxel SUV from pre- and post-treatment should represent the metabolic radiation survival fraction (SF) due to [90Y]-microspheres irradiation. In addition, to account for dose heterogeneity the SF is calculated based on the dose volume histogram (DVH) by:

$$} = \mathop \sum \limits_ \frac }} }} }\left( } \right)$$

(1)

$$}\left( } \right) = \frac}_} - }}} }}}_} - }}} }}$$

(2)

where \(V_\) is the tumour volume and \(V_\) is the sub volume corresponding to 90Y dose bin \(D_\) on the DVH. Hence using Eqs. 1 and 2, the survival fraction of cancer cells in a tumour volume with an initial volume of \(V_\) can be estimated by computing the ratio of voxel SUV from pre- and post-treatment FDG images. Next, the voxel SF (VSF) data was fitted to the GLQ model (Eq. 3, using MATLAB software) to estimate the radiobiological parameters, \(T_}}\), \(T_}}\), \(\alpha\), \(\alpha \beta\). The GLQ fit is commonly used in the field of radiobiology to derive radiobiological parameters both in vitro and in vivo studies [16].

$$} = e^ + \gamma \left( }} - T_}} } \right)} \right)}}$$

(3)

$$} = e^ - \frac }}}} }} \left( T_}} }}} \right) - \left( }} }}}} }}} \right)} \right)}}$$

(4)

$$G = \frac }} \mathop \smallint \limits_^ \dot\left( t \right)}t \mathop \smallint \limits_^ \dot^ \left( } \right) e^ } \right)}} }t^$$

(5)

For further personalisation of the GLQ, the treatment time (\(T_}}\)) was considered to be equal to the critical time (\(T_}}}\)) to include the effect of initial dose rate used for each treatment [23]:

$$T_}} = T_}}} = \frac }}\ln \left( T_}} }}} \right)$$

(6)

where T1/2 is the 90Y half-life (i.e. ≈ 2.7 days), \(\alpha\) is the cell radiosensitivity, \(R_\) is the initial does rate, and \(T_}}\) is the tumour cell proliferation (or repopulation) time. Hence by replacing the \(T_}}\) with \(T_}}}\) and also substituting the tumour cell proliferation (or repopulation) constant, \(\gamma = \ln \left( 2 \right)/T_}}\), in Eq. 4 we can obtain \(T_\), \(T_}}\), \(\alpha\), \(\alpha \beta\) parameters. Furthermore, the GLQ model includes the G-factor (or the Lea–Catcheside factor described by Eq. 5) to account for the kinetics of DNA strand break damage and repair in obtaining the true fraction of surviving cells in an irradiated cell population within the tumour [23]. In Eq. 5, the first integral represents the physical absorbed dose. The integrand of the second integral over \(t^\) refers to the first DNA single-strand break (SSB) of two SSBs needed to cause lethal DNA double-strand (DSB) damage. Also, the integral over \(t\) refers to the second SSB of remaining of two SSBs to cause a DSB. The exponential term reflects the repair and therefore reduction in induction of 2 SSB → DSB process due to decreasing dose rate [23]. Also, \(\mu\) (which is \(\ln \left( 2 \right)/T_}}}\)) is the DNA repair time constant (\(T_}}}\) is the DNA repair half-life which is ≈ 1.5 hr [23]).

The \(R_}}}\) was also calculated by Eq. 7 for further assessment of the treatment outcome.

$$R_}}} = \frac}} }}$$

(7)

Both \(T_}}}\) and \(R_}}}\) were calculated theoretically using Eqs. 6 and 7 for a range of \(\alpha , T_}}\) and \(R_\) shown in Fig. 5a–b and using the clinically driven radiobiological parameters for comparison. Furthermore, to show the impact of the proliferation acceleration on dose rate efficiency during the treatment time, the \(R_} \left( k \right)}}\) was also calculated for when the \(T_}}\) ≈ \(T_}}\) (this can represent the highest acceleration in tumour proliferation rate).

The computed lesion dose for each patient was used as an input to generate the lesion BED map. The BED calculations were performed considering the following methods:

Method1. Due to low dose rate and relatively long treatment time (90Y half-life ≈ 2.7 days), repair of sub-lethal damage might take place during the treatment duration. Additionally, highly proliferating cancers can have short \(T_}}\) (relative to the treatment time) which could impact the survival fraction and TCP. Therefore, treatment with an exponentially decaying source (integrated to fixed treatment time, T) and considering the above parameters can be computed as (Method 1):

$$}_ = D \times } - \frac}} } \right)}}}} }}$$

(8)

$$} = 1 + \left( \lambda }}} \right)\left( } \right)\left( \left( } \right) - \frac\left( \right)T}} } \right)}} }}} \right)$$

(9)

Method 2. This method is the simplified version of Method 1 where treatment time is integrated to \(T \to \infty\). This is the standard BED formulation [39] in the radionuclide therapy dosimetry to calculate the BED:

$$}_ = 1 + \left( }}} \right)\left( } \right)$$

(12)

For each tumour with N voxels, the voxel TCP was calculated first and then the overall expected TCP of tumour was obtained as the product of the expected voxel TCP [40, 41]:

$$}_}}} = e^}} V \cdot }}}$$

(13)

$$} = \mathop \prod \limits_^ \left[ }_}}} } \right]^}$$

(14)

where \(D_}}\) (i.e. \(1 \times 10^\)) and \(V\) are the density of clonogens per cm3 and the voxel volume respectively.

Finally, the \(}_}}}\) of each liver lesion was calculated by multiplying the metabolic tumour volume of that lesion with its corresponding mean SUV and the \(\Delta }_}}}\) was calculated in percentage as [22]:

$$\Delta }_}}} \left( \% \right) = \frac}_}}} - }_}}} }}}_}}} }} \times 100$$

(15)

where \(}_}}}\) and \(}_}}}\) are the \(}_}}}\) values based on pre- and post-treatment FDG PET images.