We adopt the mathematical abstraction for a predictive digital twin proposed in Kapteyn et al. (2021). We formally define a predictive digital twin in terms of the six quantities shown in Figure 2 representing the physical and digital twins: physical state, observational data, control inputs, digital state, quantities of interest, and reward, where each of these six quantities will be considered to vary with time. The physical twin refers to the specific patient and the physical state represents the physiology and anatomy of the patient, which is only partially and indirectly observable. The digital twin is characterized by a computational model or a set of coupled computational models that can represent the physical twin to the desired level. The digital state is considered to be the parameterization of the computational models comprising the digital twin. To understand the health of the patient and inform our digital state, we rely upon observational data obtained from the physical twin, such as data obtained from MRI. There is a strong relationship between the digital state and the observational data since the accuracy of the digital state depends on the type and quantity of the observational data. The control inputs in a clinical setting are the therapeutic decisions that influence the physical state of the patient, such as the dose and scheduling of medical interventions. The control inputs can also comprise scheduling decisions for observational data, such as when a patient comes in for imaging. We will use the predictive digital twin to inform our choice of control inputs. The quantities of interest (QOIs) are computational estimates of possibly unobservable patient characteristics, which are evaluated using the computational models underlying the digital twin, such as tumor cell count, time to progression and toxicity. The reward is used to quantify the performance of the patient-twin system and can encode the success or failure of the system after applying a specific control, such as the risk of under-treating a specific disease giving a therapeutic plan. These quantities describe an abstract coupled patient-twin system representing a patient physical twin and their associated computational digital twin.

The computational model underlying the predictive digital twin is the logistic growth model of tumor growth governed by the ordinary differential equation (ODE) (Benzekry et al., 2014; Jarrett et al., 2018; Hormuth et al., 2022)

where N(t) is the number of tumor cells at time t, N_initial is the initial tumor burden, ρ is the net proliferation rate of the tumor cells, and K is the carrying capacity of the tissue (i.e., the maximum number of tumor cells that can be sustained physically and biologically).

We model the effects of RT and concomitant chemotherapy as discrete treatment events that result in an instantaneous reduction of tumor cell count at treatment time, N_{post-treatment}, to a fraction of the pre-treatment cell count N_{pre-treatment}, given by

where S is the surviving fraction of tumor cells resulting from a single dose of RT and chemotherapy at time t, which is given by u_t. The surviving fraction is defined by a linear-quadratic model of cell survival (McMahon, 2019) with a multiplicative term to account for the concurrent chemotherapy effect as

where S_C is the surviving fraction resulting from a single dose of chemotherapy, S_RT is the surviving fraction resulting from a single dose of RT, and α_RT and β_RT are the radiosensitivity parameters.

The entire model formulation given by Equations (1)-(3) can be written as

where N(t; θ, u) is the number of tumor cells at time t given that θ:=[ρ,K,Ninitial,αRT]⊤∈Ω⊆ℝ+4 are the patient-specific probabilistic model parameters in the digital state and u is a vector comprising all the RT treatment doses u_t given at time t. The treatment control considered here adapts the RT dose u to specific patients as elaborated in Section 2.3. All the parameters necessary to define the ODE model are provided in Table 1. We fix the radiosensitivity parameter ratio to αRT/βRT=10 (Rockne et al., 2009) and the surviving fraction resulting from chemotherapy to S_C = 0.82 (Hormuth et al., 2021c). The system of ODEs in Equation (4) is solved via a forward Euler scheme with sufficiently small time step size of 0.2 days to ensure numerical stability. The discrete treatment events are applied at the beginning of each day.

The treatment goal considered for this digital twin illustration is to adapt the RT regimen to specific patients. The Stupp protocol (Stupp et al., 2005) is the current SOC for treatment of HGG. The SOC administers six weeks of treatment featuring five consecutive treatment days each week with a fractionated RT dose of 2 Gy/day. We consider a similar setup as the SOC dose fractionation to define the treatment control u∈U⊆ℝ+nu as a vector of weekly-fractionated RT dose with n_u = 6 being the number of weeks of RT. Under this setup, the SOC is represented by u_SOC = [2, 2, 2, 2, 2, 2] Gy/day leading to a total dose of 5∥u_SOC∥₁ = 60 Gy, where ||·||₁ denotes the ℓ₁ norm (i.e., the summation of the absolute value of the vector entries). The SOC is used as the baseline to compare with the optimized RT treatment regimens. We administer chemotherapy throughout the six weeks even when no RT is prescribed. Note that extensions of this work can consider a finer time scale, such as allowing for a varying dose on a daily basis within each week of the treatment as well as adapting the chemotherapy treatment.

To calibrate the RT-related parameters in the predictive digital twin, we need at least one observation after starting the RT treatment (see Section 2.5 for details on calibration). Thus, we consider that the patients receive the SOC RT plan during this first week and that the intra-treatment MRI is performed at the end of the first week of RT (here, simulated patients as described in Section 2.4). We fix the fractionated dose of the first week to the SOC value of 2 Gy/day and use the predictive digital twin to optimally control the RT doses for the remaining five weeks (i.e., u = [u₁ = 2, u₂, …, u₆], where u_i ∈ [0, 10] Gy/day, i = 2, …, 6). Importantly, this strategy to optimally adjust the RT plan with our predictive digital twin follows the same approach as image-driven adaptive RT methods that are currently being explored in the clinical-research setting (Nijkamp et al., 2008; Sonke et al., 2019).

The physical twin state of each patient is partially and indirectly observed through specific observational data, which are leveraged to inform the digital twin state. In the HGG setting, these observational data can be acquired non-invasively using MRI. Recently, MRI data have been used to calibrate computational models and obtain tumor forecasts (Hormuth et al., 2015, 2020, 2021c,d) after appropriate post-processing. In particular, the MRI data needs to be post-processed to extract the total tumor burden as a cell count, such that it relates to the state variable N of the ODE model in Equation (4). Since we are collapsing the spatial information of the tumor provided by the MRI data, we incur a volume-wise error in the measurement of the tumor burden. This source of error has been studied by both Mazzara et al. (2004) and Paldino et al. (2014), although it is difficult to quantify its effect on the model parameterization and ensuing forecasts.

We consider an in silico patient cohort where the physical state of the HGG tumors is simulated by generating noisy observations of the solution to the ODE model in Equation (4) using an underlying “true” parameter set θ_true, which is varied for each patient in the cohort. Note that θ_true is never seen by the predictive digital twin and the model parameters instead adopt a probabilistic formulation based on the noisy observations generated for each in silico patient. We simulate a collection of n_obs observations o=[ot1,…,otnobs] at times {ti}i=1nobs with an additive noise model as

where the noise ε_i follows a truncated normal distribution TN(0,σ2,−N(ti;θtrue,u),+∞) with a lower truncation bound of −N(t_i; θ_true, u) to avoid non-physical negative observations. The general truncated normal distribution definition TN(μ,σ2,a,b) can be read as μ and σ² specifying the mean and variance, respectively, of the general normal distribution with the truncation range as [a, b] (Cohen, 1991). For this work, we assume a constant standard deviation of σ = 2 × 10⁹ cells, corresponding to a value of 10% of the mean value from population data used to model the initial tumor burden.

We use a Bayesian formulation that combines prior knowledge with observed data to update the probabilistic distribution placed on the model parameters. The observational data from the patient are assimilated by solving an inverse problem to calibrate the parameters in the computational model (Stuart, 2010; Biros et al., 2011; Oden et al., 2016) to the specific patient. The Bayesian framework provides confidence levels for the computational model output. Prior knowledge plays a critical role in the oncology setting as it is an especially data-poor regime owing to the difficulty and expense of collecting quantitative patient information such as MRI data. Priors allow us to inject knowledge of the disease at the population level to inform the modeling process of a specific patient. Priors P(θ) for the probabilistic parameters θ are constructed from reported values of clinical data in the literature. The study by Wang et al. (2009) computed the global proliferation rate for a cohort of 31 patients. Qi et al. (2006) conducted a study to compute a plausible set of radio-sensitivity parameters for the linear-quadratic model. We define an independent truncated normal distribution based on the values from the literature as our prior distribution as shown in Table 2. A truncated normal distribution is used to enforce positivity as well as to account for physically reasonable upper and lower bounds on the model parameters.

Observational data o from the patient are assimilated to estimate the updated probability distribution of θ, otherwise known as the posterior distribution, through the Bayesian update formula given by

where P(θ|o) is the posterior distribution of the digital state conditioned upon the observed data, P(o|θ) is the likelihood of the observed data given a digital state, and P(θ) is the prior distribution that encodes our knowledge about the patient's state before any observations are made.

We consider a scenario with three available MRI observations on days 0 (post-surgery image), 20 (pre-RT image), and 27 (mid-RT image). As noted in Equation (5), we assume an additive noise model for our observational data generation. The observation on day 0 only provides knowledge about the initial tumor burden N_initial for the patient. The observation on day 20 provides information on all the probabilistic model parameters except the radiosensitivity parameter α_RT. The observation on day 27 is generated using the SOC fractionated dose and provides knowledge about all the probabilistic model parameters {ρ, K, N_initial, α_RT}. This leads us to solve the inverse problem for model calibration in two steps:

Consider an uncertain digital state defined by the posterior distributions of the model parameters generated by Bayesian model calibration as described in Section 2.5. At any point in time, we can propagate uncertainty forward in time using the model in Equation (4) to estimate our quantities of interest, which characterize the tumor control resulting from an RT regimen and the treatment toxicities. Hence, a QoI maps realizations of random model parameters θ ∈ Ω and a treatment control u∈U to a real value as M:U×Ω↦ℝ. There are a variety of clinically relevant QOIs that could be considered to characterize the tumor control resulting from a treatment regimen. These include the predicted tumor burden (e.g., defined through the tumor cell count) after a set period of time post-treatment, time to progression (TTP) of the tumor cell count beyond a given threshold, or the amount of time that the tumor cell count is kept below a specified threshold. In this work, the QoI for tumor control is based on the TTP (Saad and Katz, 2009). We consider toxicity to be proportional to the total dose of radiation delivered to the patient during RT, i.e, 5∥u∥₁. Note that the QoI for toxicity does not depend on θ and, thus, it is a deterministic quantity since we have no uncertainty in the dose administered to a patient.

We define the TTP as the amount of time that the tumor takes to proliferate to a clinically-relevant size after the conclusion of six weeks of RT, which we set as the tumor cell count right before the onset of the RT regimen (i.e., at t = 20 days). Hence, to calculate the TTP, we further define a threshold tumor cell count N_th(θ): = N(t = 20;θ, u), which does not depend on u since it is a pre-treatment quantity. Additionally, we denote the day of the conclusion of RT as t_post-RT (here, t_post-RT = 20 + 6 × 7 = 62 days). Since a second round of chemotherapy is generally prescribed to prevent or combat tumor progression three months after the conclusion of RT, we consider a finite time-horizon of t_finite = 152 days after surgery for the end of our simulations and the calculation of the TTP. We can now express the TTP as

We want to maximize the TTP, or equivalently, minimize the negative of TTP. Thus, we define our QoI for tumor control as M(u, θ) = −T_TTP(u, θ) to be compatible with the definition of a general minimization problem (see Section 2.7).

In general, the nonlinear effects of the underlying computational model and the non-parametric distributions of θ prohibit one from obtaining an analytic expression for the distribution of the QoI, hence we rely on the Monte Carlo sampling approach. We sample a realization of the parameter set θ from the posterior distribution describing a patient's digital state, solve the forward model in Equation (4) to obtain the TTP, and compute the QoI M(u, θ). The process is continued till a desired number of samples n_MC of θ are processed to obtain the realizations of TTP as shown in Figure 3A. The distribution of TTP and the connection with M(u, θ) is illustrated in Figures 3B, C.

To make decisions concerning the patient-specific design of optimal treatment regimens, there are different statistics associated with the distribution of M(u, θ) that can be used as the reward function. We use a risk measure as the statistic to characterize the risk in tumor control associated with a given treatment regimen u. The risk measure is denoted by R:U×Ω↦ℝ. Different types of risk measures can be used to quantify the risk associated with tumor control, such as the probability of exceeding a threshold TTP value and the α-quantile for a set risk level α. We use the α-superquantile risk measure (Rockafellar et al., 2000; Rockafellar and Royset, 2013, 2015; Kouri and Surowiec, 2016), which has certain desirable mathematical properties. The α-superquantile satisfies two notions of certifiability in risk-based OUU (Chaudhuri et al., 2022): (i) accounting for near-failure and catastrophic failure events, and (ii) preserving convexity of the underlying QoI to aid in convergence guarantees for risk-based optimization formulations. The α-superquantile risk measure accounts for the magnitude of the largest 100(1−α)% realizations that captures the specified portion of the worst-case scenarios through risk level α. In oncology, it is important to account for the extent to which a tumor progresses beyond a specified threshold or remains in remission but close to the threshold rather than just the frequency of exceeding the threshold. Clinical ramifications of such extreme cases could result in under-treatment of aggressive disease resulting in poor tumor control and early disease progression. This is especially important in high-grade glioma where the prognosis is already overwhelmingly poor. The α-superquantile risk measure can take into account the magnitude of the TTP exceeding a set threshold to build in a desired level of conservativeness in a data-driven way.

To define the α-superquantile risk measure, we first define the related quantity of α-quantile Q_α for risk level α∈(0, 1) as

where FM(u,θ)−1 is the inverse cumulative distribution function of M(u, θ). Then, we can define the α-superquantile Q¯α[M(u,θ)] as

where [·]⁺ = max(0, ·). When the cumulative distribution of M(u, θ) is continuous, we can view Q¯α[M(u,θ)] as the conditional expectation of the QoI conditioned on the QoI exceeding the α-quantile as Q¯α[M(u,θ)]=E[M(u,θ)∣M(u,θ)≥Qα[M(u,θ)]]. We use Algorithm 1 in Chaudhuri et al. (2022) to estimate the α-superquantile through Monte Carlo simulations. Figure 3C shows an illustration of the data-driven conservativeness inferred from the Monte Carlo samples when using α-superquantile compared to α-quantile (Q¯α[M(u,θ)]>Qα[M(u,θ)]).

The α-superquantile risk measure is used to formulate a risk-based optimization problem under uncertainty to select the optimal treatment regimen as described in Section 2.7. The α-superquantile risk measure helps in making risk-aware decisions, where the risk preference for each patient can be adjusted by changing the value of α. Hence, larger values of α lead to more risk-averse decisions. We demonstrate our approach with a value of α = 0.95.

In this work, we consider two clinical objectives: (i) minimizing the risk associated with the QoI characterizing the tumor control, and (ii) minimizing toxicity through a proportional quantity to the total RT dose. To define personalized optimal RT regimens achieving these goals, we formulate a multi-objective risk-based optimization problem under uncertainty. This method results in a suite of eligible optimal RT regimens featuring different levels of trade-off between (i) decreasing toxicity by minimizing the total dose, or (ii) minimizing the risk associated with tumor control. The multi-objective formulation is given by

where optimized therapy comes into effect after the first week of RT and we optimize the RT doses for the remaining five weeks as denoted by u = [u₁ = 2, u₂, …, u₆] (see Section 2.3). Note that we can also explore the trade-off with changing risk preferences by varying the risk level α for the α-superquantile risk measure, as discussed in Section 2.6. An illustration of the predictive digital twin for HGG that shows the specific timeline is shown in Figure 4.

To solve the multi-objective problem in Equation (10), we reformulate it into a constrained single-objective problem using the ϵ-constraint method (Haimes, 1971; Hwang and Masud, 2012) as

where D_max is the total dose to be delivered. The value of D_max is varied over a range of total doses to obtain the Pareto optimal solutions (Hwang and Masud, 2012). Each of the Pareto optimal solutions provide different balance of toxicity and tumor control. Hence, the ϵ-constraint method is suitable for our application since we have a known primary objective of minimizing the risk associated with tumor control, while reducing the dose is a secondary objective for controlling toxicity. We also have a structured way to define the range of values for D_max based on the preferred total dose, which solves the challenge of selecting meaningful levels in the ϵ-constraint method. Thus, the specific advantages of the ϵ-constraint method in the context of the multi-objective problem defined in Equation (10) are: (i) preserving properties of the underlying problem, such as convexity of QOIs and risk measures, and (ii) ease of specification of the number of Pareto optimal treatment plans by selecting the number of D_max values and solving a single-objective constrained optimization problem each time.

The solution of Equation (11) will always lead to an optimal treatment regimen with active constraint (5‖u*‖1=Dmax) since increasing dose directly leads to decreasing risk. However, after an initial analysis of the optimization problem in Equation (11), we found that it leads to sub-optimal solutions for patients whose HGG exhibits a low tumor cell proliferation rate and/or low initial tumor burden. We observed that these patients do not require the maximum total dose enforced by an active constraint at the optimal solution to achieve the same amount of tumor control after the conclusion of the RT regimen. Typically for such patients, the amount of tumor control corresponds to a maximum possible value of 132 days for the TTP realizations based on the finite time-horizon of 152 days for our simulations. In other words, the same amount of risk associated with tumor control can be obtained while reducing the total dose below D_max for these patients. To account for such cases, we add an additional penalty term using parameter λ (here, λ = 0.001) on the total delivered dose as

The single-objective constrained optimization problem in Equation (12) is solved for various selections of the total dose threshold, D_max∈{40, 50, 60, 70, 80, 100} Gy, to obtain the Pareto front of optimal solutions. The SciPy Python package (Virtanen et al., 2020) is used for solving the optimization problem. We use the basin-hopping algorithm (Wales and Doye, 1997) for multi-start local optimization with the gradient-free COBYLA (constrained optimization by linear approximation) optimizer (Powell, 1994) from the SciPy package.¹ We use 20 restarts of COBYLA with a maximum of 200 function evaluations in each restart. The risk estimation in each function evaluation during the optimization uses n_MC = 5000 forward simulations.

The risk-based multi-objective problem formulation with the specific choice of objectives used here has the following desirable properties:

To assess the tumor control achieved with the different optimized treatment regimens for varying D_max and compare them against the SOC, we analyze the time to tumor progression across the in silico patient cohort using the Kaplan-Meier estimator (Kaplan and Meier, 1958; D'Arrigo et al., 2021). In this case, a right-censor criterion is enforced to account for the finite time-horizon of HGG response simulation at t_finite = 152 days, which leads to a maximum TTP value of 132 days. For a cohort of n_p patients, we define the probability of survival to tumor progression at time t as the probability of the α-superquantile value of the TTP being longer than time t. Thus, for each treatment obtained with a different D_max, the α-superquantile value of the TTP of each patient is used to calculate the survival probability as

where d(t_i) denotes the number of patients with TTP equal to time t_i and m(t_i) are the number of patients with TTP larger than time t_i−1 for all i = 0, …, n_p. A further advantage of the probabilistic modeling approach used in the predictive digital twin is the ability to estimate the variance in survival probability numerically from the samples of the TTP instead of relying on the approximation provided by the Greenwood's formula (Greenwood, 1926).

To assess the statistical significance of different treatment plans, the logrank test (Bland and Altman, 2004) is used to compare the survival distributions between digital-twin-based treatment plans and the SOC. The logrank test is implemented through the LifeLines Python package (Davidson-Pilon, 2019). Note that the null hypothesis for this test is that there is no difference between the populations in the probability of survival, so significant p-values indicate that the considered treatment is significantly different from the SOC.

We initialize an in silico cohort of 100 patients whose true physical states θ_true are determined by sampling from the marginal distributions of the population level priors (see Section 2.4). Noisy measurements are then generated from this cohort using Equation (5) to mimic the data that would be collected in the clinic. The specific observational and control timeline used for our predictive digital twins is outlined in Figure 4. Note that we are considering a finite time-horizon of three months post-RT (t_finite = 152 days) for our simulations and optimization. Thus, the maximum possible value of TTP is 132 days [see Equation (7)].

We divide the cohort of 100 patients into three groups based on the SOC TTP α-superquantile exhibited by the patient:

We first analyze the results for Bayesian calibration and risk-aware treatment plans produced by the patient-specific predictive digital twins for three patients with varied growth and response behaviors. The underlying physical state of the three patients is provided in Table 3. Note that these patient parameters are not known to the digital twin and are only used to generate observational data for the case-study. Patients 1 and 2 are intermediate progressors and Patient 3 is an early progressor according to the patient groups defined above. We then analyze the effectiveness of our predictive digital twins over the 100 patient cohort and the three patient groups.

The patient-specific calibrated digital twin is obtained by assimilating the three observational data available for each patient as described in Section 2.5, while the priors of all patients are obtained from population clinical data as described in Table 2. The histograms in Figure 5 show the prior and the posterior distributions for the probabilistic model parameters θ in the digital state for the three case-study patients. We observe that the posterior distributions concentrate around the (unseen) true physical state θ_true given in Table 3 (shown by dashed lines in Figure 5). We can also see that the overall uncertainty in the probabilistic parameters reduces compared to the prior distribution. In the clinical context, this could be interpreted as updating our belief about the specific patient's tumor dynamics as more data is collected during the clinical management of the disease.

To forecast model uncertainty, the uncertainty in the probabilistic model parameters is propagated forward by sampling 10, 000 parameter sets from the joint posterior obtained by MCMC and solving the mechanistic model in Equation (4) to generate the posterior trajectories. Figure 6 shows the posterior trajectories for the three case-study patients with the maximum and minimum bands. The posterior trajectories simulated from the calibrated digital twins exhibit tighter bounds (as shown by the orange posterior vs blue prior shaded regions) compared to the prior trajectories, which indicate a decrease in overall uncertainty. This is a consequence of our updated belief about the patient's tumor dynamics resulting from the integration of longitudinal data into the digital twin for each individual patient. We also see that the observations for each patient are captured within the shaded region of the respective posterior trajectories, which indicates the range of posterior trajectories. The quantification of uncertainty in model outputs is crucial for estimating the risk of tumor propagation used for risk-aware clinical decision-making.

We now present the risk-aware optimal treatment regimens obtained using the calibrated digital twins. We use the α-superquantile risk measure for quantifying risk of tumor growth as described in Section 2.6 and then solve a risk-based multi-objective problem to obtain a suite of patient-specific optimal treatment regimens under uncertainty as described in Section 2.7. Figure 7 shows the results of our patient-specific treatment optimization for the three case-study patients. We plot the TTP α-superquantile against the dose threshold to show the Pareto front indicating the trade-off between the two quantities. Figure 7 shows that the SOC regimen is not the optimal dosing schedule in multiple ways. First, for the applied radiation dose of 60 Gy, the optimal radiotherapy plan can demonstrate superior tumor control in terms of longer TTP. Second, similar or better tumor control can be exerted with a reduced amount of dose compared to SOC, as shown by the TTP values at the 40 Gy and 50 Gy levels. Third, greater tumor control than the SOC regimen can be exerted for total doses ≥60 Gy as seen in Figure 7.

Figure 8 shows the entire distribution of the TTP using the optimized treatment regimens compared to the SOC treatment for the three case-study patients. For Patient 1 in Figure 8A, we observe that with total dose constraint D_max≥60 Gy we see a progressive separation between the TTP distributions obtained from the SOC and the optimized treatment plans as we increase the total dose. The TTP distribution for optimized plans move toward higher TTP values. Thus, the optimal treatment plans lead to TTP distributions that are therapeutically superior than the ones obtained from the SOC. Additionally, considering total doses lower than the 60 Gy administered by the SOC, we observe that there is not much separation in the distributions of TTP. This result highlights the possibility of reducing the total dose administered in the SOC protocol by 10-20 Gy while obtaining a similar tumor control for Patient 1. A similar trend in separation of TTP distributions is seen for Patient 3, as shown in Figure 8C. In fact, the separation is obvious even at a total dose constraint of 50 Gy, for which the TTP distribution obtained from the risk-aware optimization already indicates superior therapeutic outcomes than the one obtained from the SOC. At higher doses, the TTP distributions from the optimized results shift toward higher TTP values for patient 3 and shows separation from the TTP distribution obtained from the SOC. For Patient 2 in Figure 8B, the most likely outcome is that the disease is controlled within the finite time horizon considered in this study (i.e., t_finite = 152 days). However, the optimized treatment plans still increase the probability of treatment success. This can be seen by the higher peaks at the “end of simulation” in Figure 8B.

The SOC and optimized treatment regimens for the three case-study patients are visualized in Figure 9, which shows the dose schedule over the six weeks of RT. At lower total dose levels, treatment plans tend to feature one large dose roughly halfway through the six week treatment course. At intermediate total doses, there is a secondary dose similar to the SOC in the second week of treatment and a large dose applied toward the end of the treatment timeframe. Finally, strategies featuring a higher total dose tend to feature one large dose in the second week of treatment and then another dose toward the end of the therapy. This second dose can be larger or lower than the one delivered the second week of treatment, but it tends to be larger than the corresponding SOC dose. Overall, Figure 9 shows that optimized plans tend to leave weeks without RT with only the chemotherapy being administered during that time. This can be viewed as an additional advantage with not requiring patients to come in every week.

We now analyze the performance of the predictive digital twins at the cohort level relative to the SOC and highlight the performance for the three patient groups defined in Section 3.1. To this end, we deploy the predictive digital twins for each of the 100 in silico patients to obtain patient-specific optimal RT treatment regimens. Figure 10 shows the boxplot of TTP α-superquantile values obtained for the patient-specific optimal treatment regimens and the SOC treatment. We denote the solution of the optimization problem in Equation (12) with D_max= X Gy as “OUU: X Gy” in Figure 10. We can see that the optimal treatment regimens lead to better TTP compared to the SOC over the cohort of 100 patients with D_max≥60 Gy. We further analyze the change in TTP α-superquantile given by [(TTP α-superquantile from “OUU:X Gy”) − (TTP α-superquantile from SOC)] for each patient in the different patient groups. Figure 11 shows the boxplots for the change in TTP α-superquantile values for early and intermediate progressors. Figure 11A captures the response to optimal treatment for early progressors. For optimal solutions “OUU: 60 Gy”, which use the same total dose as the SOC, we obtain a median change in TTP α-superquantile of +2.4 days with all 16 early progressors having better TTP α-superquantile values compared to the SOC. The higher dose optimal treatment regimens provide other therapeutic options to the patient and the treating physician with up to +8.5 days median change in TTP α-superquantile. Figure 11B captures the response to optimal treatment for intermediate progressors. In this case, the optimal solution “OUU: 60 Gy” leads to a median change in TTP α-superquantile of +7 days with all 62 intermediate progressors having better TTP α-superquantile values compared to the SOC. The higher dose optimal treatment regimens provide additional therapeutic strategies with up to +21.3 days median change in TTP α-superquantile.

The late progressors are not expected to have any reduction in TTP since almost all treatment options reach the TTP α-superquantile maximum value of 132 days. One can see that the definition of late progressors coincides with the finite time horizon considered. Instead, the late progressors can have similar tumor control with reduction in total required dose. Recall that to achieve this therapeutic objective of mitigating the toxicity effects with our digital twin formulation, we use the penalty parameter λ = 0.001 in Equation (12) to drive the optimizer to lowest possible dose while obtaining the same level of tumor control, as described in Section 2.7. Figure 12 shows that this approach achieves a reduction in total radiation dose while maintaining a similar treatment efficacy as SOC for the entire cohort of 100 patients as well as the three patient groups. We can see that risk-aware optimal treatment plans can lead to a median of 10 Gy (16.7%) reduction in total dose compared to the SOC over all 100 patients, while keeping the TTP α-superquantile within ±1 day of that obtained for the SOC. Analyzing the effect on the different patient groups provides additional insight on the possible amount of dose reduction depending on the underlying dynamics of tumor response. For the early and intermediate progressors, we can see a median reduction in dose of 10 Gy and 20 Gy, respectively. However, for the late progressors we can see that the median reduction in dose is 46.7 Gy. This comparatively higher decrease in dose is a direct consequence of slower tumor growth and lower number of tumor cells in the tumors of late progressors within the post-treatment timeframe considered in this study.

We now show the results for the Kaplan-Meier survival analysis described in Section 2.8. The results on survival analysis indicate that making optimal decisions at the individual level using patient-specific predictive digital twins provide either the same survivability as SOC plans or improve it across the cohort of patients. The survival curves in Figure 13 show that all optimal therapeutic regimens generated with a total dose threshold greater than or equal to that of the SOC plan of 60 Gy have survival curves that are superior to the corresponding curve obtained with SOC plan. There is a clear gap between the 80 Gy and 100 Gy plans and the SOC with logrank p-values of 0.005 and 0.0002, respectively. For the optimized therapy with maximum allowable dose of 60 Gy, the survival curve is statistically superior to the SOC treatment option with p-value of 0.05. Figure 13 further shows that the 40 Gy plans show almost no difference in tumor control when compared to the SOC. With a p-value of 0.35≫0.05, we conclude that the optimized treatment plans with maximum allowable dose of 40 Gy are not significantly different than the SOC, which uses 60 Gy.