Predictive digital twin for optimizing patient-specific radiotherapy regimens under uncertainty in high-grade gliomas

1. Introduction

A digital twin can be defined as a mathematical model (or a collection of models) that provides a virtual representation of a specific physical object (e.g., a tumor), updates its status by assimilating object-specific data (e.g., imaging and clinical measurements of tumor growth and radiotherapy response), predicts the behavior of the object under external actions (e.g., treatments), and enables decision-making to optimize the future behavior of the object (e.g., design an optimal radiotherapy plan maximizing tumor control and minimizing toxicities) (Rasheed et al., 2020; Kapteyn et al., 2021; Laubenbacher et al., 2021; Niederer et al., 2021; Wu et al., 2022). Recent work explores the use of digital twins in healthcare and medicine to perform simulations of cardiovascular diseases (Corral-Acero et al., 2020; Peirlinck et al., 2021), enable virtual reality for surgery (Ahmed and Devoto, 2021), and enable improved decision-making in clinical oncology (Hernandez-Boussard et al., 2021; Madhavan et al., 2021; Wu et al., 2022). However, most of the work on digital twins in medicine relies on deterministic implementations. Predictive digital twins account for uncertainty through a Bayesian framework (Kapteyn et al., 2021) and provide a computational environment to support personalized risk-based management of solid tumors supported by computer forecasts of biologically-inspired mechanistic models representing these diseases and their treatments. We propose a patient-specific predictive digital twin that can address the three critical needs of: (i) accounting for uncertainty in the mechanistic model parameters by continuous integration of incoming patient data, (ii) forecasting the parameter uncertainty to estimate risk associated with the therapeutic outcomes, and (iii) supporting risk-aware clinical decision-making under uncertainty. To illustrate the methodology, we present a predictive digital twin to enable the personalized monitoring and forecast of high-grade glioma (HGG) response to radiotherapy (RT), as well as the design of optimal adaptive RT regimens for individual HGG patients (see Figure 1).

HGG is a class of brain tumors that typically exhibit an aggressive, infiltrative behavior as well as high heterogeneity in both tumor physiology and cell composition (Omuro and DeAngelis, 2013). Patients with HGG are usually treated surgical resection to reduce the tumor burden and intracranial pressure followed by adjuvant radiotherapy (RT) and chemotherapy to target residual, unresected tumor cells. The standard-of-care (SOC) RT plans account for patient-specific heterogeneity in tumor shape and location through pre-treatment anatomical or structural imaging approaches, such as T₁- and T₂-weighted magnetic resonance imaging (MRI), that can identify the residual tumor after surgery and define a surrounding 2–3 cm margin. However, the SOC RT dose and schedule generally conform to the Stupp protocol (Stupp et al., 2005), which is informed from clinical studies involving large populations, consisting of 60 Gy delivered in 30 fractions of 2 Gy. A fundamental challenge in the treatment of HGG is that response to RT is highly variable from patient to patient due to the inherent heterogeneity in both cellular architecture and tumor micro-environment (Aum et al., 2014; Hill et al., 2015), which may ultimately lead to poor treatment outcomes. This heterogeneity in tumor physiology, growth, and radiation response characteristics translates into uncertainty in the therapeutic outcomes of HGG patients receiving a standard RT protocol. This work addresses the challenge of quantifying the uncertainty in tumor characteristics from limited patient-specific data early in the course of treatment, propagating this uncertainty via mechanistic models to determine the uncertainty in treatment outcomes, and optimizing adaptive RT plans to improve overall survival through individualized predictive digital twins.

Adaptive RT has been regarded as a promising strategy to account for heterogeneity in tumor response. The varying treatment regimens aim at adjusting the dose, number, and timing of the radiation fractions in the prescribed RT plan on a patient-specific basis to account for the individual response of their HGG to radiation. Adaptive RT plans can be guided by quantitative imaging measurements to improve target delineation to deliver radiation, enable the voxelization of the RT plan (i.e., dose-painting), and facilitate the adaptation of the treatment regimen in response to observed tumor dynamics (Raaymakers et al., 2009; Jaffray, 2012; Troost et al., 2015; Kong et al., 2017). Another approach for adaptive RT planning is through a predictive framework that can be constructed by leveraging patient-specific computational forecasts of HGG growth and response to RT, which are obtained via computer simulation of biologically-inspired mechanistic models informed by routine clinical and imaging data collected before and during the course of treatment (Rockne et al., 2009; Unkelbach et al., 2014; Le et al., 2016; Poleszczuk et al., 2018; Enderling et al., 2019; Subramanian et al., 2020; Hormuth et al., 2021a; Woodall et al., 2021; Zahid et al., 2021; Lorenzo et al., 2022). However, the tumor forecasts obtained with these models are estimated using a deterministic approach. Systematically accounting for data and model uncertainties has shown promise in Bayesian calibration of tumor models (Lima et al., 2017; Oden et al., 2017; Lipková et al., 2019; Lorenzo et al., 2022). The clinical translation of current computational technologies to forecast HGG response to RT needs to accommodate a probabilistic risk assessment of pathological and therapeutic outcomes to support clinical decision-making along with progressive update of model according to incoming data. Predictive digital twins are a step toward personalized medicine that can address patient-specific clinical decision-making while accounting for the underlying uncertainty.

We construct the digital twin on an MRI-informed biology-inspired mechanistic model, which describes HGG growth and response to RT, and a probabilistic graphical model (Kapteyn et al., 2021), which accounts for uncertainty quantification and enables the risk assessment of therapeutic outcomes during tumor forecasting and RT optimization. We use prior distributions of the mechanistic model parameters informed from clinical data to initialize the digital twin. The personalization of the digital twin starts with the Bayesian model calibration for assimilating patient-specific MRI data followed by risk-aware optimized RT regimens for each patient. We solve a multi-objective risk-based optimization under uncertainty (OUU) problem to provide a suite of optimal RT regimens that trade-off maximizing tumor control and minimizing the toxicity from RT. We use the α-superquantile risk measure (Rockafellar et al., 2000) that accounts for the magnitude of exceeding a set threshold in a particular tumor characteristic. The risk-based optimization formulation allows one to account for the patient's and treating physician's risk preferences. We illustrate the predictive digital twin using an in silico population of HGG patients that is constructed by pooling clinical data of MRI measurements of HGG cellularity and mechanistic parameter values describing the dynamics of HGG growth and RT response from the literature (Qi et al., 2006; Wang et al., 2009; Hormuth et al., 2021c). We investigate varying levels of total dose to analyze the ensuing suite of therapeutic planning options. Our results show that the optimal RT plans can lengthen time to progression (TTP) with respect to the SOC regimen for the same total RT dose. Furthermore, we show that for certain patients, the proposed digital twin can provide optimal RT plans achieving similar tumor control as the SOC regimen, while significantly lower toxicities by lowering the radiation dose. We show that the range of optimal solutions also provide options for patients with aggressive cancer with increased radiation doses. We demonstrate that the optimal RT plans are non-inferior to the SOC regimen in terms of tumor control and always maintain the total radiation dose within a clinically-admissible range.

The rest of this work is organized as follows. Section 2 describes the proposed digital twin methodology to create a patient-specific predictive digital twin for HGG growth and RT response under uncertainty. Section 3 illustrates the digital twin methodology for a cohort of in silico patients. Section 4 discusses the proposed predictive digital twin framework as well as the results of our computational study, the limitations of this work, and future lines of investigation.

2. Patient-specific digital twin methodology

In this section, we describe the methodologies used to create, update, and utilize a patient-specific predictive digital twin. We begin with an overview of the components comprising the predictive digital twin and ground them in the oncology setting in Section 2.1. We then describe the mechanistic model in Section 2.2 followed by the treatment control parameters in Section 2.3 and the generation of observational data for an in silico patient in Section 2.4. The Bayesian model calibration, the propagation of uncertainty for patient-specific prognosis, and the multi-objective risk-based optimization under uncertainty problem are described in Sections 2.5, 2.6, and 2.7, respectively. We detail the survival analysis method used to assess the performance of the optimal treatment plans in Section 2.8.

2.1. Predictive digital twin formulation

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.

2.2. Mechanistic tumor growth model comprising the HGG 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 θ:=${\left[\rho ,K,N_{\mathrm{\text{initial}}},{\alpha }_{\mathrm{\text{RT}}}\right]}^{\top }\in \Omega \subseteq {ℝ}_{+}^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 ${\alpha }_{\text{RT}}/{\beta }_{\text{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.

2.3. Treatment control: patient-specific RT treatment regimen

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 $\mathit{u}\in \mathcal{U}\subseteq {ℝ}_{+}^{n_u}$ 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).

2.4. Simulated patients as proxy for observational clinical data

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 $\mathit{o}=\left[o_{t_1},\dots ,o_{t_{n_{\text{obs}}}}\right]$ at times ${\left\{t_i\right\}}_{i=1}^{n_{\mathrm{\text{obs}}}}$ with an additive noise model as

where the noise ε_i follows a truncated normal distribution $\mathcal{T}\mathcal{N}(0,{\sigma }^2,-N(t_i;{\theta }_{\text{true}},u),+\infty )$ with a lower truncation bound of −N(t_i; θ_true, u) to avoid non-physical negative observations. The general truncated normal distribution definition $\mathcal{TN}(\mu ,{\sigma }^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.

2.5. Patient-specific tumor modeling via Bayesian model calibration

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 $\mathcal{P}(\theta )$ 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 $\mathcal{P}(\theta |\mathit{o})$ is the posterior distribution of the digital state conditioned upon the observed data, $\mathcal{P}(\mathit{o}|\theta )$ is the likelihood of the observed data given a digital state, and $\mathcal{P}(\theta )$ 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:

Step 1: After obtaining the post-surgery image for the patient, the first step updates the posterior of the initial tumor burden to use the truncated normal distribution associated with the observation o_t1 at t₁ = 0 days as the mean with σ² accounting for the measurement error as $N_{\mathrm{\text{initial}}}~\mathcal{TN}(o_{t_1},{\sigma }^2,max\left\{0,o_{t_1}-2\sigma \right\},o_{t_1}+2\sigma )$. The upper and lower truncation bounds are defined by ±2σ from the observational data. The lower truncation bound is further modified to avoid non-physical negative values of initial tumor burden by using max{0, o_t1−2σ}. The updated distribution of N_initial is used as the prior distribution in Step 2.

Step 2: The Bayesian inference is performed via Markov Chain Monte Carlo (MCMC) for assimilating the remaining two observations [o_t2, o_t3] from t₂ = 20 days and t₃ = 27 days to complete the digital twin model calibration for tumor physiology and radiation response. We use the updated distribution of N_initial from step 1 and the priors for ρ, K, α_RT given in Table 2. Then the patient-specific posterior distribution for all the probabilistic model parameters are obtained by using MCMC for the remaining two observations. We use the MCMC implementation in the PyMC3 Python package (Salvatier et al., 2016) with four chains of 100, 000 samples each. The posterior distribution is characterized by 100, 000 retained samples. For a more in-depth discussion of MCMC techniques, the reader is referred to (Smith, 2013) and the references within.

2.6. Patient-specific prognosis via uncertainty propagation

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 $\mathit{u}\in \mathcal{U}$ to a real value as $M:\mathcal{U}\times \Omega \mapsto ℝ$. 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 $\mathcal{R}:\mathcal{U}\times \Omega \mapsto ℝ$. 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 $F_{M(\mathit{u},\theta )}^{-1}$ is the inverse cumulative distribution function of M(u, θ). Then, we can define the α-superquantile ${\overline{Q}}_{\alpha }[M(\mathit{u},\theta )]$ as

where [·]⁺ = max(0, ·). When the cumulative distribution of M(u, θ) is continuous, we can view ${\overline{Q}}_{\alpha }[M(\mathit{u},\theta )]$ as the conditional expectation of the QoI conditioned on the QoI exceeding the α-quantile as ${\overline{Q}}_{\alpha }[M(\mathit{u},\theta )]=\mathbb{E}\left[M(\mathit{u},\theta )\mid M(\mathit{u},\theta )\ge Q_{\alpha }\left[M(\mathit{u},\theta )\right]\right]$. 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 (${\overline{Q}}_{\alpha }[M(\mathit{u},\theta )]>Q_{\alpha }[M(\mathit{u},\theta )]$).

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.

2.7. Optimal patient-specific RT treatment regimens under uncertainty

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\Vert {\mathit{u}}^{*}{\Vert }_1=D_{\text{max}}$) 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:

• accounts for uncertainty during clinical decision-making process to provide multiple patient-specific optimal treatment regimens

• leads to a suite of patient-specific optimal treatment regimens with different levels of trade-offs between tumor control and toxicity from RT to allow more flexibility to consider the patient's and the treating physician's preferences

• provides an optimal treatment regimen option for a patient that can achieve similar tumor control as the SOC, but with reduced RT dose to mitigate toxicity effects whenever possible

• provides an optimal treatment regimen option with increased RT dose for patients with aggressive cancer, where the SOC does not provide sufficient tumor control either for patient preference or to allow for further surgical intervention

2.8. Survival analysis

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.

3. Illustrative cohort of in silico patients

This section presents an illustration of the digital twin methodology for a cohort of in silico patients. We describe the details of the computational study in Section 3.1 followed by the results for Bayesian calibration (Section 3.2), risk-based OUU (Section 3.3), and survival analysis (Section 3.4). The code used for generating the results reported for the predictive digital twin is available at https://github.com/Willcox-Research-Group/predictive-dtwin-glioma-frontiers.

3.1. Computational study format

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:

(i) Early progressors: 16 patients with SOC TTP α-superquantile ≤ 1 month after end of RT

(ii) Intermediate progressors: 62 patients with SOC TTP α-superquantile between 1-3 months after end of RT

(iii) Late progressors: 22 patients with SOC TTP α-superquantile ≥3 months after end of RT

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.

3.2. Patient-specific Bayesian calibration of the digital twin

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.

3.3. Patient-specific treatment optimization under uncertainty

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.

3.4. Survival analysis of optimal treatments

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.

4. Discussion

Our digital twin methodology provides the computational and mathematical foundation to dynamically integrate patient data with any given model of tumor growth and response. The predictive digital twin can provide risk-aware personalized treatment plans and we demonstrate it through optimized RT treatment regimens for HGG. There is a rich literature of modeling glioma growth and response across different scales (Alfonso et al., 2017; Hormuth et al., 2022); however, there are limited approaches that integrate these models with patient-specific data to optimize or adapt therapy in a dynamic fashion while accounting for underlying uncertainty. For development and demonstration of this framework, we employed an ODE model describing the total tumor cell proliferation. Response to radiotherapy and chemotherapy were reduced to an instantaneous effect at the time of treatment. While the model's biological detail is sufficient to demonstrate from end-to-end the key components of the predictive digital twin, a more biologically complex model may be needed in the clinical setting to account for inter-tumor heterogeneity in treatment delivery, proliferation rates, and resistance to therapy. Brain tumor growth and response to treatment has been described via a more complex reaction-diffusion model (Swanson et al., 2000; Alfonso et al., 2017; Hormuth et al., 2021c) that describes the spatial-temporal evolution of tumors cells throughout the brain due to invasion (i.e., the diffusion term) and proliferation (i.e., the reaction term). Our predictive digital twin can be readily adapted to other more complex models and cancer types where the prerequisite spatial-temporal or temporal data are available. We plan to extend the predictive digital twin to a more descriptive partial differential equation model that can extract the spatial information, but comes with a much higher computational cost challenge. These partial differential equation models could incorporate other components of the tumor microenvironment including the vasculature (and the blood brain barrier) which fundamentally influences the delivery and efficacy of chemotherapy and RT (Hormuth et al., 2021b). We have previously explored these components at the pre-clinical (Hormuth et al., 2020) and clinical levels (Hormuth et al., 2021c) and it would be natural to apply our digital twin framework to these models. While this digital twin focused on the delivery of RT concurrently with chemotherapy, our predictive digital twin could be adapted to other treatment modalities (e.g., immunotherapy, tumor treatment fields, convection-enhanced delivery) provided the pre-requisite data and models are available. We will also explore the critical issue of when to get the patient in to gather more imaging data.

We stress that the overall motivation of this work was the development of a digital twin framework which could be applied to any biology-based model of tumor growth. However, clinical treatment should ultimately be based upon established clinical guidelines (Horbinski et al., 2023). Indeed, one limitation of this study is the use of an in silico patient cohort whose growth and response parameters are sampled from literature values and governed by our model. We would anticipate that when applied to actual patient data, the precise gains between different RT regimen would differ to those presented in this work due to variability in the patient population and limitations of our current model. Therefore, a clinical trial would be necessary to establish the benefit of digital twin adapted RT regimens.

The optimized treatment regimens identified here tended toward a more hypofractionated (i.e., larger dose per day over a shorter period of time) or intermittent paradigm for treatment. While not widely used in the standard-of-care setting for initial treatment, hypofractionated radiotherapy is safe and tolerable (Gutin et al., 2009) and employed in the recurrent setting and for patients with poor prognoses (Trone et al., 2020). Compared to standard treatment schedules, hypofractionated radiotherapy has an increased cell kill over a shorter period of time, potential immunogenic effects, but at a potentially increased neuro-toxicity (Hingorani et al., 2012; Reznik et al., 2018). In the context of first-line treatment following surgery for HGG, the benefits of hypofractionated radiotherapy are less clear but some studies have reported similar tumor control and a palliative benefit due to its condensed treatment schedule (Floyd et al., 2004). Other modeling approaches from Leder et al. (2014), Brüningk et al. (2021), and Randles et al. (2021) have also identified optimized therapies that do not conform to the standard treatment paradigm. Recently, in the pre-clinical setting, Randles et al. (2021) predicted improved outcomes with a hyperfractioned regimen (i.e., smaller doses more than once a day) compared to standard dosing and this was validated via experiments. At the clinical level, Brüningk et al. (2021) explored intermittent radiation therapy (i.e., 6 Gy × 1 day every 6 weeks) in comparison to hyperfractionated radiotherapy for recurring disease. Brüningk et al. (2021) observed similar tumor control for both scenarios, but prolonged time to progression when additional intermittent doses of radiotherapy were delivered. In a future work, we plan to explore different clinical objectives and treatment controls, which could lead to other types of treatment regimens with further improvement in treatment efficacy.

5. Conclusions

We have developed a predictive digital twin methodology that enables end-to-end uncertainty quantification and optimization of personalized treatment regimens under uncertainty. The predictive digital twins employ a Bayesian perspective to account for uncertainty through the entirety of the clinical process to capture our knowledge (or lack thereof) of the individual dynamics. We illustrate the effectiveness of the predictive digital twin to optimize RT treatment plans for individual high-grade glioma patients. This subclass of tumors are aggressive and highly heterogeneous in nature, motivating the need for patient-specific modeling and treatment planning. We applied the predictive digital twin to a cohort of in silico patients with response and growth characteristics assigned from literature studies. We calibrated the digital twin with patient-specific data and used the calibrated digital twin to inform risk-aware clinical decision-making. We solve a multi-objective risk-based optimization under uncertainty problem to provide a suite of treatment plans balancing tumor control and toxicity. The risk-aware patient-specific optimization of the dose for the in silico study demonstrated that we were able to extend the progression-free survival for the standard total dose of 60 Gy relative to the standard-of-care. Additionally, we were able to identify potential reductions in the total dose delivered for similar tumor control response as the standard-of-care. These results represent a first step toward a practical digital twin for treatment optimization and future studies should expand this approach to data collected in the clinical setting.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation. The code used for generating the results reported for the predictive digital twin is available at https://github.com/Willcox-Research-Group/predictive-dtwin-glioma-frontiers.

Author contributions

AC, GP, MK, and KW contributed to the computational methodology. AC generated the results and plots. AC, GP, and DH wrote the first draft of manuscript. All authors contributed to conception, design of the study, manuscript revision, read, and approved the submitted version.

Funding

AC and KW acknowledge support from DARPA grant number DE-AC05-76RL01830 under the Automating Scientific Knowledge Extraction and Modeling (ASKEM) program and Department of Energy (DOE) Grant DE-SC0021239. KW acknowledges support from AFOSR MURI grant FA9550-21-1-0084. GP acknowledges support from the Office of Science, Advanced Scientific Computing Research, U.S. Department of Energy Computational Science Graduate Fellowship under Award Number DE-SC0021110. GL acknowledges the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement no. 838786. TY acknowledges support from the National Cancer Institute R01CA235800, U24CA226110, U01CA174706, and CPRIT RR160005 and is a CPRIT Scholar in Cancer Research. DH acknowledges support from CPRIT RP220225. This report was prepared as an account of work sponsored by an agency of the United States Government. Neither the United States Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Author disclaimer

The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.

Footnotes

References