Adaptive Surrogate-Assisted Algorithm for Sensitivity Analysis and Therapeutic Regime Mapping in a CAR-T Cell Therapy Model

1. Introduction

Chimeric antigen receptor T-cell (CAR-T) therapy can induce durable remissions in several refractory B-cell malignancies [1,2]. In solid tumours, however, its efficacy is limited by heterogeneous antigen expression, the tumour microenvironment and poor access to tumour cells [3]. At the same time, clinical studies of CLDN18.2-targeted CAR-T cells in gastrointestinal malignancies and GD2-targeted CAR-T cells in neuroblastoma have demonstrated the therapeutic potential of this approach [4,5] and have intensified interest in the mechanisms that determine whether treatment results in durable control, residual disease, or tumour regrowth.

The major limitations of CAR-T therapy for solid tumours include inadequate tissue infiltration, functional exhaustion of effector cells, antigen escape, and immunosuppression within the tumour microenvironment [6,7,8]. Spatial heterogeneity, hypoxia, and transport constraints can further impair cytotoxic activity and promote the persistence of resistant cell populations [9]. Because these mechanisms act jointly, computational analysis should consider not only the effects of individual parameters but also regions of parameter space associated with distinct tumour dynamics.

Mechanistic models provide a natural framework for such analyses. Ordinary differential equation (ODE) models have been used to investigate CAR-T-cell proliferation, patient preconditioning, exhaustion, and interactions between effector and tumour-cell populations [10,11,12], whereas combination-therapy models have been employed to analyse the joint effects of CAR-T cells and other anticancer interventions [13]. Models of solid tumours have addressed interactions with tumour stem cells and TGF-$\beta$-mediated immunosuppression [14], resistance mechanisms associated with antigen expression [15], and the spatial development of antigen escape in reaction–diffusion formulations [16]. Spatial models additionally link therapeutic response to cellular distributions and transport processes: following local administration, treatment outcome depends on tumour density and diffusivity, the delivery route, and the functional composition of the cellular product [17], whereas physiologically based pharmacokinetic models can be used to investigate CAR-T-cell delivery to solid tumours [18]. The importance of microenvironmental heterogeneity is also evident in diffusion-based models of tumour growth, in which spatial variation in the diffusion coefficient can alter tumour morphology [19].

A separate line of research concerns the integration of mechanistic models with data. Multiscale quantitative systems pharmacology models describe the kinetics of CAR-T-cell subpopulations [20] and support the transfer of information between experimental and clinical scales [21]. Machine-learning methods have also been used to infer latent, mechanistically interpretable parameters from clinical features [22]. However, once a model has been parameterised, its systematic exploration over a multidimensional input space still requires a large number of direct simulations. At this stage, the computational cost of a spatial model becomes a limiting factor for subsequent sensitivity analysis and response-map construction.

For spatiotemporal models, each such simulation requires the numerical solution of a nonlinear system of partial differential equations (PDEs), and the cost of a comprehensive investigation therefore increases rapidly with the dimensionality of the parameter space. The same set of simulations must, moreover, support two related tasks: quantifying how parameters contribute to variability in tumour burden and identifying qualitatively distinct dynamical regimes. Hereafter, a therapeutic regime map refers to a partition of the parameter space according to such outcomes over a specified time horizon. Global sensitivity analysis (GSA) characterises the influence of model inputs across prescribed ranges and distributions [23]; for independent inputs, first-order and total-order Sobol indices distinguish individual parameter effects from their overall contributions, including interactions [24,25]. In contrast to local perturbation around a single parameter configuration [26], global analysis covers the entire domain of interest, while the Morris method can be used for preliminary screening [27]. Variance-based analysis remains computationally demanding, and the resulting mechanism rankings depend on the chosen output, parameter ranges and input distributions. Sobol analysis has already been applied to a CAR-T model incorporating antigen heterogeneity and indirect killing of antigen-negative cells [28], further motivating the development of computationally efficient approaches for spatial models.

One way to reduce the number of evaluations of the original solver is to construct a surrogate model that approximates the input–output relationship from a limited set of simulations [29]. Polynomial chaos expansions [30], dimensionality-reduction techniques combined with kriging [31], and Gaussian processes (GPs) [32] have all been applied to complex models. Initial training designs can be generated using Latin hypercube sampling (LHS) [33], while GPs are particularly well suited to sequential refinement because they provide both predictive estimates and associated uncertainty [34]. This uncertainty can be used to direct new simulations towards insufficiently explored regions of the parameter space, although the quality of the resulting probabilistic predictions must be assessed independently [35].

Active learning has already been applied to GSA. Chauhan et al. proposed a strategy designed to refine main effects when estimating Sobol indices [36], whereas Shang et al. developed an approach based on an ensemble of polynomial chaos expansions [37]. In computational biology, surrogate models have also been used to recover global sensitivity patterns in complex agent-based models, including models of spatial tumour growth [38]. These studies indicate that surrogate quality should be assessed with respect to the downstream analytical task rather than solely through average predictive error: accurate reconstruction of individual outputs does not necessarily guarantee the stability of derived quantities, including sensitivity indices.

In parallel, sequential methods have been developed for identifying threshold regions and the surfaces that separate them. Probabilistic surrogates and sequential uncertainty reduction have been used to estimate failure probabilities and reconstruct threshold sets [39,40]; the relationship between active learning and sensitivity analysis has also been investigated in adaptive reliability analysis [41]. In the CAR-T setting, this creates a specific trade-off: GSA requires accurate approximation of the quantitative response throughout the input space, whereas regime mapping is particularly sensitive to errors near threshold boundaries. Uniform space filling is therefore not necessarily optimal under a limited computational budget, but excessive concentration of simulations near regime boundaries may compromise global approximation accuracy.

The computational foundation of the present study is a previously developed spatiotemporal CAR-T model that incorporates infiltration, functional exhaustion, antigen heterogeneity, immunosuppression, and hypoxia [42]. We developed and validated an adaptive surrogate-assisted algorithm that uses a limited number of PDE simulations to approximate model outputs, estimate global sensitivities and map therapeutic regimes. Selection of the next simulation is based on two components: epistemic uncertainty in the GP surrogate and the proximity of the predicted response to regime boundaries. Sobol indices do not enter the acquisition function directly; instead, their stabilisation is used together with independent validation and classification accuracy in the recommended stopping criterion.

The algorithm begins with a common maximin LHS design, after which GP surrogates for the final and minimum normalised tumour burdens are sequentially refined using additional PDE simulations. Adaptive enrichment is compared with a non-adaptive space-filling design under matched computational budgets and a shared initial sample. The final GP surrogates are additionally evaluated on a fully held-out blind test set and compared with linear and quadratic regression models, while the locations of predicted regime boundaries are independently verified by direct PDE simulations. Following validation, the surrogates are used to compute Sobol indices and construct pairwise regime maps for the most influential mechanisms. We assessed prediction error, the stability of Sobol indices and agreement with model-defined response classes.

The remainder of this paper is structured as follows. Section 2 introduces the spatiotemporal CAR-T therapy model and describes the surrogate-assisted computational framework, including adaptive sampling, sensitivity analysis, regime mapping, and independent validation. Section 3 presents the numerical results and evaluates surrogate accuracy, adaptive-sampling performance, global sensitivity patterns, and the structure of the response-regime maps. Section 4 discusses the methodological and biological implications of these findings and the main limitations of the study. Section 5 concludes the paper.

2. Materials and Methods

The overall computational workflow of the study is presented in Figure 1. The main stages comprise simulations of the full spatiotemporal model, construction and adaptive refinement of the GP surrogates, their independent validation, and subsequent sensitivity analysis and response-regime mapping.

2.1. Spatiotemporal CAR-T Therapy Model and Numerical Implementation

The computational basis of this study is a previously developed spatiotemporal reaction–diffusion model of CAR-T therapy for solid tumours [16,42]. In the present work, the structure of the biological interactions was left unchanged. We used the model as a computationally expensive simulator to analyse parameter effects on tumour burden and identify distinct response regimes.

Within the three-dimensional domain $\Omega ={[0,1]}^3$, we consider the densities of functionally active and exhausted CAR-T cells, C and E, tumour cells with high and reduced antigen expression, $T_A$ and $T_B$, as well as the fields of chemoattractant S, immunosuppression A, and hypoxia H. Defining $T_{\Sigma }=T_A+T_B$, the system is written as

$$\begin{array}{cc}\hfill {\partial }_{t}C=& \nabla ·(D_C\nabla C)-\nabla ·({\chi }_C\left(A\right)C\nabla S)+{\rho }_C{\displaystyle \frac{C{T}_A}{K_C+C}}-{\delta }_{CT}C{T}_A-k_e(A,H)C+{\displaystyle \frac{p_{r}E}{1+A}}-d_{C}C,\hfill \\ \hfill {\partial }_{t}E=& \nabla ·(D_E\nabla E)+k_e(A,H)C-{\displaystyle \frac{p_{r}E}{1+A}}-d_{E}E,\hfill \\ \hfill {\partial }_t{T}_A=& \nabla ·(D_T\nabla T_A)+r_T{T}_A\left(1-{\displaystyle \frac{T_{\Sigma }}{K_T}}\right)-\gamma C{T}_A-{\mu }_{AB}\left(C\right)T_A+{\mu }_{BA}\left(C\right)T_B-d_{T_A}{T}_A,\hfill \\ \hfill {\partial }_t{T}_B=& \nabla ·(D_T\nabla T_B)+r_T{T}_B\left(1-{\displaystyle \frac{T_{\Sigma }}{K_T}}\right)+{\mu }_{AB}\left(C\right)T_A-{\mu }_{BA}\left(C\right)T_B-d_{T_B}{T}_B,\hfill \\ \hfill {\partial }_{t}S=& D_S\Delta S+{\sigma }_S{T}_A-{\lambda }_{S}S,\hfill \\ \hfill {\partial }_{t}A=& D_A\Delta A+{\sigma }_A{T}_{\Sigma }+{\sigma }_M{\displaystyle \frac{H}{K_H+H}}-{\lambda }_{A}A,\hfill \\ \hfill {\partial }_{t}H=& D_H\Delta H+\eta T_{\Sigma }-{\lambda }_{H}H.\hfill \end{array}$$

(1)

The nonlinear dependencies are defined by

$$\begin{array}{cc}\hfill {\chi }_C\left(A\right)& ={\displaystyle \frac{{\chi }_0}{1+{\alpha }_{A}A}},k_e(A,H)=k_e^0(1+{\beta }_{A}A+{\beta }_{H}H),\hfill \\ \hfill {\mu }_{AB}\left(C\right)& ={\mu }_0+\nu C,{\mu }_{BA}\left(C\right)={\mu }_1{e}^{-\xi C}.\hfill \end{array}$$

(2)

A normalised computational form of the model was used. Zero-flux boundary conditions were imposed on the boundary of the domain. The initial tumour lesion was represented by a sphere of radius $0.2$ centred at $(0.5,0.5,0.5)$, with $T_A(\mathbf{r},0)=0.1$ inside the sphere and $T_B(\mathbf{r},0)=0$ throughout the domain. The initial CAR-T-cell density was set to $C=0.2$ within the layer $x\le 0.05$, restricted to a circular region of radius $0.2$ in the $(y,z)$ plane. Outside these regions, $T_A$ and C were set to zero, and zero initial conditions were used for E, S, A, and H.

The simulations used a uniform $60\times 60\times 60$ grid, a time step of $\Delta t=0.01$ and a simulation horizon of $T=30$ days. Time integration employed Strang splitting: the reaction–chemotaxis subsystem was solved explicitly using a second-order Runge–Kutta scheme, whereas the diffusion subsystem was treated implicitly using Jacobi iterations. The maximum number of iterations for the diffusion step was 60, with a stopping tolerance of ${10}^{-8}$. After each complete time step, the solution was projected onto the non-negative domain. The full solver was implemented in C++17 with OpenMP and is denoted below by $\mathcal{M}$.

The baseline coefficients define the reference parameter set around which the six mechanism multipliers were varied. This parameter set was not calibrated against clinical data. All computational experiments were performed on a workstation equipped with an Intel Core i9-9900KF 3.60 GHz processor and 64 GB of RAM.

2.2. Mechanism Parameterisation and Output Quantities

Rather than varying all model coefficients independently, we considered six dimensionless multipliers, each corresponding to a distinct biological mechanism:

$$\mathbf{s}=\left(s_{\mathrm{mig}},s_{\mathrm{cyt}},s_{\mathrm{exh}},s_{\mathrm{esc}},s_{\mathrm{env}},s_{\mathrm{hyp}}\right).$$

(3)

The chemotaxis multiplier $s_{\mathrm{mig}}$ scales ${\chi }_0$, while the random motility coefficient $D_C$ remains fixed. The remaining mechanisms are parameterised as follows

$$\begin{array}{cccc}\hfill {\chi }_0& =s_{\mathrm{mig}}{\chi }_0^{*},\hfill & \hfill D_C& =D_C^{*},\hfill \\ \hfill \gamma & =s_{\mathrm{cyt}}{\gamma }^{*},\hfill & \hfill k_e^0& =s_{\mathrm{exh}}{k}_e^{0,*},\hfill \\ \hfill {\mu }_0& =s_{\mathrm{esc}}{\mu }_0^{*},\hfill & \hfill \nu & =s_{\mathrm{esc}}{\nu }^{*},\hfill \\ \hfill {\sigma }_A& =s_{\mathrm{env}}{\sigma }_A^{*},\hfill & \hfill {\alpha }_A& =s_{\mathrm{env}}{\alpha }_A^{*},\hfill & \hfill {\beta }_A& =s_{\mathrm{env}}{\beta }_A^{*},\hfill \\ \hfill \eta & =s_{\mathrm{hyp}}{\eta }^{*},\hfill & \hfill {\beta }_H& =s_{\mathrm{hyp}}{\beta }_H^{*}.\hfill \end{array}$$

(4)

The value $s_j=1$ corresponds to the reference regime. All multipliers were assumed to be independent and uniformly distributed over the range

$$s_j\sim \mathcal{U}(0.6,1.4),j=1,\dots ,6.$$

(5)

This range represents model uncertainty around the reference parameterisation and is not interpreted as an experimentally established physiological interval.

The total tumour burden $B\left(t\right)$ and its normalised value $R\left(t\right)$ were defined as

$$B\left(t\right)={\int }_{\Omega }[T_A(\mathbf{r},t)+T_B(\mathbf{r},t)]d\mathbf{r},R\left(t\right)={\displaystyle \frac{B\left(t\right)}{B\left(0\right)}}.$$

(6)

For each simulation, we defined

$$Y_{\mathrm{end}}=R\left(T\right),Y_{\min }=\underset{0\le t\le T}{\min }R\left(t\right).$$

(7)

The quantity $Y_{\mathrm{end}}$ was used as the primary output for global sensitivity analysis, whereas $Y_{\min }$ was used together with $Y_{\mathrm{end}}$ to classify the tumour dynamics.

Four mutually exclusive computational response regimes were defined:

$$\mathcal{R}=\left\{\begin{array}{cc}\mathrm{DC},\hfill & Y_{\mathrm{end}}\le 0.70,\hfill \\ \mathrm{PR},\hfill & 0.70<Y_{\mathrm{end}}<1.00,\hfill \\ \mathrm{RL},\hfill & Y_{\mathrm{end}}\ge 1.00,Y_{\min }\le 0.95,\hfill \\ \mathrm{Failure},\hfill & Y_{\mathrm{end}}\ge 1.00,Y_{\min }>0.95.\hfill \end{array}\right.$$

(8)

Here, DC denotes durable control, PR denotes a partial response with a final tumour burden below baseline, RL denotes regrowth to or above baseline after a reduction of at least 5%, and Failure denotes the absence of a meaningful response. These categories refer exclusively to the behaviour of the mathematical model and should not be interpreted as clinical criteria for patient response.

2.3. Surrogate Model and Global Sensitivity Analysis

To reduce the number of evaluations of the full PDE solver, two independent Gaussian process (GP) models were constructed for $Y_{\mathrm{end}}$ and $Y_{\min }$ [34,35]. Before training, the six input variables were linearly rescaled to ${[0,1]}^6$, and each output was standardised. An automatic relevance determination (ARD) Matérn $5/2$ kernel was used:

$$k(\mathbf{z},{\mathbf{z}}^{\prime })={\sigma }_f^2\left(1+\sqrt{5}r+{\displaystyle \frac{5}{3}}{r}^2\right)e^{-\sqrt{5}r},r^2={\displaystyle \sum _{j=1}^6}{\displaystyle \frac{{(z_j-z_j^{\prime })}^2}{{\ell }_j^2}}.$$

(9)

The hyperparameters were estimated by maximising the marginal likelihood with three optimiser restarts; a small diagonal jitter of ${10}^{-8}$ was added for numerical stability.

Global sensitivity was evaluated for $Y_{\mathrm{end}}$ using the mean prediction of the final GP. First-order and total-order Sobol indices were calculated as

$$S_i={\displaystyle \frac{{Var}_{s_i}\left[\mathbb{E}(Y_{\mathrm{end}}\mid s_i)\right]}{Var\left(Y_{\mathrm{end}}\right)}},S_{T_i}=1-{\displaystyle \frac{{Var}_{{\mathbf{s}}_{\sim i}}\left[\mathbb{E}(Y_{\mathrm{end}}\mid {\mathbf{s}}_{\sim i})\right]}{Var\left(Y_{\mathrm{end}}\right)}}$$

(10)

using the Jansen scheme [23,24,25,43]. The base Sobol sample contained $N_{\mathrm{SA}}=16384$ points. For six inputs, direct application of the same scheme to the full PDE model would require

$$N_{\mathrm{SA}}(d+2)=16384\times 8=131072$$

(11)

full simulations. Uncertainty in the estimated indices was characterised using 95% bootstrap intervals based on 500 resamples.

2.4. Adaptive Training-Set Construction and Non-Adaptive Comparator

A common initial design consisting of 48 maximin Latin hypercube sampling (LHS) points was used [29,33]. These points were identical for the adaptive and non-adaptive branches. In addition, an independent validation LHS sample of 30 points was generated; it was not used for GP training or selection of new points and served exclusively to assess surrogate quality and the stopping criterion.

After training the GPs on the initial 48 simulations, adaptive points were selected sequentially from a fixed Sobol candidate set $\mathcal{C}\subset {[0,1]}^6$ containing 16384 points. The selection criterion combined epistemic uncertainty in the surrogate model with proximity to computational response-regime boundaries. The uncertainty component was defined as

$$U\left(\mathbf{z}\right)=\max \{{\tilde{\sigma }}_{\mathrm{end}}\left(\mathbf{z}\right),{\tilde{\sigma }}_{\min }\left(\mathbf{z}\right)\},$$

(12)

where the standard deviations were evaluated on the standardised output scale. For the three classifier thresholds, the standardised distances were defined as

$$\begin{array}{cc}\hfill d_{\mathrm{ctrl}}& ={\displaystyle \frac{|{\mu }_{\mathrm{end}}-0.70|}{{\sigma }_{\mathrm{end}}+\epsilon }},\hfill \\ \hfill d_{\mathrm{fail}}& ={\displaystyle \frac{|{\mu }_{\mathrm{end}}-1.00|}{{\sigma }_{\mathrm{end}}+\epsilon }},\hfill \\ \hfill d_{\mathrm{resp}}& ={\displaystyle \frac{|{\mu }_{\min }-0.95|}{{\sigma }_{\min }+\epsilon }}.\hfill \end{array}$$

(13)

Because the $Y_{\min }=0.95$ threshold separates RL from Failure only when $Y_{\mathrm{end}}\ge 1$, this boundary component was multiplied by a smooth logistic weight $p_F\left(\mathbf{z}\right)$, defined in Eq. (15). The resulting boundary-proximity score was

$$B\left(\mathbf{z}\right)=\max \left\{e^{-d_{\mathrm{ctrl}}^2/2},e^{-d_{\mathrm{fail}}^2/2},p_F\left(\mathbf{z}\right)e^{-d_{\mathrm{resp}}^2/2}\right\},$$

(14)

where

$$p_F\left(\mathbf{z}\right)={\left[1+\exp \left(-{\displaystyle \frac{{\mu }_{\mathrm{end}}-1}{{\sigma }_{\mathrm{end}}+\epsilon }}\right)\right]}^{-1}.$$

(15)

After min–max normalisation of U and B over the current candidate set, the acquisition function was defined as

$$A\left(\mathbf{z}\right)=0.5\overline{U}\left(\mathbf{z}\right)+0.5\overline{B}\left(\mathbf{z}\right).$$

(16)

At each iteration, the point with the largest value of A was selected, the PDE model was evaluated at the selected point, and both GPs were subsequently retrained. A maximum of 72 adaptive points was added, yielding a final adaptive training set of 120 simulations.

Diagnostic evaluations were performed after every 12 adaptive additions. The change in global sensitivity was defined as

$${\Delta }_S^{\left(k\right)}=\underset{i}{\max }\left\{|{\widehat{S}}_i^{\left(k\right)}-{\widehat{S}}_i^{(k-1)}|,|{\widehat{S}}_{T_i}^{\left(k\right)}-{\widehat{S}}_{T_i}^{(k-1)}|\right\}.$$

(17)

The recommended stopping criterion was considered satisfied after two consecutive checkpoints at which ${\Delta }_S<0.02$, the maximum NRMSE across $Y_{\mathrm{end}}$ and $Y_{\min }$ on the independent validation set was below 0.12, and regime-classification accuracy was at least 0.90. Evaluation of the stopping criterion began only after at least 24 adaptive additions. To examine algorithm behaviour beyond the point at which the criterion was first satisfied, sampling continued until the training set contained 120 points.

The non-adaptive space-filling comparator used the same initial 48 LHS points. A further 72 points were selected sequentially according to a greedy maximin criterion from a separate Sobol pool containing 32768 points, without using model outputs or GP predictions. The two branches were compared at identical training-set sizes and against the same independent validation set.

2.5. Regime Maps and Direct Verification

Following the final Sobol analysis, the mechanisms were ranked according to $S_{T_i}$. For the three mechanisms with the largest indices, all three pairwise two-dimensional sections were constructed, while the remaining multipliers were fixed at $s_j=1$. Each map was evaluated on a regular $61\times 61$ grid using the GP mean predictions for $Y_{\mathrm{end}}$ and $Y_{\min }$, followed by application of the classifier in Eq. (8).

For independent verification of the regime boundaries, 12 points were selected automatically from regions where adjacent grid points were assigned to different regimes. Additional direct PDE simulations were performed at these points; they were not added to the training set and were used solely to check regime classifications near the predicted boundaries.

In addition, a direct one-factor-at-a-time analysis was performed, comprising the reference regime and two simulations for each of the six mechanisms ($s_j=0.6$ and $s_j=1.4$), with all remaining multipliers fixed at unity. Thus, the one-factor analysis comprised 13 full-model simulations and was used to determine whether varying each mechanism increased or decreased the final tumour burden relative to the reference configuration.

2.6. Independent Blind Validation and Overall Computational Design

After completion of adaptive training-set construction, regime-map generation, and boundary-point selection, a separate fully held-out set of 30 maximin LHS points was generated using an independent random seed. These points were not used for GP training, hyperparameter tuning, adaptive point selection, evaluation of the stopping criterion, estimation of Sobol indices, or construction of the regime maps.

In addition to the GP models, two simple baseline models were trained on the final 120-point adaptive training set: multivariable linear regression and quadratic regression including all polynomial features up to second order. For each output, the three surrogate approaches were compared on the same independent test set. For $Y_{\mathrm{end}}$ and $Y_{\min }$, we calculated $R^2$, RMSE, MAE, the maximum absolute error, and the normalised error

$$\mathrm{NRMSE}={\displaystyle \frac{\mathrm{RMSE}}{s_Y}},$$

(18)

where $s_Y$ is the sample standard deviation of the corresponding PDE output on the test set. For the GP models, we additionally evaluated the proportion of true values lying within the intervals $\mu \pm 1.96\sigma$, while regime-classification accuracy was assessed for all models.

Accounting for the common initial design, the complete set of unique PDE simulations used in the analyses reported in this study comprised

$${48}_{\mathrm{initial}}+{72}_{\mathrm{adaptive}}+{72}_{\mathrm{space}-\mathrm{filling}}+{30}_{\mathrm{validation}}+{12}_{\mathrm{boundary}}+{13}_{\mathrm{one}-\mathrm{factor}}+{30}_{\mathrm{blind}}=277$$

(19)

full simulations. All intermediate results, simulation parameters, time series, trained models, and numerical data used to generate the figures were stored automatically and could be reused without rerunning the PDE solver.

The source code of the computational algorithm, including the PDE solver and the procedures for surrogate modelling, adaptive training-set construction, sensitivity analysis, and regime mapping, is available in the open GitHub repository: https://github.com/maximvpolyakov/CAR-T-surrogate-regime-mapping (accessed on 6 October 2026).

3. Results

3.1. Computational Response Regimes and One-Factor Effects

The reference configuration ($s_j=1$) corresponded to the partial response (PR) regime: the final normalised tumour burden was $Y_{\mathrm{end}}=0.9305$, while the minimum value was $Y_{\min }=0.8144$. To compare qualitatively distinct types of dynamics, representative trajectories for the four computational response regimes are shown in Figure 2. The reference trajectory showed a partial response (PR). In the durable control (DC) example, tumour burden declined and remained below baseline for the rest of the simulation. The regrowth (RL) example showed a transient reduction followed by an increase above baseline. In the Failure example, tumour burden remained at or above baseline.

Varying one mechanism at a time produced markedly different changes in final tumour burden (Figure 3). The largest variation in $Y_{\mathrm{end}}$ was associated with antigen escape and cytotoxicity. Reducing $s_{\mathrm{cyt}}$ to 0.6 increased the final burden to $Y_{\mathrm{end}}=1.1577$, shifting the system into Failure, whereas at $s_{\mathrm{cyt}}=1.4$ it decreased to 0.7749. Antigen escape exhibited the opposite trend: $s_{\mathrm{esc}}=0.6$ yielded $Y_{\mathrm{end}}=0.7226$, whereas $s_{\mathrm{esc}}=1.4$ increased it to 1.1189 and shifted the system into RL. Stronger chemotaxis improved the final outcome ($Y_{\mathrm{end}}=0.8812$ at $s_{\mathrm{mig}}=1.4$), whereas stronger exhaustion worsened it ($Y_{\mathrm{end}}=0.9815$ at $s_{\mathrm{exh}}=1.4$). Varying the immunosuppression and hypoxia multipliers over $[0.6,1.4]$ produced smaller changes in final tumour burden. Thus, the direct one-factor analysis already identified antigen escape and cytotoxicity as the most influential mechanisms, but did not characterise their contributions across the full six-dimensional parameter space.

3.2. GP Surrogate Accuracy and Independent Blind Validation

On the independent 30-point test set, the two final GP surrogates predicted $Y_{end}$ and $Y_{min}$ with low error (Figure 4). For $Y_{\mathrm{end}}$, the resulting metrics were $R^2=0.999998$, NRMSE $=0.00148$, and MAE $=2.01\times {10}^{-4}$; for $Y_{\min }$, they were $R^2=0.999884$, NRMSE $=0.01058$, and MAE $=8.28\times {10}^{-4}$. The corresponding maximum absolute errors were $6.75\times {10}^{-4}$ and $3.61\times {10}^{-3}$. Coverage of the 95% GP predictive intervals was 100% for $Y_{\mathrm{end}}$ and 96.7% for $Y_{\min }$.

All 30 test points were assigned to the same computational regimes as in the corresponding direct PDE simulations. The test set contained 3 DC, 14 PR, 10 RL, and 3 Failure cases, covering all four classes.

Linear regression was less accurate than the GP surrogates (Table 1, Figure 5). Quadratic regression also achieved high $R^2$ values; however, its NRMSE was approximately nine times higher than that of the GP for $Y_{\mathrm{end}}$ and approximately 2.5 times higher for $Y_{\min }$. The differences were substantially larger for the linear model: NRMSE reached 0.1008 and 0.0934, respectively, while regime-classification accuracy decreased to 93.3%. Like the GP, the quadratic model correctly classified all 30 blind points, but it was less accurate in its continuous predictions.

3.3. Efficiency of Adaptive Sample Enrichment

The adaptive and space-filling designs were compared at identical training-set sizes using the same 30-point validation set (Figure 6). The clearest advantage of adaptive selection was observed at intermediate computational budgets. With 72 training simulations, the adaptive GP achieved NRMSE values of 0.00269 for $Y_{\mathrm{end}}$ and 0.01571 for $Y_{\min }$, compared with 0.00316 and 0.02208, respectively, for the non-adaptive space-filling design. At 84 simulations, the adaptive strategy retained a modest advantage for $Y_{\min }$ (0.01457 versus 0.01554), although the difference between the strategies for $Y_{\mathrm{end}}$ had already become minimal.

With further expansion of the training set, the advantage of the adaptive strategy disappeared. At the maximum budget of 120 points, the space-filling variant achieved slightly lower NRMSE values: 0.00185 and 0.01029 compared with 0.00212 and 0.01467 for the adaptive variant. Regime-classification accuracy on the independent validation set remained 100% for both strategies at all checkpoint sample sizes and is therefore not shown separately in the figure. The recommended stopping criterion was first satisfied at 84 training simulations; continuation to 120 points was used to compare the strategies under the same full computational budget. Therefore, the principal advantage of the adaptive algorithm was not systematically lower error at the maximum number of simulations, but rather faster attainment of high predictive accuracy at an intermediate budget.

Projection of the training set onto the plane of the two most influential mechanisms shows that the adaptive additions did not simply provide uniform space filling (Figure 7). New points appeared both in regions of pronounced variation in $Y_{\mathrm{end}}$ and in areas that were insufficiently covered by the initial LHS design. This distribution is consistent with the combined acquisition criterion, which simultaneously accounts for GP epistemic uncertainty and proximity to computational regime boundaries. Because the figure represents a two-dimensional projection of a six-dimensional selection process, the spatial distribution of the points should be interpreted as an illustration rather than as an independent measure of algorithmic performance.

3.4. Global Sensitivity of Final Tumour Burden

Global sensitivity analysis identified two dominant mechanisms (Figure 8, Table 2). For antigen escape, the total-order index was $S_{T,\mathrm{esc}}=0.4944$, while for cytotoxicity it was $S_{T,\mathrm{cyt}}=0.4382$. The total-order indices for chemotaxis and exhaustion were substantially smaller, at 0.0340 and 0.0329, respectively, whereas those for immunosuppression and hypoxia did not exceed 0.0024.

The first-order indices were close to the corresponding total-order indices: $S_{\mathrm{esc}}=0.4913$ and $S_{\mathrm{cyt}}=0.4346$, while the sum of all first-order indices was 0.9962. This is consistent with a relatively small contribution of inter-parameter interactions to the variance of $Y_{\mathrm{end}}$ over the investigated range. The 95% bootstrap intervals for the first-order indices of immunosuppression and hypoxia included zero, reflecting uncertainty in the estimation of these small effects. The resulting global ranking was consistent with the one-factor analysis, in which the largest changes in final tumour burden were likewise observed when varying $s_{\mathrm{esc}}$ and $s_{\mathrm{cyt}}$.

3.5. Response-Regime Maps and Boundary Verification

Pairwise regime maps were constructed for the three mechanisms with the largest total-order Sobol indices: antigen escape, cytotoxicity, and chemotaxis (Figure 9). The $s_{\mathrm{esc}}$–$s_{\mathrm{cyt}}$ plane contained all four response regimes (Figure 9(a)). The combination of low antigen escape and high cytotoxicity formed the DC region; PR occupied most of the parameter plane, whereas increasing $s_{\mathrm{esc}}$ under insufficient cytotoxicity successively shifted the system into RL and, for the most unfavourable combinations, into Failure.

In the $s_{\mathrm{esc}}$–$s_{\mathrm{mig}}$ section (Figure 9(b)), DC, PR, and RL were present, whereas Failure did not arise within the investigated range. Stronger chemotaxis shifted the PR–RL boundary towards higher values of antigen escape, while a small DC region persisted at low $s_{\mathrm{esc}}$ and high $s_{\mathrm{mig}}$. In the $s_{\mathrm{cyt}}$–$s_{\mathrm{mig}}$ plane (Figure 9(c)), PR, RL, and Failure were observed: increasing cytotoxicity progressively reduced the regions associated with unfavourable regimes, whereas Failure was confined to combinations of low cytotoxicity and insufficient chemotaxis.

At all 12 verification points near the predicted boundaries, the GP-based classifier agreed with the direct PDE simulations. The maximum absolute errors relative to the PDE outputs were $2.24\times {10}^{-4}$ for $Y_{\mathrm{end}}$ and $6.43\times {10}^{-4}$ for $Y_{\min }$. Together with the independent test results, these checks support the use of the GP surrogates for sensitivity analysis and regime mapping within the investigated parameter domain.

4. Discussion

The proposed algorithm integrates global sensitivity analysis and therapeutic regime mapping using a shared, limited set of PDE simulations. The two GP surrogates trained on 120 simulations accurately predicted $Y_{\mathrm{end}}$ and $Y_{\min }$; the resulting classifier agreed with the PDE model at all 30 test points and 12 additional points near regime boundaries. These validation procedures are complementary: low average predictive error characterises approximation quality, whereas checks near threshold regions assess the correctness of regime assignment. This strategy is consistent with established approaches to GP validation and sequential refinement of threshold sets [32,35,39,40]. However, agreement at selected validation points does not guarantee accuracy over the entire transition surface.

The near-unity $R^2$ values should be interpreted in the context of a deterministic simulator and the restricted domain $\mathbf{s}\in {[0.6,1.4]}^6$. The results are consistent with a sufficiently smooth response surface that can be accurately approximated by a GP with a Matérn 5/2 kernel. Linear regression was markedly less accurate than the GP, whereas the quadratic model correctly classified all test points but exhibited larger errors in the continuous predictions. Classification accuracy alone is therefore insufficient for selecting a surrogate also used for variance-based sensitivity analysis.

The benefit of adaptive design depended on the computational budget. At intermediate stages, adaptive sampling reduced prediction error, particularly for $Y_{\min }$, and the stopping criterion was satisfied at 84 training simulations; an additional 30 validation simulations were used to monitor this criterion. At 120 training points, however, the non-adaptive space-filling design yielded lower NRMSE values. The present results therefore support the usefulness of adaptation under a constrained budget, but not its universal superiority, in agreement with studies of active learning for GSA [36,37]. The identical classification accuracy of the two strategies on the validation set likewise provides no evidence of an adaptive advantage with respect to this metric. Establishing the robustness of the observed effect will require repeated comparisons across different initial designs and separate assessment of the individual components of the acquisition criterion.

The principal computational saving arises from replacing repeated evaluations of the PDE solver with GP predictions. With $N_{\mathrm{SA}}=16384$ and six inputs, the Sobol scheme used here would require 131072 direct simulations. In the present study, the analysis was instead performed using a surrogate trained on 120 simulations; the overall budget of 277 simulations additionally included the comparator and independent validation runs. These figures quantify the reduction in the number of evaluations of the original model, but do not by themselves represent wall-clock speed-up once GP training and adaptive point selection are taken into account.

The principal mechanistic result is the dominance of antigen escape and cytotoxicity: their total-order Sobol indices for $Y_{\mathrm{end}}$ were 0.4944 and 0.4382, respectively. The directions of these effects are consistent with the established role of antigen loss or reduced antigen expression in resistance to CAR-T therapy [44,45] and with previous models of antigen heterogeneity [15,16,28]. The first-order indices summed to nearly one, indicating that parameter interactions contributed little to the variance within the investigated domain. Nevertheless, even in the presence of weak interactions, smooth variation in a quantitative response can cross classifier thresholds: all four regimes were present in the $s_{\mathrm{esc}}$–$s_{\mathrm{cyt}}$ section. Mechanism ranking and localisation of regime transitions therefore provide distinct but complementary information.

Chemotaxis illustrates this distinction particularly clearly. Despite its relatively small contribution to the variance of $Y_{\mathrm{end}}$, stronger directed infiltration shifted the PR–RL boundary towards higher values of $s_{\mathrm{esc}}$, whereas insufficient chemotaxis combined with low cytotoxicity expanded the Failure region. This dependence is consistent with the importance of CAR-T-cell delivery and functional activity in solid tumours [6,17,18,46]. Within the investigated sections, however, improved infiltration did not eliminate the adverse effect of pronounced antigen escape.

The biological interpretation depends on the reference parameter set, the independent uniform distributions of the six multipliers over $[0.6,1.4]$, the choice of $Y_{\mathrm{end}}$ and the 30-day simulation horizon. The dependence of GSA rankings on the analysis setting [26] is particularly relevant for immunosuppression and hypoxia: their low indices characterise the present computational scenario and do not negate the established importance of the tumour microenvironment and physical barriers [3,6,8,9]. The parameters and their ranges were not calibrated against clinical data; DC, PR, RL, and Failure are model-defined categories, and the two-dimensional maps represent sections through the parameter space with the remaining inputs held fixed. Thus, surrogate accuracy demonstrates faithful reproduction of the simulator, but does not establish the biological validity of its predictions.

A further limitation concerns uncertainty quantification. The bootstrap intervals for the sensitivity indices were obtained from the GP mean prediction and do not incorporate uncertainty in the surrogate itself. Future extensions should propagate surrogate uncertainty into the sensitivity estimates and consider acquisition strategies that explicitly account for uncertainty in the indices [36].

From a biological perspective, priorities for further validation include calibration against longitudinal data, refinement of parameter ranges and input dependencies, and assessment of the robustness of regime maps to the choice of outputs, thresholds, and observation horizon. A specific experimentally testable hypothesis emerging from the present analysis is that preservation of antigen expression reduces the level of cytotoxic activity required for tumour control; testing this relationship could provide a direct link between the computational regime maps and experimental design.

5. Conclusions

The adaptive algorithm combined global sensitivity analysis and therapeutic regime mapping using two GP surrogates trained on 120 PDE simulations. On an independent 30-point test set, NRMSE was 0.00148 for final and 0.01058 for minimum normalised tumour burden. Regime assignments agreed with direct simulations at all test points and 12 additional points near regime boundaries. Adaptive sampling reduced prediction error at some intermediate budgets, although the space-filling design was more accurate at 120 training points. The stopping criterion was first satisfied after 84 training simulations.

Within the investigated parameter domain, antigen escape and cytotoxicity dominated the variance in final tumour burden, while parameter interactions contributed little. Stronger chemotaxis shifted the boundary between partial response and regrowth despite its low global sensitivity index. Thus, global sensitivity ranking and regime-boundary analysis provide complementary information on the mechanisms governing therapeutic response.