PBPK-Adapted Deep Learning for Pretherapy Prediction of Voxelwise Dosimetry: In-Silico Proof of Concept

Pretherapy dosimetry prediction is a prerequisite for treatment planning and personalized optimization of the emerging radiopharmaceutical therapy (RPT). Physiologically based pharmacokinetic (PBPK) model, describing the intrinsic pharmacokinetics of radiopharmaceuticals, have been proposed for pretherapy prediction of dosimetry. However, it is restricted with organwise prediction and the customization based on pretherapy measurements is still challenging. On the other side, artificial intelligence (AI) has demonstrated the potential in pretherapy dosimetry prediction. Nevertheless, it is still challenging for pure data-driven model to achieve voxelwise prediction due to huge gap between the pretherapy imaging and post-therapy dosimetry. This study aims to integrate the PBPK model into deep learning for voxelwise pretherapy dosimetry prediction. A conditional generative adversarial network (cGAN) integrated with the PBPK model as regularization was developed. For proof of concept, 120 virtual patients with 68Ga-PSMA-11 PET imaging and 177Lu-PSMA-I&T dosimetry were generated using realistic in silico simulations. In kidneys, spleen, liver and salivary glands, the proposed method achieved better accuracy and dose volume histogram than pure deep learning. The preliminary results confirmed that the proposed PBPK-adapted deep learning can improve the pretherapy voxelwise dosimetry prediction and may provide a practical solution to support treatment planning of heterogeneous dose distribution for personalized RPT.


I. INTRODUCTION
T ARGETED radiopharmaceutical therapy (RPT) has rapidly expanded and become widely accepted as a treatment for several types of malignant tumors, including metastatic castration-resistant prostate cancer (mCRPC) [1].Encouraging results have been obtained in treating patients with mCRPC, which is known for its high morbidity and mortality [2].Several studies have demonstrated the efficacy and safety of RPT in treating mCRPC with 177 Lu-PSMA-ligands [3], [4], [5], [6].Moreover, the Phase III VISION clinical trial validated 177 Lu-PSMA-617, leading to its approval by the U.S. Food and drug administration (FDA) [7].The trial demonstrated that the progression-free survival and overall survival were prolonged compared to the best standard of care in men with advanced-stage PSMApositive mCRPC [8].However, despite the promising early results, dose personalization remains a challenge, making treatment planning vital for optimizing RPT practice [9].
In the course of RPT, PET imaging with similar targeted molecules are routinely utilized in screening patients for treatment suitability [10], [11].After the injection of therapeutic agent, the distribution of these compounds and the consequent radiation dose deposition, i.e., dosimetry, can be usually quantified based on a series of cross sectional SPECT scans [12], [13].Despite the monitoring of the dosimetry, patients are still treated with a fixed radiopharmaceutical activity and interval protocol [14].As a consequence, risks of inadequate tradeoff between therapeutic response and side effects have been a topic of concern.Treatment planning was mandated to personalize RPT similar as external beam radiotherapy.However, a critical prerequisite for treatment planning is the pretherapy prediction of the dosimetry in addition to post-therapy monitoring, which is not feasible in clinical practice.
Physiologically based pharmacokinetic (PBPK) models describe the fundamental principle of uptake procedure of pharmaceuticals, including radio-labeled ligands for PET imaging and RPT [15].They are composed of compartmentbased models capturing the absorption, distribution, metabolism, and excretion (ADME) of the pharmaceuticals within organs, which also describe interactions between the organs during the uptake procedure [16], [18].For instance, a PBPK-based model was able to describe the uptake dynamics for PSMA-directed diagnostic and therapeutic ligands within organs that are clinically hard to interpret [17], [18].
The capability of describing drug uptake procedure allows PBPK models to simulate the biokinetics of radiopharmaceuticals and to predict both organ time-integrated activity coefficients (TIACs) and the internal absorbed doses for RPT [19], [20].Pretherapy measurements, such as PET, scintigraphy, or serum measurements, can be used to customize the setting of the PBPK parameters for individual patients [17], [21], [22], [23].For instance, individualized PBPK model parameters of 90 Y-DOTATATE therapy can be obtained by fitting planar 111 In-DOTATATE scintigraphy and serum measurements to allow individualized internalization rate in the prediction of 90 Y-DOTATATE dosimetry [20].Bayesian framework was established to fit the pharmacokinetic and pharmacodynamic parameters from pretherapy PET/CT activity concentrations, planar scintigraphy, and tumor volumes for the individualization of 177 Lu-PSMA-I&T therapy [24].However, the limited measurements, such as short scanning intervals, make the customization or individualization of PBPK parameters based on pretherapy measurements as well as the extrapolation to therapy compounds still challenging.Complex relations within the ADME processes with their multilevel biochemical interactions require numerous input parameters, and are influenced by a range of individual characteristics, including underlying health conditions.Furthermore, the PBPK-based prediction is also restricted to organwise prediction, without considering the intraorgan heterogeneity.
On the other side, machine learning has been successfully applied for pretherapy prediction of dosimetry [25], [26].Despite the similar restriction in organwise prediction, the potential of artificial intelligence (AI) has been demonstrated in dosimetry prediction.Deep learning as subfield of AI, can effectively approximate complex, multidimensional, nonlinear functions in both spatial and temporal domains [27].Therefore, integrating deep learning into RPT has the potential to offer complementary insights into dose prediction and optimization [25], [26].However, the performance of deep learning models heavily relies on the quantity of available training data [28].This presents a challenge in dose prediction due to the limited availability of post-therapy dosimetry measurements.In addition, applying deep learning to PET imaging and RPT poses an extra challenge due to confounding factors caused by physical and physiological interference [29].Thus, enhancing deep learning performance with limited datasets is crucial for enabling broader adoption in the medical field.Recently, there has been intense research focused on incorporating domain knowledge into deep learning [30].Domain knowledge introduces domain-specific information, enabling the network to align more seamlessly with the intended field.Incorporating critical domain insights can help to increase interpretability and achieve better generalization, especially when data is scarce [31].Therefore, the integration of domain knowledge primarily demonstrated an impact on the results of image segmentation and denoising, as well as predicting the presence or absence of a disease of interest [32].
For the first time, we propose to consider PBPK models as domain knowledge to overcome the limitations of conventional deep learning approach for guidance of robust dosimetry prediction from pretherapy PET imaging.A PBPK-adapted deep neural network that incorporates constraint determined by PBPK models has been developed.This provides a unique opportunity to tackle the complexity of voxelwise RPT dose prediction with limited data from pretherapy PET imaging.We tested our model on realistic simulations based on physical and physiological principles, which allow systematic investigation of the accuracy and robustness with known ground truth.For the proof of concept, we investigated the accuracy and robustness of our model on critical organs that are mostly affected by PSMA-targeted RPT, namely, liver, spleen, kidneys, and salivary glands.

A. Virtual Patient Based on Realistic Physiological Simulation
A whole-body PBPK model for PSMA ligand distribution was simulated and validated to simulate the time-activity curves (TACs) of different organs [18] as shown in Fig. 1.The necessary population parameters were sampled from Gaussian distributions and are summarized in Table I.Kidney specific parameters were specifically sampled due to the importance of kidneys as one of the main dose-limiting organs in PSMA therapy and being part of the main route of radiopharmaceutical excretion [39].All required parameters that are not specified in Table I were taken from literature [18].Organs of 4-D XCAT anthropomorphic voxelized phantom [33] were associated with a specific activity calculated by the PBPK model to create the necessary pretherapy PET scan images, as well as to calculate the post-therapy dosimetry.The anatomical parameters of our phantoms were adjusted according to the body weight, height and kidney volume given by the population sampling.Phantoms consisted of 3-D cubic voxels of dimensions 4mm 3 .Each phantom was associated with different sized 3-7 ellipsoidal lesions generated using an in-house MATLAB code.The reduction in lesion volume at the end of therapy cycle was computed using the linearquadratic model along with the assumption of exponential tumor growth [24].Lesion positions of pretherapy PET and dosimetry images were randomly chosen among organs into the torso where prostate cancer metastasis has been observed in clinical practice with the organ probability being proportional to the prevalence of a specific metastatic site [34].

B. Simulation of Pretherapy PET Scan
We performed PBPK modeling over a time interval of 2h for the pretherapy PET scan simulation as shown in Fig. 1(a). 68Ga-PSMA-11 was used as a pretherapy with injected activity (115 ± 10 MBq).As described above and for each patient, an XCAT phantom with random lesions was generated where the specific organ activities were set to be the result of the PBPK model at 90 min p.i.All parameters of the PBPK model were customized with Gaussian distribution.To account for the attenuation effect of different tissues, we generated an attenuation map by mapping the linear-attenuation coefficients to corresponding phantom organ densities at 511 keV.The activity and attenuation of phantoms were then used by the in-house built forward projection simulation model integrated with attenuation, scattering and random effects to simulate a PET scan sinogram.The PET images were reconstructed from sinogram with addition of simulated Poisson noise [40] using ordered subset expectation maximization (OSEM) algorithm [41].

C. Post-Therapy Dosimetry Calculation
Similarly for pretherapy, an XCAT phantom with its lesion map was also generated for each patient, as described in Section II-A, where the specific organ activities were set in correspond to the integral of TACs of PSMA ligand as demonstrated in Fig. 1(b).A PBPK simulation was performed for 177 Lu-PSMA-I&T using injected activity (7400 ± 150 MBq) at time interval of 20 days.We calculated the dose using the dose voxel kernel (DVK) method [44], [45].The mean absorbed dose D(r T , T D ) [Gy] to target voxel r T over a defined period T D is defined as where Ã(r S , T D ) [Bq] is the simulated time-integrated activity in source voxel r S over the period T D , and S(r T ← r S ) is the radionuclide-specific quantity representing the mean absorbed dose to the target voxel r T per unit activity present in the source voxel r S .Voxel-S-values matrices simulated for soft tissues and bone tissues were obtained from the University of Bologna at www.medphys.it[49].

D. Network Architecture
We used a conditional generative adversarial network (cGAN) for the post-therapy dose prediction.The 3-D RPT dose GAN was previously developed in our group by Xue et al. [26] and optimized specifically for handling conditional pretherapy dosimetry prediction [50].Comparing to traditional GAN, cGANs incorporating contextual conditioning information can offer enhanced robustness in terms of handling specific conditions, anomalous observations, perturbed samples, and variations in the data distribution.This characteristic renders cGANs a more judicious preference for our study, wherein control over the generated outputs is requisite.
A 3-D U-net architecture [51] was used as a generator G, without batch normalization, and with strided convolution downsampling.The generator was developed to synthesize the target dosimetry image from the input.In the case of the discriminator D, a custom discriminator was created to distinguish real data from the generator output [52].In both the generator and the discriminator, Adam optimizer was used [53].Fig. 2 illustrates a schematic representation for the used network, including the PBPK constraint as a domain knowledge for robust learning.

E. Data Preparation
Overall, 100 virtual patients, including pretherapy PET images and dosimetry images with segmentation masks of Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.target organs, were created for training process, meanwhile 20 patients were kept aside for testing to assess the model ability to accurately predict dosimetry values.PET images were in standardized uptake value (SUV) units, while the dosimetry images were in Gray (Gy).In both cases, we performed min-max normalization.The images were split into 32 × 32 × 32 patches for network training by applying a convolution of stride 16.This step was done to increase the overall network robustness for spatial perturbations and minimize the computational resources needed.

F. PBPK Loss Term
To create a PBPK-based loss term, the PBPK model was simulated with the same parameters used previously for virtual patient simulations in Sections II-B and II-C to generate PBPK-based dose map.The segmented organs (liver, spleen, kidneys and salivary glands) were taken, and their values were set to the corresponding cumulated activity value given by the model.Since the lesions have a significantly higher activity than the surrounding tissue, they were also segmented out from the corresponding organs if present.The organ doses were calculated with the DVK method as previously described in Section II-C.The PBPK equations describing organwise activity were used to constrain the average dose within each organ.Moreover, instead of fixed dose values, specified ranges were determined to allow spatial distribution of heterogeneous doses within organs.Namely, the PBPK-based boundaries for a subset of organ j within a patch k of size of N are defined as where b l,j is the lower-constraint boundary, delta l defines the width of the lower constraint and is a real number in the interval [0, 1], b u,j is the upper constraint boundary, delta u defines the width of the upper constraint and is a real number in the interval [0, 2].Values of both delta l and delta u were determined empirically in our study.G i is the voxel generated by the generator, and C i is the corresponding voxel of the PBPK predicted dose value, and N is the total number of voxels.The conditional operator I(C i , G i ) is added to control generating non-negative dose of voxel G i for C i voxel in organ j, and is defined as Finally, the PBPK loss term can be defined as where G k,j is defined as Due to the fact that dose values are only non-negative, an additional loss term was added

G. Training Procedure
The final Generator loss function L G is given by where L rec is the L1 loss used to compute the voxelwise L1 distance between the generated output and dosimetry Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
image.Comparing to L2 loss after investigation, L1 loss is less sensitive to data outliers and large deviations.Hence, it forces the network to produce a less blurry output, and guides the network to roughly prediction of low-frequency details, providing stable voxelwise estimations [26].L adv is the binary cross-entropy loss (BCE) between the generator and discriminator outputs where the generated output has been labeled as real, L PBPK is the PBPK-based loss term, L nonneg is the loss term of dose non-negativity, and λ [1,4] represent the weights of each term, respectively.In the case of discriminator training, we used BCE loss function.The loss weights' values were determined empirically, relying on experimental evidence.we set λ 1 = 0.025, λ 2 = 250, λ 3 = 200, and λ 4 = 200.For a fair comparison, the nonconstrained and constrained networks were trained on equal epoch numbers with excluding and including L PBPK terms.

A. Synthetic Dataset Validation
An example of coronal pretherapy PET image as well as the calculated post-therapy dosimetry on virtual patients are shown in Fig. 3.
The synthetic pretherapy PET images visually appear quite similar to actual PET data.In addition, the observed SUV values are within the clinically observed values.The doses calculated for individual organs and lesions are within the range reported in clinical trials [54] and possess a broad variability needed to adequately cover the sample space.However, the dosimetry images are much sharper, but they still exhibit heterogeneity similar to the usual dosimetry images observed in practice.The most likely reasons are that dosimetry images are usually based on interpolating activities from SPECT images sampled at several time intervals.SPECT images appear quite noisy and blurred due to the nature of imaging technique and physical phenomena involved [29].

B. Influence of the PBPK Constraint on Generated Dosimetry Image Quality
Fig. 4 illustrates typical examples of coronal dosimetry image synthesized by our 3-D RPT dose GAN with and without the added PBPK loss function.There is no considerable image quality variation that could be observed within the organs in the case of the generated images by the network with the added loss term.The added term that acts as a regularizer to constrain any major deviations from the targeted organs.Hence, it does not dramatically impact any image portion outside of the given organs.

C. Influence of the PBPK Constraint on Dose Prediction Accuracy
Dose-volume histograms (DVHs) were constructed for the test set in order to explore the impact of equation-consistent behavior on the organ-specific generalization.The PBPK constrained model produces a DVH that resembles more the ground truth than the unconstrained model as shown for a sample patient in Fig. 5.To quantitatively assess the variation in the dose accuracy by adding the PBPK constraint, mean absolute percentage error (MAPE) was calculated for each of the organs present in the DVH individually.The results were compared between the unconstrained and constrained models (Fig. 6).Overall, the general MAPE obtained for liver is 19.21±8.97%without and 6.16±2.77%with the PBPK constraint.In the case of spleen, MAPE of 28.33±7.18%without and 15.95±10.7%with the constraint was achieved.MAPE of kidneys is 52.63±21.02%without and 29.93±10.09%with the constrained model.Whilst MAPE of salivary glands is 33.1±11.63%without and 19.72±7.91%with the constrained PBPK model.Moreover, the robustness of the PBPK constraint was assessed using four distinct datasets comprising 80, 60, 40, and 20 patients, each accompanied by corresponding test Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.datasets.The result revealed that our methodology maintained robustness even with limited training data, as evidenced in Fig. 7.

D. Dose Distribution Assessment Within the Organs
Table II gives the calculated values of half-maximal dose (D50) and equivalent uniform dose (EUD) of tested subject.The EUD was used as a dose-response index for evaluating the cumulative dose inhomogeneity and distribution within the three organs of interest.Meanwhile, D50 is a treatment parameter used to estimate the accumulative dose that produces 50% of the maximum response.All EUDs and D50s were calculated directly from the corresponding DVH [55].Both values of D50 and EUD for each organ were closer to the ground truth values in case of the PBPK constrained model, indicating that the introduced PBPK constraint makes the overall dose distribution closer to the ground truth.

IV. DISCUSSION
A critical challenge for the prediction of post-therapy dose distribution based on pretherapy imaging is the huge information gap.The pretherapy imaging measures the uptake at a single time point, while the post-therapy dose distribution is the overall effect of radiation-matter interaction over a long dynamic course.In clinical practice, 177 Lu-PSMA-I&T as a therapeutic tracer is injected at several-fold higher activity than the imaging tracer, 68 Ga-PSMA-11, and their PSMA ligands

TABLE II SUMMARY OF HALF-MAXIMAL DOSE (D50) AND EUD VALUES
WITHOUT AND WITH THE PBPK CONSTRAINT are similar but not identical.In addition, their biodistribution may be influenced by tumor burden [18], injected activity and cold ligand amounts [16].Xue et al. [25] demonstrated that the variation of their biological half-lives between patients is although considerable small, but can not be neglected (11.2 ± 6.0 h).According to previous studies [12], [13], post-therapy dose distribution correlates with SUV values of pretherapy imaging.Such relation demonstrates that taking the pretreatment information into account could assist the estimation of the post-therapy dosimetry [17], [21], [22].
As PBPK model describes the fundamental principles of RPT dose distribution, the prediction based on PBPK model brings the advanced domain knowledge of the spatiotemporal biodistribution of the therapeutic tracers.Therefore, it assists to integrate pretherapy measurements into individualized prediction of post-therapy dosimetry [56].The combination of population-based parameters and individually inferred parameters from pretherapy measurements can reduce the complexity and improves the robustness in the pretherapy dosimetry prediction.However, the hand-crafted determination and inference of individualized parameters may not fully explore the potential of pretherapy measurements in the customization of the PBPK models and limits the accuracy of individualization.
Evidence has demonstrated the spatial heterogeneity of dose distribution influences the response of RPT [42], [43].Although voxelwise dosimetry monitoring has been established for RPT [44], [45], the pretherapy prediction of voxelwise dosimetry is technically more challenging due the increased complexity of the inference.The PBPK-based prediction is in principle only organwise.On the other side, pretherapy imaging, such as PET, has already captured the intraorgan heterogeneity of tracer biodistribution, which can potentially support inference of voxelwise distribution.However, conventional inference scheme, such as the applied Bayesian estimation, is unable to accommodate the complexity of voxelwise parameter customization or dose estimation.
Moreover, prediction of post-therapy dose distribution of conventional PBPK model is only organwise dosimetry, which assume that all micro-compartments related to a specific compartment contribute to therapeutic tracer concentrations observed on pretherapy images.Therefore, an organ-based approach is unable to reveal the heterogeneity of dose distribution, and hence is not sufficient for accurate post-therapy dose prediction or realization of an optimal treatment planning.Xue et al. [26] demonstrated the potential of deep learning to estimate voxelwise post-therapy dosimetry quantitatively based on pretherapy PET imaging.However, complex mechanism with multiple confounding factors may hamper the performance of deep learning models in RPT dose prediction from cross sectional measurements, especially with limited data set size.Therefore, it is concerned about the interoperability and generalizability of the pure data-driven deep learning method [46].This may hamper the clinical translation into the treatment planning for RPT.
Defining a domain-based acceptable range of values could help improve the model quality, and prevent over-fitting to the underlying noise [31].Therefore, the combination of both PBPK constraint and deep learning can complement the advantages of both approaches and overcome their limitations.The results of our experiments confirmed our hypothesis that PBPK-integrated deep learning can enhance dose prediction, as opposed to an unconstrained model, indicating its feasibility, as demonstrated in Fig. 6.A novel loss term was devised relying on PBPK-based calculations to constrain the network output into acceptable ranges.The advantage of full PBPK model is that it takes well-established patient physiological characteristics into account, as well as the general drug distribution inference based on pretherapy measurements [16].As a result, this has contributed to improve the robustness and generalizability of the developed model to adapt new data that it has not previously encountered.Therefore, the term combines PBPK predictions with voxelwise prediction to boost dose prediction accuracy.The advantage of the method presented here is that it is modular and can be easily extended to other organs or tumor lesions, given that proper segmentation is done.The newly developed model significantly increased the network's accuracy in the PBPK constrained organs under this study, as shown in Fig. 7.
Castration-resistant prostate cancer (CRPC) can metastasize in all organs examined in this study [34].Lesions located in organs of interest are a known cause of prediction inaccuracy.Hence, a dose prediction method should attempt to minimize the perturbation caused by the lesion presence.The performance of the PBPK guided model seems to be agnostic to the presence or absence of lesions within the organs of Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
interest.Probably even better accuracy can be obtained in organs with metastasis if the PBPK constraint is included for the lesion as well.This will be explored in a future study.Moreover, a prerequisite for the application of the proposed PBPK-integrated deep learning model is the availability of automatic organ and tumor augmentations, which is not accurately available at the moment.However, several attempts has been made for automatic segmentation on PSMA PET imaging [47], [48].The proved potential confirmed that such kind of segmentation will be feasible in near future.
This study has several limitations.The first is that, a specific PBPK model was used for XCAT phantom creation.Apparently, a model which supports our aim.Therefore, the most important anatomical and physiological parameters are included in the model.However, the current limitations of the PBPK models and the individualization of the parameters may restrict the performance of our model.For example, the current PBPK model contains some fixed parameters, which are assumed to be true for all patients.However, according to Hardiansyah et al. [23], a marginal effect of changed fixed parameters could be proved to the TIACs, although the impact might be still minimal on the estimation of organ TIACs in PRRT using the PBPK model.For the second limitation, although the PBPK soft constraint has demonstrated its merit in both noisy and noise-free data-sets, the synthetic dataset was created based on theoretical PBPK settings where almost all kinetic parameters of the model are set to predefined ranges which are not supposed to be strictly Gaussian in nature, and consequently also the TIACs.The intraindividual variability from the individual dosimetry uncertainty has further been reported to be about 10-20% in several studies [23], [57], which may further challenge the potential performance of the proposed model in a real situation.Another limitation of our study is that, the S-values used in data simulation are not patient specific and applied by organ-level dose calculation tools were obtained based on reference phantoms, which are not intended for individualized dose calculation.However, our approach demonstrated good potential in dose prediction.In future work, the potential of this approach in dose distribution prediction should be to attempt to utilize this constraint on the clinics-acquired patient data and in different clinical contexts.

V. CONCLUSION
This study proposed a novel PBPK-integrated deep learning to combine the PBPK-based domain knowledge in the pretherapy prediction of voxelwise dose distribution using deep learning.A novel loss function term was devised relying on PBPK-based estimations to constrain the network output within acceptable ranges.The term combines PBPK predictions with identified organs to boost the dose prediction accuracy of most affected vital organs.Overall, the results presented here are encouraging for improving deep learningbased dosimetry prediction.They demonstrate that integrating pharmacological and physiological knowledge might help to overcome the limitations of AI as a tool for personalized RPT treatment planning.

Fig. 3 .
Fig. 3. Visualization of coronal image for pretherapy PET image (a) and post-therapy dosimetry calculation (b) of a simulated patient.

Fig. 4 .
Fig. 4. Visualization of coronal image for cGAN predicted dosimetry without (a) and with the PBPK constraint (b).

Fig. 5 .
Fig. 5. Sample of a DVH for dosimetry predicted by unconstrained model (a) and PBPK constrained model (b).

Fig. 6 .
Fig. 6.Comparison of the prediction accuracy between the unconstrained and constrained PBPK model.

Fig. 7 .
Fig. 7. Summary plot depicting dosimetry prediction performance of the unconstrained and constrained PBPK model with different amount of data.