# Modeling and Parameter Estimation of Tumor Spheroids from Bottom-View Imaging: Beyond the Spherical Assumption

Tiphaine Delaunay, Faiza Laanani, Léa Pechtimaldjian, Sandra Oucherif, Christèle Etchegaray, François Moisan, Annabelle Collin

► To cite this version:

Tiphaine Delaunay, Faiza Laanani, Léa Pechtimaldjian, Sandra Oucherif, Christèle Etchegaray, et al. Modeling and Parameter Estimation of Tumor Spheroids from Bottom-View Imaging: Beyond the Spherical Assumption. 2026. [hal-05449367](https://hal.science/hal-05449367v1)

HAL Id: hal-05449367

https://hal.science/hal-05449367v1

Preprint submitted on 8 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Modeling and Parameter Estimation of Tumor
Spheroids from Bottom-View Imaging: Beyond the
Spherical Assumption

Tiphaine Delaunay¹, Faiza Laanani², Léa Pechtimaldjian², Sandra Oucherif², Christèle Etchegaray³, François Moisan², and Annabelle Collin⁴

¹LAGA UMR CNRS 7539, Université Sorbonne Paris Nord, France

²Univ. Bordeaux, Inserm, Bordeaux Institute of Oncology (BRIC), UMR 1312, F-33076 Bordeaux, France

³Univ. Bordeaux, CNRS, Inria, Bordeaux INP, IMB, UMR 5251, 33400 Talence, France

⁴Laboratoire de Mathématiques Jean Leray, Nantes Université, Nantes, F-44322, France

Abstract

Tumor spheroids are widely used *in vitro* models for studying tumor growth and treatment response, yet their analysis is often limited by incomplete geometric measurements. In practice, only a two-dimensional bottom view is available, which can lead to strong biases when spheroids deviate from an ideal spherical shape, in particular for flattened aggregates. We propose a mechanistic modeling and data assimilation framework that accounts for both poor data quality and model uncertainty. A biologically motivated partial differential equations model with axial symmetry is introduced and reduced to a low-dimensional system compatible with sparse observations and ellipsoidal shapes. Uncertainties in initial conditions and parameters are handled using a Luenberger-type observer coupled with a reduced-order unscented Kalman filter. The approach is validated on synthetic data and applied to experimental *in vitro* spheroid data under propranolol treatment, illustrating how model-based integration of limited measurements improves the interpretation of spheroid growth dynamics and treatment-induced morphological changes.

# 1 Introduction

Mathematical modeling using ordinary and partial differential equations has become a central tool for understanding tumor growth across both clinical and biological contexts. These approaches are increasingly used not only to interpret *in*
*vitro* experiments, but also to support clinical decision-making and treatment design as illustrated by selected studies [1, 11, 22, 3, 9, 17, 4]; this list is by no means exhaustive.

In this work, we focus specifically on the mathematical modeling of tumor spheroids, a widely used *in vitro* system meant to emulate the three-dimensional architecture of solid tumors [7, 21]. However, despite their apparent geometric simplicity, spheroids often deviate substantially from idealized spherical shapes. Furthermore, a major experimental limitation is that only their bottom view – essentially a 2D projection – is observed, while their height, internal structure, and true geometry are severely under-measured. Even though the height can either be approximated through crude laws based on the observed 2D intensity, or inferred by sacrificing a few spheroids from the cohort (hereby allowing only very limited temporal information to be collected) it is unreasonable to assume that the quality of this measurement is of the same order as that of the bottom-view radius measurement.

This incompleteness becomes critical for certain cell lines, such as spheroids generated from skin-tumor cells, which tend to form flattened, pancake-like structures rather than spheres [16]. In such cases, approximating the spheroid volume under a spherical assumption using only the bottom-view area may introduce significant modeling biases. In particular, this may lead to the erroneous conclusion that a drug accelerates tumor growth, whereas the observed increase in the bottom area may in fact result from spheroid flattening rather than from an actual increase in volume [6].

Our objective is to develop a mathematical approach that integrates mechanical modeling with parameter estimation within a framework that accounts for both model error and poor data quality. We introduce a simple yet biologically expressive model able to reproduce key features of tumor spheroids, including the coexistence of proliferative and quiescent cell populations, an axial symmetry, and distinct dynamical behaviors along the vertical axis compared to the bottom view radius. We design an Luenberger type observer [10] capable of compensating for both uncertainties in initial conditions and model error in a setting where data are limited. This observer is then coupled with a reduced-order unscented Kalman filter [8, 13] to estimate model parameters robustly. Finally, we apply the methodology to real experimental data in the context of propranolol treatment for angiosarcomas [23], illustrating how model-based integration of sparse measurements can refine biological interpretation.

Section 2 introduces the PDE model and its order reduction, which is primarily based on axial symmetry. Section 3 presents the data assimilation framework, illustrated using synthetic data. Finally, Section 4 reports the results obtained on real experimental data.

# 2 Modeling

## 2.1 General model

By denoting the tumor spheroid domain at time $t$ by $\Omega(t)$ assumed to be smooth, the evolution of the proliferative density $P$ and the quiescent density $Q$ are supposed
to satisfy the following advection-diffusion system for $t > 0$ and $\vec{x} \in \Omega(t)$,

$$ \left\{ \begin{aligned} \partial_t P(t, \vec{x}) + \nabla \cdot (\vec{v}(t, \vec{x}) P(t, \vec{x})) &= \tau^G(t)(P(t, \vec{x}) + Q(t, \vec{x})) - \tau^{P2Q}(t, \vec{x})P, \\ \partial_t Q(t, \vec{x}) + \nabla \cdot (\vec{v}(t, \vec{x}) Q(t, \vec{x})) &= \tau^{P2Q}(t, \vec{x})P(t, \vec{x}), \end{aligned} \right. \quad (1) $$

where $\tau^G$ represents the growth of tumor cells – proliferating and quiescent – and $\tau^{P2Q}$ the rate of transition from proliferating to quiescent cells. Using the saturation hypothesis

$$ P + Q = 1, \text{ in } \Omega(t), \quad (2) $$

we have

$$ \nabla \cdot \vec{v}(t, \vec{x}) = \tau^G(t)(P(t, \vec{x}) + Q(t, \vec{x})) = \tau^G(t), \quad \Omega(t). \quad (3) $$

One can show using Reynolds theorem that the tumoral volume defined as

$$ V(t) = \int_{\Omega(t)} d\vec{x}. \quad (4) $$

satisfies $V'(t) = \tau^G(t)V(t)$. In the case where we assume that $\tau^G(t) = ae^{-bt}$, the volume follows a Gompertz model with parameters $(a, b, V_0) \in (\mathbb{R}^{+*})^3$ with $V_0 = V(0)$ i.e. $V(t) = V_0 e^{\frac{a}{b}(1-e^{-bt})}$, for more details see for example [3].

Drawing from the work of [12] and assuming that the spheroid is fully symmetrical, this model can be rewritten in radial and relative coordinates. Interestingly, this implies that the velocity field $\vec{v}$ is entirely determined by its divergence, therefore no assumption on the rheology of the tumor spheroid is needed. However, in our context where we assume that the spheroids are not symmetrical, a Darcy's law – involving a pressure denoted by $\pi$ corresponding to the pressure that tumor cells exert on the surrounding medium – must be considered to close the system. To integrate the fact that we would like to be able to consider more complex shapes than spherical ones, we define the pressure $\pi$ by:

$$ \pi(t, \vec{x}) := \frac{1}{2} \vec{x}^T M(t) \vec{x}, \quad (5) $$

the matrix $M$ being defined by

$$ M(t) = \begin{pmatrix} M_{x_1}(t) & 0 & 0 \\ 0 & M_{x_2}(t) & 0 \\ 0 & 0 & M_{x_3}(t) \end{pmatrix} \quad \text{such that} \quad \left\{ \begin{array}{l} M_{\alpha} \in \mathbb{C}^0([0, T]), \alpha = x_1, x_2, x_3, \\ \displaystyle\sum_{\alpha=x_1, x_2, x_3} k_{\alpha} M_{\alpha}(t) = -\tau^G(t), \end{array} \right. \quad (6) $$

where the coefficients $M_{\alpha}$, $\alpha = x_1, x_2, x_3$ will be chosen later. Then, Darcy's law writes

$$ \vec{v}(t, \vec{x}) = -K \nabla \pi(t, \vec{x}) = -KM(t)\vec{x}, \quad (7) $$

where $K$ is the constraint tensor defined by

$$ K = \begin{pmatrix} k_{x_1} & 0 & 0 \\ 0 & k_{x_2} & 0 \\ 0 & 0 & k_{x_3} \end{pmatrix} \quad (8) $$

with $k_{x_1}, k_{x_2}, k_{x_3} \in \mathbb{R}$. The tensor $K$ describes the deformability of the tumor over time.
**Remark 1.** The conditions given in (6) ensure the compatibility of the Darcy law (7) with the divergence constraint (3). More precisely, at t fixed, if M(t) satisfies (6), then the pressure π given by (5) is the unique solution to the following Poisson problem

$$ \begin{cases} -\nabla \cdot (K \nabla \pi)(t, \vec{x}) = \tau^G(t), & \Omega(t), \\ \pi(t, \vec{x}) = (\frac{1}{2} \vec{x}^T M(t) \vec{x})|_{\partial \Omega(t)}, & \partial \Omega(t) \end{cases} \quad (9) $$

Indeed, from the definition of π (5) and using (6), we have

$$ \nabla \cdot (K \nabla \pi)(t, \vec{x}) = \sum_{\alpha=x_1,x_2,x_3} k_{\alpha} M_{\alpha}(t) = -\tau^G(t), $$

which corresponds with the first equation of (9). Trivially, π defined by (5) is the solution to the second equation of (9). By uniqueness of the Poisson problem with non-homogeneous Dirichlet condition, π is the unique solution of (9).

**Remark 2.** No explicit boundary conditions on P and Q are prescribed here. For purely advective dynamics on a closed and smooth domain Ω(t), the system is fully determined by the divergence condition and initial data. If diffusion terms were added, appropriate boundary conditions would be required.

## 2.2 Ellipsoidal solutions

We now aim at selecting the matrix M appearing in Darcy's law (7) so that the tumor domain, denoted Ω(t), remains an ellipsoid over time. We assume that the ellipsoid is centered at the origin and described by its three radii L<sub>x1</sub>(t), L<sub>x2</sub>(t), L<sub>x3</sub>(t):

$$ \Omega(t) = \{ \vec{x} = (x_1, x_2, x_3) \in \mathbb{R}^3 : \frac{x_1^2}{L_{x_1}(t)^2} + \frac{x_2^2}{L_{x_2}(t)^2} + \frac{x_3^2}{L_{x_3}(t)^2} \le 1 \}. $$

Because the velocity is linear (Darcy's law with a quadratic pressure),

$$ \vec{v}(t, \vec{x}) = -K M(t) \vec{x}, $$

the boundary evolution is governed by a linear deformation. A point located on the boundary remains on the boundary ∂Ω(t) (hence the domain remains ellipsoid) if and only if each radius satisfies an ODE of the form

$$ L'_{\alpha}(t) = -k_{\alpha} M_{\alpha}(t) L_{\alpha}(t), \quad \alpha = x_1, x_2, x_3, $$

with the constraint $\sum_{\alpha} k_{\alpha} M_{\alpha}(t) = -\tau^G(t)$. In order to impose controlled (and possibly anisotropic) growth along each axis, we select parameters g<sub>α</sub> with

$$ g_{\alpha} \ge 0, \quad \sum_{\alpha=x_1,x_2,x_3} g_{\alpha} = 3, $$

and enforce

$$ L'_{\alpha}(t) = \frac{g_{\alpha}}{3} \tau^G(t) L_{\alpha}(t). $$

Injecting this expression into the previous relation yields

$$ M_{\alpha}(t) = -\frac{\tau^G(t) g_{\alpha}}{3 k_{\alpha}}. $$

This choice guarantees that:
• the radii follow Gompertz-type laws whenever $\tau^G(t)$ does,

• the volume growth is consistent with the divergence constraint,

• the domain remains an ellipsoid at all times.

A simple choice consists in taking $g_\alpha = k_\alpha$, in which case

$$M_{\alpha}(t) = -\frac{\tau^{G}(t)}{3}, \quad L'_{\alpha}(t) = \frac{k_{\alpha}}{3}\tau^{G}(t)L_{\alpha}(t), \qquad (10)$$

with the normalization $\sum_{\alpha} k_{\alpha} = 3$. This keeps the model interpretable while allowing anisotropic ellipsoidal deformation of the tumor spheroid.

## 2.3 Radial equations for the 3D model

The idea now is to rewrite the expression for the dynamics of $P$ and $Q$ given by (1) in the context of an ellipsoid domain, and assuming that their dynamics are invariant by rotation. To do this, we introduce a change of coordinates

$$\vec{x}_L(t) = (x_{1,L}(t), x_{2,L}(t), x_{3,L}(t)) := \left( \frac{x_1}{L_{x_1}(t)}, \frac{x_2}{L_{x_2}(t)}, \frac{x_3}{L_{x_3}(t)} \right) \qquad (11)$$

that corresponds to normalization in each direction. We now introduce the radius

$$r(\vec{x}_L(t)) = \|\vec{x}_L(t)\| \in [0, 1].$$

Under the assumption of invariance by rotation, we introduce, for any $(t, \vec{x}) \in [0, T] \times \Omega(t)$ and any associated function $u$, the notation :

$$u(t, \vec{x}) = \bar{u}(t, r(\vec{x}_L(t))) \text{ for any } P, Q \text{ and } \tau^{P2Q}.$$

The idea is now to express the dynamics of proliferating cells in these new coordinates in order to simplify their expressions. One can show the following theorem, see Appendix A for the proof.

**Theorem 1.** In the ellipsoidal domain $\Omega(t)$, under the assumption of invariance by rotation, the dynamics of proliferating cells $P(t)$ and quiescent cells $Q(t)$ given by (1) are simplified in $\Omega(t)$ by

$$\left\{ \begin{array}{l} \partial_t \bar{P} = \tau^G(t)(1 - \bar{P}) - \bar{\tau}^{P2Q} \bar{P}, \\ \bar{Q} = 1 - \bar{P}. \end{array} \right. \qquad (12)$$

To close the system, the evolution of the rate $\tau^{P2Q}$ has to be described. Since it represents the lack of oxygen, one can consider the following logistic function

$$\bar{\tau}^{P2Q}(t, (L_{\alpha})_{\alpha}, r) = \tau^0 - \frac{\tau^0 - \tau^1}{1 + \exp\left(\sqrt{L_{x_1}^2(t) + L_{x_2}^2(t) + L_{x_3}^2(t)}(1-r)/s - d/s\right)}, \qquad (13)$$

with $\tau^0 > \tau^1 \ge 0$, $r \in [0, 1]$ et $d, s \ge 0$. The parameter $d$ corresponds to the distance from the tumor front at which the oxygen is easily accessible. The parameter $s$ corresponds to the slope of the function. The parameters $\tau^0$ (resp. $\tau^1$) corresponds to the value of $\tau^{P2Q}$ far from (resp. close to) the tumor boundary.
### 3 Data assimilation strategy

Combining the information from the previous section and accounting for model disturbances and uncertainties in the initial conditions, we assume that the growth of the spheroid, represented by its three axes ($L_{x_1}$, $L_{x_2}$, $L_{x_3}$) and the radial proliferative density $\bar{P}$, can be described by the following PDE/ODE system for $t > 0$ and $0 \le r = \|x_L(t)\| \le 1$:

$$ 
\left\{
\begin{aligned}
\partial_t \bar{P}(t,r) &= \tau^G(t)(1 - \bar{P}(t,r)) - \bar{\tau}^{P2Q}(t, (L_\alpha)_\alpha, r) \bar{P}(t,r) + \nu_P(t,r), \\
\bar{P}(0,r) &= P^0(r) + \xi_P(r), \\
L'_\alpha(t) &= \frac{k_\alpha}{3} \tau^G(t) L_\alpha(t) + \nu_\alpha(t), \quad \alpha = x_1, x_2, x_3, \\
L_\alpha(0) &= L_\alpha^0 + \xi_{L_\alpha},
\end{aligned}
\right.
\qquad (14) $$

where $L_\alpha^0$ and $P^0$ are the prior values of the initial conditions, assumed known, while $\xi_{L_\alpha}$ and $\xi_P$ represent the associated uncertainties. The terms $\nu_P$ and $\nu_\alpha$ correspond to model disturbances, accounting for errors arising from model approximations. In the following, we consider that $\tau^G(t) = ae^{-bt}$, so that the spheroid volume follows a Gompertz law.

Our objective is to estimate the model parameters $a, b, k_x, k_y, k_z$ using the available observations (all the radii but potentially with very poor quality or only two of them). One can easily show that the parameters of $\tau^{P2Q}$ are not identifiable and their values are fixed to reproduce heterogeneity similar to the data available in [12] ($d = 0.4, s = 0.05, \tau^0 = 0.02, \tau^1 = 0.005$).

To do so, we follow these steps:

**Step 1** First, assuming that the true values of the parameters are known, we construct state observers whose role is to correct over time the errors due to model disturbances and uncertainties in the initial conditions. Three cases will be considered:

(a) when observations are $t \mapsto L_\alpha(t)$, $\alpha = x_1, x_2, x_3$,

(b) when observations are $t \mapsto L_\alpha(t)$, $\alpha = x_1, x_2$, together with the initial volume $V_0 = \frac{4}{3}\pi L_{x_1}(0)L_{x_2}(0)L_{x_3}(0)$,

(c) when observations are $t \mapsto L_\alpha(t)$, $\alpha = x_1, x_2$.

We will theoretically show that in cases (1) and (2) and without disturbance ($\nu_P = \nu_\alpha = 0$, $\alpha = x_1, x_2, x_3$), all uncertainties (full axes and proliferative density) can be corrected, whereas in case (3), only the uncertainties on the first two axes can be corrected.

**Step 2** Second, the state observers will be coupled with a Reduced-order Unscented Kalman Filter (RoUKF), allowing not only to correct for model disturbances and uncertainties in initial conditions but also to estimate the model parameters.

Numerical simulations on synthetic data, including model disturbances and measurement noise, will be used to illustrate these results<fnref n="1" />.

<footnote n="1">¹Once the paper is accepted, the code and the pre-processed real data will be made publicly available online.</footnote>
## 3.1 State observer

**Design of the observer and convergence** We consider the following state observer

$$ 
\left\{
\begin{aligned}
\partial_t \hat{P}(t,r) &= \tau^G(t)(1 - \hat{P}(t,r)) - \bar{\tau}^{P2Q}(t, (\hat{L}_\alpha)_\alpha, r) \hat{P}(t,r), \\
\hat{P}(0,r) &= P^0(r), \\
\hat{L}_\alpha'(t) &= \frac{k_\alpha}{3} ae^{-bt} \hat{L}_\alpha(t) + G_\alpha(\hat{L}_\alpha, (L_\beta)_{\beta=x_1,x_2,x_3}), \quad \alpha = x_1, x_2, x_3, \\
\hat{L}_\alpha(0) &= \hat{L}_\alpha^0,
\end{aligned}
\right.
\quad (15) $$

where the gain observers $G_\alpha$ are defined for $\gamma_{x_1}, \gamma_{x_2}, \gamma_{x_3} > 0$ by:

$$ G_{x_1}(\hat{L}_{x_1}, (L_\beta)_\beta) = -\gamma_{x_1}(\hat{L}_{x_1} - L_{x_1}), \quad (16) $$

$$ G_{x_1}(\hat{L}_{x_2}, (L_\beta)_\beta) = -\gamma_{x_2}(\hat{L}_{x_2} - L_{x_2}), \quad (17) $$

$$ G_{x_3}(\hat{L}_{x_3}, (L_\beta)_\beta) = \begin{cases} -\gamma_{x_3}(\hat{L}_{x_3} - L_{x_3}), & \text{Case (a)}, \\ \hat{L}_{x_3} \left( \gamma_{x_1}(\hat{L}_{x_1} - L_{x_1}) / \hat{L}_{x_1} + \gamma_{x_2}(\hat{L}_{x_2} - L_{x_2}) / \hat{L}_{x_2} \right), & \text{Case (b)}, \\ 0, & \text{Case (c)}. \end{cases} \quad (18) $$

We assume that there are unique solutions of System (14) and (15) and that these solutions are smooth (this can be proved by a Banach contraction mapping theorem). In order to analyze the convergence of the observer we will study the estimation errors defined by $\tilde{L}_\alpha = L_\alpha - \hat{L}_\alpha$ and $\tilde{P} = \bar{P} - \hat{P}$. One can prove the following theorem, see Appendix B for the proof.

**Theorem 2.** Let $\tau^0 \neq b$ and $\tau^G(t) = ae^{-bt}$. When there is no disturbance (i.e. $\nu_P = 0$ and $\nu_{L_\alpha} = 0$ for all $\alpha$) and non-null priors $L_\alpha^0$, $\alpha = x_1, x_2, x_3$, we have that $t \to L_{x_1}(t)$ and $t \to \tilde{L}_{x_2}(t)$ converge exponentially to 0 when $t$ goes to $\infty$. Furthermore, in Cases (a) and (b), $t \to \tilde{L}_{x_3}(t)$ converges exponentially to 0 when $t$ goes to $\infty$. Moreover, in Case (a) if $\min(\gamma_{x_1}, \gamma_{x_2}, \gamma_{x_3}) - \tau^0 > 0$ and $\min(\gamma_{x_1}, \gamma_{x_2}, \gamma_{x_3}) + b - 2\tau^0 > 0$ and in Case (b) if $\min(\gamma_{x_1}, \gamma_{x_2}) - \tau^0 > 0$ and $\min(\gamma_{x_1}, \gamma_{x_2}) + b - 2\tau^0 > 0$, then $t \to \|\tilde{P}\|_{L^2([0,1])}$ converges to 0 when $t$ goes to $\infty$.

**Numerical illustrations with synthetic data** We numerically illustrate the observer model on the three cases. We consider the following initial conditions:

$$ 
\begin{aligned}
L_{x_1}(0) = 0.35 \text{ mm}, \quad \hat{L}_{x_1}(0) = 0.15, \quad L_{x_2}(0) = 0.35, \quad \hat{L}_{x_2}(0) = 0.5, \\
L_{x_3}(0) = 0.175, \quad \hat{L}_{x_3}(0) = \frac{L_{x_1}(0)L_{x_2}(0)L_{x_3}(0)}{\hat{L}_{x_1}(0)\hat{L}_{x_2}(0)} - 0.05, \\
\bar{P}(0, r) = 0.8, \quad \hat{P}(0, r) = 0.7,
\end{aligned}
$$

except for Case (b) where we consider $\hat{L}_{x_3}(0) = \frac{L_{x_1}(0)L_{x_2}(0)L_{x_3}(0)}{\hat{L}_{x_1}(0)\hat{L}_{x_2}(0)}$.

A temporally colored model disturbance on $L_{x_3}$ modeled by an Ornstein-Uhlenbeck process is considered with a mean-reversion rate of 10 and a noise amplitude of 0.002. The numerical approximation is obtained using an explicit Euler scheme ($T = 30$ h, $\Delta t = 0.1$ h) for the time discretization. For the spatial discretization, we use a first-order finite-difference scheme for $P$ with spatial step $\Delta x = 0.01$. The parameter values are set to $a = 0.04, b = 0.13, k_x = 2.4, k_y = 2.4, k_z = -1.8$.
<figure><img src="image_2.png" /><figcaption>Figure 1: Temporal evolution of the three axes over time (left: $L_{x1}$, middle: $L_{x2}$ and right: $L_{x3}$). First row: Case (a). Second row: Case (b). Third row: Case (c).</figcaption></figure>

To handle the time sampling of the data (here each hour), two strategies may be considered: either using the data only at observation times or relying on time interpolation. Here, we consider one observation every hour and we linearly interpolate the data. Additive white noise of realistic intensity (0.01 for $L_{x1}$ and $L_{x2}$, 0.005 for $L_{x3}$) is added. Finally, the following gain values are considered: $\gamma_{x1} = \gamma_{x2} = \gamma_{x3} = 0.5$.

Figure 1 illustrates the time evolution of the three axes in the three scenarios (Top: Case (a), Middle: Case (b), Bottom: Case (c)) for both the target and observer solutions. Figure 2 shows the corresponding evolution of the proliferative density. For comparison for the third axis, the simulation obtained with the same parameter set but without the added disturbance is also displayed, highlighting both the nature and the impact of the introduced model noise.

In Case (a), the errors on the initial condition and on the model operator are well corrected for the three axes as well as for *P*, although the convergence toward the target proliferative density remains slow. Decreasing the value of *b* or increasing the gain accelerate it. Moreover, as the observation time window is generally short in real datasets, it seems difficult to expect a perfect convergence for the proliferative
<figure><img src="image_3.png" /><figcaption>Figure 2: Snapshots of the proliferative density *P*. First column: Case (a). Second column: Case (b). Third column: Case (c).</figcaption></figure>

density. The main difference with Case (b) is that the model disturbance cannot be corrected: $\hat{L}_{x_3}$ converges to a *noisy* version of $L_{x_3}$ without disturbance (the noise being inherited from that of $L_{x_1}$ and $L_{x_2}$). In Case (c), only $\hat{L}_{x_1}$ and $\hat{L}_{x_2}$ converge to their target values. These numerical results are consistent with the theoretical predictions.

## 3.2 Joint state and parameter observer

As previously mentioned, errors in the initial conditions and model disturbances are addressed using the Luenberger observer introduced earlier, while the parameters are estimated via a population reduced-order Kalman-based filter restricted to the parametric space. To combine these approaches, we adopt the joint state-parameter estimation framework presented in [14] and further developed in [13]. This strategy

<table><caption>Table 1: Mean squared error (MSE) and bias for the three parameters *a*, *b*, and $r_{x_3}$ across the four observation scenarios. Disturbance indicates whether a model disturbance was added.</caption><thead><tr><th>Observation</th><th>Disturbance</th><th>MSE <i>a</i></th><th>Bias <i>a</i></th><th>MSE <i>b</i></th><th>Bias <i>b</i></th><th>MSE <i>r</i><sub>x<sub>3</sub></sub></th><th>Bias <i>r</i><sub>x<sub>3</sub></sub></th></tr></thead><tbody><tr><td>3 axes</td><td>No</td><td>0.054</td><td>0.145</td><td>0.022</td><td>0.063</td><td>0.015</td><td>0.090</td></tr><tr><td>3 axes</td><td>Yes</td><td>0.082</td><td>0.15</td><td>0.022</td><td>0.062</td><td>0.055</td><td>0.085</td></tr><tr><td>2 axes and <i>V</i><sub>0</sub></td><td>No</td><td>0.11</td><td>0.092</td><td>0.013</td><td>0.033</td><td>0.14</td><td>0.0065</td></tr><tr><td>2 axes</td><td>No</td><td>0.12</td><td>0.12</td><td>0.025</td><td>0.052</td><td>0.13</td><td>0.011</td></tr></tbody></table>
<figure><img src="image_4.png" /><figcaption>Figure 3: Estimated parameters a, b, and $r_{x3}$ versus target values for 500 replicates. First column: all three axes observed (case (a)) with no model disturbance. Second column: all three axes observed (case (a)) with a temporally colored disturbance on $L_{x3}$. Third column: $L_{x1}$, $L_{x2}$, and $V_0$ observed; nudging on $L_{x3}$ (case (b)). Fourth column: only $L_{x1}$ and $L_{x2}$ observed; no nudging on $L_{x3}$ (case (c)).</figcaption></figure>

assumes that a suitable Luenberger observer has already been designed to control the trajectory, which is the case in our study.

Regarding the parameters $k_{x1}$, $k_{x2}$, $k_{x3}$, only two need to be estimated since $k_{x1} + k_{x2} + k_{x3} = 3$. Assuming axial symmetry along z, we have $k_{x1} = k_{x2}$, leaving only one parameter to estimate. To facilitate interpretation and estimation, instead of directly estimating $k_{x1}$, we estimate the ratio between $k_{x1}$ (or $k_{x2}$) and $k_{x3}$, denoted by $r_{x3} = \frac{k_{x3}}{k_{x1}}$, assuming $k_{x1} \neq 0$. Using the equality

$$k_{x1} + k_{x2} + k_{x3} = (2 + r_{x3})k_{x1} = 3,$$

we must ensure that $r_{x3} \neq -2$. In our case, we restrict $r_{x3} \in ]-2, 2[$. Biologically, this constraint reflects that the spheroid cannot collapse or excessively elongate along the z-axis relatively to the horizontal axes; in other words, the deformation along z remains within realistic bounds compared to radial expansion in the horizontal plane. To incorporate this constraint along with positivity constraints for a and b, we consider the following prior distributions:

$$ \operatorname{logit}\left(\frac{r_{x3}+2}{4}\right) \sim \mathcal{N}\left(\operatorname{logit}\left(\frac{m_{r_{x3}}+2}{4}\right), \sigma_{r_{x3}}\right), \\ \log(a) \sim \mathcal{N}(\log(m_a), \sigma_a), \quad \log(b) \sim \mathcal{N}(\log(m_b), \sigma_b). $$
<figure><img src="image_5.png" /><figcaption>Figure 4: Temporal evolution for a replicate in Scenario 2. Top: state variables ($L_{x_1}$, $L_{x_2}$, $L_{x_3}$) over time. Bottom: parameter estimates ($a$, $b$, $r_{x_3}$) over time.</figcaption></figure>

To validate our methodology and implementation, we first consider a synthetic data benchmark. The goal is to simulate scenarios close to reality to guide estimation strategies for real applications. For all cases, $N_R = 500$ replicates are generated. Target initial conditions are drawn from Gaussian distributions: $L_{x_1}(0), L_{x_2}(0) \sim N(0.35, 0.01)$ and $L_{x_3}(0) \sim N(0.175, 0.01)$. As previously, additive white noise of realistic intensity (0.005 for $L_{x_1}$ and $L_{x_2}$, 0.001 for $L_{x_3}$) is added, and hourly time linear interpolation is applied. The initial condition is set to the noisy first observation. Target parameters are drawn according to the distributions defined above with $r_{x_3} = -0.75$, $\sigma_{r_{x_3}} = 0.5$, $m_a = 0.04$, $\sigma_a = 0.3$, $m_b = 0.13$ and $\sigma_b = 0.3$, while prior values are fixed to $p_{r_{x_3}} = 0$, $\sigma_{p_{r_{x_3}}} = 1$, $p_a = 0.01$, $\sigma_p = 1$, $p_b = 0.08$ and $\sigma_p = 1$. The observer state gains are set to $\gamma_{x_1} = \gamma_{x_2} = \gamma_{x_3} = 0.5$. Four simulation scenarios are considered:

* All three axes are observed (case (a)) with no model disturbance.

* All three axes are observed (case (a)) with a temporally colored disturbance on $L_{x_3}$ modeled by an Ornstein-Uhlenbeck process.

* $L_{x_1}$, $L_{x_2}$, and $V_0$ are observed; nudging is applied to $L_{x_3}$ (case (b)).

* Only $L_{x_1}$ and $L_{x_2}$ are observed; no nudging on $L_{x_3}$ (case (c)) and the prior of this axis is fixed to 0.2.

The covariance matrix of the observation errors is assumed to be time-independent and is fixed at $0.005^2 I_3$, except in Scenario 2, where a larger variance of $0.01^2$ is assigned to the third axis to account for the lower confidence in its measurements. For each parameter $\theta \in \{a, b, r_{x_3}\}$, we compute

$$ \text{BIAS}(\theta) = \frac{1}{N_R} \sum_{i=1}^{N_R} (\theta_i^{\text{target}} - \theta_i^{\text{estim}}), \quad \text{MSE}(\theta) = \frac{1}{N_R} \sum_{i=1}^{N_R} (\theta_i^{\text{target}} - \theta_i^{\text{estim}})^2, $$
see Table 1 for the results. Figure 3 shows estimated versus target values for each parameter (row) and scenario (column). Figure 4 illustrates the temporal evolution of the 3 axes and parameters for a replicate in the second scenario. One can see that parameters now evolve over time. For comparison, the simulation obtained with the same parameter set but without the added disturbance is also plotted, highlighting the nature and effect of the model noise introduced.

The results indicate that when no information is available on $L_{x_3}$, parameter estimates deteriorate significantly, as reflected in both the MSE and the plots, particularly for $a$ and $r_{x_3}$, even when $V_0$ is well known and the state observer incorporates this information. This highlights a theoretical or practical identifiability issue. However, when some partial knowledge of $L_{x_3}$ dynamics is available - for instance, through a correlation between pixel intensity and spheroid height, modeled here as a model error – the parameters can still be reasonably estimated. The good quality of the parameter estimation is due to the use of a Luenberger observer, which corrects not only the initial conditions but also the model disturbances. Importantly, even if the model is not fully accurate, the biological interpretation of the parameters remains meaningful, especially for comparative analyses, such as between untreated and treated conditions. Therefore, in our real application case, we will still aim to acquire data on $L_{x_3}$, even if of very poor quality, in order to mimic Scenario 2.

# 4 Application to a real case

To illustrate the capabilities of the 3D growth model and the data assimilation procedure, we examine the growth of angiosarcoma tumor spheroids with and without propranolol treatment.

## 4.1 Problematic and data presentation

**Angiosarcoma and propranolol** Angiosarcomas (AS) is a soft tissue tumor of endothelial origin. In fact, AS is a rare, clinically aggressive, difficult-to-treat malignant tumor with a 5-year survival rate of about 20-35% [23]. Chemotherapy remains the primary treatment for metastatic AS, but its use is limited due to comorbidities in older patients and the toxicity of the drugs used [5]. Among the targeted therapies under consideration, VEGF targeting and the use of tyrosine kinase inhibitors have shown promise for the treatment of AS [20].

However, as there is no satisfactory therapy for AS, the development of new treatment methods is required. In this context, propranolol – a non-selective beta-blocker commonly used to treat cardiovascular diseases – has emerged as a novel therapeutic solution for vascular tumors, especially angiosarcomas [2]. Indeed, due to its antiproliferative and antiangiogenic properties, it has been shown to be potentially suitable for cancer therapy. A therapeutic effect has been demonstrated in patients with angiosarcomas, but it remains unclear which patients will benefit and at what dose. Although research on infantile hemangioma has uncovered part of the non-cytotoxic mechanism, current study models for angiosarcoma are not advanced enough to decipher complex antitumor mechanisms [19, 15].

To illustrate the impact of propranolol on angiosarcoma growth, we analyze the 3D dynamics of tumor spheroids treated or not with propranolol, using the 3D
growth model and its reduced formulation presented in Section 2, together with
the data assimilation procedure described in Section 3.

<figure><img src="image_6.png" /><figcaption>Figure 5: Top: Available data with bottom-view imaging at time $t=20h$ for the three cases (control, propranolol at 10 µM and 50 µM) on the left. Front and side views of a control spheroid at time $t=0h$ on the right. Middle and bottom: Directly extracted data (mean intensity and $L_\alpha$, $\alpha = x_1, x_2$) from the bottom-view images on the left. Linear regression between height $H(0, \cdot)$ and population-averaged intensity $I(0, \cdot)$ among control spheroids using a corona-based partition</figcaption></figure>

and reconstruction of $L_{x_3}$ on the right.

**Data presentation** The biological experiments consist of culturing spheroids in multi-well plates. Three experimental setups are considered: a control experiment corresponding to free growth, and two experiments in which propranolol is added at two different concentrations (10 µM and 50 µM). For each experimental condition, the temporal evolution of approximately ten spheroids is monitored. We have at our disposal: (1) bottom-view images of the spheroids acquired every hour over a 30-hour period for the three experimental conditions (control, propranolol 10 µM and 50 µM), see Figure 5-Top-Left; and (2) a 3D acquisition including bottom and side views of a spheroid in the control case at initial time, see Figure 5-Top-Right. The
3D acquisition, although outside the scope of this study, was not obtained under exactly the same biological conditions as the bottom-view image series. Further details on the experimental protocols for bottom-view images can be found in [18].

Using the OpenCV Python library, all bottom-view images are automatically segmented and two relevant quantitative features are extracted: (1) the spheroid radius over the 30-hour period, allowing for the direct estimation of the time-dependent data $L_{x_1}$ and $L_{x_2}$, see Figure 5-Bottom-Left; and (2) the mean pixel intensity over the same period, see Figure 5-Middle-Left.

It can be observed in Figure 5-Bottom-Left that, at $t = 30$ h, the radii $L_{x_1}$ and $L_{x_2}$ reach an average value of approximately 0.47 mm in the control case, whereas they exceed 0.5 mm when propranolol is added. Focusing solely on these measurements could therefore suggest that propranolol enhances tumor growth. However, in all cases, the mean intensity decreases over time (see Figure 5-Middle-Left), indicating a progressive flattening of the spheroids. More specifically, the mean intensity is approximately 0.89 in the control case and decreases to about 0.84 in the presence of propranolol. These observations suggest that the flattening of the spheroids is amplified when propranolol is present, indicating that propranolol may induce alterations in spheroid cohesion. Our objective is to estimate the parameters of the model introduced later in order to quantify the impact of propranolol.

The tests performed on synthetic data in Section 3.2 motivate the use of data along all spatial axes. Even when the data associated with $L_{x_3}$ are highly noisy, these tests show that incorporating such information is preferable to having no data at all, as it still improves parameter estimation despite model uncertainty and measurement noise. To this end, we aim to exploit the available measurements to reconstruct even very noisy time-dependent data for $L_{x_3}$, as described in the next paragraph.

Regarding the side-view acquisition, it should be noted that the spheroid is sacrificed in order to obtain images from multiple angles, which explains why such data are available for only one spheroid. Consequently, instead of considering the ten spheroids in each cohort individually, the measured quantities (radii and mean intensity) are averaged within each cohort. This results in three sets of measurements corresponding to a representative mean spheroid for each experimental condition.

**$L_{x_3}$ construction with intensity** Direct observations of $L_{x_3}(t)$ over time are not accessible, since only a single 3D view is available at the initial time. To overcome this limitation, we propose an indirect estimation based on the relationship between image intensity and spheroid height.

We first introduce some geometric preliminaries. Along the $x_1$-axis, let $l_{x_1} \in ]-L_{x_1}(t), L_{x_1}(t)[$ denote the radial coordinate measured from the spheroid center. We define $H(t, l_{x_1})$ as the height (diameter) of the spheroid at time $t$ and position $l_{x_1}$. In particular, at the spheroid center ($l_{x_1} = 0$), we have $H(t, 0) = 2L_{x_3}(t)$. Assuming an ellipsoidal geometry (see Figure 6), the height profile reads

$$H(t, l_{x_1}) = 2\sqrt{1 - \frac{l_{x_1}^2}{L_{x_1}(t)^2}} L_{x_3}(t).$$

At the initial time $t = 0$, the 3D view suggests that $L_{x_3}(0) \approx \frac{1}{2}L_{x_1}(0)$ (see
In the plane $x_2 = 0$, elliptical contour: $\frac{x_1^2}{L_{x_1}(t)^2} + \frac{x_3^2}{L_{x_3}(t)^2} = 1$.

<figure><img src="image_7.png" /><figcaption>Figure 6: Illustration of the definition of $(t, l_{x_1}) \rightarrow H(t, l_{x_1})$.</figcaption></figure>

Figure 5-Top-Right). Hence, the initial height profile is

$$H(0, l_{x_1}) = \sqrt{1 - \frac{l_{x_1}^2}{L_{x_1}(0)^2}} L_{x_1}(0). \qquad (19)$$

We denote by $I(t, l_{x_1})$ the normalized intensity at time $t$ and position $l_{x_1}$. The control bottom-view images at initial time are partitioned into 25 concentric coronas with radii $(l_{x_1})_{i=1,...,25} \subset (0, L_{x_1}(0))$. For each corona, the mean intensity is computed and then averaged over all control spheroids, yielding the dataset $\{I_{\text{data}}(0, (l_{x_1})_i)\}_{i=1,...,25}$.

Using (19), we associate to each radius $(l_{x_1})_i$ the corresponding height $H(0, (l_{x_1})_i)$. The pairs $(I_{\text{data}}(0, (l_{x_1})_i), H(0, (l_{x_1})_i))$ are then used to fit a linear relation between image intensity and spheroid height (see Figure 5-Middle-Right).

Finally, assuming that this intensity-height relation is independent of time, for each time point $t$ and each experiment (control and propranolol), the intensity at the spheroid center $I(t, 0)$ is converted into a height using the fitted law, providing an estimate of the central height

$$H(t, 0) = 2L_{x_3}(t),$$

and thus of $L_{x_3}(t)$ (see Figure 5-Bottom-Right).

## 4.2 Estimation results

To mimic the previously presented Scenario 2, we perform joint state-parameter estimation using the state observer on all three axes. Specifically, we consider the case where the different state observers are defined by (16), (17), and (18) (Case- (a)). We keep the same prior distributions for the parameters and the same state
<figure><img src="image_8.png" /><figcaption>Figure 7: Control case. Temporal evolutions of the 3D spheroid (with estimated pa- rameters). The 2D slice shows the proliferation profile P ( r ) in grayscale, with black corresponding to high proliferation and white to low proliferation at times 0, 10, 20 and 30h.</figcaption></figure>

gain values. Since the real data are less noisy, we only adjust the observation covari- ance, setting its diagonal to [0.001$^{2}$, 0.001$^{2}$, 0.002$^{2}$] and assuming zero covariance between the axes.

Table 2 shows the estimated parameters in the real dataset. The parameter a represents the initial growth rate of the spheroid, b − 1 corresponds to the character- istic time over which the growth capacity decreases, and r$_{x}$$_{3}$ describes the growth ratio along the x , y , and z axes. A negative (resp. positive) value of r$_{x}$$_{3}$ indicates a decrease (resp. an increase) along the z axis, while values closer to − 2 correspond to spheroids becoming increasingly flattened.

The results clearly indicate that propranolol strongly affects spheroid growth: both doses lead to a reduction in the initial growth rate a . Moreover, the decrease in r$_{x}$$_{3}$ suggests that the spheroids flatten more rapidly under propranolol treat- ment, reflecting a faster loss of cell adhesion and structural integrity, see Figures 7 and 8 for a 3D reconstruction including the density of proliferative cells. Our esti- mation explains the counter-intuitive observation that the spheroid’s bottom-view area increases under treatment, as the faster flattening along the z axis causes the spheroid to lose structural integrity and collapse, spreading cells over a larger
<figure><img src="image_9.png" /><figcaption>Figure 8: Propranolol (10 µM) case. Temporal evolutions of the 3D spheroid (with estimated parameters) The 2D slice shows the proliferation profile $P(r)$ in grayscale, with black corresponding to high proliferation and white to low proliferation at times 0, 10, 20 and 30h.</figcaption></figure>

# 5 Conclusion

In this work, we addressed the challenge of quantitatively analyzing tumor spheroids from partial and noisy data, especially when their shape deviates from perfect sphericity. We developed a 3D tumor model incorporating three key features: the coexistence of proliferative and quiescent cell populations, axial symmetry, and distinct dynamics along the vertical axis compared to the radial direction. Under the assumptions that spheroid evolution can be approximated by ellipsoidal shapes and that the cells dynamics' have rotational invariance, the model was reduced to a system of three coupled ODEs for the three axes, combined with a 1D PDE describing the dynamics of proliferative cells. Although this model does not yet include cell death, it can be readily extended, and the proposed methodology remains fully applicable in more complex versions.

We further proposed state observers for different observation scenarios: full 3-axis data, only 2 axes with known initial volume, and only 2 axes without initial volume. These observers were theoretically analyzed and tested numerically in a context where modeling errors affect the third axis. The state observer was then coupled with a Kalman filter to estimate model parameters. Simulations on
<table><caption>Table 2: Parameter estimation in the real context.</caption><thead><tr><th>Parameter</th><th>Control</th><th>Propranolol 10µM</th><th>Propranolol 50µM</th></tr></thead><tbody><tr><td>a</td><td>0.033</td><td>0.0024</td><td>0.0037</td></tr><tr><td>b</td><td>0.12</td><td>0.11</td><td>0.10</td></tr><tr><td>r<sub>x3</sub></td><td>-0.98</td><td>-1.94</td><td>-1.90</td></tr></tbody></table>

synthetic data highlighted that including even very noisy measurements of the third axis is preferable to having no data at all for robust parameter estimation.

This approach was applied to experimental spheroid data from angiosarcoma cells treated with propranolol. The third axis was reconstructed using intensity data from the bottom view. Results show that, although the bottom-view area may appear to increase under propranolol, the spheroid is actually collapsing along the vertical axis, and cellular proliferation is strongly reduced compared to control spheroids.

We are confident that this work provides a promising framework for studying 3D spheroids even when they are not perfectly spherical and data are limited, allowing the interpretation of phenomena linked to cell adhesion and spatial organization.

## References

[1] S. Benzekry, C. Lamont, A. Beheshti, A. Tracz, J. M. Ebos, L. Hlatky, and P. Hahnfeldt. Classical mathematical models for description and prediction of experimental tumor growth. *PLoS computational biology*, 10(8):e1003800, 2014.

[2] W. Chow, C. N. Amaya, S. Rains, M. Chow, E. B. Dickerson, and B. A. Bryan. Growth attenuation of cutaneous angiosarcoma with propranolol-mediated β-blockade. *JAMA dermatology*, 151(11):1226-1229, 2015.

[3] A. Collin, H. Bruhier, J. Kolosnjaj, M. Golzio, M.-P. Rols, and C. Poignard. Spatial mechanistic modeling for prediction of 3d multicellular spheroids behavior upon exposure to high intensity pulsed electric fields. *AIMS bioengineering*, 9(2):102-122, 2022.

[4] A. Collin, C. Copol, V. Pianet, T. Colin, J. Engelhardt, G. Kantor, H. Loiseau, O. Saut, and B. Taton. Spatial mechanistic modeling for prediction of the growth of asymptomatic meningiomas. *Computer Methods and Programs in Biomedicine*, 199:105829, 2021.

[5] V. Florou and B. A. Wilky. Current and future directions for angiosarcoma therapy. *Current treatment options in oncology*, 19(3):14, 2018.

[6] L. Guillaume, L. Rigal, J. Fehrenbach, C. Severac, B. Ducommun, and V. Lobjois. Characterization of the physical properties of tumor-derived spheroids reveals critical insights for pre-clinical studies. *Scientific reports*, 9(1):6597, 2019.

[7] F. Hirschhaeuser, H. Menne, C. Dittfeld, J. West, W. Mueller-Klieser, and L. A. Kunz-Schughart. Multicellular tumor spheroids: an underestimated tool is catching up again. *Journal of biotechnology*, 148(1):3-15, 2010.
[8] S. Julier and J. Uhlmann. A new extension of the Kalman filter to nonlinear systems. In *Proc. of AeroSense: The 11th Int. Symp. on Aerospace/Defence Sensing, Simulation and Controls*, 1997.

[9] J. P. Ko, E. J. Berman, M. Kaur, J. S. Babb, E. Bomsztyk, A. K. Greenberg, D. P. Naidich, and H. Rusinek. Pulmonary nodules: growth rate assessment in patients by using serial ct and three-dimensional volumetry. *Radiology*, 262(2):662-671, 2012.

[10] D. Luenberger. An introduction to observers. *IEEE Transactions on Automatic Control*, 16:596-602, 1971.

[11] X. Mao, S. McManaway, J. K. Jaiswal, P. B. Patel, W. R. Wilson, K. O. Hicks, and G. Bogle. An agent-based model for drug-radiation interactions in the tumour microenvironment: Hypoxia-activated prodrug SN30000 in multicellular tumour spheroids. *PLoS computational biology*, 14(10):e1006469, 2018.

[12] T. Michel, J. Fehrenbach, V. Lobjois, J. Laurent, A. Gomes, T. Colin, and C. Poignard. Mathematical modeling of the proliferation gradient in multicellular tumor spheroids. *Journal of theoretical biology*, 458:133-147, 2018.

[13] P. Moireau and D. Chapelle. Reduced-order Unscented Kalman Filtering with application to parameter identification in large-dimensional systems. *ESAIM: Control, Optimisation and Calculus of Variations*, 17(2):380-405, 2011.

[14] P. Moireau, D. Chapelle, and P. Le Tallec. Joint state and parameter estimation for distributed mechanical systems. *Computer Methods in Applied Mechanics and Engineering*, 1987(6-8):659-677, 2008.

[15] F. Moisan, S. Oucherif, P. Kaulanjan-Checkmodine, S. Prey, B. Rousseau, M. Bonneu, S. Claverol, E. Gontier, S. Lacomme, L. Dousset, et al. Critical role of aquaporin-1 and telocytes in infantile hemangioma response to propranolol beta blockade. *Proceedings of the national academy of sciences*, 118(7):e2018690118, 2021.

[16] J. Munoz-Garcia, C. Jubelin, A. Loussouarn, M. Goumard, L. Griscom, A. Renodon-Cornière, M.-F. Heymann, and D. Heymann. In vitro three-dimensional cell cultures for bone sarcomas. *Journal of Bone Oncology*, 30:100379, 2021.

[17] B. Ocaña-Tienda, A. Eroles-Simó, J. Pérez-Beteta, E. Arana, and V. M. Pérez-García. Growth dynamics of lung nodules: implications for classification in lung cancer screening. *Cancer Imaging*, 24(1):113, 2024.

[18] S. Oucherif. *Étude du ciblage de la voie bêta-adrénergique par le propranolol dans le traitement des tumeurs vasculaires cutanées*. PhD thesis, Université de Bordeaux, 2021.

[19] E. Pasquier, N. André, J. Street, A. Chougule, B. Rekhi, J. Ghosh, D. S. Philip, M. Meurer, K. L. MacKenzie, M. Kavallaris, et al. Effective management of advanced angiosarcoma by the synergistic combination of propranolol and vinblastine-based metronomic chemotherapy: a bench to bedside study. *EBioMedicine*, 6:87-95, 2016.

[20] V. Ravi, E. M. Sanford, W.-L. Wang, J. S. Ross, N. Ramesh, A. Futreal, S. Patel, P. J. Stephens, V. A. Miller, and S. M. Ali. Antitumor response of vegfr2-and vegfr3-amplified angiosarcoma to pazopanib. *Journal of the National Comprehensive Cancer Network*, 14(5):499-502, 2016.
[21] S. Shah and G. G. DSouza. Modeling tumor microenvironment complexity in vitro: Spheroids as physiologically relevant tumor models and strategies for their analysis. *Cells*, 14(10):732, 2025.

[22] A. Stéphanou, P. Ballet, and G. Powathil. Hybrid data-based modelling in oncology: successes, challenges and hopes. *Mathematical Modelling of Natural Phenomena*, 15:21, 2020.

[23] R. J. Young, N. J. Brown, M. W. Reed, D. Hughes, and P. J. Woll. Angiosarcoma. *The lancet oncology*, 11(10):983–991, 2010.

A Proof of Theorem 1

In this section, we give the proof of Theorem 1.

*Proof.* First, using the expression of the velocity divergence (3) in the first equation of (1), the dynamics of P reads

$$ \partial_t P + \vec{v} \cdot \nabla P = \tau^G(t)(P + Q)(1 - P) - \tau^{P2Q}P. \quad (20) $$

Let us calculate the time derivative of P with the new variables.

$$ \begin{aligned} \partial_t P(t, \vec{x}) &= \partial_t \bar{P} + \partial_r \bar{P} \left( \frac{\vec{x}_L}{\|\vec{x}_L\|} \cdot \left( -x_1 \frac{L'_{x_1}}{L_{x_1}^2}, -x_2 \frac{L'_{x_2}}{L_{x_2}^2}, -x_3 \frac{L'_{x_3}}{L_{x_3}^2} \right)^T \right), \\ &= \partial_t \bar{P} - \partial_r \bar{P} \left( \frac{\vec{x}_L}{\|\vec{x}_L\|} \cdot \left( x_{1,L} \frac{L'_{x_1}}{L_{x_1}}, x_{2,L} \frac{L'_{x_2}}{L_{x_2}}, x_{3,L} \frac{L'_{x_3}}{L_{x_3}} \right)^T \right), \\ &= \partial_t \bar{P} - \partial_r \bar{P} \left( \frac{\vec{x}_L}{\|\vec{x}_L\|} \cdot \frac{\tau^G(t)}{3} (x_{1,L}k_{x_1}, x_{2,L}k_{x_2}, x_{3,L}k_{x_3})^T \right), \\ &= \partial_t \bar{P} - \partial_r \bar{P} \frac{\tau^G(t)}{3} \frac{\vec{x}_L \cdot (K\vec{x}_L)}{\|\vec{x}_L\|}. \end{aligned} \quad (21) $$

Let us now calculate the space derivative of P with these new variables.

$$ \nabla_{\vec{x}} P(t, \vec{x}) = \partial_r \bar{P} \frac{1}{r(\vec{x}_L)} \left( \frac{x_1}{L_{x_1}^2}, \frac{x_2}{L_{x_2}^2}, \frac{x_3}{L_{x_3}^2} \right) = \partial_r \bar{P} \frac{1}{r(\vec{x}_L)} \left( \frac{x_{1,L}}{L_{x_1}}, \frac{x_{2,L}}{L_{x_2}}, \frac{x_{3,L}}{L_{x_3}} \right). $$

We express the dynamics of P given by (20) with the new variables.

$$ \begin{aligned} \partial_t \bar{P} + \left( \frac{1}{r(\vec{x}_L)} \left( v_{x_1} \frac{x_{1,L}}{L_{x_1}} + v_{x_2} \frac{x_{2,L}}{L_{x_2}} + v_{x_3} \frac{x_{3,L}}{L_{x_3}} \right) - \frac{\tau^G(t)}{3} \frac{\vec{x}_L \cdot (K\vec{x}_L)}{\|\vec{x}_L\|} \right) \partial_r \bar{P} \\ = \tau^G(t)(\bar{P} + \bar{Q})(1 - \bar{P}) - \bar{\tau}^{P2Q} \bar{P} \quad (22) \end{aligned} $$

On the other hand, from the explicit expression for $M(t)$ given by (10), the velocity defined by (7) becomes

$$ \vec{v} = \frac{\tau^G(t)}{3} KX \iff \begin{pmatrix} v_{x_1} \\ v_{x_2} \\ v_{x_3} \end{pmatrix} = \begin{pmatrix} \frac{\tau^G(t)}{3} k_{x_1} x_1 \\ \frac{\tau^G(t)}{3} k_{x_2} x_2 \\ \frac{\tau^G(t)}{3} k_{x_3} x_3 \end{pmatrix}. $$
Using this result in (22), we then have

$$ \partial_t \bar{P} + \left( \frac{1}{r(\vec{x}_L)} \frac{\tau^G(t)}{3} \left( k_{x_1} x \frac{x_{1,L}}{L_{x_1}} + k_{x_2} y \frac{x_{2,L}}{L_{x_2}} + k_{x_3} z \frac{x_{3,L}}{L_{x_3}} \right) - \frac{\tau^G(t)}{3} \frac{\vec{x}_L \cdot (K\vec{x}_L)}{r(\vec{x}_L)} \right) \partial_r \bar{P} \\ = \tau^G(t)(\bar{P} + \bar{Q})(1 - \bar{P}) - \bar{\tau}^{P2Q} \bar{P}. \quad (23) $$

Using the expression of the new variables (11), we get

$$ \partial_t \bar{P} + \left( \frac{1}{r(\vec{x}_L)} \frac{\tau^G(t)}{3} \left( k_{x_1} x_{1,L}^2 + k_{x_2} x_{2,L}^2 + k_{x_3} x_{3,L}^2 \right) - \frac{\tau^G(t)}{3} \frac{\vec{x}_L \cdot (K\vec{x}_L)}{r(\vec{x}_L)} \right) \partial_r \bar{P} \\ = \tau^G(t)(\bar{P} + \bar{Q})(1 - \bar{P}) - \bar{\tau}^{P2Q} \bar{P}. $$

From the saturation hypothesis (2), in the tumoral domain $\Omega(t)$, we have $\bar{P}+\bar{Q}=1$, then we can conclude that the dynamics of $\bar{P}$ is as follows:

$$ \partial_t \bar{P} = \tau^G(t)(1 - \bar{P}) - \bar{\tau}^{P2Q} \bar{P}. $$

Finally, using the saturation hypothesis (2), the dynamics of $\bar{Q}$ in the tumoral domain $\Omega(t)$ reads

$$ \bar{Q} = 1 - \bar{P}. $$

<figure><img src="image_10.png" /></figure>

# B Proof of Theorem 2

In this section, we present the proof of Theorem 2.

## Proof. Step 1: Convergence of $\check{L}_{x_1}$ and $\check{L}_{x_2}$

We first establish the convergence of the errors $\check{L}_{x_1}$ and $\check{L}_{x_2}$. Using (14), (15) with gains defined by (16) and (17), the error dynamics $\check{L}_\alpha(t)$ for $\alpha = x_1, x_2$ satisfy

$$ \check{L}'_\alpha(t) = \left( \frac{k_\alpha}{3} ae^{-bt} - \gamma_\alpha \right) \check{L}_\alpha(t). \quad (24) $$

Integrating this expression yields

$$ \check{L}_{\alpha}(t) = \check{L}_{\alpha}(0)e^{\frac{k_{\alpha}}{3}\frac{a}{b}(1-e^{-bt})-\gamma_{\alpha}t}, \quad \alpha = x_1, x_2, \quad (25) $$

which establishes exponential convergence to zero as $t \to \infty$.

## Step 2: Convergence of $\check{L}_{x_3}$

**Case (a):** With the gain from (18)-Case (a), we obtain

$$ \check{L}_{x_3}(t) = \check{L}_{x_3}(0) \exp \left( \frac{k_{x_3}}{3} \frac{a}{b} (1 - e^{-bt}) - \gamma_{x_3} t \right), \quad (26) $$

which proves exponential convergence to zero.
**Case (b):** To deal with the gain from (18)-Case (b), we introduce the volume of the ellipsoid defined by the observers radii:

$$ \hat{V}(t) := \frac{4}{3}\pi \hat{L}_{x_1}(t) \hat{L}_{x_2}(t) \hat{L}_{x_3}(t). \qquad (27) $$

The corresponding volume error is

$$ \check{V}(t) := \hat{V}(t) - V(t) = \frac{4}{3}\pi \left( \hat{L}_{x_1}(t) \hat{L}_{x_2}(t) \hat{L}_{x_3}(t) - L_{x_1}(t) L_{x_2}(t) L_{x_3}(t) \right). $$

Its time derivative satisfies

$$ \begin{aligned} \check{V}' &= \frac{4}{3}\pi \left( \hat{L}'_{x_1} \hat{L}_{x_2} \hat{L}_{x_3} + \hat{L}_{x_1} \hat{L}'_{x_2} \hat{L}_{x_3} + \hat{L}_{x_1} \hat{L}_{x_2} \hat{L}'_{x_3} \right. \\ & \qquad \left. - L'_{x_1} L_{x_2} L_{x_3} - L_{x_1} L'_{x_2} L_{x_3} - L_{x_1} L_{x_2} L'_{x_3} \right), \\ &= (k_{x_1} + k_{x_2} + k_{x_3}) \frac{a}{3} e^{-bt} \check{V}, \\ &= ae^{-bt} \check{V}, \end{aligned} $$

where we used $k_{x_1} + k_{x_2} + k_{x_3} = 3$, (18) and (15). Since the initial condition is $\check{V}(0) = 0$, we conclude $\check{V}(t) = 0$, and hence $\hat{V}(t) = V(t)$.

Using the ellipsoid volume formula and the fact that $\hat{V} = V = \frac{4\pi}{3} L_{x_1} L_{x_2} L_{x_3}$, we express the error as

$$ \check{L}_{x_3} := \hat{L}_{x_3} - L_{x_3} = \frac{3}{4\pi} \frac{\hat{V}}{\hat{L}_{x_1} \hat{L}_{x_2}} - \frac{3}{4\pi} \frac{V}{L_{x_1} L_{x_2}} = L_{x_1} L_{x_2} L_{x_3} \left( \frac{1}{\hat{L}_{x_1} \hat{L}_{x_2}} - \frac{1}{L_{x_1} L_{x_2}} \right). $$

Further algebraic expansions give

$$ \begin{aligned} \check{L}_{x_3} &= -\frac{L_{x_3}}{\hat{L}_{x_1} \hat{L}_{x_2}} (\hat{L}_{x_1} \hat{L}_{x_2} - L_{x_1} L_{x_2}) \\ &= -\frac{1}{2} \frac{L_{x_3}}{\hat{L}_{x_1} \hat{L}_{x_2}} (\hat{L}_{x_2} + L_{x_2}) (\hat{L}_{x_1} - L_{x_1}) - \frac{1}{2} \frac{L_{x_3}}{\hat{L}_{x_1} \hat{L}_{x_2}} (\hat{L}_{x_1} + L_{x_1}) (\hat{L}_{x_2} - L_{x_2}). \end{aligned} $$

Substituting $\hat{L}_{\alpha} = \check{L}_{\alpha} + L_{\alpha}$ and using the following bounds

$$ |L_{\alpha}(0)| \le |L_{\alpha}(t)| \le |L_{\alpha}(0)| e^{\frac{k_{\alpha}}{3} \frac{a}{b}}, \quad \alpha = x_1, x_2, x_3, \qquad (28) $$

$$ 0 \le |\check{L}_{\alpha}(t)| \le |\check{L}_{\alpha}(0)| e^{\frac{k_{\alpha}}{3} \frac{a}{b}}, \quad \alpha = x_1, x_2, \qquad (29) $$

we can estimate $|\check{L}_{x_3}(t)|$ in terms of $|\check{L}_{x_1}(t)|$ and $|\check{L}_{x_2}(t)|$. Combined with the exponential decay of $\check{L}_{x_1}$ and $\check{L}_{x_2}$ from (25), one can show that there exist constants $C_1, C_2 > 0$ such that

$$ |\check{L}_{x_3}(t)| \le C_1 e^{\frac{k_{x_1}}{3} \frac{a}{b} (1 - e^{-bt}) - \gamma_{x_1} t} + C_2 e^{\frac{k_{x_2}}{3} \frac{a}{b} (1 - e^{-bt}) - \gamma_{x_2} t}, \qquad (30) $$

which proves exponential convergence of $\check{L}_{x_3}$ to zero as $t \to \infty$.

### Step 3: Convergence of $\bar{\tau}^{P2Q}$

We now prove that $\|\bar{\tau}^{P2Q}(t, \hat{L}_{\alpha}(t), \cdot) - \bar{\tau}^{P2Q}(t, L_{\alpha}(t), \cdot)\|_{L^2([0,1])} \to 0$ where $\alpha = x_1, x_2, x_3$, as $t \to \infty$. For the remainder of the proof, we omit the dependence to $L_{\alpha}$ and $\hat{L}_{\alpha}$ in the rates of transition from proliferating to quiescent cells meaning that
we will denote $\bar{\tau}^{P2Q}(t,r) := \bar{\tau}^{P2Q}(t, L_{\alpha}(t), r)$ and $\hat{\tau}^{P2Q}(t,r) := \bar{\tau}^{P2Q}(t, \hat{L}_{\alpha}(t), r)$, $\alpha = x_1, x_2, x_3$.
For brevity, we define

$$ \beta(t,r) := \exp\left(\sqrt{L_{x_1}^2(t) + L_{x_2}^2(t) + L_{x_3}^2(t)} \frac{1-r}{s} - \frac{d}{s}\right), \quad (31) $$

and similarly

$$ \hat{\beta}(t,r) := \exp\left(\sqrt{\hat{L}_{x_1}^2(t) + \hat{L}_{x_2}^2(t) + \hat{L}_{x_3}^2(t)} \frac{1-r}{s} - \frac{d}{s}\right). \quad (32) $$

By the definitions of $\bar{\tau}^{P2Q}$ and $\hat{\tau}^{P2Q}$ (we recall that by hypothesis $\tau^0 > \tau^1 \ge 0$), we have

$$ \begin{aligned} \|\hat{\tau}^{P2Q}(t, \cdot) - \bar{\tau}^{P2Q}(t, \cdot)\|_{L^2([0,1[)} &= (\tau^0 - \tau^1) \left\| \frac{1}{1 + \hat{\beta}(t, \cdot)} - \frac{1}{1 + \beta(t, \cdot)} \right\|_{L^2([0,1[)} \\ &= (\tau^0 - \tau^1) \left\| \frac{\beta(t, \cdot) - \hat{\beta}(t, \cdot)}{(1 + \hat{\beta}(t, \cdot))(1 + \beta(t, \cdot))} \right\|_{L^2([0,1[)} \\ &\le (\tau^0 - \tau^1) \|\beta(t, \cdot) - \hat{\beta}(t, \cdot)\|, \end{aligned} $$

as for all $t, \sup_{r \in [0,1]} [\beta(t,r) \ge 0$ and $\sup_{r \in [0,1]} [\hat{\beta}(t,r) \ge 0$. To bound $\|\beta(t, \cdot) - \hat{\beta}(t, \cdot)\|$, we use the following mean value theorem to the exponential function : for any $p_0, p_2 \in \mathbb{R}$,

$$ |\exp(p_2) - \exp(p_0)| \le |p_2 - p_0| \exp(\max(p_0, p_2)). $$

For $t \ge 0, r \in ]0, 1[$, this implies that

$$ |\beta(t,r) - \hat{\beta}(t,r)| \le \frac{|1-r|}{|s|} \left| \sqrt{L_{x_1}^2(t) + L_{x_2}^2(t) + L_{x_3}^2(t)} - \sqrt{\hat{L}_{x_1}^2(t) + \hat{L}_{x_2}^2(t) + \hat{L}_{x_3}^2(t)} \right| \\ \exp \left( \frac{1-r}{s} \max(\sqrt{L_{x_1}^2(t) + L_{x_2}^2(t) + L_{x_3}^2(t)}, \sqrt{\hat{L}_{x_1}^2(t) + \hat{L}_{x_2}^2(t) + \hat{L}_{x_3}^2(t)}) \right). $$

Using $0 \le 1 - r \le 1$ and the growth of the exponential function, one get

$$ |\beta(t,r) - \hat{\beta}(t,r)| \le \frac{1}{|s|} \left| \sqrt{L_{x_1}^2(t) + L_{x_2}^2(t) + L_{x_3}^2(t)} - \sqrt{\hat{L}_{x_1}^2(t) + \hat{L}_{x_2}^2(t) + \hat{L}_{x_3}^2(t)} \right| \\ \exp \left( \frac{1}{s} \max(\sqrt{L_{x_1}^2(t) + L_{x_2}^2(t) + L_{x_3}^2(t)}, \sqrt{\hat{L}_{x_1}^2(t) + \hat{L}_{x_2}^2(t) + \hat{L}_{x_3}^2(t)}) \right). $$

As the right term does not depend on $r$, we get the same result for $\|\beta(t, \cdot) - \hat{\beta}(t, \cdot)\|$. We have previously shown that $(L_{\alpha})_{\alpha=x_1,x_2,x_3}$, and $(\hat{L}_{\alpha})_{\alpha=x_1,x_2}$ are bounded on $[0, \infty)$ (see (28) and (29)). This directly implies that $(\hat{L}_{\alpha})_{\alpha=x_1,x_2}$ are bounded on $[0, \infty)$ and using the same argument as in Step 2, $\hat{L}_{x_3}$ is also bounded. Thus, there exists a constant $C^{st} > 0$ such that

$$ \|\hat{\tau}^{P2Q}(t, \cdot) - \bar{\tau}^{P2Q}(t, \cdot)\|_{L^2([0,1[)} \le C^{st} \left| \sqrt{L_{x_1}^2(t) + L_{x_2}^2(t) + L_{x_3}^2(t)} - \sqrt{\hat{L}_{x_1}^2(t) + \hat{L}_{x_2}^2(t) + \hat{L}_{x_3}^2(t)} \right| \quad (33) $$

Then by reverse triangle inequality, we have

$$ \|\hat{\tau}^{P2Q}(t, \cdot) - \bar{\tau}^{P2Q}(t, \cdot)\|_{L^2([0,1[)} \le C^{st} \sqrt{(\hat{L}_{x_1} - L_{x_1})^2(t) + (\hat{L}_{x_2} - L_{x_2})^2(t) + (\hat{L}_{x_3} - L_{x_3})^2(t)}. \quad (34) $$
**Case (a):** With the gain defined in (18)-Case (a), the exponential convergence of $\tilde{\mathcal{L}}_{x_1}, \tilde{\mathcal{L}}_{x_2}, \tilde{\mathcal{L}}_{x_3}$ (from (25) and (26)) implies that there exists $C^{st} > 0$ such that

$$ \|\hat{\tau}^{P2Q}(t, \cdot) - \bar{\tau}^{P2Q}(t, \cdot)\|_{L^2(]0,1[)} \le C^{st} \exp(-\min(\gamma_{x_1}, \gamma_{x_2}, \gamma_{x_3})t). \quad (35) $$

**Case (b):** For the gain defined in (18)-Case (b), using (25) and (30) in (34), gives that there exists $C^{st} > 0$ such that

$$ \|\hat{\tau}^{P2Q}(t, \cdot) - \bar{\tau}^{P2Q}(t, \cdot)\|_{L^2(]0,1[)} \le C^{st} \exp(-\min(\gamma_{x_1}, \gamma_{x_2})t). \quad (36) $$

### Step 4: Convergence of $\tilde{\mathcal{P}}$

The error dynamics for $\tilde{\mathcal{P}}$ satisfy for all $r \in ]0, 1[$,

$$ \partial_t \tilde{\mathcal{P}}(t, r) = \partial_t \hat{\mathcal{P}}(t, r) - \partial_t \bar{\mathcal{P}}(t, r) = -(\tau^G(t) + \bar{\tau}^{P2Q}(t, r))\tilde{\mathcal{P}} + (\bar{\tau}^{P2Q}(t, r) - \hat{\tau}^{P2Q}(t, r))\hat{\mathcal{P}}(t, r). $$

The integration of this equation yields

$$ \tilde{\mathcal{P}}(t, r) = (\tilde{\mathcal{P}}(0, r) - \mathcal{Q}(t, r))\mathcal{G}(t, r), \quad (37) $$

where

$$ G(t, r) := \exp \left( - \int_{0}^{t} \tau^{G}(s) + \bar{\tau}^{P2Q}(s, r) ds \right), \quad (38) $$

and

$$ \begin{aligned} Q(t, r) &:= \int_{0}^{t} \exp \left( \int_{0}^{s} \tau^{G}(u) + \bar{\tau}^{P2Q}(u, r) du \right) (\hat{\tau}^{P2Q}(s, r) - \bar{\tau}^{P2Q}(s, r)) \hat{\mathcal{P}}(s, r) ds \\ &=: \int_{0}^{t} q(s, r) ds. \end{aligned} \quad (39) $$

We will show that for all $r \in ]0, 1[$, $G(t, r)$ converges to 0 when $t \to \infty$ and that $Q(t, r)$ remains bounded. From the definition of $\tau^G$ and $\bar{\tau}^{P2Q}$, we have for all $r \in ]0, 1[$,

$$ \tau^G(t) + \bar{\tau}^{P2Q}(t, r) = ae^{-bt} + \tau^0 - \frac{\tau^0 - \tau^1}{1 + \beta(t, r)} \ge ae^{-bt} + \tau^1, \text{ with } \beta \text{ defined by (31).} $$

By the comparison theorem, we have for all $r \in ]0, 1[$,

$$ \lim_{t \to +\infty} \int_0^t \tau^G(s) + \bar{\tau}^{P2Q}(s, r) ds \ge \lim_{t \to +\infty} \int_0^t ae^{-bs} + \tau^1 ds = \lim_{t \to +\infty} \frac{a}{b}(1 - e^{-bt}) + \tau^1 t = +\infty, $$

which implies

$$ \lim_{t \to +\infty} G(t, r) = \lim_{t \to +\infty} \exp \left( - \int_{0}^{t} \tau^{G}(s) + \bar{\tau}^{P2Q}(s, r) ds \right) = 0. \quad (40) $$

Next, we study $Q$. We first show that $\hat{\mathcal{P}}$ is bounded. From the dynamics of $\hat{\mathcal{P}}$ given (15), we get

$$ \begin{aligned} \hat{\mathcal{P}}(t, r) = \exp \left( - \int_{0}^{t} \tau^{G}(s) + \hat{\tau}^{P2Q}(s, r) \right) \hat{\mathcal{P}}(0, r) \\ + \exp \left( - \int_{0}^{t} \tau^{G}(s) + \hat{\tau}^{P2Q}(s, r) \right) \int_{0}^{t} \left[ \tau^{G}(s) \exp \left( \int_{0}^{s} \tau^{G}(u) + \hat{\tau}^{P2Q}(u, r) \right) \right]. \end{aligned} \quad (41) $$
From the definition of $\tau^G$ and $\hat{\tau}^{P2Q}$, we can prove that for all $r \in]0, 1[, \text{ for all } t \ge 0$,

$$0 \le \tau^G(t) + \hat{\tau}^{P2Q}(t, r) = ae^{-bt} + \tau^0 - \frac{\tau^0 - \tau^1}{1 + \hat{\beta}(t, r)} \le ae^{-bt} + \tau^0, \text{ where } \hat{\beta} \text{ is given by (31).} \quad (42)$$

Using the lower bound, one can deduce that for all $r \in]0, 1[$ and for all $t \ge 0$,

$$\int_0^t \tau^G(s) + \hat{\tau}^{P2Q}(s, r) \ge 0 \implies \exp\left(-\int_0^t \tau^G(s) + \hat{\tau}^{P2Q}(s, r)\right) \le 1,$$

thus, the first term of (41) is bounded since $\hat{P}(0, \cdot)$ is also bounded. For the second term, using the upper bound in (42), we have for all $r \in]0, 1[$ and for all $t \ge 0$,

$$\int_0^t \tau^G(s) + \hat{\tau}^{P2Q}(s, r) \le \int_0^t ae^{-bs} + \tau^0 = \frac{a}{b}(1 - e^{-bt}) + \tau^0 t. \quad (43)$$

Hence, for all $r \in]0, 1[$ and for all $t \ge 0$,

$$\tau^G(t) \exp\left(\int_0^t \tau^G(s) + \hat{\tau}^{P2Q}(s, r)\right) \le ae^{\frac{a}{b}} \exp\left(-e^{-bt} + \tau^0 t - bt\right) \le ae^{\frac{a}{b}} e^{-(b-\tau^0)t},$$

which leads to

$$\int_0^t \tau^G(s) \exp \left( \int_0^s \tau^G(u) + \hat{\tau}^{P2Q}(u, r) \right) \le ae^{\frac{a}{b}} \int_0^t e^{-(b-\tau^0)s} = \frac{a}{\tau^0 - b} e^{\frac{a}{b}} (1 - e^{-(b-\tau^0)t}).$$

Thus, there exists $C^{st} > 0$ such that

$$\|\hat{P}(t, \cdot)\|_{L^2([0,1])} \le C^{st}(1 - e^{-(b-\tau^0)t}). \quad (44)$$

We now focus on $q$ defined in (39) by:

$$q(t, r) = \exp\left(\int_0^t \tau^G(s) + \bar{\tau}^{P2Q}(s, r)\right) \left(\hat{\tau}^{P2Q}(t, r) - \bar{\tau}^{P2Q}(t, r)\right) \hat{P}(t, r), \quad t \ge 0, r \in]0, 1[.$$

From the definition of $\tau^G$ and $\bar{\tau}^{P2Q}$, we have

$$\int_0^t \tau^G(s) + \tau^{P2Q}(s, r) \le \int_0^t (ae^{-bs} + \tau^0) = \frac{a}{b}(1 - e^{-bt}) + \tau^0 t \le \frac{a}{b} + \tau^0 t.$$

Using this result and (44), we obtain that there exists $C^{st} > 0$ such that

$$\|q(t, \cdot)\|_{L^2([0,1])} \le C^{st} \|\hat{\tau}^{P2Q} - \bar{\tau}^{P2Q}\|_{L^2([0,1])} (e^{\tau^0 t} + e^{-(b-2\tau^0)t}). \quad (45)$$

**Case (a):** Using (35) in (45), we obtain that there exists a positive constant $C^{st}$ such that

$$\|q(t, \cdot)\|_{L^2([0,1])} \le C^{st} (e^{-(\min(\gamma_x, \gamma_y, \gamma_z) - \tau^0)t} + e^{-(\min(\gamma_x, \gamma_y, \gamma_z) + b - 2\tau^0)t}). \quad (46)$$

Integrating in time, it follows for all $t \ge 0$,

$$\|Q(t, \cdot)\|_{L^2([0,1])} \le \int_0^t \|q(t, \cdot)\|_{L^2([0,1])} \\ \le C^{st} \left( \frac{1 - e^{-(\min(\gamma_x, \gamma_y, \gamma_z) - \tau^0)t}}{\min(\gamma_x, \gamma_y, \gamma_z) - \tau^0} + \frac{1 - e^{-(\min(\gamma_x, \gamma_y, \gamma_z) + b - 2\tau^0)t}}{\min(\gamma_x, \gamma_y, \gamma_z) + b - 2\tau^0} \right).$$
Since $min(\gamma_x, \gamma_y, \gamma_z) - \tau^0 > 0$ and $min(\gamma_x, \gamma_y, \gamma_z) + b - 2\tau^0 > 0$, both exponents are negative which proves the desired result.

**Case (b):** Using (36) in (45), we obtain that there exists a positive constant $C^{st}$ such that

$$ \|q(t, \cdot)\|_{L^2(]0,1[)} \le C^{st}(e^{-(\min(\gamma_x, \gamma_y)-\tau^0)t} + e^{-(\min(\gamma_x, \gamma_y)+b-2\tau^0)t}). \quad (47) $$

Since $min(\gamma_x, \gamma_y) - \tau^0 > 0$ and $min(\gamma_x, \gamma_y) + b - 2\tau^0 > 0$, again both exponents are negative, hence we obtain that the integral of $q$ does not diverge and implies the convergence of $\tilde{P}$ to zero. $\square$

<figure><img src="image_11.png" /></figure>