Multicomponent mathematical model for tumor volume calculation with setup error using single-isocenter stereotactic radiotherapy for multiple brain metastases

We evaluated the tumor residual volumes considering six degrees-of-freedom (6DoF) patient setup errors in stereotactic radiotherapy (SRT) with multicomponent mathematical model using single-isocenter irradiation for brain metastases. Simulated spherical gross tumor volumes (GTVs) with 1.0 (GTV 1), 2.0 (GTV 2), and 3.0 (GTV 3)-cm diameters were used. The distance between the GTV center and isocenter (d) was set at 0–10 cm. The GTV was simultaneously translated within 0–1.0 mm (T) and rotated within 0°–1.0° (R) in the three axis directions using affine transformation. We optimized the tumor growth model parameters using measurements of non-small cell lung cancer cell lines’ (A549 and NCI-H460) growth. We calculated the GTV residual volume at the irradiation’s end using the physical dose to the GTV when the GTV size, d, and 6DoF setup error varied. The d-values that satisfy tolerance values (10%, 35%, and 50%) of the GTV residual volume rate based on the pre-irradiation GTV volume were determined. The larger the tolerance value set for both cell lines, the longer the distance to satisfy the tolerance value. In GTV residual volume evaluations based on the multicomponent mathematical model on SRT with single-isocenter irradiation, the smaller the GTV size and the larger the distance and 6DoF setup error, the shorter the distance that satisfies the tolerance value might need to be.


Introduction
Brain metastases of non-small-cell lung cancer (NSCLC) are the most common intracranial brain metastasis, occurring in roughly 40-50% of cancer patients [1]. Stereotactic radiosurgery (SRS) and stereotactic radiotherapy (SRT) have become increasingly important in the management of brain metastases, because the use of these treatments has improved the systemic disease control and reduced the effects of radiation on normal tissue [2][3][4]. A singleisocenter irradiation technique for SRS and SRT has been used for multiple brain metastases [5][6][7]. The singleisocenter technique irradiates multiple targets simultaneously in a single isocenter, which has the advantage of 1 3 shortening the dose delivery time compared to conventional irradiation in which one isocenter is set for each target [8,9]. It has also been reported that the dose distribution to the target and normal tissues provided by the single-isocenter technique is the same or better compared to conventional irradiation [6,10]. When multiple targets such as brain metastases are to be irradiated, the singleisocenter technique is thus effective.
In many single-isocenter irradiation cases, the isocenter is not located at the center of the target, and single-isocenter irradiation has a disadvantage in that the effect of patient setup errors, including translational and rotational errors, may be greater than in conventional irradiation methods where the isocenter is set at the center of the target. Several research groups have examined the uncertainty in the target dose due to the setup error of single-isocenter irradiation [11][12][13][14]. We reported the optimal distances from the isocenter where targets from 1.0 to 3.0 cm satisfy several tolerances for a set of patient setup errors [11,12]. Rojas et al. derived the optimal planning target volume (PTV) margin required to meet the tolerances [13]. Although there have been evaluations of the target's tumor control probability (TCP) and normal tissue complication probability (NTCP) in the presence of some setup errors, most of those studies performed only physical evaluations [14,15].
Mathematical models have long been used to elucidate physiological and pathological processes in humans, and the models have been extended to many areas of medicine [16,17]. Mathematical models have been applied for calculations of complex responses in physiology and pathology, and many studies have been reported the use of mathematical models for medical care, the effects of vaccines, and the prediction of anticancer drug effects [18][19][20]. The ordinary differential equation (ODE) mathematical model is widely used in applications such as the computation of tumor growth models because of its simplicity and stability [21,22]. Evaluations using mathematical models based on ODEs for the radiotherapy of tumors have also been reported [23][24][25][26]. By using a mathematical model based on an ODE, it is possible to evaluate a tumor's volume as an output in order to determine the effects of radiotherapy on the tumor and to derive the optimal irradiation schedule and dose. Therefore, radiobiological calculations such as those for a tumor response can be made from physical dose indices of radiotherapy, and the effects of the radiotherapy on the tumor volume can be evaluated. However, we have found no reports of an evaluation of the radiobiological effects of a physical dose reduction to the tumor due to setup errors that use a mathematical model based on an ODE.
In this study, we evaluated the treatment effect by using the calculated dose to the tumor based on a mathematical ODE model when a patient setup error occurred in singleisocenter irradiation for multiple brain metastases. We also derived the distance from the isocenter that satisfied the tumor tolerance and the set residual rate.

Phantom design and dose distribution generation
The diameters of the target spheres, which were simulated gross tumor volumes (GTVs), were set as follows: 1.0 cm (GTV 1), 2.0 cm (GTV 2), and 3.0 cm (GTV 3). MATLAB ver. 2021b software (MathWorks, Natick, MA, USA) was used. The coordinates (unit: cm) of the GTVs were set so that the distance (d) from the center of the GTV to the set at the origin of the coordinate axes was 0-10 cm [14,27]. The dose distribution was adjusted so that D95, which is the prescribed dose administered to 95% of the GTV volume of each GTV, is covered by the 27 Gy in three fractions prescribed dose.

Calculation of the GTV dose considering six degrees-of-freedom (6DoF) setup errors
We used affine transformation to calculate the translation and rotation of the GTV in the axis through the isocenter. The center coordinates of the GTV were set as GTV x , GTV y , and GTV z . To account for the six degrees-of-freedom (6DoF) setup error for the GTV, the three-dimensional (3D) coordinate points of the GTV were transformed to the isocenter. The 6DoF within approximately 1.0° and 1.0 mm accuracy were used because it has been reported that correction can be made in radiotherapy to within 1.0° and 1.0 mm in brain SRT [28][29][30].

Mathematical model for tumor growth with radiotherapy
The multicomponent mathematical model (MCM), which models growth by distinguishing between tumor cell states and growth rates (such as active and quiescent tumors) was used. Combining the GTV growth shown by the MCM and the cell lethality represented by the radiation calculation by using the linear quadratic (LQ) model, the volume calculation of a GTV in radiation therapy can be expressed. In order to simplify the evaluation using the model in this study, we optimized the MCM with M = 2 used in a previous study [25]. Active tumors were defined as the types T 1 , T 2 , T m , where V 1 , V 2 , V m are their volumes. V Q is the volume of quiescent cells (T Q ), and V ND is the volume of nondividing cells (T ND ). K 1 and K 2 are constants that are related to the growth rate and carrying capacity of the tumor. The dose rate D i was set as a constant D of D95 (Gy), and the sum of V 1 , V 2 , and V Q was defined as the GTV volume.

3
V 1i , V Qi , and D i are the i-th GTV volumes and the absorbed dose with V 1 , V Q , and D, respectively; N is the total number of fractions in radiotherapy. Additionally, D was the value of dose per fraction to the GTV when a 6DoF setup error occurred (Fig. 1). The parameter t intra is the irradiation time, and t inter is the time between irradiations in fractional dose irradiation [31]. The tumor volumes in fractionated irradiation for V 2 , V Q , and V ND were calculated in the same way using Eq. (1).
We used the measurement data of tumor growth of two non-small cell lung cancer (NSCLC) cell lines, A549 and NCI-H460 (H460), to evaluate the effect of SRT on brain metastasis [32][33][34]. MATLAB using the SimBiology toolbox was used to optimize the parameter values of the MCM to match the measured values of A549 and H460 cells, respectively. The α/β values of the parameters for A549 and H460 cells representing the sensitivity to radiation, respectively, were obtained from other studies [25,34,35]. Table 1 shows the MCM calculation parameters obtained using A549 and NCI-H460 cells.
(1)  Fig. 1 The multicomponent mathematical model (MCM) simulates tumor growth. In this study, the MCM was used to calculate the GTV volume by dividing tumor cells into active tumor cells with M = 2, quiescent tumor cells, and nondividing tumor cells. The GTV volume was calculated by inputting the D95 (Gy) in dose per fraction (RT) using a DVH for each GTV into the MCM when 6DoF setup errors occurred. In addition, to simulate fractionated irradiation in radiotherapy, the following equation was defined:

Distance from the isocenter that satisfied the GTV tolerance value
The pre-irradiation GTV volume is used as the reference, and the volume at the end of the irradiation is used to calculate the GTV residual volume rate [25]. We defined the tolerance values of the GTV residual volume rate as 10%, 35%, and 50% using the GTV residual volume at the end of irradiation, referring to the values of other studies [36][37][38], and we derived the distance from the isocenter at which the GTV satisfied the allowable value when each 6DoF error occurred. The status "no derived distance" was defined as a case in which the distance derived to satisfy the tolerance value exceeded 10 cm, since the maximum distance between the GTV and the isocenter evaluated in this study was 10 cm.  Fig. 3 The MCM-based GTV 1 (1.0-cm) volume of A549 and H460 cells with a 6DoF setup error and varying distance from the isocenter in SRT 1 (1.0 cm). The first irradiation (n = 1) was performed at t = 0, the second irradiation (n = 2) at t = 1, and the third irradiation (n = 3) at t = 2. A 6DoF setup error for the GTV occurred, and the greater the reduction in the dose to the GTV was, the greater the residual tumor volume was.

Evaluation of GTV volumes with 6DoF setup errors based on the MCM
The GTV residual volume rate was calculated at the end of the third irradiation by using the GTV volume (Fig. 4). The GTV residual volume rate was increased with the increases in the distance from the isocenter and the 6DoF error for both A549 and H460 cells.
We used the GTV residual volume rate as the reference volume of GTVs before irradiation that GTV 1 a diameter of 1.0 cm, 523 mm 3 , GTV 2 has a diameter of 2.0 cm, 4187 mm 3 , and GTV 3 has a diameter of 3.0 cm, 14,130 mm 3 . In the A549 cells, the residual volumes rates were 5.4%

The distance between the isocenter and the GTV required to satisfy the tolerance values for the residual volume for A549 and H460 cells
We derived the distances between the isocenter and the GTV required to satisfy the tolerance values (10%, 35%, and 50%) using the GTV residual volume rate calculated by the MCM. Table 2 shows the distances necessary to satisfy the tolerance values for each GTV with varying GTV sizes and 6DoF setup errors for A549 cells. When the tolerance value was 10%, the distances were 2.0 cm and 0.3 cm for GTV 1 and In the (T, R) = (0.5 mm, 0.5°) condition of GTV 3, the distance was > 10 cm, and thus no distance was calculated ( Table 2).
As shown in Table 3

Discussion
In single-isocenter irradiation, the physical dose to the GTV is expected to decrease due to various factors such as the GTV size, the distance from the isocenter, and the 6DoF setup error. It is necessary to evaluate the radiobiological effect of the reduced physical dose on the GTV since the therapeutic effect of radiotherapy depends on whether or not the tumor is lethal [39]. Mathematical ODE models have proven to be a concise and powerful tool in for the analyses of various phenomena in tumor cells [16][17][18][19][20][21][22][23][24][25][26]. It is possible to analyze the status and development of tumors by using appropriate mathematical models and arbitrary biological data. The use of tumor data related to radiation would make it possible to simultaneously evaluate the radiobiology of tumors, as in the present study, in addition to the conventional physical evaluation for tumors [40,41]. Figure 4 shows that there was difference in the residual volume between the A549 and H460 cells because the recovery and lethal effects of radiation on GTVs depend on the MCM parameters α and β. The α and β values represent the intrinsic radiation sensitivity of irradiated cells, with cells with higher α and β values being more sensitive to radiation. The ratio of α and β values (α/β) is a measure of the fractionation sensitivity of a cell; cells with a high α/β ratio are less sensitive to the sparing effects of fractionation. The α/β value can also be used to compare the effects of fractionated radiotherapy and the equivalent total dose for different fractionation schedules [25,31,32]. We observed that the α and β values for the radiosensitivity of A549 cells were greater compared to those of the H460 cells, resulting in a higher GTV residual volume rate. The difference in the residual volume rates of both the A549 and H460 cells were smaller in GTV 3 than in GTV 1 and GTV 2 (Fig. 4). As the absolute volume of the GTV increased, the difference between the residual volume rate of the A549 and H460 cells became smaller when it was evaluated based on the residual volume of the GTV, since the effect on the dose reduction to the GTV due to a 6DoF error is small.
We used the tolerance values 10%, 35%, and 50% to calculate the distance from the isocenter. The larger the tolerance value was, the more conditions were found to satisfy the tolerance value even at a greater distance from the isocenter (Tables 2, 3). The contents of Tables 2 and 3 show that when the tolerance values were set to 35% and 50%, the tolerance rate was met for almost all conditions (except for GTV 1) regardless of the distance from the isocenter. When the tolerance value was 35%, the distances that satisfied the tolerance value under the (T, R) = (0.5 mm, 0.5°) condition were 8.0 cm for the A549 cells and 9.0 cm for the H460 cells in GTV 1. Only the condition of (T, R) = (1.0 mm, 1.0°) did not meet the tolerance value when the tolerance value was 50% for GTV 1. The distance that was necessary to satisfy the tolerance value varied depending on the tolerance value to be set and on the type of tumor cells. The relationship between the residual tumor volume rate and local control   [42]. Our present findings using 50% as the tolerance value are similar to the clinical results reported by Kraft et al. We also compared the distance reported in a study that used a physical tolerance value of 10% with a 6DoF setup error (T, R) = (0.5 mm, 0.5°) for a 1.0-cm GTV with the distances that met the tolerance value obtained with the biological tolerance value of 10% in the present investigation [11]. We observed that shorter distances of 2.0 cm and 2.4 cm from the isocenter were required to satisfy the tolerances for A549 and H460 cells, compared to the 9.3 cm distance using physical tolerances. This result shows that the effect of a setup error at 6DoF for a GTV using the biological tolerance value was larger than that of physical tolerance value. A further accumulation of clinical data clarifying the relationship between the tolerance value of the GTV residual volume rate and therapeutic effects in single-isocenter irradiation for brain metastases is necessary. The radiobiology modeling was reported to have uncertainties. The K involved in the growth rate and carrying capacity of the tumor used in the MCM model has a significant impact on tumor growth. Enderling et al. evaluated the uncertainty inherent in K and calculated its impact on tumor growth [43]. The model in this study was optimized to match the measured values by referring to the values of K in NSCLC and constraining the range of values [44]. Although the tumor growth of A549 and H460 in this study was modeled using a single measured data set, tumor growth may differ even for the same tumor cells [45]. Therefore, the same tumor cells may have different values of K, which could affect the results of this study. The some studies have been made to overcome the uncertainties associated with radiobiological models. Vaghi et al. reported a dramatic improvement in both accuracy and precision by using Bayesian inference to resolve uncertainties in the radiobiological model [24]. In addition, Kosinsky et al. addressed uncertainty by simulating a total of 2000 virtual experiments, including 10 animals in each experiment, to match the typical size of the actual experimental group, while accounting for uncertainty and variability in the model parameters for each treatment scenario [46]. Further evaluation of this study is also considered necessary to account for radiobiological model uncertainties in the future. The reliability of the α/β ratio in the LQ model affecting tumor lethality has a strong influence on the estimation of radiotherapy effects. Nevertheless, it has been reported that α/β values have both uncertainty and a certain range of values [47]. In radiotherapy, a more detailed evaluation could be expected by using a mathematical model, it is necessary to obtain a further accumulation of cell-specific parameters such as α and β values that vary among each tumor cell; in addition, improvements of the model are considered necessary to evaluate the tumor survival rate and the degree of treatment effect using a mathematical model [48,49]. There have been studies of an improved LQ model that considers sublethal damage recovery, repopulation, and oxygen effects, and it was proposed that such a model may be used to calculate the effects of radiation on tumors more accurately [50][51][52].
Our study has the limitations to address. First, the 6DoF setup erros in this study were calculated unidirectional and systematic. Second, the use of the 6DoF setup error for each fraction means that the same value was repeatedly given to the GTV for the evaluation of the GTV dose, although setup errors are a combination of random and systematic errors. It is necessary to simulate clinical conditions by using the values of the 6DoF setup error during irradiation for further evaluations, since the GTV dose was evaluated by repeatedly giving the same 6DoF setup error value to GTV for each fraction. Third, only the D95 for the GTV was evaluated, and a D100 and other dose prescriptions were not examined. Forth, tumor repair and lethal effects after the end of irradiation were not evaluated, since the GTV residual rate was evaluated using the tumor residual volume at the end of SRT irradiation. Finally, we modeled the tumor volume of H460 cells, and the volume after the irradiation exceeded the tumor volume modeled under the same condition on GTV 3 (3.0 cm). We evaluated the tumor volume of H460 cells by extrapolation of the MCM under those conditions.

Conclusion
The smaller the GTV size and the larger the isocenter distance and 6DoF setup error, the shorter the distance might need to be set in evaluations of the residual GTV volume based on the multicomponent mathematical model for SRT with single-isocenter irradiation applied to multiple brain metastases.