A Hybrid Machine Learning Framework for Mechanistically Interpretable Latent Parameter Inference in a Spatiotemporal CAR-T Therapy Model for Solid Tumours

1. Introduction

Chimeric antigen receptor T-cell (CAR-T) therapy has become a major advance in modern cellular immunotherapy. In B cell haematological malignancies, it yields high rates of complete remission and durable disease control in patients with refractory or relapsed disease [1,2,3,4]. These outcomes have been made possible by the rational design of CAR constructs, the optimisation of co-stimulatory domains, and the refinement of cell manufacturing protocols [4,5]. However, translating these advances to solid tumours has proved considerably more challenging, necessitating the development of new biological and computational strategies to improve therapeutic efficacy [6,7,8,9,10,11].

Solid tumours are characterised by marked spatial heterogeneity, complex tissue architecture, limited availability of tumour-specific antigens, and temporal variability in antigenic profiles [6,7,8,11,12]. One of the principal mechanisms underlying failure of CAR-T therapy in this group of diseases is antigen escape, whereby selective pressure leads to loss or downregulation of target expression and the emergence of resistant clones [11,13,14]. An additional limitation is the risk of on-target off-tumour toxicity associated with the expression of potential targets in normal tissues [6,10,11,15]. Therefore, the development of CAR-T approaches for solid tumours requires the simultaneous consideration of antigen heterogeneity, safety, durability of therapeutic effect, and the potential for personalisation [12,16,17].

An equally important limitation is insufficient infiltration of CAR-T cells into tumour tissue. To exert their cytotoxic effects, effector cells must survive after administration, reach the tumour site, overcome vascular and stromal barriers, distribute throughout the tissue volume, and retain functional activity within an adverse microenvironment [6,7,10,15]. A major contributor to this problem is the mismatch between the chemokine profile produced by the tumour and the receptor repertoire of CAR-T cells, which reduces directed migration and accumulation of effector cells [18]. Accordingly, engineering strategies aimed at enhancing the recruitment and retention of immune cells within tumours are being actively developed, including the modification of CAR-T cells to secrete interleukin-7 (IL-7) and C-C motif chemokine ligand 19 (CCL19), which improves infiltration, survival, and antitumour activity in preclinical models [19]. More broadly, it is precisely these spatial constraints on CAR-T cell delivery and migration that make models capable of capturing tissue heterogeneity and local interactions between cellular populations particularly valuable [14,17]. In the wider context of spatial oncological modelling, diffusion-based approaches have already been used to describe tumour growth in heterogeneous environments, making it possible to account for the influence of tissue heterogeneity on tumour population dynamics [20].

A key determinant of therapeutic failure in solid tumours is the immunosuppressive tumour microenvironment. It encompasses hypoxia, nutrient deprivation, immunoregulatory cytokines, suppressive myeloid cells, regulatory T cells, and the expression of inhibitory ligands, which collectively impair the proliferation, cytotoxicity, and long-term persistence of CAR-T cells [6,7,15,16]. Under these conditions, CAR-T cells enter a state of functional impairment or exhaustion, characterised by reduced production of effector molecules, metabolic impairment, and sustained expression of inhibitory receptors [21,22,23,24,25]. In solid tumours, the hypoxic and metabolically adverse microenvironment is particularly important, as it limits the energetic capacity of T cells and promotes the development of an exhausted phenotype [26]. Thus, the spatiotemporal dynamics of CAR-T therapy are determined not only by direct cytotoxicity, but also by the balance among cell delivery, activation, exhaustion, and local mechanisms of immune suppression [11,16,17].

Despite these limitations, recent years have also brought encouraging clinical evidence supporting the fundamental applicability of CAR-T approaches to solid tumours. In recurrent high-grade gliomas, the safety and clinical activity of locoregional delivery of IL-13R$\alpha$2-targeted CAR-T cells have been demonstrated [27]. In gastrointestinal tumours, promising results have been reported for CLDN18.2-specific CAR-T cell therapy [28]. In hepatocellular carcinoma, clinical activity of GPC3-targeted CAR-T cells has been reported [29], while in neuroblastoma GD2 CAR-T cells have shown antitumour activity without pronounced on-target off-tumour toxicity [30]. These studies confirm the genuine clinical potential of CAR-T therapy for solid tumours; however, treatment response remains highly variable and depends on the biological characteristics of both the tumour and the patient, thereby necessitating a shift from population-level analytical frameworks to personalised models [10,12].

Against this backdrop, the role of quantitative and computational approaches capable of formalising CAR-T cell kinetics, tumour burden dynamics, and the influence of the microenvironment on therapeutic response is becoming increasingly important. Models of the cellular kinetics and pharmacodynamics of commercial CAR-T products have already been developed to describe cell expansion, contraction, and persistence following infusion [31,32,33]. To analyse interactions between tumour cells and CAR-T cells under controlled conditions, mathematical models based on real-time killing assay data have been proposed, enabling the respective contributions of CAR-T cell proliferation and exhaustion to be deconvoluted [34]. More sophisticated mechanistic approaches are also being developed that account for multiple CAR-T cell interactions with tumour cells, tumour spatial structure, antigen escape, and bystander killing effects [35,36,37]. In parallel, models of intracellular signal transduction and methods for predicting the in vitro potency of cellular products are being established, and these may be used for the parameterisation and validation of broader models of immune response [38,39]. Taken together, these studies show that mathematical modelling has already become an important tool for interpreting CAR-T therapy; however, most existing models are either focused on haematological malignancies or remain insufficiently personalised to describe solid tumours.

In parallel, computational oncology has seen the emergence of a research direction focused on the development of patient-specific mechanistic models, digital twins, and hybrid frameworks that combine mechanism-oriented modelling with machine learning methods [40,41,42,43]. Such approaches are particularly valuable in settings where experimental and clinical data are fragmentary and model parameters are only partially observable. Under these conditions, machine learning methods may be used not as a “black box” for direct outcome prediction, but rather as a tool for inferring the latent parameters of a mechanistic model from clinical, molecular, and experimental features [42,43]. This strategy has already been applied to the personalisation of reaction–diffusion models of tumour dynamics using machine learning methods [44]. In the context of CAR-T therapy for solid tumours, this is particularly important because therapeutic response is determined by a combination of difficult-to-measure characteristics, including the intensity of infiltration, the rate of exhaustion, the level of immunosuppression, antigen heterogeneity, and the efficiency of local cytolysis [13,18,25,26].

Thus, the contemporary literature points to a need for hybrid, mechanistically interpretable models of CAR-T therapy for solid tumours that account for the spatiotemporal structure of the tumour and enable limited patient-specific personalisation. Particularly promising is the use of published experimental and clinical data for population-level calibration of the mechanistic core of the model, followed by the inference of a limited set of latent patient-specific parameters using machine learning methods. This approach should be regarded as a feasibility study and a hypothesis-generating framework that supports the stratification of mechanisms of response and resistance, rather than as a clinically validated tool for precise individual prediction [34,40,42,43].

In this study, we propose a hybrid computational approach to the personalisation of a spatiotemporal model of CAR-T therapy for solid tumours. It is based on a mechanistic description of the interactions among CAR-T cells, tumour populations, and microenvironmental factors, complemented by machine learning algorithms for inferring patient-specific parameters from published data and clinical features.

The principal aim of this study is to establish a hybrid framework in which machine learning is used to infer a reduced set of patient-specific latent parameters of a spatiotemporal model of CAR-T therapy from available clinical data, while the mechanistic core is used to provide a biological interpretation of these parameters in terms of therapeutic response and mechanisms of resistance in solid tumours.

The remainder of this paper is organised as follows. Section 2 presents the materials and methods. Section 2.1 formulates the spatiotemporal reaction-diffusion-chemotaxis model of CAR-T therapy for solid tumours. Section 2.2 describes the numerical solution of the system and the software implementation of the computational algorithm. Section 2.3 defines the problem of patient-specific model personalisation. Section 2.4 characterises the source of the clinical data and the procedure used for variable preparation. Section 2.5 considers machine learning methods for inferring the hidden patient-specific parameters of the model. Section 2.6 specifies the validation protocol and the performance metrics used. Section 3 presents the results of the study. Section 3.1 analyses the accuracy of recovering patient-specific latent parameters. Section 3.2 considers an auxiliary direct empirical baseline for binary clinical outcomes. Section 3.3 provides the main mechanistic analysis of the latent factors underlying response and resistance. Section 4 discusses the findings, the limitations of the proposed approach, and directions for its further development. Section 5 sets out the principal conclusions of the study.

2. Materials and Methods

2.1. Mechanistic Spatiotemporal Model of CAR-T Therapy for Solid Tumours

As the mechanistic core of the personalised framework, we use a previously developed spatiotemporal reaction-diffusion-chemotaxis model describing the coupled dynamics of CAR-T cells, tumour subpopulations, and microenvironmental factors [45]. In this study, machine learning methods do not replace the underlying biology of the process, but are instead used to infer hidden patient-specific parameters within an interpretable system of equations.

We consider a bounded domain $\Omega \subset {\mathbb{R}}^3$ corresponding to the volume of the tumour tissue and its immediate microenvironment. The model variables depend on the spatial coordinate $\mathbf{x}\in \Omega$ and time $t\ge 0$. Let $C(\mathbf{x},t)$ denote the density of functional CAR-T cells, and let $E(\mathbf{x},t)$ denote the density of exhausted CAR-T cells. Furthermore, let $T_A(\mathbf{x},t)$ denote the density of tumour cells with preserved expression of the target antigen, whereas $T_B(\mathbf{x},t)$ denotes the density of tumour cells with reduced or lost antigen expression. In addition, we introduce the following microenvironmental fields: $S(\mathbf{x},t)$, the concentration of a chemoattractant mediating directed migration of CAR-T cells; $A(\mathbf{x},t)$, the concentration of a generalised soluble immunosuppressive factor (such as transforming growth factor-$\beta$ (TGF-$\beta$), interleukin-10 (IL-10), adenosine, among others); and $H(\mathbf{x},t)$, the normalised degree of hypoxia.

The governing system of equations takes the following form

Here, the auxiliary relationships are defined as follows

$${\chi }_C\left(A\right)={\displaystyle \frac{{\chi }_0}{1+{\alpha }_{A}A}},k_e(A,H)=k_e^0\left(1+{\beta }_{A}A+{\beta }_{H}H\right),{\mu }_{AB}\left(C\right)={\mu }_0+\nu C,$$

(8)

$${\mu }_{BA}\left(C\right)={\mu }_{1}exp(-\xi C),f\left(H\right)={\displaystyle \frac{H}{K_H+H}}.$$

(9)

Within the present model, infiltration is represented by a combination of diffusion and chemotaxis of CAR-T cells in heterogeneous tissue. CAR-T cell proliferation is induced by contact with antigen-positive tumour cells and is limited by saturation. Exhaustion is described as a microenvironment-regulated transition of functional cells to a dysfunctional state. Cytotoxic activity is implemented through the direct elimination of the $T_A$ subpopulation. Antigen escape is represented as phenotypic switching between $T_A$ and $T_B$, with both the direction and the intensity of switching depending on therapeutic pressure. The immunosuppressive microenvironment and hypoxia not only alter CAR-T cell motility, but also accelerate exhaustion, thereby establishing conditions conducive to resistance.

The system of Equations ()–() is supplemented by the following initial conditions

$$C(\mathbf{x},0)=C_0\left(\mathbf{x}\right),E(\mathbf{x},0)=E_0\left(\mathbf{x}\right),T_A(\mathbf{x},0)=T_{A,0}\left(\mathbf{x}\right),T_B(\mathbf{x},0)=T_{B,0}\left(\mathbf{x}\right),$$

(10)

$$S(\mathbf{x},0)=S_0\left(\mathbf{x}\right),A(\mathbf{x},0)=A_0\left(\mathbf{x}\right),H(\mathbf{x},0)=H_0\left(\mathbf{x}\right),\mathbf{x}\in \Omega ,$$

(11)

where the functions $T_{A,0}$ and $T_{B,0}$ describe the initial spatial structure of the tumour and its antigenic heterogeneity, while $C_0$ specifies the initial distribution of the infused cellular product. No-flux boundary conditions are imposed on the boundary of the domain $\Omega$

$${\displaystyle \frac{\partial T_A}{\partial n}}={\displaystyle \frac{\partial T_B}{\partial n}}={\displaystyle \frac{\partial E}{\partial n}}={\displaystyle \frac{\partial S}{\partial n}}={\displaystyle \frac{\partial A}{\partial n}}={\displaystyle \frac{\partial H}{\partial n}}=0,\mathbf{x}\in \Omega ,$$

(12)

$$\left(D_C\left(\mathbf{x}\right)\nabla C-{\chi }_C\left(A\right)C\nabla S\right)·\mathbf{n}=0,\mathbf{x}\in \partial \Omega ,$$

(13)

which corresponds to the absence of flux of cells and soluble factors across the outer boundary of the modelled volume.

The principal quantities used in the model and their biological interpretation are summarised in Table 1. Detailed parameter ranges are provided in [45].

Thus, in the proposed approach, personalisation is based not on direct empirical prediction of clinical outcome, but on the adaptation of biologically interpretable parameters of a spatiotemporal mechanistic model. This makes it possible to preserve the causal structure of the description, support mechanistically interpretable stratification of specific resistance mechanisms, and use the model as a foundation for further analysis of individual CAR-T therapy scenarios.

2.2. Numerical Solution of the Model and Software Implementation

The system of Equations ()–() is solved numerically using the finite difference method on a uniform three-dimensional rectangular grid with operator splitting by physical processes. This construction makes it possible to separate the diffusive component from the nonlinear reaction-chemotaxis terms, thereby providing an acceptable compromise among stability, accuracy, and computational cost.

The computational domain is defined as the rectangular parallelepiped

$$\Omega =[0,L_x]\times [0,L_y]\times [0,L_z],$$

(14)

which is discretised by a uniform grid

$$x_i=i\Delta x,i=0,\dots ,N_x-1,y_j=j\Delta y,j=0,\dots ,N_y-1,z_k=k\Delta z,k=0,\dots ,N_z-1,$$

(15)

where

$$\Delta x={\displaystyle \frac{L_x}{N_x-1}},\Delta y={\displaystyle \frac{L_y}{N_y-1}},\Delta z={\displaystyle \frac{L_z}{N_z-1}}.$$

(16)

In the implementation, the domain dimensions and grid resolution are specified by the parameters Lx, Ly, Lz and nx, ny, nz.

For the diffusion terms, a standard second-order central finite difference approximation in space is used. For any variable $u\in \{C,E,T_A,T_B,S,A,H\}$, the Laplacian operator is approximated by a seven-point stencil

$$\begin{array}{cc}\hfill {\left({\nabla }_h^{2}u\right)}_{i,j,k}\approx & {\displaystyle \frac{u_{i+1,j,k}-2u_{i,j,k}+u_{i-1,j,k}}{\Delta x^2}}+{\displaystyle \frac{u_{i,j+1,k}-2u_{i,j,k}+u_{i,j-1,k}}{\Delta y^2}}\hfill \\ \hfill & +{\displaystyle \frac{u_{i,j,k+1}-2u_{i,j,k}+u_{i,j,k-1}}{\Delta z^2}}.\hfill \end{array}$$

(17)

The chemotactic term in the equation for functional CAR-T cells is discretised in conservative flux form, thereby reducing numerical artefacts and remaining consistent with the physical interpretation of directed migration.

Time integration is based on Strang splitting. Denoting by

$$U=\left(\begin{array}{c}C\\ E\\ T_A\\ T_B\\ S\\ A\\ H\end{array}\right)$$

(18)

the vector of state fields, we represent the system in the form

$${\displaystyle \frac{dU}{dt}}=\mathcal{R}\left(U\right)+\mathcal{D}\left(U\right),$$

(19)

where $\mathcal{R}$ includes the reaction and chemotactic terms, whereas $\mathcal{D}$ includes diffusion. One time step of size $\Delta t$ is then computed as

$$U^{n+1}={\Phi }_{\mathcal{R}}\left({\displaystyle \frac{\Delta t}{2}}\right){\Phi }_{\mathcal{D}}(\Delta t){\Phi }_{\mathcal{R}}\left({\displaystyle \frac{\Delta t}{2}}\right)U^n.$$

(20)

The reaction-chemotaxis sub-system is integrated using an explicit two-stage second-order Runge-Kutta scheme

$$K_1=\mathcal{R}\left(U^m\right),K_2=\mathcal{R}(U^m+h{K}_1),$$

(21)

$$U^{m+1}=U^m+{\displaystyle \frac{h}{2}}(K_1+K_2).$$

(22)

To enhance stability for stiff parameter sets, the implementation allows the explicit step to be subdivided into substeps using the parameter exp_substeps. After each substep, a projection onto the physically admissible domain is performed to prevent the occurrence of negative concentration and density values.

The diffusion sub-system for each variable is solved implicitly. For an arbitrary field u, the following relation is used

$${\displaystyle \frac{u^{*}-u^{\mathrm{rhs}}}{\tau }}=D_u{\nabla }_h^2{u}^{*},\mathrm{or}\left(I-\tau D_u{\nabla }_h^2\right)u^{*}=u^{\mathrm{rhs}},$$

(23)

where $u^{\mathrm{rhs}}$ is the known value of the corresponding component of the solution, obtained after the preceding reaction-chemotaxis substep and used as the right-hand side in the implicit treatment of diffusion. The resulting linear system is solved using the Jacobi iterative method. The number of iterations and the stopping criterion are specified by the parameters diff_iters and diff_tol.

Reflecting boundary conditions are imposed for all variables. In the software implementation, this is achieved by referencing the nearest admissible grid node whenever an index falls outside the domain, which is equivalent to a zero normal derivative and, consequently, to zero flux across the boundary. The initial conditions are generated by a separate procedure. In the test configuration, a spherical tumour nodule located at the centre of the domain and a layer of CAR-T cells positioned near one of the boundaries are considered.

The numerical solver is implemented in C++, using the Params, Grid3D, Field, State, and Deriv structures to store parameters, grid geometry, and state fields. To accelerate computation, loops over grid nodes are parallelised using OpenMP. The results are saved in CSV format, with each record containing the time, node indices, coordinates, and the values of all seven variables $(t,i,j,k,x,y,z,C,E,T_A,T_B,S,A,H)$. This output organisation is convenient for subsequent visualisation, statistical analysis, and integration with external computational modules.

Figure 1 schematically shows the structure of the computational solver. At the highest level, the implementation comprises parameter and grid initialisation, generation of the initial conditions, the time-stepping loop, and result storage. Within each time step, Strang splitting is applied: the reaction–chemotaxis sub-system is integrated using an explicit RK2 scheme, whereas the diffusion sub-system is solved implicitly for each component of the state vector using the Jacobi iterative method.

To ensure reproducibility of the results, the full source code of the numerical solver has been made available in an open GitHub repository: https://github.com/maximvpolyakov/CAR-T-model (accessed on 6 April 2026).

2.3. Formulation of the Limited Inverse Problem of Patient-Specific Model Personalisation

In this work, personalisation is understood as a limited parametric adaptation of a spatiotemporal mechanistic model of CAR-T therapy to an individual patient. The objective is not to predict clinical outcome directly from observable features, but to estimate a set of latent parameters that enables the mechanistic description to be aligned with the available observations.

The full set of model parameters is represented as

$$\Theta ={\Theta }^{\mathrm{pop}}\cup {\Theta }^{\mathrm{pat}}\cup {\Theta }^{\mathrm{tr}},$$

(24)

where ${\Theta }^{\mathrm{pop}}$ denotes the population-level parameters, pre-calibrated using literature and experimental data, ${\Theta }^{\mathrm{pat}}$ denotes the patient-specific parameters, and ${\Theta }^{\mathrm{tr}}$ denotes the treatment-regimen parameters. Here, ${\Theta }^{\mathrm{tr}}$ includes the observed characteristics of CT041 administration and the associated initial conditions, whereas from the set ${\Theta }^{\mathrm{pat}}$ we identify a limited subset for personalisation

$${\tilde{\Theta }}^{\mathrm{pat}}\subseteq {\Theta }^{\mathrm{pat}}.$$

(25)

The reduced vector of patient-specific parameters is defined as

$${\tilde{\theta }}^{\mathrm{pat}}=\left({\theta }_{\inf },{\theta }_{\mathrm{exh}},{\theta }_{\mathrm{cyt}},{\theta }_{\mathrm{esc}},{\theta }_{\mathrm{env}},{\theta }_{\mathrm{init}}\right),$$

(26)

where ${\theta }_{\inf }$ characterises the effective infiltration of CAR-T cells into the tumour, ${\theta }_{\mathrm{exh}}$ denotes the propensity for functional exhaustion, ${\theta }_{\mathrm{cyt}}$ denotes the integral antitumour potency, ${\theta }_{\mathrm{esc}}$ denotes the intensity of antigen escape, ${\theta }_{\mathrm{env}}$ denotes the severity of the immunosuppressive microenvironment, and ${\theta }_{\mathrm{init}}$ denotes the individualised initial conditions associated with tumour volume and baseline antigenic heterogeneity. The exact composition of this vector may subsequently be refined in accordance with data availability and the results of identifiability analysis.

For patient p, two groups of observable quantities are distinguished. The first forms the feature vector

$$x^{\left(p\right)}\in \mathcal{X},$$

(27)

which may include tumour type, pretreatment characteristics, the level of CLDN18.2 expression, the CT041 dose, and other variables available before or at the early stages of therapy. The second group forms the vector of observed responses

$$y_{\mathrm{obs}}^{\left(p\right)}\in \mathcal{Y},$$

(28)

which comprises quantities comparable to the outputs of the mechanistic model, including the response category according to RECIST 1.1, characteristics of peripheral CAR-T cell kinetics, the grade of cytokine release syndrome (CRS), peak inflammatory markers, and time to event measures represented in a quantitatively interpretable form.

The personalisation problem is formulated as a limited inverse problem. The aim is to determine a vector ${\tilde{\theta }}_{\mathrm{pat}}^{\left(p\right)}$ that, for fixed ${\theta }^{\mathrm{pop}}$ and known ${\theta }_{\mathrm{tr}}^{\left(p\right)}$, provides the best possible agreement between the model and the patient’s observed data. To this end, the following Cauchy problem is solved:

$$U_t=\mathcal{M}\left(U;{\theta }^{\mathrm{pop}},{\tilde{\theta }}_{\mathrm{pat}}^{\left(p\right)},{\theta }_{\mathrm{tr}}^{\left(p\right)}\right),U(\mathbf{x},0)=U_0\left(\mathbf{x};{\tilde{\theta }}_{\mathrm{pat}}^{\left(p\right)},{\theta }_{\mathrm{tr}}^{\left(p\right)}\right),$$

(29)

where $U={(C,E,T_A,T_B,S,A,H)}^{\top }$ is the vector of state fields, and the operator $\mathcal{M}$ is defined by the system in Section 2.1. On the basis of the numerical solution, the observable functionals

$${\widehat{y}}^{\left(p\right)}=\mathcal{G}\left(U^{\left(p\right)}\right),$$

(30)

are computed and compared with $y_{\mathrm{obs}}^{\left(p\right)}$.

In general form, the parameter-estimation problem is written as

$${\tilde{\theta }}_{\mathrm{pat}}^{\left(p\right)}=arg\underset{\theta \in {\tilde{\Theta }}^{\mathrm{pat}}}{min}\left[\mathcal{L}\left(\mathcal{G}\left(U\left(\theta \right)\right),y_{\mathrm{obs}}^{\left(p\right)}\right)+\mathcal{R}\left(\theta \right)\right],$$

(31)

where $\mathcal{L}$ is the mismatch function between the model-derived and observed quantities, whereas $\mathcal{R}$ is a regularisation term that constrains biologically implausible parameter values.

Because the clinical data available from the literature are incomplete, personalisation in the present study is necessarily partial. Not all components of ${\tilde{\theta }}^{\mathrm{pat}}$ can be reliably inferred for every patient. Therefore, the primary focus is placed on those effective parameters that are both influential for the model and supported by the available variables described in Section 2.4.

In this formulation, machine learning methods do not replace the mechanistic model, but instead approximate the mapping

$${\mathcal{F}}_{\varphi }:\mathcal{X}\to {\tilde{\Theta }}^{\mathrm{pat}},$$

(32)

which maps the patient’s observable features to the reduced set of latent parameters. Accordingly, the overall study pipeline takes the form

$$x^{\left(p\right)}\stackrel{{\mathcal{F}}_{\varphi }}{\to }{\tilde{\theta }}_{\mathrm{pat}}^{\left(p\right)}\stackrel{\mathcal{M}}{\to }{U}^{\left(p\right)}\stackrel{\mathcal{G}}{\to }{\widehat{y}}^{\left(p\right)}.$$

At the first stage, the latent parameters are estimated; at the second, the individualised dynamics of the system are computed; and at the third, clinically interpretable output characteristics are extracted. The overall logic of the proposed personalisation framework is schematically shown in Figure 2.

2.4. Clinical Data Source and Variable Preparation

As the primary literature source, we used a clinically coherent study [46], which reported interim results of a study of CT041, a CLDN18.2-specific CAR-T therapy, in patients with previously treated tumours of the digestive system. This choice was driven not by an intention to confine the analysis to a single study, but by the need to minimise inter-study heterogeneity in reporting formats while at the same time preserving a minimally sufficient set of observations for the partial constraint of latent parameters in the mechanistic model and their subsequent clinical interpretation.

Only those variables that could be unambiguously matched to model variables or used as personalisation features were included in the working dataset. These comprised tumour type, the level of CLDN18.2 expression, information on prior therapy, the CT041 dose, response category according to RECIST 1.1, measures of objective response rate (ORR), disease control rate (DCR), progression-free survival (PFS), and overall survival (OS), the grade of cytokine release syndrome (CRS), characteristics of peripheral CAR-T cell expansion and persistence, as well as peak post-infusion inflammatory markers. Variables lacking a clear quantitative interpretation, temporal reference, or specification of measurement units were excluded from the analysis.

Information extraction was performed from the main text of the article, the tables, and the supplementary materials. Priority was given to data presented in tabular form or as explicit numerical values. Graph digitisation was used only when no tabulated representation was available. For each variable, the measurement units, temporal reference point, and level of data aggregation were recorded. After extraction, temporal characteristics were normalised relative to the first CT041 infusion, response categories were coded according to RECIST 1.1, survival endpoints were expressed in months, and the drug dose was retained as the absolute number of administered CAR-T cells.

A separate harmonisation procedure was applied to variables related to antigen expression and cellular kinetics. Because CLDN18.2 positivity in the original study was defined immunohistochemically as a staining intensity of at least $2+$ in no fewer than 40% of tumour cells, this variable was encoded as an ordinal categorical variable. For the quantitative analysis of CAR-T cell kinetics, the primary metrics used were $C_{max}$, $T_{max}$, and persistence duration. Inflammatory markers, including CRP, IL-6, and ferritin, were considered in terms of their peak post-infusion values.

The use of a single clinically coherent data source at this stage reduced variability associated with inter-study differences in design, endpoint composition, and completeness of reporting, thereby providing a more homogeneous basis for the initial feasibility assessment of the hybrid pipeline. Missing values were handled conservatively. When exact numerical values were unavailable, neither model-based nor multiple imputation was applied; instead, the corresponding variables were excluded from quantitative calibration and used only for qualitative interpretation.

As a result, a harmonised set of variables was assembled, including tumour type, pretreatment characteristics, CT041 dose, the level of CLDN18.2 expression, response category, ORR/DCR, PFS and OS, CRS grade, peak inflammatory markers, and characteristics of peripheral CAR-T cell kinetics. This dataset provides a reproducible foundation for population-level calibration of the mechanistic core, limited inference of latent patient-specific parameters, and their subsequent mechanistic interpretation.

After variable selection, harmonisation, and rescaling to a comparable range, the resulting patient-oriented feature set was further analysed by means of a heat map (Figure 3). This visualisation provides an intuitive overview of inter-patient variability and the structure of the harmonised set of clinical features across the cohort as a whole, and introduces the subsequent machine-learning stage, in which this same feature set is used to infer the reduced patient-specific parameters of the model.

2.5. Machine Learning for Inference of Patient-Specific Model Parameters

In this work, the machine learning module is used to approximate the mapping

$${\mathcal{F}}_{\varphi }:\mathcal{X}\to {\tilde{\Theta }}_{\mathrm{pat}},$$

which reconstructs the reduced vector of mechanistic model parameters from the patient’s clinical features. The target values are the reference estimates

$${\tilde{\theta }}_{\mathrm{ref}}^{\left(p\right)}=\left(\begin{array}{c}{\theta }_{\inf }^{\left(p\right)}\\ {\theta }_{\mathrm{exh}}^{\left(p\right)}\\ {\theta }_{\mathrm{cyt}}^{\left(p\right)}\\ {\theta }_{\mathrm{esc}}^{\left(p\right)}\\ {\theta }_{\mathrm{env}}^{\left(p\right)}\\ {\theta }_{\mathrm{init}}^{\left(p\right)}\end{array}\right),$$

obtained by solving the limited inverse problem for each patient.

Gradient boosting on decision trees was used as the primary method. This approach is well suited to small datasets, accommodates mixed tabular data, and is capable of capturing non-linear relationships among features. In addition, tree-based ensembles are typically more robust and interpretable than more complex neural network architectures under conditions of limited data availability. The model hyperparameters were tuned in the inner loop of nested cross-validation using GridSearchCV, with the mean absolute error as the optimisation criterion. The following hyperparameter grid was considered: number of trees $n_{\mathrm{estimators}}\in \{100,200\}$, maximum tree depth $d_{max}\in \{2,3\}$, learning rate $\eta \in \{0.05,0.1\}$, and the minimum number of samples per leaf $\in \{1,2\}$, with the subsampling parameter fixed at $\mathtt{subsample}=0.8$. Restricting tree depth and learning rate in this way was intended to reduce the risk of overfitting in the setting of a small sample size and to preserve the stability of out-of-fold estimates.

Separate regression models were trained for each component of the reduced vector. A distinct model was constructed for each component of the reduced vector

$$f_j(·;{\varphi }_j):\mathcal{X}\to \mathbb{R},j=1,\dots ,d,$$

so that

$${\mathcal{F}}_{\varphi }\left(x\right)=\left(\begin{array}{c}{f}_1(x;{\varphi }_1)\\ ⋮\\ f_d(x;{\varphi }_d)\end{array}\right).$$

This construction makes it possible to tune the models in accordance with differences in the identifiability and statistical properties of individual biological parameters.

To ensure compliance with biological constraints, each parameter ${\theta }_j$ is transformed into a continuous variable before model training. Let the admissible range of the parameter be given by the interval $[a_j,b_j]$. The following transformation is then used

$$z_j^{\left(p\right)}=log{\displaystyle \frac{{\theta }_{\mathrm{ref},j}^{\left(p\right)}-a_j+\epsilon }{b_j-{\theta }_{\mathrm{ref},j}^{\left(p\right)}+\epsilon }},$$

where $\epsilon >0$ is a small constant. After prediction, the inverse transformation is applied:

$${\widehat{\theta }}_j^{\left(p\right)}=a_j+(b_j-a_j)\sigma \left({\widehat{z}}_j^{\left(p\right)}\right),$$

where $\sigma (·)$ is the logistic function. This guarantees that the predicted values remain within the admissible biological domain.

To improve robustness to noise in the pseudo-labels, training is performed using the robust Huber loss function

$${\mathcal{J}}_j\left({\varphi }_j\right)={\displaystyle \frac{1}{|P_j|}}\sum _{p\in P_j}{w}_p{\ell }_{\delta }\left(f_j(x^{\left(p\right)};{\varphi }_j)-z_j^{\left(p\right)}\right)+\lambda \Omega \left({\varphi }_j\right),$$

where $P_j$ denotes the set of patients for whom an estimate of the parameter ${\theta }_j$ is available, $w_p$ denotes the observation weights, $\Omega \left({\varphi }_j\right)$ is the regularisation term, and ${\ell }_{\delta }$ is the Huber loss. This formulation reduces the influence of errors associated with incomplete clinical data and the inaccuracy of the inverse problem solution.

Continuous features are log-transformed, whereas categorical and ordinal variables are retained in their original form. Hyperparameter tuning is performed exclusively on the training subsets within nested cross-validation.

After training, the predicted parameter vector

$${\widehat{\tilde{\theta }}}_{\mathrm{pat}}={\mathcal{F}}_{\varphi }\left(x\right)$$

is passed to the mechanistic model in order to compute the system dynamics and the resulting prediction of therapeutic response. Thus, machine learning serves as a tool for personalising the latent parameters of the mechanistic model.

2.6. Validation Protocol and Performance Metrics for Latent Parameter Recovery and Baseline Prediction of Clinical Outcomes

Validation of the proposed approach in the present study is structured around two tasks, both directly reflected in Section 3. The first task is to assess the quality of out-of-fold recovery of the reduced vector of patient-specific latent parameters,

$${\tilde{\Theta }}^{\mathrm{pat}}=({\theta }_{\inf },{\theta }_{\mathrm{exh}},{\theta }_{\mathrm{cyt}},{\theta }_{\mathrm{esc}},{\theta }_{\mathrm{env}},{\theta }_{\mathrm{init}})$$

from the available patient clinical features. The second task is an auxiliary evaluation of a direct empirical baseline for binary clinical outcomes, aimed at determining what prognostic signal is contained in the observed features themselves, without the stage of inferring the latent parameters of the mechanistic model. The main mechanistic analysis of the latent factors underlying response and resistance, presented later in the Results section, is regarded as an interpretative stage and is not accompanied by a separate formal performance metric.

Because of the small sample size, evaluation is performed at the patient level using nested cross-validation. In the outer loop, K-fold cross-validation with patient-specific splitting is applied. All data from a given patient are assigned entirely either to the training set or to the test set, thereby preventing information leakage. Within the training portion of each outer split, all stages of model development are carried out, including feature transformation, handling of missing values, hyperparameter tuning, and training. Any preprocessing statistics are estimated exclusively on the training subset of the corresponding fold and are then applied to the test patients without further adjustment.

For the outer split m, the mapping

$${\mathcal{F}}_{\varphi }^{\left(m\right)}:\mathcal{X}\to {\tilde{\Theta }}^{\mathrm{pat}}$$

(33)

is trained. For a test patient $p\in {\mathcal{D}}_{\mathrm{test}}^{\left(m\right)}$, an out-of-fold estimate of the latent parameters is computed from the feature vector $x^{\left(p\right)}$ as

$${\widehat{\tilde{\theta }}}_{\mathrm{pat}}^{\left(p\right)}={\Pi }_{\mathcal{A}}\left({\mathcal{F}}_{\varphi }^{\left(m\right)}\left(x^{\left(p\right)}\right)\right),$$

(34)

where ${\Pi }_{\mathcal{A}}$ denotes the projection onto the set of biologically admissible values. The predicted values are compared with the reference estimates ${\tilde{\theta }}_{\mathrm{ref}}^{\left(p\right)}$, obtained independently by solving the limited inverse problem using the patient’s available observable data. These reference values are interpreted not as true biological parameters, but as pseudo-labels consistent with the mechanistic model and the available measurements.

For component-wise evaluation of the recovery quality of each latent component ${\theta }_j$, four metrics are used. The first metric is the mean absolute error

$${\mathrm{MAE}}_{{\theta }_j}={\displaystyle \frac{1}{N_{\mathrm{test}}}}\sum _p\left|{\widehat{\theta }}_j^{\left(p\right)}-{\theta }_{\mathrm{ref},j}^{\left(p\right)}\right|.$$

(35)

The second metric is the root mean square error

$${\mathrm{RMSE}}_{{\theta }_j}=\sqrt{{\displaystyle \frac{1}{N_{\mathrm{test}}}}{\sum }_p{\left({\widehat{\theta }}_j^{\left(p\right)}-{\theta }_{\mathrm{ref},j}^{\left(p\right)}\right)}^2}.$$

(36)

The third metric is the Spearman rank correlation coefficient ${\rho }_{{\theta }_j}$, which characterises how well the inter-patient ranking of the corresponding latent component is preserved. The fourth metric is the coefficient of determination

$$R_{{\theta }_j}^2=1-{\displaystyle \frac{{\sum }_p{\left({\widehat{\theta }}_j^{\left(p\right)}-{\theta }_{\mathrm{ref},j}^{\left(p\right)}\right)}^2}{{\sum }_p{\left({\theta }_{\mathrm{ref},j}^{\left(p\right)}-{\overline{\theta }}_{\mathrm{ref},j}\right)}^2}},$$

(37)

where ${\overline{\theta }}_{\mathrm{ref},j}$ denotes the mean value of the reference component across the test patients. The joint use of $\rho$, MAE, RMSE, and $R^2$ makes it possible to assess not only the absolute error, but also the preservation of rank structure and the proportion of explained variability.

In the second validation block, a direct empirical baseline is considered for two binary endpoints: achievement of an objective response and disease control. Its role is limited to that of an external reference point and does not replace the main mechanism-oriented study pipeline. The quality of binary prediction is evaluated on the basis of out-of-fold predictions using balanced accuracy, the $F_1$ score, AUROC, and AUPRC. Balanced accuracy is used as a more robust measure under conditions of potential class imbalance; the $F_1$ score characterises the trade-off between recall and precision for the positive class; AUROC assesses ranking quality across the full range of thresholds; and AUPRC is particularly informative for rare positive outcomes.

In addition, for each binary endpoint, agreement between the observed and predicted cohort-level frequency of the positive outcome is analysed. For this purpose, the observed proportion of patients with a positive outcome is compared with the corresponding proportion reconstructed from the out-of-fold predictions of the baseline model. This assessment does not replace individual classification metrics, but allows evaluation of how well the model reproduces the aggregate event frequency at the level of the entire cohort.

All final estimates are calculated at the patient level from out-of-fold predictions. For the metrics of the direct empirical baseline, 95% confidence intervals obtained by patient-level bootstrap estimation are additionally reported. Under conditions of a small sample size, the main emphasis is placed on comparing metric magnitudes across tasks and latent components, rather than on formal testing of differences among a large number of competing models.

3. Results

3.1. Inference of Patient-Specific Latent Model Parameters

At the first stage, the quality of the patient-specific latent representation, expressed as a reduced vector, was evaluated ${\tilde{\Theta }}_{pat}=({\theta }_{\inf },{\theta }_{\mathrm{exh}},{\theta }_{\mathrm{cyt}},{\theta }_{\mathrm{esc}},{\theta }_{\mathrm{env}},{\theta }_{\mathrm{init}})$.

The reference estimates of the latent parameters demonstrated pronounced inter-patient heterogeneity (Figure 4). The latent profiles differed substantially across patients and could not be reduced to a single averaged pattern, indicating the heterogeneity of the resulting latent space.

The quality of out-of-fold recovery of the latent parameters from the observed clinical features proved to be heterogeneous across the components of the latent vector (Figure 5 and Table 2). The best results were obtained for the parameter ${\theta }_{\mathrm{env}}$ ($\rho =0.858$, MAE $=0.021$, RMSE $=0.031$, $R^2=0.831$). Good agreement between the reference and predicted values was also observed for ${\theta }_{\mathrm{exh}}$ ($\rho =0.585$, $R^2=0.388$). For ${\theta }_{\inf }$ and ${\theta }_{\mathrm{cyt}}$, recovery quality was moderate. Despite positive rank correlations ($\rho =0.569$ and $0.563$, respectively), the coefficient of determination remained limited ($R^2=0.182$ and $0.230$). For ${\theta }_{\mathrm{esc}}$, predictive performance remained low ($R^2=-0.027$), whereas ${\theta }_{\mathrm{init}}$ did not show reliable agreement between the reference and predicted estimates ($\rho =-0.023$, $R^2=-0.357$).

The error distributions shown in Figure 5 further confirm differences in the reproducibility of individual latent components. The lowest MAE and RMSE values were observed for ${\theta }_{\mathrm{env}}$, whereas ${\theta }_{\mathrm{esc}}$ exhibited the highest median errors and the widest spread. For ${\theta }_{\mathrm{exh}}$, the error distributions remained comparatively narrow, while for ${\theta }_{\inf }$ and ${\theta }_{\mathrm{cyt}}$ they occupied an intermediate position.

The results show that recovery accuracy differs substantially across latent components. In the current configuration, the most reproducible parameters are those associated with the microenvironment and exhaustion, whereas reliable recovery was not achieved for the parameters related to antigen escape and, in particular, the individualised initial conditions.

3.2. Auxiliary Empirical Baseline for Binary Clinical Outcomes

We next consider an auxiliary empirical baseline that makes it possible to assess the extent to which the available observable features themselves contain signal for distinguishing binary clinical outcomes at the patient level. The analysis was performed for two endpoints – achievement of an objective response and disease control – and served as an external reference point for the subsequent interpretation of the contribution of mechanistic latent-parameter inference. Baseline performance was evaluated in an out-of-fold setting using balanced accuracy, $F_1$ score, AUROC, and AUPRC (Figure 6).

For the endpoint related to achievement of an objective response, classification performance remained limited: balanced accuracy was 0.62, the $F_1$ score was 0.63, AUROC was 0.57, and AUPRC was 0.60. For the disease-control endpoint, the metrics were higher: balanced accuracy was 0.64, the $F_1$ score was 0.74, AUROC was 0.66, and AUPRC was 0.87. Thus, disease control was predicted from the observed features substantially better than objective response.

Agreement between the observed and predicted cohort-level frequencies was additionally evaluated (Figure 7). For the endpoint of achieving an objective response, the observed proportion of patients with a positive outcome was 0.49, whereas the predicted frequency was 0.51. For the endpoint of achieving disease control, the corresponding values were 0.73 and 0.76, respectively.

Overall, the direct empirical approach showed that the observed clinical features carry only limited signal for individual prediction of the more stringent endpoint of objective response, whereas for disease control this signal is somewhat more robust. Thus, the auxiliary baseline should not be regarded as a standalone final result of the study, but rather as an external reference point against which the principal value of the proposed approach lies in translating clinical heterogeneity into a space of interpretable latent mechanisms.

3.3. Mechanistic Interpretation of Latent Factors Underlying Response and Resistance

After evaluating the quality of recovery of the patient-specific latent parameters and considering the auxiliary empirical baseline, the latent space was then interpreted in terms of the mechanisms underlying response and resistance. To this end, the values of the components of the reduced vector ${\tilde{\Theta }}_{pat}=({\theta }_{\inf },{\theta }_{\mathrm{exh}},{\theta }_{\mathrm{cyt}},{\theta }_{\mathrm{esc}},{\theta }_{\mathrm{env}},{\theta }_{\mathrm{init}})$ were compared across patients with positive and negative clinical outcomes. Effect size was assessed using Cliff’s delta separately for objective response and for disease control (Figure 8).

For objective response, positive effects were observed primarily for the parameters related to CAR-T cell cytotoxic potency and infiltration, whereas negative effects were found for the parameters associated with exhaustion, antigen escape, microenvironmental suppressiveness, and initial conditions. Accordingly, the probability of response increases when preserved effector function is combined with effective tumour infiltration by CAR-T cells, whereas resistance is associated with intensified immunosuppression, loss of the antigenic target, and dysfunction of the effector pool.

For disease control, the direction of the effects was preserved, although the contrast between the groups was less pronounced. The parameters most consistently associated with lack of disease control were those related to antigen escape and unfavourable initial conditions, whereas the contributions of infiltration and cytotoxic potency remained positive. This is consistent with the notion that achieving disease control requires a less stringent combination of favourable mechanisms than does the development of an objective response.

For a two-dimensional interpretation of the latent space, the landscape of the predicted response probability was examined in the $({\theta }_{\inf },{\theta }_{\mathrm{env}})$ plane, reflecting the balance between effective CAR-T cell infiltration and microenvironmental suppressiveness (Figure 9).

The probability of response is determined not by any single parameter, but by the patient’s position within the latent space (Figure 9). More favourable regions correspond to the combination of high infiltration and a less pronounced suppressive microenvironment, whereas displacement towards the opposite region is accompanied by a lower probability of response. The separation between responders and non-responders is continuous rather than discrete, indicating the existence of intermediate states between sensitivity and resistance.

Thus, the latent space of the model reflects mechanistically interpretable differences among patients. Response is associated with higher CAR-T cell infiltration and cytotoxic potency, whereas resistance is associated with exhaustion, antigen escape, and a suppressive microenvironment.

4. Discussion

The results obtained indicate that hybrid personalisation of a spatiotemporal model of CAR-T therapy in solid tumours is feasible, although the informativeness of clinical features is not uniform across different latent mechanisms. The parameters most robustly recovered were those associated with the immunosuppressive microenvironment and CAR-T cell exhaustion, whereas antigen escape and individualised initial conditions remained substantially less identifiable. The integral effects of the microenvironment and dysfunction are reflected across several observable channels, including CAR-T cell kinetics, inflammatory markers, and toxicity, whereas spatial antigen heterogeneity and the initial intratumoural architecture manifest themselves only indirectly in aggregated clinical data [6,7,11,21,26,47]. Accordingly, differences in recovery accuracy reflect not only algorithmic limitations, but also the fundamental incompleteness of the observations.

The latent space in the proposed model does not impose a rigid division of patients into responders and non-responders, but rather defines a continuous landscape of sensitivity and resistance. This indicates that response to CAR-T therapy in solid tumours is determined by the combined effects of infiltration, preservation of effector function, antigenic stability, and the level of local suppression, rather than by a single dominant mechanism [10,12,17]. Within this framework, a higher probability of response corresponds to the combination of effective infiltration and a less pronounced suppressive microenvironment, whereas resistance is associated with intensified exhaustion, antigen escape, and unfavourable initial conditions. This interpretation is consistent with the literature, in which inadequate CAR-T cell trafficking, chemokine mismatch, hypoxia, and inhibitory signalling are regarded as key barriers to therapy in solid tumours [16,18,19,47].

The lower predictability of objective response than of disease control is unsurprising. Disease control may be achieved under conditions of partial suppression of tumour growth, whereas objective response requires a deeper degree of cytoreduction and, consequently, more precise information on the spatial distribution of the target, the local density of infiltration, and the selection of antigen-negative clones. This is consistent both with the clinical data on CLDN18.2-specific CAR-T therapy [28,46] and with current understanding of the role of antigen escape and functional T-cell dysregulation in the development of resistance [13,14,23,25,48]. From the perspective of mathematical modelling, the results extend the line of work in which interpretable parameters are used to decompose the contributions of proliferation, exhaustion, spatial constraints, and antigen escape, while advancing it to the level of limited patient-specific adaptation [31,32,34,35,36,37].

From a theoretical perspective, the study supports the rationale for hybrid frameworks in which machine learning does not replace the mechanistic model, but is instead used to infer latent parameters from incomplete clinical data [40,42,43]. From a practical perspective, this makes it possible to regard latent parameters as a language for stratifying mechanisms of resistance: insufficient infiltration points to the potential value of strategies that enhance CAR-T cell trafficking and retention, dominance of microenvironmental suppression and exhaustion indicates the need for armoured constructs or combined immunomodulation, and a high risk of antigen escape highlights the relevance of multi-antigen platforms [16,17,19,47,48]. In a broader context, this approach is consistent with the development of patient-specific digital twins in oncology and digital medicine [49,50,51].

The present study should be interpreted in the context of an early stage in the development of a hybrid patient-specific model of CAR-T therapy for solid tumours. The cohort size and the predominant use of a single clinically coherent data source limit the statistical generalisability of the results; however, at this stage, such a design made it possible to minimise inter-study heterogeneity and to assess the internal consistency of the proposed personalisation pipeline. Importantly, the latent patient-specific parameters in the problem considered here are not directly observable; therefore, the study used surrogate mechanistic target variables obtained within the framework of a limited inverse problem, and these should be regarded as an operational approximation rather than as independent ground-truth labels. In addition, spatial intratumoural processes are reflected only indirectly in the available clinical data. Consequently, the results presented here primarily support the feasibility and mechanistic interpretability of the hybrid approach, rather than its readiness for clinical application. In this sense, the study constitutes a hypothesis-generating proof-of-concept framework that provides a foundation for subsequent external validation in independent cohorts and for expansion of the set of multimodal observations [27,28,49,50,51].

5. Conclusions

The aim of the present study was to develop and evaluate a hybrid framework for patient-specific calibration of a spatiotemporal model of CAR-T therapy in solid tumours, in which a reaction-diffusion-chemotaxis mechanistic core is complemented by machine learning to infer latent parameters from clinical data.

The results show that the proposed hybrid framework makes it possible to partially infer biologically interpretable latent parameters of a spatiotemporal model of CAR-T therapy from a limited set of clinical observations and to use them for the mechanistic interpretation of response and resistance scenarios. The parameters most reliably identified were those associated with the immunosuppressive microenvironment and CAR-T cell exhaustion, whereas antigen escape and individualised initial conditions remained substantially less observable. The auxiliary direct baseline model based on clinical features showed that, in the absence of mechanistic structure, the more stringent clinical endpoint exhibits only limited predictability. Thus, the current version of the approach demonstrates the feasibility and interpretability of limited patient-specific personalisation, but requires further independent external validation, expansion of the spectrum of tumour targets, and integration of more spatially informative data sources.