Analysing the Structural Identifiability and Observability of Mechanistic Models of Tumour Growth

TL;DR

This paper analyzes whether parameters and states in tumor growth models can be reliably inferred from data, using 20 published models and an open-source tool.

Contribution

The paper provides a systematic analysis of structural identifiability and observability for 20 tumor growth models, including those with interventions.

Findings

Many tumor growth models are not structurally identifiable or observable under standard measurement setups.

The identifiability and observability of models depend heavily on the choice of experimental design and measured variables.

The paper provides computational implementations to reproduce the analysis and guide model selection.

Abstract

Mechanistic cancer models can encapsulate beliefs about the main factors influencing tumour growth. In recent decades, many different types of dynamic models have been used for this purpose. The integration of a model’s differential equations yields a simulation of the behaviour of the system over time, thus enabling tumour progression to be predicted. A requisite for the reliability of these quantitative predictions is that the model is structurally identifiable and observable, i.e., that it is theoretically possible to infer the correct values of its parameters and state variables from time course data. In this paper, we show how to analyse these properties of tumour growth models using a well-established methodology, which we implemented previously in an open-source software tool. To this end, we provide an account of 20 published models described by ordinary differential equations, some of which incorporate the effect of interventions including chemotherapy, radiotherapy, and immunotherapy. For each model, we describe its equations and analyse their structural identifiability and observability, discussing how they are affected by the experimental design. We provide computational implementations of these models, which enable readily reproducing results. Our results inform about the possibility of inferring the parameters and state variables of a given model using a specific measurement setup, and, together with the corresponding methodology and implementation, they can be used as a blueprint for analysing other models not included here. Thus, this paper serves as a guide to select the most appropriate model for each application.

Figures

Tables

Table 1. Summary of the models considered in this paper.

Model Name	Acronym	Treatment	Reference
Exponential	EXP	None	[5]
Power Law	POW	None	[29]
Lotka–Volterra	L–V	None	[30]
Gompertz	GOM	None	[29]
Logistic	LOG	None	[29]
Von Bertalanffy	BERT	None	[29]
Radio Base	RAD	Radiotherapy	[31]
Carrying Capacity Radio	RCAP	Radiotherapy	[32]
Radio Necrotic	NECR	Radiotherapy	[33]
Immuno CAR-T cells CRS	CRS	Immunotherapy	[34]
Immuno-Hematological CAR-T cells	HCART	Immunotherapy	[35]
Immuno-Radio	IMRAD	Immunotherapy	[36]
CML Tumour-Immune Int	LEUK	Immunotherapy	[37]
Cancer–Immunity Cycle	CYCLE	Chemotherapy	[38]
Cytostatic and Cytotoxic Effects	CYTO	Chemotherapy	[39]
Cancer Immune Chemotherapy Vitamins	CICV	Chemotherapy	[40]

Full text

1. Introduction

In recent decades, ever-growing progress has been achieved regarding the understanding of cancer biology. The result is an increasingly complex picture, well illustrated by the enlargement of those traits considered to be hallmarks of cancer [1,2,3]. Likewise, the use of mathematical modelling—and, more specifically, dynamic modelling—has become increasingly common in the biomedical sciences. As a result, mathematical oncology has become an established area of mathematical biology, with its own set of hallmarks [4,5]. There are currently a wealth of computational tools that support the biomedical modelling process on the one hand and that enable fully exploiting the capabilities of mathematical models on the other [6]. The integration of mathematical and computational tools—sometimes known as the (dynamic) systems biology approach [7,8,9]—has facilitated the modelling of tumour dynamics and the exploration of the therapeutic efficiency of different treatments.

In this paper, we focus on mathematical models that are mechanistic: i.e., they provide explanations of the mechanism by which a phenomenon takes place; and dynamic: i.e., they are based on (differential) equations that describe the key aspects of their dynamics. These models can be used for many purposes [10]. First of all, they encapsulate the available knowledge (or beliefs) about a process, enabling results to be produced and to test if a hypothesis is consistent with the current knowledge. Furthermore, they can make quantitative predictions [11], ideally with an estimate of the associated uncertainties. This predictive capability provides many possibilities; crucially, it allows not only to anticipate the result of a treatment but also to optimise its design.

Thus, in the context of mathematical or computational oncology, these models can be used to understand complex phenomena such as the impact of hypoxia on tumour growth, the interactions between tumour and immune cells, and the effects of therapeutic regimens. There are currently a wide variety of tumour growth models of varying mathematical formalisms, purposes, and levels of detail (see, e.g., the reviews [5,12,13,14,15] and the references therein). In this context, Ref. [16] provided a conceptual review of tumour growth dynamics, showing that the growth patterns of the tumours can be grouped into a limited set of categories (exponential, surface, sigmoidal, atypical, and multistep), each reflecting the insights that drive tumour progression. Complementarily, Ref. [17] carried out a systematic comparison of intrinsic growth laws, fitting classical models such as exponential, logistic, Gompertz, von Bertalanffy, and power-law dynamics to experimental data from diverse cancer types. Most previous comparative studies focused either on classifying intrinsic tumour growth laws or fitting them to experimental datasets. Our study differs in scope as we systematically analyse the structural identifiability and observability of mechanistic ODE-based models, including those that integrate therapeutic intervention. One aim of this paper is to provide a compact account of some of the most relevant ones, facilitating appropriate model choice for a given purpose. To this end, we outline a set of mathematical models that describe tumour growth, with and without therapeutic interventions, including radiotherapy, immunotherapy, and chemotherapy. We consider models based on ordinary differential equations (ODEs), which can describe processes that are deterministic and homogeneous. For some of them, stochastic extensions have been published, which enable studying the influence of random fluctuations. We assign an acronym to each model, and we provide their equations in a systematic notation. Table 1 lists the 16 basic models considered in this paper; since we also discuss variants of some of them, the total number of analysed models is 20.

The application of these models goes beyond understanding the biological aspects of cancer; their ultimate goal would arguably be to guide the design and optimisation of treatments, maximising their efficacy while minimising adverse effects for the patient. The use of optimal control in this context has a long history that includes the design of treatments based on chemotherapy [18,19], immunotherapy [20], and radiotherapy [21]. In order to perform successfully in those applications, a model must be properly calibrated [6]. The key tasks of the calibration process include identifiability and observability analysis [22], parameter estimation [23], sensitivity analysis [24], uncertainty quantification [25], experimental design [26], and model selection [27].

In this work, we focus on the first of the tasks in the above list. These analyses assess whether it is theoretically possible to infer the model parameters (identifiability) and internal dynamics or state variables (observability) from measurable outputs. The relevance of identifiability in clinical applications was recently emphasised by the authors of [28], who illustrated how a lack of identifiability can undermine predictive accuracy even in carefully constructed cancer models. To help avoid these pitfalls, we perform a detailed analysis of the structural identifiability of the parameters of the collected models, and of the observability of their state variables. We provide all the files that enable reproducing the results in the models folder of the STRIKE-GOLDD toolbox, which we used to perform the analyses (https://github.com/afvillaverde/strike-goldd, accessed on 1 August 2025).

The remainder of this paper is organised as follows: Section 2 introduces the methodological framework, including the mathematical formulation of ODE models and the tools used for structural identifiability and observability analysis. Section 3 presents a catalogue of tumour growth models, classified according to the treatments that they include: no therapy, radiotherapy, immunotherapy, or chemotherapy. Section 4 reports and discusses the results of the analyses for each group of models and explains the implications for parameter estimation and model exploitation. Finally, Section 5 summarises the main conclusions and outlines directions for future research.

2. Methods

In this section, we present the mathematical framework used to describe tumour dynamics based on ordinary differential equations (ODEs). We then introduce the notions of structural identifiability and observability, which are key to determine whether model parameters and internal states can be inferred from measurable outputs. Finally, we describe the computational tools used to perform these analyses, focusing on the STRIKE-GOLDD Matlab toolbox 4.2.2.

2.1. Modelling Framework

Ordinary differential equations are the most widely used framework for modelling tumour growth. In their general form, they can be written as

$$[eqn]$$

where f and h are analytic functions, which are in general nonlinear; $[eqn]$ is the vector of state variables at time t; $[eqn]$ is the input vector; $[eqn]$ , the output vector; and $[eqn]$ is the parameter vector. In the following, we will sometimes omit the time dependence for ease of notation; i.e., we may simply write $[eqn]$ . A model’s outputs are the measurable functions, which are often a subset of the state variables, or a simple function of them (e.g., a sum). For the models in which we are interested, the inputs are the treatments or therapeutic actions.

2.2. Structural Identifiability and Observability

In this study, we followed a systematic approach to evaluate the structural properties of ODE-based models of the form (1). We focused on two fundamental aspects: structural identifiability and observability, which are essential to ensure the model’s reliability for parameter estimation and model prediction. Structural identifiability refers to the theoretical possibility of uniquely determining the model parameters from ideal noise-free input–output data. When a parameter is not identifiable, there are infinite different parameter sets that result in identical outputs, hindering model calibration and predictive use [22]. A related concept, observability, refers to the ability to infer the internal state variables of the system based on output measurements over time. We refer to structural identifiability and observability jointly with the acronym SIO [41]. Figure 1 illustrates these two concepts.

It should be noted that the SIO properties can be analysed before collecting any data (which is why they are sometimes called a priori properties) since they do not depend on the specific measurement values but only on the output definition. In other words, they are fully determined by Equation (1), and their analysis implicitly assumes that it is possible to measure the output $[eqn]$ at any possible instant $[eqn]$ and with perfect precision. Thus, SIO are necessary conditions for a model’s successful identification; sufficient conditions are provided by their practical counterparts [22,42,43]. Practical identifiability analysis requires numerical—not symbolic—methodologies, and it is highly application-dependent, so it is out of the scope of the present study.

Therefore, to assess the way in which the availability of measured outputs influences SIO, in our analyses, we considered in some cases various measurement scenarios, including individual and combined outputs. We explored different measurement configurations, such as isolated versus combined output observations, to understand their impact on the structural properties of the system. For example, if a model describes the evolution over time of the tumour volumes of two different cell populations—one susceptible to radiotherapy and the other one resistant to it—as in the example of Figure 1, the model’s SIO properties may be different if we measure the volume of the susceptible population, of the resistant population, of both of them separately, or of their sum. We would also like to note that in this paper we analyse local versions of structural identifiability and observability. Even though they are local properties, since they are structural they are also generally valid. That is, if a parameter is structurally locally identifiable (SLI), its true value can be distinguished from any other values in a neighbourhood around it, and this result holds for every possible value in parameter space. The only possible exceptions to this generic character, if they exist, belong to a set of measure zero; in other words, even if a parameter is SLI, it may become unidentifiable for a particular value, or at most for a finite number of values. Furthermore, empirical studies [44] suggest that the distinction between locally and globally identifiable is typically of little importance in practice since most SLI parameters are also globally identifiable, and for those that are not there is often no possibility of confusion.

2.3. Methodology for Analysing Structural Identifiability and Observability

There are several approaches to the analysis of the structural local identifiability and observability (SIO) of a nonlinear model. Here, we use the one based on differential geometry, which was originally proposed in [45], among other works. This approach includes the state variables and parameters in an augmented state vector, $[eqn]$ and derives an observability–identifiability matrix $[eqn]$ of the form [46]

$$[eqn]$$

The above matrix is built from Lie derivatives of the output. The first-order Lie derivative, $[eqn]$ is defined as

$$[eqn]$$

Higher-order Lie derivatives are recursively calculated as

$$[eqn]$$

Once the $[eqn]$ matrix has been obtained, we can determine the SIO of a model as follows: if $[eqn]$ , all the parameters are structurally locally identifiable, and all the state variables are observable. If $[eqn]$ , the model is not fully identifiable or observable. In this case, we determine the identifiability of each parameter and the observability of each state variable by removing the corresponding matrix column and recalculating the matrix rank, as described, e.g., in [41].

We performed these analyses using the STRIKE-GOLDD toolbox, version 4.2.2 [47], an open-source Matlab toolbox based on differential geometry methods for structural system analysis. This tool enables the examination of model properties independently of specific parameter values or experimental data and is freely available on GitHub. An alternative method is the Matlab toolbox StrucID [48], which is based on the calculation of sensitivities. It is usually faster than STRIKE-GOLDD and can hence be convenient for some analyses. Here, we used it for one model (IMRAD) that was particularly computationally expensive.

Some models (e.g., those that include radiotherapy) involve time-dependent treatment inputs such as impulses, which are often modelled with Dirac delta functions. We treated these impulses as constants during the analysis. This representation enables the study of identifiability and observability with the approach described above without compromising the mathematical consistency of the analysis.

3. A Catalogue of Models

In this section, we provide an overview of some representative ODE models of tumour growth found in the literature. We provide their equations, and we classify the variables of each model as inputs (which can be manipulated, and usually correspond to treatments), outputs (i.e., the quantities or functions that can be measured), states (whose variation over time is provided by the model’s differential equations), and parameters (constant quantities, which may be either known or unknown). This classification is a requirement for performing SIO analyses. Table 2 lists the parameters, state variables, inputs, and outputs of the ODE models considered in this paper, along with their acronyms and equations.

3.1. Tumour Growth Models Without Therapy

In the following subsection, we explore models that describe tumoural growth in the absence of any type of treatment.

3.1.1. Exponential Model (EXP)

The most basic growth model (EXP1) is the one provided by the following equation [5]:

$$[eqn]$$

Its solution, which is an exponential over time, is

$$[eqn]$$

where $[eqn]$ is the tumour volume over time, $[eqn]$ the net growth rate of the tumour, and $[eqn]$ the initial volume of the tumour.

While this model can typically accurately describe early growth stages, it fails to describe the reduced growth rates and eventual saturation due to the lack of nutrients and oxygen that vascular tumours develop in vivo, and also in avascular tumours cultured in vitro. To account for these aspects, a slightly more complex model can be used (EXP2):

$$[eqn]$$

Its solution is

$$[eqn]$$

with $[eqn]$ , and $[eqn]$ represents the population carrying capacity, which is the value approached by $[eqn]$ as time tends to infinity.

The model can be made more flexible by introducing an additional parameter $[eqn]$ , resulting in (EXP3) the following:

$$[eqn]$$

which leads to the following solution:

$$[eqn]$$

The possible output of these models is obviously the tumour volume, which is often measured over time to assess the growth dynamics and response to treatments. A common limitation of these models is the difficulty of relating their parameters to the behaviour of single cells.

3.1.2. Power Law Model (POW)

The Power Law model [29] adopts a different approach for describing how cells actively divide and grow in a lower spatial dimension than the full organoid. It is provided by the equation

$$[eqn]$$

with $[eqn]$ , which has the following explicit solution:

$$[eqn]$$

where $[eqn]$ is the number of cells of an organoid, and a is the rate at which the cells divide. This model is useful, for example, when the growth of the tumour is restricted to its surface. Its primary measurable variable—i.e., its possible output—is the population size, which can be interpreted as the tumour volume in certain biological contexts.

3.1.3. Lotka–Volterra Model (L–V)

The classical Lotka–Volterra model [30] can be used to describe the dynamics of several populations, e.g., of radio-sensitive and radio-resistant tumour cells. Without modelling the effect of radiotherapy explicitly, the governing equations of this model can be written as

$$[eqn]$$ $$[eqn]$$

where $[eqn]$ and $[eqn]$ represent the volumes of the control and resistant populations, $[eqn]$ and $[eqn]$ their growth rates, and $[eqn]$ and $[eqn]$ their carrying capacities. Additionally, $[eqn]$ and $[eqn]$ describe the effect that radiosensitive cells have on resistant cells and vice versa.

The possible outputs for this model are the tumour volumes of the two interacting populations, $[eqn]$ and $[eqn]$ . These volumes represent the measurable variables, capturing the growth dynamics of each tumour population over time. It may be possible to measure each of these volumes individually, although it is easier to obtain measurements of their sum, $[eqn]$ . To this end, one promising approach is the use of dual-reporter systems in preclinical models for in vivo systems. Here, the sensitive and resistant clones are genetically labelled with distinct fluorescent proteins, enabling the tracking of each subpopulation through intravital microscopy or fluorescence molecular tomography [49].

3.1.4. Gompertz Model (GOM)

This well-known model [29] assumes that an initial growth rate $[eqn]$ decreases exponentially within time, according to a decay parameter $[eqn]$ . The differential equation is (GOM1) as follows:

$$[eqn]$$

which, for $[eqn]$ , has the solution

$$[eqn]$$

There is an alternative formulation of this model (GOM2) that considers an initial growth rate $[eqn]$ and a carrying capacity K as follows:

$$[eqn]$$

Its solution is

$$[eqn]$$

The possible output for both versions of the Gompertz model is the population size, N.

3.1.5. Logistic Model (LOG)

A model similar to the previous one is the so-called logistic model, where the growth rate decays linearly with the size of the population until the population reaches the carrying capacity K; it is provided by the equation

$$[eqn]$$

where $[eqn]$ , and with solution

$$[eqn]$$

The main output is obviously the population size, N, which can be measured over time as the population grows and approaches its limit, K.

3.1.6. Von Bertalanffy Model (BERT)

The model proposed by von Bertalanffy [29] assumes growth as a three-dimensional ball; cell division happens only on the surface, and the deaths are uniformly distributed across the organoid at a rate $[eqn]$ . This model is provided by the equation

$$[eqn]$$

with $[eqn]$ . It has another more general version that replaces the exponent $[eqn]$ with a parameter $[eqn]$ . Thus, it considers that the subset of actively dividing cells occupies a lower spatial dimension. This variant has the equation

$$[eqn]$$

which, for $[eqn]$ , has the explicit solution

$$[eqn]$$

reaching the carrying capacity $[eqn]$ as $[eqn]$ . Reparameterising in terms of b, $[eqn]$ , and K, it can be rewritten as

$$[eqn]$$

3.2. Tumour Growth Models with Radiotherapy

In this subsection, we list models that describe the way in which radiotherapy affects tumour growth.

3.2.1. Radio Base Model (RAD)

This model was presented in [31] and is provided by (RAD1)

$$[eqn]$$

with solution for $[eqn]$ :

$$[eqn]$$

where $[eqn]$ , n is the number of radiotherapy fractions, d the dose administered at times $[eqn]$ ( $[eqn]$ ), and $[eqn]$ denotes the Dirac delta function, defined such that, for an arbitrary smooth function $[eqn]$ ,

$$[eqn]$$

$[eqn]$ and $[eqn]$ represent the tumour volumes immediately before and after radiotherapy, and the survival fraction (SF) is defined as $[eqn]$ , so

$$[eqn]$$

In order to analyse the model with the approach described in Section 2, the input must be infinitely differentiable. Hence, in our analyses, we replace the delta functions with a constant input that quantifies the average effect of a radiotherapy dose over a short interval. Note that, while this is an approximation, it is arguably closer to reality than a Dirac delta function itself, which cannot be physically realised in practice.

Cell death resulting from radiotherapy can be integrated with the logistic and generalised logistic growth models in the following way (RAD2):

$$[eqn]$$

with $[eqn]$ , and $[eqn]$ for the logistic growth law. The possible output for both models is the tumour volume V, which can be measured over time as the tumour grows and responds to the doses of radiation.

3.2.2. Carrying Capacity Radio Model (RCAP)

Other authors [32] have proposed a model in which the doses of radiotherapy instantaneously reduce the carrying capacity of the tumour in the logistic model instead of affecting the tumour volume:

$$[eqn]$$

where $[eqn]$ describes how effective a given radiation dose is at killing the proliferative portion of the tumour. It quantifies the effect of radiotherapy on the tumour’s carrying capacity, rather than on its volume, because radiation damage affects both tumour vasculature and its microenvironment. In this way, after each dose, the maximum tumour size is effectively being reduced. The possible output for this model is the tumour volume V, which can be measured over time as the tumour grows and responds to radiotherapy.

3.2.3. Radio Necrotic Model (NECR)

This ODE model [33] distinguishes between viable tumour volume, $[eqn]$ , and dead or necrotic tissue, $[eqn]$ , and it involves the equations

$$[eqn]$$ $$[eqn]$$

where $[eqn]$ is the constant rate under which cells undergo necrosis, and $[eqn]$ is the rate at which the necrotic material undergoes exponential decay. In addition, $[eqn]$ is the death rate of tumour volume V. The possible outputs for this model are the viable tumour volume $[eqn]$ and the necrotic material $[eqn]$ .

3.3. Tumour Growth Models with Immunotherapy

Here, we describe models of immunotherapy, a family of treatments based on the idea of administering drugs to help the immune system fight cancer.

3.3.1. Immuno CAR-T Cells Cytokine Release Syndrome Model (CRS)

In [34], a multi-layer mathematical model that describes the dynamics of CAR-T cell therapy and the associated cytokine release syndrome (CRS) was presented.

The first variable represents the CAR-T cells that were administered to the patient (injected CAR-T cells, $[eqn]$ ) to combat cancer. Their dynamics is provided by

$$[eqn]$$

where $[eqn]$ represents the antigen–receptor binding interaction that is described by

$$[eqn]$$

The second compartment includes the CAR-T cells that have been activated and are in the process of proliferating (expanding CAR-T cells, $[eqn]$ ). These cells multiply in response to the presence of tumour cells that express the specific antigen that they recognise. Their equation is

$$[eqn]$$

The third compartment represents CAR-T cells that survived the expansion phase and differentiated into memory cells (persistent CAR-T cells, $[eqn]$ ). They are important for the long-term immune response as they can be reactivated if needed. Their equation is as follows:

$$[eqn]$$

The fourth compartment represents the tumour cells that express the specific antigen that the CAR-T cells are designed to recognise and attack (antigen-positive tumour cells, $[eqn]$ ); this variable describes how CAR-T cells affect the tumour burden:

$$[eqn]$$

Likewise, the dynamics of antigen-negative tumour cells ( $[eqn]$ ) is provided by

$$[eqn]$$

In this case, $[eqn]$ is a constant fraction that represents a reduced cytotoxicity due to antigen absence.

The model parameters are A, the half-saturation constant; $[eqn]$ , $[eqn]$ , and $[eqn]$ , the death rates of injected, expanding, and persistent CAR-T cells, respectively; $[eqn]$ , the tumour growth rate; and K, the tumour cell carrying capacity.

An additional layer describes cytokine dynamics and cytokine release syndrome (CRS). It assumes that cytokine release is mediated by macrophage activation, driven by three main mechanisms: the release of inflammatory mediators by activated CAR-T cells, the release of DAMPs (Damage-Associated Molecular Patterns, which are molecules released by stressed or dying cells that trigger immune activation) from tumour cell death, and contact-dependent activation between macrophages and CAR-T cells via CD40/CD40–ligand interaction. Under these assumptions, the model for macrophage activation and cytokine release reads as

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

where $[eqn]$ and $[eqn]$ are the number of naive (monocytes) and activated macrophages at time t, $[eqn]$ is the $[eqn]$ concentration, $[eqn]$ is the natural production rate of naive macrophages, $[eqn]$ is the death rate of macrophages, $[eqn]$ the endogenous production rate of $[eqn]$ , $[eqn]$ the $[eqn]$ natural decay rate, and $[eqn]$ the rate of $[eqn]$ release by activated macrophages, whose activation rate is provided by

$$[eqn]$$

where $[eqn]$ is the contribution of the antigen-binding-release of inflammatory signals by CAR-T cells, $[eqn]$ is the contribution of the tumour-killing-mediated release of DAMPs, and $[eqn]$ the contribution of CAR-T cell–macrophage contact.

Some of the eight state variables in the model are easier to measure than others. We have considered the set $[eqn]$ , $[eqn]$ , $[eqn]$ , and $[eqn]$ as reasonable outputs; we justify this decision as follows.

Regarding the expanding CAR-T cells ( $[eqn]$ ), the measurement of the expansion of these cells is essential to check the immediate effectiveness of the treatment and to determine how well the immune system is adapting to fight the tumour; this can be conducted through blood tests using flow cytometry, a relatively common procedure in immunotherapy. However, this process still requires specialised equipment and technical personnel.

One important variable in the model is the concentration of the cytokine $[eqn]$ , which plays a major role in the immune and inflammatory responses. $[eqn]$ is commonly measured to track immune activity in the context of CAR-T therapy, where excessive immune activation can lead to CRS [34]. $[eqn]$ levels are relatively easy to measure through blood tests (ELISA). This makes $[eqn]$ one of the more straightforward and low-effort variables to monitor while providing critical insight into the immune response and potential toxicities associated with treatment.

Another output could be the activated macrophages ( $[eqn]$ ) and monocytes ( $[eqn]$ ), whose activity helps to assess the broader immune response beyond the CAR-T cells themselves. These types of cells can be measured through blood tests and flow cytometry, but identifying and distinguishing between them requires additional markers, so the process would be more complex.

Lastly, tumour cells with antigen-positive ( $[eqn]$ ) and antigen-negative ( $[eqn]$ ) expression can indicate whether resistance to the treatment is developing. Unfortunately, measuring these variables often involves more invasive techniques, such as biopsies or advanced imaging, to detect changes in tumour antigen expression. Hence, even though it could be interesting to know these levels, it may require too much effort.

3.3.2. Immuno-Hematological CAR-T Cell Model (HCART)

This model [35] describes the interaction between tumour cells and CAR-T cells (effector and memory populations) in immunodeficient mouse models. The model incorporates tumour growth, CAR-T cell dynamics, and the formation of immunological memory. It uses three ODEs to explore the factors influencing therapy outcomes, including cytotoxicity, tumour-induced immunosuppression, and long-term memory formation.

Effector CAR-T cells ( $[eqn]$ ) proliferate at a certain rate ( $[eqn]$ ) in response to tumour burden, die or differentiate into memory cells ( $[eqn]$ ), are activated by interaction with tumour cells ( $[eqn]$ ), and are suppressed by tumour-induced immunosuppressive mechanisms ( $[eqn]$ ):

$$[eqn]$$

Memory CAR-T cells ( $[eqn]$ ) are generated from effector CAR-T cells ( $[eqn]$ ), activated back into effector cells when interacting with tumour cells ( $[eqn]$ ), and decay naturally at a rate $[eqn]$ :

$$[eqn]$$

Tumour cells (T) grow logistically with a maximum growth rate r and carrying capacity $[eqn]$ . They are eliminated by effector CAR-T cells through a cytotoxic effect ( $[eqn]$ ):

$$[eqn]$$

3.3.3. Immuno-Radio Model (IMRAD)

This model [36] describes the tumour response to radio-immunotherapy, i.e., the combination of radiotherapy (RT) with monoclonal antibodies such as $[eqn]$ PDL1 and $[eqn]$ CTLA4. This model describes the dynamics of tumour cell populations, T-cells, and antibody release.

The equations of the dynamics of tumour cells (C), damaged tumour cells ( $[eqn]$ ), T-cells ( $[eqn]$ ), and antibody release (A) are, respectively,

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

In addition, these three equations are needed to fully define the model:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

3.3.4. CML Tumour-Immune Interaction Model (LEUK)

This model [37] analyses the interactions between Chronic Myeloid Leukemia (CML) cells and immune effector cells, particularly in predicting outcomes after tyrosine kinase inhibitor (TKI) treatment cessation. The model integrates CML cell dynamics with immune responses, incorporating parameters for immune cell functionality and tumour load.

The main equations describing the model are as follows.

Growth of active tumour cells, T:

$$[eqn]$$

Quiescent tumour cells, $[eqn]$ :

$$[eqn]$$

Immune effector cells, E:

$$[eqn]$$

Functional tumour-immune interaction as a function of tumour load:

$$[eqn]$$

Anti-tumour activity of immune cells:

$$[eqn]$$

Suppression of immune function by tumour burden:

$$[eqn]$$

Tumour-dependent modulation of immune recruitment:

$$[eqn]$$

Immune effector cells as a function of tumour load:

$$[eqn]$$

The model parameters are $[eqn]$ , the net proliferation rate of active tumour cells; $[eqn]$ , the carrying capacity of tumour cells; $[eqn]$ , the rate at which active tumour cells enter quiescence; $[eqn]$ , the rate at which quiescent tumour cells re-enter the active state; $[eqn]$ , the TKI-induced cell kill rate (dependent on treatment timing); $[eqn]$ , the immune cell activity per tumour cell; $[eqn]$ , the suppression factor of immune function by tumour load; $[eqn]$ , the production rate of immune effector cells; $[eqn]$ , the death rate of immune effector cells; and $[eqn]$ , the tumour-mediated suppression of immune cell recruitment.

3.4. Tumour Growth Models with Chemotherapy

3.4.1. Cancer–Immunity Cycle Model (CYCLE)

A recent model [38] considers the synchronisation of chemotherapy and immunotherapy to enhance treatment efficacy within the cancer–immunity cycle. It adopts an evolutionary game theory approach where healthy cells, cancer cells, and T-cells compete in a non-transitive dynamic, similar to a rock–paper–scissors game, modulated by the tumour microenvironment (TME). The aim is to explore how synchronising dosing intervals with the fundamental period of the cancer–immunity cycle ( $[eqn]$ ) can improve treatment effectiveness, allowing for lower doses with reduced toxic effects. The main equations describing the model are as follows:

Replicator equations for healthy cell ( $[eqn]$ ) and cancer cell ( $[eqn]$ ) fractions:

$$[eqn]$$ $$[eqn]$$

Growth of T-cells (n):

$$[eqn]$$

where

$$[eqn]$$

Fitness function for each subpopulation (healthy cells and cancer cells):

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

Tumour volume growth ( $[eqn]$ ):

$$[eqn]$$

Instantaneous growth rate ( $[eqn]$ ):

$$[eqn]$$

Interpolated payoff matrix:

$$[eqn]$$

where

$$[eqn]$$

The variables in these equations are as follows: $[eqn]$ and $[eqn]$ are the fractions of healthy and cancer cells in the total population, respectively. The variable n represents the T-cell response level, which indicates immune activity. The chemotherapy control variable, $[eqn]$ , ranges from 0 (no chemotherapy) to 1 (maximum chemotherapy), while the immunotherapy control variable $[eqn]$ takes values between 2 (no immunotherapy) and higher values to indicate active immunotherapy. The timescale parameter $[eqn]$ determines the relative speed of the immune response compared to cancer cell dynamics, and $[eqn]$ represents the natural immunity’s limitation in counteracting tumour growth entirely. The parameter $[eqn]$ is the instantaneous growth rate of the cancer cell fraction, and K represents the carrying capacity or upper limit of tumour volume. The matrices $[eqn]$ and $[eqn]$ are the payoff matrices for tumour growth and regression, respectively, while $[eqn]$ is an interpolation function that modulates the transition between tumour growth and regression based on immune response level.

3.4.2. Tumour Evolution with Cytostatic and Cytotoxic Effects (CYTO)

This recent model presented by [39] includes the dynamics of tumour volume $[eqn]$ , which are governed by two key parameters: a and b.

The parameter $[eqn]$ represents the cytostatic effect; it inhibits tumour growth without necessarily causing its reduction, so it acts like a mechanism of slowing down the proliferation of the tumour cells.

The parameter $[eqn]$ , on the other hand, corresponds to the cytotoxic effect, modelling the active destruction of tumour cells. When $[eqn]$ , the model accounts for tumour regression, typically due to the action of chemotherapy or the effect of the immune system. The equation of the model is

$$[eqn]$$

Both parameters (a and b) interact with a time-scaling function $[eqn]$ , whose derivatives $[eqn]$ and $[eqn]$ control how fast the tumour grows or shrinks over time. Since $[eqn]$ is assumed to be a known and smooth function, for the purpose of structural identifiability analysis, we treated it as a known input. The interplay between a, b, and the shape of $[eqn]$ determines whether the tumour volume $[eqn]$ grows, stabilises, or decreases over time.

To assess the SIO of this model, it is necessary to rewrite it as a system of first-order derivatives, resulting in the following two equations (CYTO2):

$$[eqn]$$ $$[eqn]$$

3.4.3. Cancer Growth Model with Chemotherapy and Boosting of the Immune System (CICV)

The CICV model presented by [40] describes the interaction between tumour cells $[eqn]$ and immune cells $[eqn]$ , incorporating the effects of chemotherapy and vitamin intake. It is defined as the following system of ODEs:

$$[eqn]$$ $$[eqn]$$

The parameter $[eqn]$ is the tumour’s growth rate, $[eqn]$ defines the inverse of its carrying capacity, $[eqn]$ is the immune killing rate, and $[eqn]$ is a Michaelis–Menten-type saturation constant. Chemotherapy is represented by $[eqn]$ , the rate at which it kills tumour cells directly. The parameter $[eqn]$ represents the baseline production rate of immune cells, $[eqn]$ quantifies the suppression that tumour cells conduct on the immune system, and $[eqn]$ models the stimulation of immune cells due to tumour presence, following a Holling type III response, which is modulated by the constant $[eqn]$ . The parameter $[eqn]$ is the death of immune cells due to natural and treatment causes. Finally, $[eqn]$ represents the rate at which regular vitamin intake enhances immune cell production.

4. Results and Discussion

In this section, we report the results of the analyses of structural identifiability and observability of the models described in the previous section. The results are summarised in Table 3, Table 4, Table 5 and Table 6. Each model is analysed at least once, and possibly more (for multiple output configurations, indicated in the columns ‘Measured Outputs’) if there are several feasible options.

Table 3 shows the results for models without treatment. All of them are fully identifiable and observable, although it is important to note that, in the case of the Lotka–Volterra model, it is necessary to measure the sum of both volumes (or, naturally, both of them separately, if possible) to achieve identifiability.

Table 4 shows the results for the radiotherapy models. Here, the most remarkable outcome is the realisation that neither the RAD1 nor RAD2 variants are identifiable. The remaining models do not exhibit identifiability issues.

Table 5 shows the results of the immunotherapy models. These models are typically complex, incorporating multiple immune components and a larger number of parameters. This feature is prone to induce unidentifiability issues due to overparameterisation. In the CRS model, only 13 out of the 23 parameters are identifiable, and none of the state variables are observable. HCART has slightly better structural properties, with T being observable and measured, but with $[eqn]$ and $[eqn]$ unobservable and several parameters remaining unidentifiable. IMRAD is unobservable, with several parameters (p, $[eqn]$ , b, $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , and $[eqn]$ ) and the initial conditions $[eqn]$ being non-identifiable from outputs C and A. To improve identifiability, we would need to measure more system outputs or fix some parameter values. In the LEUK model, only three out of ten parameters ( $[eqn]$ , $[eqn]$ , and $[eqn]$ ) are identifiable, even if all its states (T, $[eqn]$ , and E) are directly measured. It should be noted that, for liquid cancers such as leukemia, the underlying biology often involves nonlocal interactions between spatially separated compartments, for instance between bone marrow niches and circulating blood cells. In these cases, nonlocal or delay-based models have been proposed to capture the distributed nature of hematopoietic regulation [50,51].

Lastly, Table 6 summarises the results of the chemotherapy models. In the CYCLE model, key system parameters such as the payoff matrices ( $[eqn]$ and $[eqn]$ ), the interpolation factor a, the timescale parameter $[eqn]$ , and the carrying capacity K are identifiable, while $[eqn]$ and $[eqn]$ remain unidentifiable; all the state variables are observable if V is the measured output. The CYTO model (or, more precisely, its first-order version, CYTO2) is fully identifiable and observable under the assumed measurements. The identifiability of the CICV depends on the measured variable: with the output C, all the parameters are identifiable; however, when I is the measured output, parameters like $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , and $[eqn]$ become unidentifiable. It can be seen that the chemotherapy models show structurally better identifiability properties compared to the immunotherapy models.

Overall, the comparative analysis of Table 3, Table 4, Table 5 and Table 6 reveals clear differences among the model families. The classical growth models without treatments (Table 3) are fully identifiable and observable under natural measurement assumptions. In contrast, the radiotherapy models (Table 4) display more heterogeneity: while some formulations remain identifiable, the models RAD1 and RAD2 are not. The immunotherapy models (Table 5) show some structural limitations since many parameters are unidentifiable and several states are unobservable (especially CRS and IMRAD), which might be due to their increased complexity. Finally, the chemotherapy models (Table 6) tend to present favourable structural features, with most parameters and states identifiable from feasibly measurable outputs. These results suggest that models of increasing biological realism, especially in immunotherapy, are also those that pose more challenges for their identification.

From a practical point of view, unidentifiability means that there are infinite parameter sets that can reproduce the same measured outputs. As a result, it is impossible to determine the correct parameter estimates from experimental data with any degree of certainty; the results of a parameter estimation procedure are very likely to be wrong. In turn, due to the strong relation between (un)identifiability and (un)observability, it may also be impossible to infer the value of certain state variables, and their predictions may be unreliable. This limits the use of the model for therapy optimisation since treatment simulations could vary depending on the assumed parameter values.

5. Conclusions and Future Directions

In this paper, we have assembled a set of 20 dynamic models of tumour growth found in the literature. We have focused on mechanistic models that describe the physical and biochemical details of the processes involved in tumour growth using ODEs. Naturally, it would be unrealistic to attempt to include every model ever presented; instead, we have selected a number of them that include a wide array of mechanisms. Other models not included in this report can be analysed using the approach described here. Our selection includes some ODE models without therapeutic interventions, and others that describe the effect of treatments including radiotherapy, immunotherapy, and chemotherapy.

Our comparative results highlight clear differences between families of models: while classical growth laws and most chemotherapy models show favourable structural properties in regard to identifiability and observability, several immunotherapy models suffer from severe unidentifiability, at least for certain experimental configurations. These findings emphasise the importance of assessing structural properties before model calibration.

Overall, the abundance of models and the variety of types naturally lead to the question of how to choose the most appropriate one for a particular purpose. The decision should take into account, among other factors, the level of detail required and the availability of the data that is necessary to constrain the model. To facilitate this choice, for each model, we have discussed the measurements that are typically available, or feasible, and how they relate to the model variables. Sometimes, there may be more than one possible model that is compatible with the requirements; in this case, one may use a model selection technique to determine the one that provides the best balance between complexity and ability to fit the data [27].

Another important aspect to consider is whether the model allows us to estimate its unknown components reliably. This depends on two properties: identifiability, which refers to the ability to determine the model parameters from the available data, and observability, which refers to the ability to infer the internal variables from the measurements. If a model is not identifiable or observable, its predictions may be uncertain or misleading. Therefore, these properties should be analysed before applying a model, especially when it is used for designing treatments or making patient-specific predictions. In this paper, we have contributed to this endeavour by analysing the structural identifiability and observability (SIO) of the models. Since the SIO properties depend on the output functions, i.e., on the variables that are available for measurement, we have considered those possible combinations of measurements that are experimentally feasible. For some models, more than one output configuration is possible; in those cases, we have performed several analyses for each model. Our results, summarised in Table 3, Table 4, Table 5 and Table 6, shed light on our theoretical ability to obtain mechanistic insights from such models and can be used as a guide for choosing an appropriate model for a particular application.

The use of mechanistic models in oncology has been increasingly common since the 1990s. In recent years, these models have been instrumental in the development of personalised medicine. It is foreseeable that their use will only become more widespread in the near future as the key component of two emerging concepts: virtual patients and digital twins [52]. A virtual patient is basically a mathematical model calibrated with data obtained from a patient cohort. Thanks to this alignment with measurements, the simulations of a properly built virtual patient are expected to yield physiologically plausible outcomes. Hence, they can be used for simulating clinical trials for the purpose of drug development or therapy optimisation. The concept of digital twins [53] resembles that of virtual patients but with two key differences. First, a digital twin is typically built with data obtained from a single patient and is therefore a model of an individual, not of a target population. Second, digital twins are regularly fed with new data, thus enabling their use for cancer monitoring on a patient-specific basis. Digital twins are thus an evolution of virtual patient cohorts that are increasingly becoming an important part of precision medicine, enabling personalised predictions and in silico testing of therapeutic strategies [54]. For example, in a classic paper, Swiernak et al. [55] described virtual patient cohorts calibrated with clinical data to optimise chemotherapy protocols in breast cancer, demonstrating their potential to simulate clinical trials in silico. More recently, Ref. [56] discussed oncology digital twins that integrate longitudinal patient data to provide individualised treatment guidance and continuous monitoring.

Whereas mechanistic models are grounded in biological and physical understanding, recent research points to their strategic integration with data-driven approaches as a promising direction for the future. Unlike mechanistic models, data-driven models are not based on first principles and cannot be used for understanding a phenomenon, but they are capable of making predictions; however, the results lack explainability or interpretability. The emerging concept of mechanistic learning [57], which encompasses various ways of combining mechanistic models with machine and deep learning, offers a path to overcome some of the limitations that are inherent to each approach separately. Such hybrid methods promise to combine the interpretability of mechanistic models with the flexibility and data assimilation capabilities of learning-based algorithms. The incorporation of artificial intelligence (AI) techniques and big data to mechanistic models was recently discussed in [15]. Data-driven techniques can assist in estimating model parameters from high-dimensional data, correct residual errors, or serve as surrogates for complex models, thereby reducing computational costs and enabling real-time applications. The limitations of AI include the need for large datasets for training, which can typically be obtained only for patient cohorts and hence may fail to capture individual patient features. Crucially, AI predictions are difficult to generalise outside the conditions in which the training data was generated (overfitting to the training set). On the other hand, AI also offers important strengths for oncology modelling. Machine learning methods can detect complex patterns in high-dimensional clinical or omics data [58] and complement, as mentioned, mechanistic models in many different ways, such as predicting treatment responses [59]. Combined with mechanistic frameworks, AI can reduce computational costs by acting as surrogate models, enable real-time updating of digital twins, and facilitate the integration of heterogeneous data sources.

Moving forward, the field will benefit from a systematic exploration of different integration strategies, each suited to specific modelling goals, data availability, and clinical contexts. As data complexity and computational costs continue to grow, the combination of mechanistic and data-driven modelling has great potential to improve the predictive accuracy, adaptability, and relevance of tumour growth models.