Can the Kuznetsov Model Replicate and Predict Cancer Growth in Humans?

Several mathematical models to predict tumor growth over time have been developed in the last decades. A central aspect of such models is the interaction of tumor cells with immune effector cells. The Kuznetsov model (Kuznetsov et al. in Bull Math Biol 56(2):295–321, 1994) is the most prominent of these models and has been used as a basis for many other related models and theoretical studies. However, none of these models have been validated with large-scale real-world data of human patients treated with cancer immunotherapy. In addition, parameter estimation of these models remains a major bottleneck on the way to model-based and data-driven medical treatment. In this study, we quantitatively fit Kuznetsov’s model to a large dataset of 1472 patients, of which 210 patients have more than six data points, by estimating the model parameters of each patient individually. We also conduct a global practical identifiability analysis for the estimated parameters. We thus demonstrate that several combinations of parameter values could lead to accurate data fitting. This opens the potential for global parameter estimation of the model, in which the values of all or some parameters are fixed for all patients. Furthermore, by omitting the last two or three data points, we show that the model can be extrapolated and predict future tumor dynamics. This paves the way for a more clinically relevant application of mathematical tumor modeling, in which the treatment strategy could be adjusted in advance according to the model’s future predictions. Supplementary Information The online version contains supplementary material available at 10.1007/s11538-022-01075-7.


Introduction
Cancer immunotherapy with immune checkpoint inhibitors has revolutionized the treatment of patients with solid tumors in the last ten years. In addition to chemotherapy and molecularly-targeted therapy, immunotherapy provides a new set of tools for the oncology toolkit (Wheeler et al. 2004). In several tumor types such as melanoma, non-small cell lung cancer (NSCLC), and genito-urinary tumors, immunotherapy has markedly improved the average life expectancy of patients with advanced disease. Both laboratory and clinical experiments have verified the importance of the immune system in fighting cancer (de Pillis et al. 2006;Farrar et al. 1999;O'Byrne et al. 2000). Patients who suffer from acquired immunodeficiency syndrome (AIDS) are very susceptible to having some rare forms of cancer (de Pillis et al. 2006;Dalgleish and O'Byrne 2002). This also shows the significant role the immune system plays against cancer.
One of the fundamental problems in treating patients with cancer immunotherapy is the lack of predictive biomarkers. Ideally, before the treatment begins, patients could be selected for immunotherapy, but existing biomarkers fail to deliver a high predictive value in most tumor types (Chatterjee and Zetter 2005). In addition, most patients who initially respond to immunotherapy experience a relapse: the tumor later on develops immune escape mechanisms due to evolutionary pressure. Forecasting the time of relapse or treatment resistance is of high practical relevance (Anderson and Quaranta 2008;Rockne and Scott 2019). However, predictions of such changes in the tumor behavior are currently not possible in clinical routine. The main problem is that most biomarkers such as tumor mutational burden (TMB) are static, i.e., they are measured at a given time point but are not dynamically updated as the tumor evolves.
In other complex systems such as financial markets (Ledoit et al. 2001), climate systems (Manabe 1983) or complex industrial processes (Thompson and Kramer 1994), differential equation models can provide a prediction of the behavior of the system over time. By analogy, in oncology, a number of mathematical models to predict tumor growth over time have been developed in the last decades (Norton et al. 1976). Most notably, multiple of these models explicitly include the interactions of tumors with the immune system and are therefore in principle suited to model response and resistance to cancer immunotherapy (Kogan et al. 2012). de Boer and Hogeweg (1986) modeled the cellular immune reaction to tumors. They demonstrated that small doses of antigens lead to tumor dormancy (de Pillis et al. 2006). Kirschner and Panetta (1998) linked the oscillations in the tumor size and the long-term tumor regression to the dynamics among immune cells, tumor cells, and Interleukin-2 (de Pillis et al. 2006). The most prominent of these models was presented by Kuznetsov et al. (1994). Kuznetsov's model has served as a blueprint for many other related models (de Pillis et al. 2006;Rhodes and Hillen 2019;Makhlouf et al. 2020) and has been investigated in several theoretical studies (Bellomo and Preziosi 2000;Kolev 2003;de Pillis et al. 2006;Owens and Bozic 2021).
However, none of these established oncological models are currently being used in the clinic. What is more, very few of these models have been systematically fitted to actual clinical data. While some studies have fitted models to murine tumors on a small scale (de Pillis et al. 2006;Benzekry et al. 2014;Vaghi et al. 2020), the pronounced differences between mice and humans preclude the transfer of such insights to realworld cancer patients (Ruggeri et al. 2014).
The structures of these mathematical models are well defined (Kuznetsov et al. 1994;Kolch et al. 2015;Tyson et al. 2011;Fröhlich et al. 2018). However, in the complex biological environment of cells, little is known about the associated parameters and kinetic constants. The parameter values are essential for quantitative modeling and prediction of cancer progression. In mechanistic models, one can integrate the data from various experimental procedures and sources, and design in silico experiments to generate hypotheses for underlying mechanisms (Clegg and Gabhann 2015;Baker et al. 2018). By fitting the model to the experimental data, we reverse-engineer the parameters of the system. Parameter estimation of mathematical cancer models remains a major bottleneck on the way to model-based and data-driven medical treatment of the future.
In this study, we use a mathematical model based on Kuznetsov's model to characterize the interactions between the growing tumor and the immune system, and aim to fill this conceptual gap in the literature. We use a large dataset of thousands of cancer patients who underwent cancer immunotherapy as part of clinical trials. We then investigate how well the model can represent the actual tumor volume changes over time in these patients. After estimating the parameters of the model, we conduct an identifiability analysis to examine the uniqueness of the estimated parameters (i.e., whether we have over-fitting). Finally, we investigate if the model can be used to forecast treatment response or relapse under immunotherapy.
The remainder of the article is structured as follows. First, we briefly present the acquisition of patients data and its pre-processing in Sect. 2. In Sect. 3, we then provide the model and all the methods used to estimate model parameters and fit the data, conduct parameter identifiability analysis and extrapolate the model for tumor growth prediction. Finally, before drawing conclusions in Sect. 5, we present and discuss the obtained results in Sect. 4.

Data Acquisition
We briefly provide here the declaration and sources of the experimental data. For more details, please refer to Ghaffari Laleh et al. (2022).

Declarations and Data Sharing
We followed the Declaration of Helsinki and International Ethical Guidelines for Biomedical Research Involving Human Subjects developed by the Council for International Organizations of Medical Sciences (CIOMS). In this study, we used a publicly available set of anonymized patient data shared by Ghaffari Laleh et al. (2022), which is originally derived from five large clinical trials, as we describe below. Patients gave their informed consent for data analyses as part of the original clinical trials. No specific ethical approval was sought or required for this retrospective analysis of a publicly available dataset.

Data Measurement and Pre-processing
The original data for this study is from five clinical trials which were designed to evaluate the efficiency of Atezolizumab (an immune checkpoint inhibitor). Table 1 shows the original number of patients and their treatment arms. Four out of these five trials evaluated the effect of Atezolizumab on NSCLC and cohort GO29293 reported this efficiency on bladder cancer. In two of the cohorts (GO28753, GO28915), patients' responses to Atezolizumab treatment were compared to the outcome of the second treatment arm who received Docetaxel (a chemotherapy drug) as treatment. In all the trials, the longest diameter and shortest diameter (LD and SD) of the target and non-target lesions (measured manually based on the CT scans) alongside the time intervals are reported. Several patients have only one or two data points because of tumor progression and potentially the death of patients. In this study, we use the anonymized and publicly available subset of data from Ghaffari Laleh et al. (2022), which was created by selecting the patients with three or more measurement points. For each of these patients, only the LD measurement for one target lesion has been selected. For this reason, the total number of original patients of clinical trials has been decreased from 2693 to 1472. Considering patients that have very few data points will lead to an over-parameterized problem. The version of the Kuznetsov model we used, which we discuss in Sect. 3.1, has six parameters. The model in this case would easily fit any combination of data points with many (essentially arbitrary) values for the parameters. On the other hand, disregarding patients with few data points (patients with tumor progression) introduces potential bias to the estimation problem and might limit the model application in a clinical setting. All the pre-processing details for the data generation have been described in more details in Ghaffari Laleh et al. (2022). Moreover, before using the data we pre-process it by first removing repetitive and null inputs. Following Ghaffari Laleh et al. (2022), Faustino-Rocha et al. (2013 and Shevtsov et al. (2019), we converted the measured LD in mm to the number of tumor cells (TC) by considering 8 × 10 −6 mm 3 /TC, where we consider spherical shapes of the lesions and that the TC have 3/4 of the lesion volumes. Moreover, we consider only patients with more than six net measurements because we are conducting model extrapolation in Sect. 3.4 in which the last two or three data points are omitted when estimating the model parameters. Thus, we preserve a minimum number of data points  for parameter estimation also in this case (at least three). Figure 1 shows the number of patients per study and arm before and after data pre-processing. Finally, the total number of patients has been reduced from 1472 to 210.

Methods
In this section, we introduce the mathematical model of Kuznetsov et al. (1994) with slight modifications. Afterward, we show the optimization formulations for: the estimation of the model parameters by fitting the clinical data and for the parameter identifiability analysis. We finally investigate the extrapolation capabilities of the model.

Mathematical Tumor Model
To predict the frequently observed phenomena in clinics like tumor dormancy and tumor size oscillation, the tumor mathematical model has to include terms related to the response of the immune system. The inclusion of the entire immune system in the mathematical model can be very difficult (Perelson and Weisbuch 1997). The antitumor immune response has highly nonlinear dynamics which are complicated and not well understood. Therefore, models that describe the immune system response to tumor presence should necessarily focus on those elements of the immune system that have the highest effects on tumor dynamics (de Pillis et al. 2006). Kuznetsov's model (Kuznetsov et al. 1994) describes the response of the cytotoxic T lymphocyte (CTL) to the growth of an immunogenic tumor. Usually, a cell-mediated immune response to a tumor takes place. The cytotoxic T lymphocytes and natural killer (NK) cells play the main role. The model considers immunogenic TC that are attacked by cytotoxic effector cells (EC). The EC can be, for example, CTL or NK cells. The model takes into account the possibility of EC inactivation as well as the infiltration of the TC by EC. TC and EC interaction is described through the following reactions: where EC-TC denotes conjugates of effector and tumor cells and EC * and TC * are the inactivated effector and lethally-hit tumor cells, respectively. We define E, T , C, E * , T * as the number of EC, TC, EC-TC conjugates, EC * , and TC * , respectively. The non-negative kinetic parameters, k 1 , k −1 , k 2 and k 3 , describe the rates of the interactions. EC-TC conjugates can reversibly decompose without damaging the cells with the kinetic rate k −1 . However, they can also irreversibly result in EC * or TC * with respective kinetic rates k 2 and k 3 . The following system of nonlinear differentialalgebraic equations describes those interactions, which is a slightly simplified version of the model of Kuznetsov (Kuznetsov et al. 1994).
The rate of flow of mature EC to TC localization area is characterized by the generation term s. This rate is unaffected by the presence of TC. The destruction or migration of EC from the localization region of TC is represented by the elimination rate h. The model does not take into account any TC or EC-TC conjugates migration. Both multiplication and death of TC are included in parameter a that characterizes the maximum growth rate of TC population. The function F(C, T ) represents the accumulation rate of the cytotoxic EC in the TC localization area due to tumor existence (stimulated accumulation), where f and g are positive constants. The EC accumulation, F(C, T ), is due to signals, like released cytokines, generated by the EC in EC-TC conjugates. Thus this stimulated accumulation has some maximum value when T becomes large. The equations describing the rate of change of E * and T * are not included in the system because they are irreversibly formed and thus have no effect on the other variables, and our target is to model T and E only. In Kuznetsov et al. (1994), the model includes a sink term in the rate of change equation of T that represents TC growth limitation due to biological environment conditions. It considers, for example, resources competition like oxygen and substrates. We do not consider this term because growth limitations of even high initial T are associated with high rates of cytotoxic EC accumulation as well as the absence of their activity suppression by TC (de Boer and Boerlijst 1994).
Following the suggestion of Kuznetsov et al. (1994), we consider a quasi-steadystate assumption for (1c), i.e., dC /dt t ≈ 0, because C is formed and dissociated at much faster rates compared to the multiplication and influx of the EC, as well as the lysis of the lethally-hit TC. Thus, C ≈ K ET , where K = k l /(k 2 + k 3 − k −1 ). As a result, (1a) and (1b) become: For further analysis and use of the model in parameter estimation and identifiability analysis, we use the same strategy of Kuznetsov et al. (1994) for non-dimensionalizing model equations. We non-dimensionalize (2a) and (2b) by considering concentration scales E 0 = 10 7 cells and T 0 = 10 9 cells for EC and TC, respectively (Kuznetsov et al. 1994). We non-dimensionalize t by relating it to the deactivation rate of TC and introducing τ = k 2 K T 0 t/100. The final model formulation is: where Thus, the final model is composed of two ordinary differential equations (ODEs) with two variables, x and y, with their respective initial values, x 1 and y 1 , at the initial normalized time, τ 1 , and six unknown parameters, σ , ρ, η, μ, δ, and α.

Data Fitting and Parameter Estimation
To determine the parameter values for the nonlinear system (3a) and (3b) that best describe the experimental data, we conduct a regression analysis in the nonlinear least-squares sense by minimizing the sum of the squared residuals. The considered residuals are the differences between the measured values of tumor lesion longest diameters (converted to tumor number of cells as previously discussed) and the ones calculated from the model. In the present contribution, we do not seek global values of parameters, i.e., the same values for all patients. Instead, we solve the optimization problem for each patient individually to identify the parameter values that best describe the data of that patient. This constitutes a first step to check whether the model can describe the experimental data at all and analyze the ranges of parameter values. We also estimate the initial normalized value of EC number, x 1 , because it is unknown. In contrast, the initial normalized value of TC number, y 1 , is provided experimentally and does not need to be estimated. The nonlinear least square problem for each patient j ∈ J = {1, 2, . . . , J } is expressed as follows: subject to model equations (3a) and (3b), where J is the total number of patients considered. After data pre-processing, J = 210 patients. The index i ∈ I j is for the observed experimental values, where I j = {1, 2, . . . , N j } with N j being the total number of observed values for patient j.
The non-dimensionalized model-predicted and observed values of TC number at the normalized time τ i j are y i j (τ i j , p j , x 1 j ) andỹ i j , respectively. The initial value of the non-dimensionalized EC number at τ 1 j , x 1 j , has lower and upper bounds x l 1 j and x u 1 j , respectively. The vector p j contains the non-dimensionalized parameter values of the model equations, (3a) and (3b), for patient j, p j = σ j , μ j , δ j , α j , ρ j , η j . The lower and upper bounds of the components of p j are the components of vectors p l j and p u j , respectively. The decision variables of the optimization problem are thus the components of p j and x 1 j . The dynamic optimization problem (OP1) is nonconvex and nonlinear, and can thus have multiple (suboptimal) local solutions. Hence, global optimization techniques are required to guarantee the global optimal solution, p opt j and x opt 1 j .

Parameter Identifiability Analysis
Parameter identifiability analysis determines if model parameters can be uniquely estimated (Walter and Pronzato 1997). Different definitions of identifiability analysis are available in the literature. Miao et al. (2011) reviewed several methods of parameter identifiability analysis for nonlinear ODE models and distinguished between different methodologies including structural and practical identifiability analyses. In the former analysis, one determines if a given structure of a model allows the realization of unique parameters when certain measured variables are provided (Walter and Pronzato 1997). However, it only provides necessary conditions for identifiability because it does not take into consideration parameters precision (Jung et al. 2019;Raue et al. 2011). In contrast, practical identifiability aims to predict confidence intervals for the estimated parameters (Jung et al. 2019;Gábor et al. 2017). It can be conducted locally (in the neighborhood of the estimated parameter values) or globally over the entire range of values. We here carry out the latter analysis and evaluate it globally to improve the confidence in the parameter values that are determined by solving (OP1).
We conduct the global practical identifiability analysis by determining the smallest box that contains the so-called feasible parameter set P e j for each patient j, as suggested in Jung et al. (2019). This set includes all values of parameters in p j for which the differences between model predictions y i j (τ i j , p j , x opt 1 j ) and optimal model predictions y opt i j (determined by solving (OP1)) fall within certain defined bounds ∀i ∈ I j , that is where y , ε is the percentage of deviation, and P j is the set of parameter values in p j bounded by the components of p l j and p u j defined in (OP1). The set P e j is depicted in dark gray as shown in Fig. 2 for the case of a two-dimensional vector p j . We approximate the nonconvex set P e j by a rectangular box (light gray color), whose edges are formed by the extreme values of the elements of p j ( p min j 1 , p max j 1 , p min j 2 , p max j 2 ). We determine these extreme values by solving a series of constrained dynamic optimization problems. For a vector p j that consists of K parameters, the optimization problem for parameter number k is formulated ∀ j ∈ J, and ∀k ∈ K as follows (Jung et al. 2019;Paulen et al. 2016): (3a) and (3b), with K = {1, 2, . . . , K } and K = 6. We set ε to 20%, which means that the differences between model predictions for p j = p max j and p j = p min j , and optimal model predictions for p j = p opt j fall within 20% of the optimal prediction values ∀i ∈ I j . For each parameter p j k in p j , its minimum value p min k j in p min j and its maximum value p max j k in p max j are determined by solving (OP2). Therefore, (OP2) is solved K times for each patient j ∈ J. As a result, the approximation of P e j is determined, and its edges are the elements of p min j and p max j . The set P e j allows for the determination of confidence regions of the estimated parameter values in p opt j . When P e j covers a large space in the direction of p j k , the estimated parameter p opt j k is not identifiable, in the sense that there is a large range of values for p j k that could lead to a good fit. When the set covers a small space, p opt j k is identifiable because the parameter is determined to a sufficient accuracy. To decide on identifiability a certain threshold is thus needed. We here do not define a cutoff, we rather analyze the identifiability qualitatively.
If estimates of the errors in tumor length measurements were available, then depending on the confidence intervals of the measurements, the threshold value of 20% in the practical identifiability analysis could be changed accordingly. For instance, for wide error margins in the measurements, the threshold value should be increased to account for those wide margins. However, we do not have estimates of measurement errors. The longest and shortest diameters were measured manually by radiologists in the original clinical trial. Such manual measurement is currently the state of the art, although our previous experience indicates that there can be a 1-2 mm margin of error in these measurements (McNitt-Gray et al. 2015).

Tumor Growth Prediction
Beyond being able to reproduce experimental data a posteriori, a more clinically relevant application of mathematical tumor modeling would be if the model was able to predict tumor growth. This could lead to model-based tumor treatment, as the supplied doses to patients could be adjusted in advance according to model predictions. In Sect. 3.2, we fitted the model to all experimental data points when estimating the parameters. In order to compare model extrapolation capabilities and future predictions to the clinical data, we now do not include the last ζ points of the data when fitting the model and solving (OP1). For each patient j, we solve an optimization problem The optimal values of the decision variables obtained from this problem are called p ext j and x ext 1 j . We then integrate the model equations for p j = p ext j and x 1 j = x ext 1 j from τ 1 j to τ (N j ) j , thus extrapolating beyond the data used for fitting. The results of the integrated y ext i j = y i j (τ i j , p ext j , x ext 1 j ) between τ (N j −ζ ) j and τ (N j ) j are the extrapolated part of the model that can be compared with the remaining ζ data points to gauge the extrapolation capabilities. We consider two model extrapolation cases in which ζ is equal to two and three.
Moreover, we formulate another optimization problem to investigate how far model extrapolation could deviate from the actual values. We aim to find two "extreme-case" lines that are designed to be as far away from each other at the final time point (τ (N j ) j ) while both being within some θ tolerance of the found optimal fit (y ext i j ) for the fitted time before extrapolation starts. For this, we solve the following optimization problem ∀ j ∈ J: model equations (3a) and (3b), where p upp j and p low j are the vectors that contain the parameter values for the upper and lower "extreme-case" model extrapolation deviations, respectively. The initial values of the non-dimensionalized EC number at τ 1 j for these upper and lower deviation cases are x upp 1 j and x low 1 j , respectively. We set θ to 10%. By integrating (3a) and (3b) from τ 1 j to τ (N j ) j , we get y , which allow the comparison of these two "extreme-case" extrapolation deviations with the remaining ζ data points after τ (N j −ζ ) j (start of extrapolation). If estimates of the measurement errors were available, then the deviations of the "extreme case" predictions from the experimental measurements could be assessed while taking into consideration those error estimates. One could then conclude if the "extreme case" predictions still lie within the measurement uncertainty or not. However, as aforementioned, we do not have estimates of the measurement errors.

Implementation
We implement the model, (3a) and (3b), in MATLAB R2019b (MATLAB 2019). All optimization problems are solved in the MATLAB version of the global optimization toolbox MEIGO using the enhanced scatter search metaheuristic (eSS) method (Egea et al. 2014). The eSS is stochastic and employs some elements of the scatter search and path re-linking methodologies (Egea et al. 2010). We set the maximum number of function evaluations, the maximum CPU time and the maximum absolute violation of the constraints to 10 5 , 100 s and 10 −5 , respectively. We use a 50% probability of biasing the search toward bounds and the dynamic hill climbing (DHC) (Yuret and de La Maza 1993) as a local search method. For all aforementioned optimization problems, we set x l 1 j and x u 1 j to 10 −2 and 10 2 , respectively. All elements of p l j and p u j are set to 10 −2 and 10 2 , respectively. Because of the non-dimensionalization described in Sect. 3.1, we expect that the parameter values are close to one. Therefore, we arbitrarily chose those bounds. Wider bounds might result in lower optimal objective values and hence better fits, but the computational cost would increase too. We also tried for several patients to use wider bounds (e.g., 10 −3 and 10 3 ). However, it did not yield significant improvements for the majority of the cases.
The ODE (3a) and (3b) are solved using the variable-step, variable-order (VSVO) solver based on the numerical differentiation formulas (NDFs) of orders one to five (ode15s) (Shampine and Reichelt 1997). We set the relative and absolute error tolerances to 10 −3 and 10 −6 , respectively. The solution refinement factor is one, and the maximum step size is 0.1(τ (N j ) j − τ 1 j ).
Although MEIGO is a global solver, since it is stochastic, the solution depends on the initial guesses. Thus, the global optimum is not guaranteed. For all patients, we performed multiple optimization runs from different initial guesses and chose the decision variables that resulted in the lowest objective functions. In general, the improvements were very slight. However, for some patients whose profiles are difficult to predict (e.g., wavy profiles), the multi-start optimization did improve the optimal solution that was found.

Results and Discussion
We now show the results of estimation, identifiability, and predictions. We show the results of six selected patients in detail and give performance measures for all 210 patients. We selected those six patients in a way to provide the different profiles of TC dynamics. Data fitting and growth prediction results of all 210 patients are provided in the supplementary material.

Data Fitting and Parameter Estimation
We fitted the parameters of a modified Kuznetsov model on a dataset of solid tumors in human patients under immunotherapy treatment. The model predicted the different tumor growth profiles represented by a selection of representative patients, as well as in the total (unselected) cohort. As shown in Fig. 3 and Table 2, model prediction and experimental data profiles are qualitatively and quantitatively very close for these patients. The mean absolute error (MAE), the root-mean-square error (RMSE), and the coefficient of determination (R 2 ) for the selected six patients are given in Table 2. We found that the model gave very high goodness of fit as measured by R 2 . Across 210 patients in all studies, an average R 2 of 0.784 was achieved with median and range values of 0.896 and 1.594, respectively. According to the statistical Lilliefors normality test, this R 2 distribution does not follow a Normal distribution. Moreover, Table 4 provides the number of patients and the R 2 values per study and per arm. Study 1 and Study 5 have higher R 2 than the remaining studies. Arm 1 has the highest R 2 in all studies except for Study 2. Although direct comparison of this performance with the previous work in Ghaffari Laleh et al. (2022) is not possible, comparing the MAE of the selected six patients with the reported average MAE in the previous study indicates good fitting performance of the developed Kuznetsov model. Furthermore, by analyzing the goodness of fit in individual patients, we found that the modified Kuznetsov model was able to fit clinically interesting patterns. In particular, the modified Kuznetsov model was able to predict relapse after initial tumor response (patient #207 in Fig. 3) and other types of fluctuating behavior, solving a key limitation of previously used simpler models as in Ghaffari Laleh et al. (2022).
We fitted the parameters of the modified Kuznetsov model to a large clinical dataset obtained from five clinical trials. In general, the distributions of the resulting parameter values were similar between the studies (Fig. 4). These ranges can be useful for further studies because they enable other researchers to determine plausible ranges and optimal boundaries when fitting the same model to other datasets, thereby simplifying the optimization procedure. A global parameter estimation, however, was not performed in this study and could be attempted in future studies. Figure 4 shows the estimated values (black dashes) of model parameters of all 210 patients. All parameter values vary per patient. The values are distributed all over the bounds, except for α, which has a maximum of 6.331. The parameter α is the normalized parameter for a, which represents the maximum growth rate of TC population. In addition, most of the values of μ and ρ are close to the upper bound. Although parameter values are quite distributed between the bounds, the aforementioned findings can help in narrowing the expected ranges of values of parameters when global parameter estimation is targeted. Moreover, the distribution of parameter values is compared among the considered five studies. As we can see in Fig. 4, the distribution densities between the bounds of parameter values are the same for all studies.
We also provide box plots for the results of each parameter in Fig. 4. The medians of the parameter values in all studies and arms have close values. In addition, we performed statistical Lilliefors normality tests, for the parameter value distributions in each study and arm. The parameter σ follows Normal distributions in Study 1, Study 2, Study 5 and all the Arms of Study 3 except Arm 6. The parameter ρ follows Normal distributions in Study 1, Study 2, Study 3, and Arm 1 of Study 4. Figure 5 provides the identifiability analysis results of the estimated model parameter values of the selected patients. Depending on the patient and the parameter, the estimated parameter values can be unique or take other values. For patient #22, the approximated feasible parameter set covers a small space in all parameters directions. Their estimated parameter values are close to being unique and thus identifiable. On the other hand, the approximated feasible parameter set for patient #203 covers a large space in all parameter directions. Therefore, the found values are not unique and other combinations of values could lead to good fits and predictions. For the remaining selected patients, the feasible parameter sets can be small or large depending on the patient and the parameter.  Figure 6 shows the averages of the indicator functions for each parameter range (found by the identifiability analysis) for all patients. These averages provide the most frequent ranges of parameter values within their bounds. They are equivalent to histograms where the output is between zero and one (indicator function output). The ranges are quite distributed all over the parameter bounds for all parameters. The only exception is α, for which there is a maximum bound after which no values could be found. The most frequent ranges are close to the median values in Fig. 4, except for μ, for which it is at the lower bound of the parameter value.

Parameter Identifiability Analysis
When the target is to estimate parameter values of the tumor dynamics model of a certain patient, identifiable values are what we need. However, when the aim is to find global parameters values for all patients, large spaces of the feasible parameter sets can be desirable for finding those global values, in which they are independent of a considered patient. In summary, the results of Figs. 5 and 6 show that for several combinations of parameter values, good data fitting and model predictions can be achieved. This opens the potential for global parameter estimation, in which the estimated values of all or some model parameters are the same for all patients. However, using data from patients diagnosed with different cancer types (e.g., NSCLC and bladder cancer) might make it unlikely to find global values for all parameters. Nevertheless, we are not certain about that. Also, immunotherapy is tumor-agnostic and we wanted to pursue this tumor-agnostic approach. The applied version of the Kuznetsov model has six parameters. The model could describe different tumor growth profiles of both the NSCLC and bladder cancer when estimating the parameter values individually. It Global parameter estimation would have the advantage of fixing part or all of the Kuznetsov model parameters so that the estimation problem is simplified when applied to each patient in clinical practice. However, fixing all the model parameters and only estimating the initial values of the EC contradicts personalized modeling. In personalized modeling, patient-specific characteristics, treatment types, and metadata are included in the model as combined or additional terms. One could thus fix the parameters that are not specific to the patient or treatment types and estimate the others. The parameters specific to treatment types could be also fixed if only a certain Fig. 6 The results of the averages of the indicator functions for each parameter range (found by the identifiability analysis) for all patients. The graphs are equivalent to histograms indicating the frequency of the estimated ranges among the patients. Ordinates: indicator functions values, between zero and one (estimated ranges frequencies). Abscissas: non-dimensionalized model parameters. The most frequent ranges of parameter values are close to the median values in Fig. 4, except for μ, the range is most frequent at the lower bound of the parameter value. These median values are here the vertical lines (Color figure online) treatment type or strategy is applied. Parameters for metadata terms like gender, age, etc., could be also fixed within each group. Therefore, the approach of global parameter estimation should be applied to the parameters involved in the model terms describing common phenomena (which are tumor-agnostic) among tumor types, independent of treatment type or metadata, such as the terms in the Kuznetsov model.

Tumor Growth Prediction
In clinical decision making, a possible role of mechanistic models is forecasting tumor growth during treatment, potentially enabling physicians to adjust the treatment strategy earlier. We found that the modified Kuznetsov model indeed was able to extrapolate beyond the initial time points when the last two or three data points are not included when fitting the model. The solid black lines in Fig. 7 show model extrapolation results of the selected patients when the last two data points are not considered for fitting. For the six patients, the model quantitatively forecasts tumor dynamics, except for patient #186. The last data point of patient #186 is almost impossible to forecast because it suddenly shifts upward after a mild and constant increase in tumor growth. However, the model can still qualitatively predict the growth. For the other patients, the predictions are very close to the experimental data. In Table 3, the first sub-table "optimal Extra." provides the MAE, RMSE, and R 2 of the selected patients, as well as the mean absolute percentage error (MAPE) for the extrapolated part, defined as: The values of R 2 for the six patients are close to one except for patient #186 due to the aforementioned explanation. Compared to the previous work in Ghaffari Laleh et al. (2022), model extrapolation results here have higher R 2 values, specifically, an R 2 of 0.979 was reached for patient #207. The average values of R 2 and MAPE of all 210 patients are also provided in the table. The R 2 value is 0.419 with median and range values of 0.685 and 6.717, respectively. According to the statistical Lilliefors normality test, this R 2 distribution does not follow a Normal distribution. Model extrapolation results when omitting the last two data points for all 210 patients are provided in the supplementary material. By analyzing the performance of model extrapolation results of the 210 patients, we could not generally classify the tumor dynamic profiles into different subtypes. However, the model could not predict the future growth for several patients where tumor dynamics suddenly increase after steady profiles (e.g., patient #186). In addition, the performance of model extrapolation for patients receiving chemotherapy is better than for those receiving immunotherapy. Also, the model performed in general well when extrapolating curves that have slight slope changes after starting future predictions. In addition, we provide model extrapolation results while omitting the last three data points for the selected six patients in the supplementary material. For the latter results, the model could also predict the growth very well. Table 4 provides the R 2 and MAPE values per study and per arm for model extrapolation when omitting the last two data points. In general, and as expected, the calculated average R 2 for the extrapolation experiment is lower than the R 2 when using all the data points in all the studies. Similar to data fitting, Study 1 and Study 5 have higher R 2 than the remaining studies. In Study 3, Arm 2 and Arm 3 have the lowest R 2 . Interestingly, in Study 2 and Study 4, the extrapolation performance of the model is markedly better for the Docetaxel group (Arm 2) than for the Atezolizumab group (Arm 1). This indicates that patients receiving docetaxel have a more predictable response trajectory compared to patients receiving immunotherapy in whom unexpected patterns of tumor response can occur even at later time points. Moreover, in Study 3 where all patients received immunotherapy, Arm 1 (MPDL3280A-1a) and Arm 3 (MPDL3280A-1a) have better model performance in both fitting and extrapolation than the rest of the arms.
Moreover, Fig. 7 shows how far model extrapolation could deviate from the actual values. The results of the two "extreme-case" (dash-dotted blue and the dashed magenta) lines are designed to be as far away from each other at the final time point (τ (N j ) j ) while both being within some 10% tolerance of the found optimal fit for the fitted time before extrapolation starts. For some patients, the "extreme-case" lines can significantly deviate from the actual values, especially for the upper case (dashdotted blue lines), as for patients #22 and #83. In contrast, the "extreme-case" lines for patients #163 and #186 are very close to the optimal extrapolated ones. Sub-tables "Upper Extra." and "Upper Extra." in Table 3 provide the MAE, RMSE, R 2 and MAPE for the selected patients for the two "extreme-case" extrapolations.
To sum up, the model can forecast tumor dynamics of the patients and the "extremecase" extrapolation scenarios were conducted to check the worst model predictions. For some patients the "extreme-case" extrapolation results deviate from the actual values, for others they remain close to the optimal extrapolation results.

Conclusion
In this study, we show that quantitative mathematical models can be used to describe and forecast the behavior of cancer. Previous studies have used the same datasets to fit very simple ODE models to the tumor volume measurements over time (Ghaffari Laleh et al. 2022). However, in Ghaffari Laleh et al. (2022), it was observed that all established ODE models were not able to fit "U"-shaped trajectories well. In clinical terms, patients who relapsed after an initial response, or patients who showed a delayed response, were not adequately represented in these previous models. Compared to this, the present study evaluates a more complex model which has the benefit of being able to describe a larger variety of real-world time series. This specific model is a slight simplification of the Kuznetsov model (Kuznetsov et al. 1994), which has not been linked with or validated in large amounts of quantitative real-world human data. Specifically, it could quantitatively fit the Kuznetsov model to a large dataset of 1472 patients. Data are collected from patients undergoing immunotherapy or chemotherapy treatments.
In the parameter estimation for each patient, we found that some parameters for some patients are not unique (identifiability analysis). This means that many combinations of parameter values could lead to good fitting and predictions. This opens the potential for global parameter estimation, in which parameter values are the same for all patients. However, since we did not consider patients with few data points (patients with early tumor progression and potentially their death), potential bias to the estimation problem is introduced that might limit the clinical applicability of the model with  Fig. 7). The Upper Extra. sub-table relates values to the upper "extreme-case" extrapolation deviation results (the dash-dotted (blue) line in Fig. 7). The Lower Extra. sub-table relates values to the lower "extreme-case" extrapolation deviation results (the dashed (magenta) line in Fig. 7) global parameters. Still, the model could predict tumor growth (2-3 omitted measurements) which could indicate practical usefulness as a predictive biomarker. This is a more clinically relevant application of mathematical tumor modeling. Specifically, the model fitting and prediction could potentially describe and forecast the behavior of cancer, improve the understanding of underlying biological mechanisms, and provide model approaches for cancer treatments, as the supplied doses to patients could be adjusted in advance according to model predictions. Future studies should attempt a global parameter estimation, use deterministic and less computationally demanding solution methods. Another possibility is the reduction of the number of model parameters in order to have less number of parameters to estimate and thus enhance the possibility to reach global values. In addition, the mixed-effects modeling (population approach), which allows the simultaneous modeling of tumor dynamics and inter-individual variability within a statistical framework, could be considered as a potential work for future investigation.  The model complexity is mostly limited by the availability of suitable data. Therefore, if more data become available, it would be possible to include other terms in the model, e.g., metadata, explicitly representing different types of immune cell populations or multiple cancer cell clones.
A general limitation of our approach is that mechanistic models of tumor growth are competing against statistical models of tumor growth, including machine learning models. In order to forecast clinical response to immunotherapy and to guide clinical decisions, both types of models could be applied. However, it is still unclear which approach will ultimately translate to clinical application. While mechanistic models are attractive because they integrate expert knowledge into their predictions, they are also constrained by these assumptions. Our study contributes empirical evidence that mechanistic models can make clinically relevant predictions in a relevant use case, but further work is required to validate this in other disease contexts and compare the predictions to statistical models.