Experimental validation of an innovative approach in biokinetics study for personalised dosimetry of molecular radiation therapy treatments

One of today’s main challenges in molecular radiation therapy is to assess an individual dosimetry that allows treatment to be tailored to the specific patient, in accordance with the current paradigm of ‘personalized medicine’. The evaluation of the absorbed doses for tumor and organs at risk in molecular radiotherapy is typically based on MIRD schema acquiring few experimental points for the assessement of biokinetic parameters. WIDMApp, the wearable individual dose monitoring apparatus, is an innovative approach for internal dosimetry based on a wearable radiation detecting system for individual biokinetics sampling, a Monte Carlo simulation for particle interaction, and an unfolding algorithm for data analysis and integrated activity determination at organ level. A prototype of a WIDMApp detector element was used to record the photon emissions in a body phantom containing 3 spheres with liquid sources (18F, 64Cu and 99m Tc) to simulate organs having different washout. Modelling the phantom geometry on the basis of a CT scan imaging, the Monte Carlo simulation computed the contribution of each emitting sphere to the signal detected in 3 positions on the phantoms surface. Combining the simulated results with the data acquired for 120 h, the unfolding algorithm deconvolved the detected signal and assessed the decay half-life (T 1/2) and initial activity values (A(0)) that best reproduces the observed exponential decays. A 3%–18% level of agreement is found between the actual A(0) and T 1/2 values and those obtained by means of the minimization procedure based on the Monte Carlo simulation. That resulted in an estimation of the cumulated activity <15%. Moreover, WIDMApp data redundancy has been used to mitigate some experimental occurrences that happened during data taking. A first experimental test of the WIDMApp approach to internal radiation dosimetry is presented. Studies with patients are foreseen to validate the technique in a real environment.


Introduction
In molecular radiotherapy (MRT), radioactive tracers with specificity for the tumor to treat are administered to selectively damage diseased cells via the radiation emitted by the radionuclide.However, given their systemic administration, these molecules not only accumulate on the tumor, but also diffuse to healthy tissues, which thus receive an unavoidable radiation dose.
For this reason, the accurate determination of absorbed dose by healthy organs in MRT plays a crucial role in maximizing the effectiveness and lowering toxicity (Stabin 2008).Furthermore, current European regulation (EU Directive 59/2013) considers mandatory the treatment planning and verification in all patients undergoing radiotherapy procedures, including MRT.
In MRT the absorbed dose to lesion and organs at risk (OAR) is traditionally calculated using the MIRD schema (Bolch et al 2009): the absorbed dose is calculated as the product between a time-dependent variable, the cumulated activity Ã, and the S-value.Ã is the area under the biokinetic curve and its assessment remains today one of the most challenging aspect in MRT.The S-value is the mean absorbed dose per cumulated activity.The methods to determine the biokinetic parameters for adminstration of 177 Lu and 131 I, in which a sequence of SPECT scans (usually 4 or 5) is acquired after the administration of the radiopharmaceutical, are summarized in pamphlet MIRD 24 and 26 (Dewaraja et al 2013, Ljungberg et al 2016).Less resource-demanding approaches, as for example the use of an external probe or integration dosimeters (e.g.thermoluminescent dosimeters, TLDs), are much less powerful (Strigari et al 2023).This approach results in an uncertainty on internal dosimetry calculations larger than 30% (Flux et al 2006).
WIDMApp (the Wearable Individual Dose Monitoring Apparatus) is the recent proposal of a new approach in internal radiation dosimetry (IRD) for MRT with both 177 Lu and 131 I due to a presence of gamma emissions.This approach exploits a system with minimal impact on clinical practice and on personnel work and allows for a personalized assessment of biokinetics, taking advantage of a in vivo sampling with a rate higher than the one of the sequential nuclear imaging sessions.
A feasibility study of WIDMApp based on a Monte Carlo simulation of a simplified thyroid treatment with 131 I is described in Morganti et al (2021).
The WIDMApp system is based on 3 components: 1.An array of light-weight particle detectors that can be comfortably worn even continuously over 24 h per day by the patient for several days after the radio tracer administration.These detectors, having at least one of them for each organ of interest, continuously monitor the in vivo photon emission rate of the source volumes and reconstruct the time-counts curves (TCCs, one for each detector).
2. A Monte Carlo simulation that models the system composed by patient plus detectors and the activity distribution on the basis of a single SPECT/CT scan acquired some hours after the injection of the radiophamaceuticals.Then, it computes the probability that an isotope decaying in the given organ eventually yields a signal in one of the detectors, providing the 'probability matrix'.
3. An unfolding algorithm that uses this probability matrix and the experimental data (TCC) to obtain the time activity curve (TAC) of each organ of interests and computes the cumulated activity, from which the absorbed mean dose is estimated.
In this paper, the first experimental validation of the core functioning principle of the WIDMApp approach is presented, using a NEMA IEC Body Phantom Set TM filled with 3 different radionuclide to simulate biological excretion of three different organs.Aim of the validation is to address some of the practical aspects that will characterize the applicative case in real patients, such as detector efficiency and positioning and the use of real SPECT/CT images as input to the Monte Carlo simulation.These are the main issues to be investigated before the future perspective of in vivo tests for the WIDMApp system.

The experimental setup
A sketch of the experimental setup is in figure 1.

The phantom
A NEMA IEC PET Body Phantom filled with water was used for mimicking a human waist.Three internal spheres of diameters 17, 22 and 28 mm were filled with saline solutions of 18 F, 99m Tc and 64 Cu respectively.The three different radionuclides with different half-lives were used to simulate emitting volumes with different washout times inside the body.Three acquisition points on the phantom surface were identified, each facing one of the radioactive spheres to maximize the exposure to that emitting volume.

The radiation detecting element
The detector used in this test is a first prototype of the simple and compact sensor element of WIDMApp.It is composed by a 15 × 15 × 3 mm 3 p-terhpenyl scintillating crystal (Angelone et al 2014), coupled with a 2 × 2 matrix of 6 × 6 mm 2 Silicon Photomultipliers (C-Series SiPM by ONSEMI).The assembly is then encapsulated in a light tight 3D printed ABS (Acrylonitrile butadiene styrene) cylinder, 25 mm diameter, 10 mm thick.
In this detector prototype, specifically designed not for clinical use but only for laboratory tests, the SiPMs are controlled and readout via cable connection using the ArduSiPM board, which is a portable battery-powered low cost board, designed to acquire and elaborate the SiPM waveforms and used to calculate the counting rate (in-depth description of electronics and signal processing can be found in Bocci et al (2014)).The ArduSIPM board communicates with a PC with a cable connection for data recording and displaying.
This prototype respects the basic requirements in terms of performance (i.e. light yield and photon sensitivity) and of lightness and portability, even though it is not yet optimized for the WIDMApp scope.
When necessary, appropriate solutions were introduced in this test to overcome the availability of only one detector element and some limitations in the read-out electronics as the lack of a control system for temperature and humidity which influencing the gain of the SIPMs and represent an important issue.Since the sensor prototype does not provide the information necessary to correct the response of the SIPMs, the experiment was carried out in a room with controlled temperature and humidity to avoid correlations between the detector response and external conditions.
The main limitation was the nonlinearity of the detector response on the entire range of counting rate.The prototype was tuned on low activity laboratory sources (within 1 kBq).In the experimental test, as in a real treatment, a wide range of counts per second (CPS) was expected, depending on the initial activity of the source faced by the sensor and the acquisition time.A correction factor for counting rate >1 kHz was applied (see section 2.4.1).

The data acquisition
The acquisitions of 3 detecting elements in the three positions facing the spheres (figure 1) was realized by cyclically alternating the position of the unique available prototype.Sensor placement was driven by a robotic arm with 6 degrees of freedom (Pololu electronics (R)) to make it reproducible.On each position the arm stopped for 180 s for data acquisition before moving to the next location.The counting rate registered in each position was the data point of the TCC of that position in that moment.The use of a single moving detector is due to the availability of only one single-channel readout system.This setup increases the complexity and adds errors to the measurement, however this allows us to simulate the inevitable and unpredictable movement of the detectors due to patient movements during an in vivo measurement.
Two radioactive set-ups were realized: 1-Isotope: only one sphere at a time was filled with the radioactive solution.This configuration was realized to evaluate the contribution of the filled sphere to each TCC over an acquisition time of ∼50 h and was repeated by changing the filled sphere (one dataset for each radionuclide).Data collected with this setup were used to measure the detector efficiency in this specific validation and will not be necessary for the WIDMApp application in the clinical practice as the detector efficiency and calibration will be evaluated with a laboratory measurements on the isotopes used in MRT.3-Isotopes: all the three phantom spheres were filled with the radioactive solutions and data were stored for ∼120 h to determine the complete TCC.
The initial activities (A(0)) of radioactive solutions used in the two experimental configurations together with physical half-times are in table 1.

The time-counts curves
As already described in Morganti et al (2021), the emissions of each emitting volume contribute to the detected signal.Therefore, in configuration 3-Isotopes the counting rate registered at the instant t by the detector in position j (with j ranging from 1 to 3 as in figure 1) was where CPS i,j (0) is the counting rate at the test start time t = 0, λ i is the decay constant of the ith radionuclide and i = 1, 2, 3 indicates 18 F , 99m Tc, 64 Cu respectively.The counting rate contribution CPS ij (t) is the product of the radionuclide activity A i (t) times the probability that a photon yielded in the sphere S i interacts in the detector in position j, i.e. the 'probability matrix' element ij MC  (see section 2.2) and the detecting efficiency of the jth sensor for the emissions of the ith radionuclide, ij Det  .Therefore in this study the ij Det  values were evaluated with the 1-Isotope experiments, using the espression: rate variation as a function of time (CPS ij (t)) for the 1-isotope acquisitions is shown in figure 2. Single exponential fits, fixing the experimentally measured dark current and the T 1/2 to be the physical one, allow to estimate the CPS ij (0) values and, hence, to calculate the detector efficiency.

Monte Carlo simulation
The probability matrix used to unfold the cumulative TCC data into single organ TAC was calculated starting from a MC simulation of the experimental setup.The radiation tracking in matter was performed using Geant4 (Agostinelli et al 2003, Allison et al 2016), a toolkit widely used in particle physics, and constantly benchmarked for medical applications (Arce 2021).In the simulation, the structure of the NEMA phantom was included via its own CT scan.This is indeed relevant not only to allow a very precise reconstruction of its geometry, but also to test the use of real imaging data of the patient, that is a core step in the complete WIDMApp approach.CT DICOM files were processed in order to reproduce the actual geometry voxel by voxel, whose material and density were evaluated according to a calibration table provided for the CT scanner.Active volumes of the detectors positioned on the surface of the CT imported geometry were also included in the simulation.The three detector positions were obtained by inserting makers in the NEMA phantom during the CT scan in correspondence with the three positions assumed by the detector prototype during data taking.The three volumes, as reported in figure 1, were then manually drawn within the CT images using markers as reference positions.
As far as primary particles generation is concerned, regions of interest (ROIs) outlines were drawn on the CT scan, coinciding with the emitting sphere volumes.The simulation program then uniformly samples the proper primary particles inside each of these ROIs separately, and propagates the generated particles into the phantom.The particle propagation has been simulated using the most precise electromagnetic physics list available among those developed by the Geant4 collaboration, i.e. mStandardPhysics_option4.The secondary production cuts for electrons, positrons and photons were settled to 1 mm and 250 eV.No variance reduction techniques were applied at this level.For each simulation run, the number of primary particles that eventually interact in the given detector was stored.Based on the results of the MC simulations, a 3 × 3 matrix in which for each of the 3 source spheres (rows) the probability ( i j MC ,  ) to have a signal in one of the 3 detectors (columns) was obtained.

Unfolding algorithm
The unfolding algorithm that allows to estimate TACs from TCCs were first described in detail and validated with the feasibility study in Morganti et al (2021).
In brief, a minimization procedure was performed to identify the set of 6 free TAC parameters (3 pairs of A (0) and T 1/2 ) that best reproduces all the TCCs at the same time.In this study, this is achieved by define merit functions ( ( ) 1 2 (one for each radionuclide source) that consider all points from the three TCCs: Counting rate as a function of time for the 1-isotope mode acquisitions: 18 F sphere (top), 99m Tc sphere (center), 64 Cu sphere (bottom).In each plot, the colors represent counts registered in the three detector positions.The CPS(0) value for each dataset is also shown.
where S j (t k ) is the experimental TCC value in the detector position j at the time t k and ( ( ) 1 2 is the test value the signal will be compared with.Finally, σ j (t k ) is the standard deviation estimation for each measure.The aforementioned previous study suggested that good minimization of the intercorrelated merit functions are provided by the Nelder Mead algorithm (Nelder and Mead 1965), that was thus used also for this study.

Data conditioning
The continuous data acquisition enabled by the wearable detector foreseen for WIDMApp provides a great abundance and redundancy of data allowing for example the detection of changes in the functional shape of an organ clearance throughout the observation period.It is also a powerful tool for verifying data consistency and integrity.On these occasions, redundancy allows eliminating inconsistent data without spoiling the precise measurement of the TCC.Systematic errors can be individuated and corrected as well, as shown in this first experimental test.
During the 3-Isotopes mode acquisition, an acquisition interruption and a change in operating condition compatible with a systematic, position dependent, mis-position of the sensor occured due to an accidental power outage in the building and the consequent restart of the robotic arm.The systematic shift of the data before and after the power cut at 40 h is clearly visible in the 3 exponential fits of the position 3 TCC (figure 3).
This effect, as well as other possible mis-positiong errors, can be corrected thanks to the large amount of data collected by applying a scaling factor to experimental points after the gap.This factor was calculated as the one that minimizes the Pearson's χ 2 of the fit on the single 3-exponential function (e.g. one exponential for each radionuclide that contribute to the detector response), and was found to be 1 (i.e.negligible correction), 0.95 and 0.90 for detector position 1, 2 and 3 respectively.

Saturation correction
The detector readout system implemented in the experimental setup has a dead time (i.e.time it takes for the system to register another event) of 1.2 μs.Consequently, the probability of missing some particles crossing starts already at about 1 kCPS, with a nonlinear increase as the particle flow increases.To account for this effect, particularly relevant in the first hours of measurements, a data calibration tailored on a Monte Carlo simulation was implemented.In this simulation, events are generated with a time difference extracted from a poissonian distribution centered on the expected signal rate.Considering a windows of 1.2 μs it is possible to obtain the number of events lost due to the dead time for different signal rates.The number of expected counts, corresponding to all the particles which interact in the detector (including the ones lost due to the dead time) was calculated starting from the counts recorded by each detector and an calibration function obtained from an exponential fit of the data obtained from the aforementioned Monte Carlo simulation with the exponential function y = A(1 − e − x/ B ). Systematic errors on expected counts are calculated from the errors propagation of the calibration function parameters A and B obtained from the fit.

Results
The detector efficiency matrix obtained by means of the 1-Isotope mode measurements is reported in table 2 and allows exploiting equation (2) to correlate detector counting and sources activity.
This correlation was performed by the unfolding algorithm using the probability matrix provided from the Monte Carlo simulation and the TCCs experimentally collected.The probability matrix reconstructed from the simulation of the NEMA structure is shown in figure 4.
The TCCs recorded in the three detector positions in the mode 3-Isotopes is shown in figure 5.
The results are shown in table 3, in terms of sources initial activities and decay half life.The largest source of uncertainty was the detector efficiency, Det  , defined in equation (2) and showed in table 2. Therefore, to estimate the results uncertainties, the unfolding algorithm was runned 10 000 times, sampling the efficiencies from a normal distribution as obtained from the 1-Isotope mode calculation.
Results of the WIDMApp approach in table 3 show an excellent agreement (considering the state of the art standards Gear et al 2018) between the measured and real values both for the decay times and for the sources initial activities.In particular, goal of WIDMApp is to evaluate the effective half lives, which are estimated with an error of ∼3%−8% in all 3 source cases. that a given primary particle originating in the given ith active sphere (row) eventually gives an energy deposition in the jth detector (columns).
Figure 5. Counting rate as a function of time in the 3 detector positions after data conditioning in 3-Isotopes setup: detector position 1 (left), detector position 2 (center), detector position 3 (right).Data (blue dots) are compared with the shape (red line) of the superimposition of three exponential decays with the physical half-lives of the three isotopes and the dark counts set to the measured value.
Table 2. Detector efficiencies calculated from the 1-isotope mode.The energy of the detected photons is also reported.Uncertainties on the Det  values were calculated by propagating the statistical error on counting rate and the one on initial activity.The uncertainty on Monte Carlo results were neglected due to the good statistics, being smaller than 1%.As expected, the detecting efficiency depends on the energy of the detected photons.

Discussion
After a first, purely Monte Carlo based feasibility study (Morganti et al 2021), the aim of the present study is to show an experimental test of the WIDMApp approach, using an antrophomorphic phantom to mimic a patient, and three different decaying radionuclides ( 18 F, 64 Cu and 99m Tc) to simulate organs having different washout behavior.
The strength of WIDMApp is the extensive and frequent sampling of the individual biokinetics, which allows to obtain the TAC evolution with better time resolution and less resources required over the entire treatment period: from the uptake to the early washout and the long tail.In fact, a better knowledge of the TAC behaviour in time for each emitting organ would largely improve the accuracy of the cumulated activity estimations, and, therefore, of the absorbed doses, particularly in the case of long term retain.For that reason the accuracy of the T 1/2 estimation for the TAC determination is of crucial importance for WIDMApp.
The results in terms of TAC determination are extremely encouraging.Already in the current prototype phase, WIDMApp was able to provide an estimate of the cumulated activity ˜( ) is the decay constant) for the three sources within 15%, considering the propagation of errors of the parameters obtained from the minimization algorithm.This is a remarkable result, which suggests that the accuracy in measuring the mean absorbed dose with WIDMApp could actually be competitive compared to the current Standard Of Care for MRT treatment with 177 Lu and 131 I.
Indeed, it has to be noted that in this study both A(0) and λ were estimated by fitting data acquired by WIDMApp.In the clinical case, instead, the activity distribution in the body will be quantified by a SPECT/CT scan (the same one used as input to the simulation for modelilng the patient's structures) acquired in the first 24-48 h after injection.
Furthermore, the artifice to experimentally simulate organs with different washouts using isotopes with different decay ways (2 β + and 1 γ) and, therefore, different energy of the photons to be detected is a more complex scenario compared to the clinical application where it will not be necessary to carry out acquisitions with the 1-Isotope mode.As far as sources initial activities are concerned, the agreement is very good in the case of 18 F and 64 Cu, while the measured 99m Tc half-life shows a 18% deviation from the expected value.This difference is expected to be ascribed to the lower efficiencies to photons emitted by this isotope shown by the detector.In fact, as shown in table 2, detector efficiency on 99m Tc is about one half the one on 18 F and 64 Cu, due to the differences in these photons energies (140 keV versus 511 keV).In the actual treatment the photons emitted by the administered radionuclide will be uniform in energy at any point in the body.The response of the detector will be fully characterized with dedicated laboratory measurements with respect to the emissions of that radionuclide and the detection efficiency quantified before the treatment.
Moreover, the test was performed with a very early prototype of the hardware system, making the measurement less straightforward than it could have been.First and foremost, the lack of a real multi-channel detector: the need to use a single sensor operated by a robotic arm obviously implies a certain loss in reproducibility of detector positioning, that is expected to reflect in all the analyzed results.In principle, the technology or the shape and dimensions of the crystal could be appropriately designed for each anatomical district and tuned on the radionuclide and on the treatment specificity.However, aim of this study was to investigate the validity of the concepts underlying the WIDMApp approach to the IRD.Hence, the detector prototype used was not specifically designed for this setup and a data correction that will not be needed when the final detector will be available was required.However, the data redundancy that is at the base of the WIDMApp approach allowed the successful controlling and correction of these experimental occurrences, thus confirming its role in mitigating possible similar eventualities in the real application case.
The detector limitations and the approximations introduced to simulate the presence of three organs with different clearance do not alter significantly the relevance of this test.Its real, main limitation is the lack of a representation of the background that is due to diffusion of the radiopharmaceutical throughout the body, and particularly in the blood compartment.Besides being difficult to practically accurately represent such a background in the used setup, its characterization is out of the scope of this first experimental validation of the WIDMApp approach, whose aim was to provide the first verification of the efficacy of the proposed technique via laboratory measurements.Tests with patients will be needed to properly assess the impact of this kind of background using dedicated sensors that will be placed in anatomical region of the patient with high vascularization and far from the OARs (e.g.thigh and/or calf).The background due to the hematic compartment could be also extrapolated with the measurement of the activity in the blood extracted from the patient at different intervals after the radiopharmaceutical injection, that is already carried out today in many medical centers.Another limitation of this study is the use of radioactive sources whose TACs depend only on the physical half-life of the three considered radionuclides and parametrizable by a single monoexponential decrease.In a real clinical scenario it will be necessary to take into account the biological uptake and clearance for each organ, which complicate the functional form of the TAC The number and type of functions that will be used for the description of the TACs in patients must be assumed based on the data present in the literature and clinical practice.The use of more complicated functions potentially increases the uncertainties of the estimate of Ã.
Furthermore, in the present work the possible limitations due to the wearability of the final device were not taken into consideration.The design of the new miniaturized sensors and the ways in which patient data will be taken will be studied, taking into consideration patient comfort and the implementation in the clinical practice.
The next challenge to assess the robusteness of the WIDMApp approach will be to test it on patients undergoing MRT.In this case, the in vivo measurement will permit to evaluate the realistic impact of the background from the hematic compartment and, in general, from the body.The actual biokinetics expressed by the patient will be a real benchmark to validate the dosimetric approach and to confirm the results obtained with laboratory tests.

Conclusions
The accurate determination of dose absorbed by OARs in MRT treatments and, mostly, the possibility to perform such measurements on all patients are a challenging task.European legislation (EU Directive 59/2013) considers this subject mandatory but, in fact, only a minority of the medical centers have the resources (both human and instrumental) to provide this service basing their data on nuclear medicine imaging.
Extending the serial SPECT/CT approach to all patients, even limiting it to a few time points, requires huge resources in term of manpower and hospital availability to adapt the instrumentation and management of medical staff and patients for this purpose.
WIDMApp aims to minimize the impact on the public health system by limiting to only one SPECT/CT acquisition, leaving to an economical system (managed by a properly instructed patient) the collection of data for the precise and individual valuation of the dose accumulated by tumors and OARs.Furthermore, data can be collected in any time frame that is considered useful for an accurate characterization of each individual's own pharmacokinetics.
The present study is the first experimental validation of WIDMApp approach for an internal dosimetry assessment of MRT patients based on personalised biokinetic measurements.Using a body phantom with three emitting volumes mimicking organs with different biokinetic behaviors, the estimated time decay constants of the emitting volumes resulted in agreement with the true values within 3%-8%.These results confirmed the WIDMApp's high potentiality of assessing the cumulated activity at organ level with high margin of precision.

Figure 1 .
Figure 1.A scheme of the experimental setup.The three spheres represent the emitting points filled with the radionuclide.The three orange squares shows the position where the detector has been placed to acquire the datasets.

Figure 2 .
Figure2.Counting rate as a function of time for the 1-isotope mode acquisitions: 18 F sphere (top), 99m Tc sphere (center),64 Cu sphere (bottom).In each plot, the colors represent counts registered in the three detector positions.The CPS(0) value for each dataset is also shown.

Figure 3 .
Figure 3. Raw counting rates as a function of time for detector Position 3 in 3-Isotopes mode.The two lines represent two separate 3-exponential fits, whose different slope suggest the appearance of some offset after the power outage at ∼40 h.

Figure 4 .
Figure 4. Monte Carlo probability matrix.Each elements represent the the probability ijMC

Table 3 .
Final Results: initial source activity and decay half life, as reconstructed by the WIDMApp algorithm and the relative uncertainty as obtained from the minimization procedure.Their respective deviation from the real values in table 1 is also shown.