# Spatial mechanistic modeling for prediction of 3D multicellular spheroids behavior upon exposure to high intensity pulsed electric fields

Annabelle Collin, Hadrien Bruhier, Jelena Kolosnjaj, Muriel Golzio, Marie-Pierre Rols, Clair Poignard

► To cite this version:

Annabelle Collin, Hadrien Bruhier, Jelena Kolosnjaj, Muriel Golzio, Marie-Pierre Rols, et al. Spatial mechanistic modeling for prediction of 3D multicellular spheroids behavior upon exposure to high intensity pulsed electric fields. AIMS bioengineering, 2022, 9 (2), pp.102-122. 10.3934/bioeng.2022009. hal-03700689

HAL Id: hal-03700689

https://inria.hal.science/hal-03700689

Submitted on 21 Jun 2022

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 Attribution 4.0 International License
Spatial mechanistic modeling for prediction of 3D multicellular
spheroids behavior upon exposure to high intensity pulsed
electric fields

Annabelle Collin¹,*, Hadrien Bruhier¹, Jelena Kolosnjaj², Muriel Golzio², Marie-Pierre Rols², and Clair Poignard¹

¹Univ. Bordeaux, CNRS, INRIA, Bordeaux INP, IMB, UMR 5251, F-33400 Talence,
France

²Institut de Pharmacologie et de Biologie Structurale, Université de Toulouse, CNRS,
UPS, 31077 Toulouse, France

*Corresponding author

June 21, 2022

Abstract

The objective of this work was to investigate the growth specificities of cancer cells spheroids subjected to pulsed electric field. Multicellular HCT-116-GFP spheroids were exposed to different electric field intensities and the volume of multicellular spheroids was monitored by fluorescence and bright field microscopy. Thanks to an advanced mathematical model, based on differential equations and well-adapted estimation strategies, our modeling enables us to characterize the multicellular spheroids growth after permeabilizing pulsed electric field. In particular, we identify the percentage of cells which are destroyed and the percentage of cells which exhibit an altered growth pattern for different magnitudes of the electric field. We also quantify the growth resumption upon reversible and partially irreversible electroporation. Our preliminary results provide a first quantification of the impact of electroporation on multicellular spheroids growth, and suggest a booming growth of partially irreversible electric pulses, leading to an accelerated regrowth.

**Keywords:** PDE modeling; parameters estimation; electroporation; spheroids

# 1 Introduction

Electroporation or electropermeabilization is a physical phenomenon that refers to the oc-
currence of permeant structures within the cell membranes upon cells exposure to well defined
pulses of electric fields. Depending on the pulse parameters, such as the intensity of the elec-
tric field, their duration and the number (and the frequency) of the train of pulses, defects
of different size may be induced in cell membrane, which make it permeable to non permeant
molecules. This permeabilization can be either reversible (does not affect the cell viability) or
irreversible (meaning that it can lead to cell death). In reversible electroporation, cell mem-
branes reseal after several minutes and the cells continue their lifecycle. This process has led to
a combined cancer therapy, referred to as electrochemotherapy, where electric pulses are used to
permeabilize the cell and highly improve the local penetration of hydrophilic anticancer drugs
[14]. Irreversible electroporation is the biophysical phenomenon that leads to cell death. This
process is currently under development in clinical oncology for nonsurgical and non thermal
tumor ablation [16]. It is worth noting that electric field induced membrane permeabilization
is just one among many other phenomena [10], which can occur upon cells/tissue exposure to
high electric field pulses.

While irreversible electroporation is efficient in most of the cases, with a complete ablation
rate of about 75% of treated tumor nodules after one irreversible electroporation procedure [3],
it is well-known that for some patients the ablation is unsuccessful and event may lead to an
accelerate growth of the tumor. In order to improve our understanding on the impact of the
pulsed electric fields on the tumor response, we exploited multicellular tumor spheroids as an
*in vitro* micro-tumor model. Multicellular tumor spheroids are three-dimensional structures
composed of cancer cells that can accurately reproduce the behavior of small solid tumors in
their preliminary avascular stage. Multicellular spheroids are becoming an essential tool in
cancer research as they provide an intermediate complexity between 2D monolayer cell models
and *in vivo*-grown solid tumors. Spheroids closely resemble solid micro-tumors in many aspects,
such as the heterogeneous architecture, internal gradients of signaling factors, nutrients, and
oxygenation. Moreover, the cells within the multicellular spheroids interact between each-other
via cell-to cell interactions and their growth kinetics are similar to those observed in tumors *in
vivo*. These similarities provide great potential for studying the biological properties of tumors,
thus they are often used for drug screening as well as used to determine the therapeutic efficacy.
Such models already demonstrated their relevance to address the electroporation phenomena [8,
7].

Standard clinical protocols of irreversible electroporation use about one hundred electric train of monopolar pulses of 100 µs duration of 1500 V/cm of nominal amplitude at the frequency of 1 Hz. Our objective here was to quantify, thanks to mathematical modeling, the impact of different pulsed electric field of different nominal magnitudes on the growth of the spheroids. To mimick the pulses of the clinical protocols of IRE ablation, the train of pulses consist of 80 monopolar pulses of 100 µs duration at 1 Hz with an amplitude ranging from 500 to 2000 V·cm⁻¹. The aim was to characterize the multicellular spheroids regrowth after the pulse delivery with respect to the pulse magnitude. Throughout the paper the indicated field magnitudes are the *nominal magnitudes*, that is the voltage between the electrodes divided by their distance. It is worth noting it is not the electric field seen by the cells.

Thanks to mathematical models and calibration strategy, we fit the model with the exper-
imental data to assess global behaviors among the spheroid population. Two approaches were
considered. The first one consisted in studying only the resumption of the growth after the
exposure to pulsed electric field with the Gompertz model, which is a classical ODE model for
tumor growth [2, 4, 9, 12, 13, 17]. Gompertz models are well-known to properly track the
free growth of multicellular spheroids and tumors. The model consisted in a nonlinear ordinary
differential equation written on the volume, where the multicellular spheroids are considered
as a homogeneous entity. The second model is a more general model and is based on partial
differential equations representing the evolution of populations of tumor cells controlled by a
“carrying capacity”. Three populations of tumor cells have been considered: the proliferative
cells, the quiescent cells and the cells whose functioning is altered by the effect of the electric
pulse. This more complex approach involves additional parameters but enables representing the
heterogeneity of the response of the tumor cells within the multicellular spheroids.

The models were confronted to biological data corresponding to the volume evolution of
multicellular spheroids submitted to electric pulses with different intensities. To fit the bio-
logical data with our models, we considered an estimation strategy of the parameters based
on a population approach allowing to compensate for sparse sampling times and measurement
uncertainties by constraining the variability of the parameters in the population.
The manuscript is divided as follows: in Section 2, the data containing various cohorts depending on the intensity of the electric field are presented. In Section 3, the mathematical models and the strategy to estimate the parameters using the biological data are introduced. Section 4 gives all the obtained results which are then discussed in Section 5.

# 2 Materials and method

## 2.1 Biological experiments

Biological experiments were conducted at the Institute of Pharmacology and structural biology (IPBS), Toulouse, France, and involved a human colorectal carcinoma cell line (HCT116) (ATCC® CCL-247TM), expressing the green fluorescent protein (GFP) (HCT116-GFP). Cells were grown under standard conditions (5% CO₂, 37 °C) in the Dulbecco's Modified Eagle Medium (DMEM, Gibco-Invitrogen, Carlsbad, CA, USA) containing 4.5 g/L glucose, L-Glutamine and pyruvate, 1% of penicillin/streptomycin, and 10% of fetal bovine serum. Multicellular cell spheroids were made by seeding 500 cells per well in Costar® Corning® Ultra-low attachment 96-well plates (Fisher Scientific, Illkirch, France). Plates were kept in 5% CO₂ humidified atmosphere at 37 °C. Multicellular spheroids acquired a cohesive structure and were submitted to pulsed electric field 3 days following the seeding. After the electric exposure (detailed below), the spheroids were repositioned in ultra-low attachment plates, in which they were kept for 10 days. Fresh medium was added to the spheroids at a 3 to 4 days interval. Multicellular spheroids growth was monitored by fluorescence and bright field videomicroscopy using the IncuCyte Live Cell Analysis System Microscope (Essen BioScience IncuCyte™, Herts, Welwyn Garden City, UK) at a magnification multiplied by × 10.

Electroporation was achieved following the delivery of 80 unipolar pulses at 0 V·cm⁻¹ (control) or 500 V·cm⁻¹ or 1000 V·cm⁻¹ or 2000 V·cm⁻¹, the duration of a pulse was 100 µs, and the pulses were applied at a frequency of 1 Hz. The pulses were applied to multicellular spheroids in a low conduction pulsing medium (see [5, 18]) with a parallel plate stainless steel electrode (the plate distance was 0.4 cm) connected to the Electrocell S 20 generator (Leroy Biotech, St Orens, France), following which the spheroids were rinsed in phosphate saline buffer and transferred to the cell growth medium. Once again let us emphasize that the values 500 V·cm⁻¹ or 1000 V·cm⁻¹ or 2000 V·cm⁻¹ corresponds to the amplitudes of the *nominal electric field*. Presented results were obtained in two independent experiments. For each experiment, 6 to 8 spheroids were pulsed. Two control groups were used and are denoted either as *free growth* or 0 V·cm⁻¹. In free growth condition, the spheroids were kept in growth medium. For 0 V·cm⁻¹ condition, the spheroids were transferred to the low conduction buffer for 10 minutes, following which the spheroids were rinsed in phosphate saline buffer and re-transferred to growth medium.

Images were analyzed with the ImageJ software (U.S. National Institute of Health, Bethesda, MD, USA). The software was used to determine the mean fluorescence intensity of spheroids, and to measure spheroids area in the equatorial plane. The representative micrographs of multicellular spheroids followed by fluorescence microscopy of control (0 V·cm⁻¹), 500 V·cm⁻¹, 1000 V·cm⁻¹ and 2000 V·cm⁻¹ treated spheroids, immediately after the application of electric pulses are shown on Figure 1.
<figure><img src="image_2.png" /><figcaption>Figure 1: Bright field micrographs of multicellular spheroids followed by fluorescence microscopy of control (0 V·cm⁻¹), 500 V·cm⁻¹, 1000 V·cm⁻¹ and 2000 V·cm⁻¹-treated spheroids, after the application of electric pulses (N = 80, t = 100 µs, ν = 1 Hz). Green color: green fluorescing protein in living cells and red color: propidium iodide. Top: immediately after pulses application. Bottom: 1 h, 4 h, 1 day and 10 days after pulses application.</figcaption></figure>

The green color reflects the green fluorescing protein in living cells and the red color corresponds to propidium iodide, a red fluorescent probe that penetrates the cells when permeabilization occurs on cell membrane. Note that at 0 V·cm⁻¹ no red fluorescence is visible, while at higher electric field the degree of red color increases, as increasingly more propidium iodide penetrates the cells illustrating the cell permeabilization. Figure 1-Top (resp. Bottom) shows the spheroids immediately (resp. 1 h, 4 h, 1 day and 10 days) after the shock. Note that 500 V·cm⁻¹ pulses do not affect cell viability (as indicated by green fluorescence and the growth of the cellular mass), while at 1000 V·cm⁻¹ we have a transient disappearance of green fluorescence, and at 2000 V·cm⁻¹ we have a permanent disappearance of green fluorescence. Figure 2 shows the representative micrographs of treated spheroids (1000 V·cm⁻¹), that illustrates the localization of spheroid regrowth following the pulse delivery. At 5 h following electroporation the green fluorescence generally begins to fade (complete extinction of fluorescence is only observed in one spheroid). Twenty four hours after treatment, the fluorescence is either extinct, or still observed inside the spheroid or in the periphery.
<figure><img src="image_3.png" /></figure>

Figure 2: Localization of spheroid regrowth following electroporation at 1000 V·cm⁻¹. Top: micrographs of representative spheroids 5 h post treatment. The green fluorescence generally begins to fade (complete extinction of fluorescence is only observed in one spheroid-top right panel). Middle: same spheroids as in top panels at 24 h after treatment. The fluorescence is either still observed inside the spheroid or in the periphery (first and second middle micrographs, respectively, or is completely extinct). Bottom: micrographs of same spheroids as on top panels, at 7 days after treatment, all exhibiting green fluorescence that indicate spheroids viability and growth. The regrown zones can be either attached to the treated spheroids or scattered in space.

## 2.2 Data presentation

Figure 3 shows the evolution of the volume of the multicellular spheroids in control (Left panel) and electroporation (Right panel) configurations. In total 83 multicellular spheroids were tracked over 250 hours. Measurements, as plotted on Figure 3 were obtained from fluorescence micrographs as described in Section 2.1.

The cohort is divided as follows:

* 24 multicellular spheroids correspond to the control spheroids (free growth). We exclude the spheroids 3 and 24 as dynamics strongly differs from the others.

* 59 multicellular spheroids correspond to electroporated spheroids divided into 4 categories depending on the electric field strength.

Different groups are denoted EF0, EF500, EF1000 and EF2000, respectively, which corresponds to measures obtained in groups exposed to 0 V·cm⁻¹, 500 V·cm⁻¹, 1000 V·cm⁻¹ and 2000 V·cm⁻¹, respectively.
<figure><img src="image_4.png" /><figcaption>Figure 3: Volume evolutions of multicellular spheroids. Left-Control experiment (24 cases). Right-Electroporation experiment (59 cases divided into 4 categories: 13 for EF0, 16 for EF500, 14 for EF1000 and 16 for EF2000).</figcaption></figure>

## 2.3 Elimination of the outliers

The cohort EF0 contains 13 multicellular spheroids. The cohort EF500 contains 16 multicellular spheroids among which 4 have been excluded: 2 because they did not see the electric field (spheroids 18 and 29) and 2 because their dynamics strongly differed from the others (spheroids 16 and 17). The cohort EF1000 contains 14 multicellular spheroids among which 3 have been excluded: 1 because there was no data after $t > 90$ h (spheroid 43), 1 because the time 90 h (which is a fundamental time for the estimation process) was missing (spheroid 37) and 1 because its dynamics strongly differed from the others (spheroid 38). The cohort EF2000 contains 16 multicellular spheroids among which 2 have been excluded: 1 because it did not see the electric field (spheroid 59) and 1 because its dynamics strongly differed from the others (spheroid 47).

## 2.4 Division of the cohort EF1000 into 2 subgroups

Data presented in Figure 4, indicates that according to multicellular spheroids response to the electric field, the cohort EF1000 could be split into 2 subgroups. As two kinds of behaviors appeared in cohort EF1000, we study them separately:

* Group A (6 spheroids): the multicellular spheroids were very damaged by the electric shock, the volume of green fluorescing cells reached almost 0 mm³ after 48 hours,

* Group B (5 spheroids): the dynamics of the multicellular spheroids were very similar to the dynamics observed in cohort EF500.

Precisely, multicellular spheroids with a volume at 43h below (resp. above) 0.01 mm³ could be grouped into Group A (resp. Group B), see Figure 4.
<figure><img src="image_5.png" /><figcaption>Figure 4: Volumes after exposure to electric pulses (t = 43 h) for cohort EF1000 (logarithmic scale).</figcaption></figure>

# 3 Methods

## 3.1 Modeling of volume evolution

**Gompertz model for free growing volume** Gompertz model is a well-known ordinary differential equation (ODE) for tumor growth, which fits well with the free growth of *in vitro* and *in vivo* tumors. The reader may refer to [2] for a comparison of different ODE models on *in vivo* data and [9, 4, 12, 13, 17, 6] for validations. The tumor volume arising from the Gompertz equation reads

$$V(t) = V_0 e^{\frac{a}{b}(1 - e^{-bt})}, \qquad (1)$$

where $V_0$ corresponds to the initial volume at time $t = 0$, $b^{-1}$ corresponds to the characteristic time at which the tumor growth capacity decreases and $a$ is a constant corresponding to the initial growth rate.

**Population approach to fit the data** To estimate the three parameters $V_0$, $a$ and $b$, a population approach has been chosen. The mixed-effect approach [19] consists of pooling all the multicellular spheroids together and estimating a global distribution of the model parameters in the population. More precisely, the individual parameters $a^j$ and $b^j$ (where $j$ denotes the individual) are assumed to be realizations of a random variable decomposed into two parts:

$$\{V_0, a, b\}^j = \{V_0, a, b\}^{pop} + \{\tilde{V}_0, \tilde{a}, \tilde{b}\}^j, \forall j \in [1, N_s],$$

where $N_s$ is the number of considered spheroids, $V_0^{pop}$, $a^{pop}$ and $b^{pop}$ correspond to the fixed effects and $\tilde{V}_0^j$, $\tilde{a}^j$ and $\tilde{b}^j$ correspond to the random effects and have been assumed with mean zero. The parameters that we want to estimate are not physiological and there are no priors for their means and standard deviations. Then, we assume that

$$a^{pop}, b^{pop} \sim \mathcal{N}(0, 0.1^2) \text{ and } V_0^{pop} \sim \mathcal{N}(0, 0.05^2),$$

emphasizing that these assumptions do not seem too restrictive. The estimation strategy is not sensitive to these priors. We then assume that the measurement error for all volumes is proportional

$$V^{j,MRI}(t) = V^{j,estim}(t)(1 + e),$$
where $V_k^{j,estim}$ corresponds to the estimated volume using ODE systems and e the measurement error. We assume that it follows a Gaussian law $N(0, 0.1^2)$. The standard deviation is estimated.

This estimation strategy is implemented in the software Monolix [1], which maximizes the likelihood using the stochastic approximation expectation maximization (SAEM) algorithm [11]. For all the results obtained with Monolix in this paper, the convergence of the algorithm has been reached.

The population estimation strategy for Gompertz model will be applied on the 22 multicellular spheroids of the control experiment and the 13 multicellular spheroids EF0 of the electroporation experiment in order to verify that the cohorts are similar. The objective is to validate the assumption that there is no impact of the electrodes on the growth.

A correlation between the parameters a and b in a population of tumors has been observed frequently in Gompertz model. In [21], a reduced Gompertz model is introduced consisting in estimating $V_0$ and b as previously and $k = \frac{a}{b}$ as a constant parameter in the population i.e.

$$ \{V_0, b\}^j = \{V_0, b\}^{pop} + \{\tilde{V}_0, \tilde{b}\}^j, \forall j \in [1, N_s], \\ \{k\}^j = \{k\}^{pop}, \forall j \in [1, N_s]. $$

The use of the reduced Gompertz model improves the stability of the estimation process and will be used in what follows.

There were fewer measurements in the electroporation experiment than for the control experiment. It is therefore necessary to check whether the SAEM algorithm works correctly with a low number of spheroids and a low number of measurements. To do this, we create a cohort with a certain number of individuals denoted by y, drawn at random from the control group. Then, we randomly select x consecutive times to estimate the parameters. This operation is performed 50 times in order to minimize errors caused by chance. The results are given in Appendix.

**Gompertz model for multicellular spheroids after electroporation pulses** In this paragraph, we want to study the growth of multicellular spheroids after the electric shock in order to answer the following questions.

* After the electric shock, is there a resumption of growth?

* Is it still a Gompertz growth?

* What are the impacts on the parameters of the model?

First of all, for most of the spheroids which have been electroporated (excepted of spheroids for cohort EF2000), a resumption of the growth occurs, see Figure 3. Letting $t_{init}$, the time of resumption and $t_{last}$ the time of the last measurement, we rewrite the reduced Gompertz model in an inverted form

$$ V(t) = V_{last}e^{ke^{-bt_{last}(1-e^{b(t_{last}-t)})}}. $$

For the 3 cohorts EF500, EF1000-A and EF1000-B, we will estimate the parameters $V_{last}$, b and k for each fixed $t_{last}$ using the population approach presented previously (in which $V_0$ plays the role of $V_{last}$). The parameter $t_{last}$ will be estimated by considering the minimal value given reasonable errors.

This approach based only on the volume evolution will show some limitations motivating us to consider 3D evolution of the multicellular spheroids.
## 3.2 Modeling the 3D evolution of spheroid

In this section, we complexify the growth model to introduce a heterogeneity in the cell population. The model has been designed to describe the 3D evolution of the spheroids with two constraints: (1) to be compatible with the Gompertz model at the volume scale, (2) to depend on few parameters which are identifiable with the available data.

**General equations** We will consider 3 compartments of cells: the proliferative ones whose the density denoted by $P$, the quiescent ones whose the density denoted by $Q$ and the cells with a functioning altered by the impact of the electrical shock whose the density denoted by $F$. By denoting the tumor spheroid domain at time $t$ by $\Omega(t)$, the evolution of the density $P$ is supposed to satisfy the following equation

$$ \partial_t P + \nabla \cdot (\vec{v}P) = f(t, P, Q), \quad \Omega(t), \qquad (2) $$

and the density of $Q$ follows

$$ \partial_t Q + \nabla \cdot (\vec{v}Q) = g(t, P, Q), \quad \Omega(t), \qquad (3) $$

where $\vec{v}$ denotes the velocity field that describes the motion of tumor cells. The evolution of the density $F$ is supposed to satisfy the following equation

$$ \partial_t F + \nabla \cdot (\vec{v}F) = 0, \quad \Omega(t). \qquad (4) $$

**Remark 1.** In the case of free growth (when there is no electroporation), the density $F$ equals to 0 for all $t$ and for all $x$.

Using the saturation hypothesis $P+Q+F = 1$ in $\Omega(t)$, we have $\nabla \cdot v = f(t, P, Q)+g(t, P, Q)$. This implies that the volume defined by

$$ V(t) = \int_{\Omega(t)} dX $$

verifies

$$ \begin{aligned} V'(t) &= \frac{d}{dt} \int_{\Omega(t)} dX = \int_{\partial\Omega(t)} (v \cdot n) dX = \int_{\Omega(t)} \nabla \cdot v dX \\ &= \int_{\Omega(t)} [f(t, P, Q) + g(t, P, Q)] dX. \end{aligned} $$

Our objective is to have the volume following a Gompertz law when there is no electroporation (then $P + Q = 1, \Omega(t)$). One can prove – following the proof given in [6] – that a way to do so consists in imposing

$$ f(t, P, Q) + g(t, P, Q) = \tau_G(t)(P(t, x) + Q(t, x)), \forall t, \qquad (5) $$

where $\tau_G(t) = ae^{-bt}$ corresponds to the growth rate.

**Modeling of the tumor heterogeneity** The function $g$ has to be defined in order to obtain realistic behavior. Assuming that proliferative cells become quiescent cells, we assume that $g(t, P, Q) = \tau_{PtoQ}P$. The function $\tau_{PtoQ}$ depending on $(t, x)$ has to be chosen to describe the appearance of a high quiescent proportion of cells in the center of the spheroid.
**Electroporation modeling** When electroporation pulses occur, three phenomena are considered:

(1) a part of cells is destroyed,

(2) the functioning of a part of the cells is modified i.e. a part of proliferative and quiescent cells changes their functioning,

(3) the initial growth rate is modified.

Denoting by $t_{as}$ the time just after the electrical shock, these three phenomena can be mathematically modeled by:

(1) the death of proliferative and quiescent cells leads to $P(t_{as}, x) = 0$, for $x \in \omega_P \subset \Omega(t < t_{as})$ and $Q(t_{as}, x) = 0$, for $x \in \omega_Q \subset \Omega(t < t_{as})$ then $\Omega(t_{as})$ decreases,

(2) the functioning modification of proliferative and quiescent cells leads to $F(t_{as}, x) = \lambda(P(t_{as}, x) + Q(t_{as}, x))$, for $x \in \Omega(t_{as})$,

(3) the value of the parameter is decomposed into two parts:

$$a = \begin{cases} a_{bs}, & \forall t < t_{as}, \\ ma_{bs}, & \forall t \ge t_{as}, \end{cases} \qquad (6)$$

where $m$ is a constant (always estimated superior to 1) which can be interpreted as a boost of the growth after the electrical shock.

**Radial equation for free growth** Inspiring from the work of [15], the spherical symmetry of spheroid can be used to rewrite the model 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.

Since the volume $V = \frac{4}{3}\pi R^3$ follows a Gompertz law when free growth is considered as (5) is verified then simple calculations shows that the radius of the spheroid $R$ follows also a Gompertz law: $R(t) = R_0 e^{\frac{a}{3b}(1-e^{-bt})}$ with $R_0 = (\frac{3}{4\pi}V_0)^{\frac{1}{3}}$ since it satisfies

$$R'(t) = \frac{\tau_G(t)}{3} R(t). \qquad (7)$$

For all $x \in \Omega(t)$, we denote the normalized radial coordinate ($r = 0$ at the center of the spheroid, $r = 1$ at the surface)

$$r(t, x) = \frac{\|x\|}{R(t)}.$$

Under the hypothesis of invariance by rotation, we can then define, for all $(t, x)$:

$$\tilde{u}(t, r(t, x)) = u(t, x), \text{ for unknowns } P, Q, \tau_{PtoQ} \text{ and } \vec{v}.$$

We denote by $\tilde{v}$ the radial component of the velocity $\tilde{\vec{v}}$ and the invariance by rotation implies that $\tilde{\vec{v}} = \tilde{v}\tilde{e}_r$. As the explicit formula of $R$ is known, it is sufficient to compute the unknown $\tilde{P}$ and then $\tilde{Q}$ can be deduced by $\tilde{Q} = 0$ outside the tumor and $\tilde{Q} = (1 - \tilde{P})$ inside the tumor. One can show, see Appendix for the proof-that $\tilde{P}$ follows for each $r \in [0, R(t)]$, the following ordinary differential equation

$$\tilde{P}' = \tau_G(1 - \tilde{P}) - \tilde{\tau}_{PtoQ}\tilde{P}.$$
To close the system, the evolution of the rate $\tilde{\tau}_{PtoQ}$ has to be described. Seen as representing the lack of oxygen, one can consider the following logistic function

$$ \tilde{\tau}_{PtoQ}(t, r) = \tilde{\tau}_{PtoQ}^{begin} - \frac{\tilde{\tau}_{PtoQ}^{begin} - \tilde{\tau}_{PtoQ}^{end}}{1 + e^{\frac{(R(t)(1-r)-d)}{s}}}. $$

The parameter $d$ corresponds to the distance from the tumor front for which the oxygen is easily accessible. The parameter $s$ corresponds to the slope of the function. The parameters $\tilde{\tau}_{PtoQ}^{begin}$ (resp. $\tilde{\tau}_{PtoQ}^{end}$) corresponds to the value of $\tilde{\tau}_{PtoQ}$ far from (resp. close to) the tumor boundary.

**Remark 2.** When the hypothesis of spherical symmetry is not valid anymore, the equation on the velocity $\nabla \cdot v = f + g$ is not sufficient enough to close the system. For example, in [6], a Darcy law $v = -\nabla\pi$ where $\pi$ is the pressure exerted by tumor cells to healthy cells is considered. Furthermore, the evolution of the rate $\tau_{PtoQ}$ is not valid anymore as it directly depends on $R$ and $r$. One could consider the distance to the tumor front instead of the radius. The study of the stability of the results for almost radial situations is referred to future work.

**Radial equation for electroporation process** When considering electric pulses, the volume $V$ (resp. the radius $R$) does not follow Gompertz laws as we have:

$$ V'(t) = \tau_G(t) \int_{\Omega(t)} (P(t, x) + Q(t, x)) dX. $$

We denote by $I$ the following function of $t$ and $r$

$$ I(t, r) = \int_{0}^{r} (\tilde{P}(t, \underline{r}) + \tilde{Q}(t, \underline{r})) \underline{r}^{2} d\underline{r} $$

Using radial coordinate and the fact that $V = \frac{4}{3\pi}R^3$, one can show that

$$ R'(t) = R(t)\tau_G(t)I(t, 1). $$

Following the proof given in Appendix, one can show that

$$ \partial_t \tilde{P}(t, r(t, x)) + \tau_G(r^{-2}I(t, r) - rI(t, 1))\partial_r \tilde{P}(t, r(t, x)) = \tau_G(1 - \tilde{F})(1 - \tilde{P}) - \tilde{\tau}_{PtoQ} \tilde{P}. \quad (8) $$

and

$$ \partial_t \tilde{F}(t, r(t, x)) + \tau_G(r^{-2}I(t, r) - rI(t, 1))\partial_r \tilde{F}(t, r(t, x)) = -\tau_G(1 - \tilde{F})\tilde{F}, \quad (9) $$

then $\tilde{Q}$ can be determined using the fact that $\tilde{Q} = 1 - (\tilde{P} + \tilde{F})$. Six parameters have to be estimated: $V_0, k, b, p, m$ and $\lambda$ using volume evolution. As previously, the electroporation data are separated in four groups: EF500, EF1000-A, EF1000-B and EF2000. To properly estimate the values of the parameters, we do not consider the 4 times just after pulses: $t = 18, 24, 30, 36$ h. Indeed, due to the osmotic shock in cells constituting the spheroids after the electrical shock, very large values of volumes have been measured. As our model does not integrate this phenomenon, only the values after the rebuilding of the multicellular spheroids ($t > 36$ h) will be used to estimate the parameters.

The software Monolix [1] cannot be used here to estimate the parameters as our model is a PDE system. We then will use a Matlab library called nlmefitsa in which the SAEM algorithm is implemented.
**Considered errors** To validate all our estimation strategy, we will consider the mean squared error (MSE) defined by

$$ \text{MSE} = \frac{1}{N_s} \sum_{s=1}^{N_s} \frac{1}{N_t^s} \sum_{i=1}^{N_t^s} |V_{ds}^i - V(t_s^i)|^2, \quad (10) $$

and the $l^1$-norm defined by

$$ \text{err}_1 = \frac{1}{N_s} \sum_{s=1}^{N_s} \frac{1}{N_t^s} \sum_{i=1}^{N_t^s} |V_{ds}^i - V(t_s^i)|, \quad (11) $$

where $N_t^s$ corresponds to the number of observations for each spheroid $s \in [1, N_s]$, $V_{ds}^i$ corresponds to the data at time $t_s^i$ for each spheroid $s$ and each observation $i \in [1, N_t^s]$ and $V(t_s^i)$ to the predicted volume value at time $t_s^i$.

# 4 Results

## 4.1 Volume fits

**Free growth** Table 1 gives the mean and the standard deviation of the estimated values of $V_0$, $a$ and $b$ for the cohort of control experiment (second column) and for the cohort EF0 of electroporation experiment (third column). As expected, the values are very close meaning that the electrodes do not perturb the growth. Figure 5-Top shows the corresponded fits. The mean and the standard deviation of the estimated values of $V_0$, $a$ and $b$ coupling the 2 cohorts are given in the last column of Table 1. Figure 5-Bottom-Left shows the strong correlation between parameters $a$ and $b$ ($p = 6.02 \times 10^{-9}$) validating the use of the reduced Gompertz model hence the fits using the coupling of the 2 cohorts are given in Figure 5-Bottom-Right and the values of the estimated parameters in Table 2.

<table><caption>Table 1: Estimated parameters of Gompertz model in case of free growth.</caption><thead><tr><th>Parameters</th><th>Control Experiment</th><th>Electroporation experiment EF0</th><th>Cohorts coupling</th></tr></thead><tbody><tr><td>V<sub>0</sub></td><td>0.052 (mean) - 0.0024 (std)</td><td>0.055 (mean) - 0.00069 (std)</td><td>0.050 (mean) - 0.0024 (std)</td></tr><tr><td>a</td><td>0.024 (mean) - 0.0023 (std)</td><td>0.019 (mean) - 0.0017 (std)</td><td>0.023 (mean) - 0.0027 (std)</td></tr><tr><td>b</td><td>0.011 (mean) - 0.00093 (std)</td><td>0.0091 (mean) - 0.00055 (std)</td><td>0.011 (mean) - 0.0012 (std)</td></tr><tr><td>MSE</td><td>0.00016</td><td>0.00011</td><td>0.00014</td></tr><tr><td>err<sub>1</sub></td><td>0.0089</td><td>0.0079</td><td>0.0086</td></tr></tbody></table>

<table><caption>Table 2: Estimated parameters of reduced Gompertz model in case of free growth.</caption><thead><tr><th>Parameters</th><th>Free growth (control case and EF0)</th></tr></thead><tbody><tr><td>V<sub>0</sub></td><td>0.051 (mean) - 0.0029 (std)</td></tr><tr><td>k</td><td>2.18 (mean) - 0.0 (std)</td></tr><tr><td>b</td><td>0.010 (mean) - 0.0017 (std)</td></tr><tr><td>MSE</td><td>0.00023</td></tr><tr><td>err<sub>1</sub></td><td>0.010</td></tr></tbody></table>
<figure><img src="image_6.png" /><figcaption>Figure 5: Volume evolution in case of free growth. Left-Top: fits using classical Gompertz model for the cohort of control experiment. Right-Top: fits using classical Gompertz model for the cohort EF0 of electroporation experiment. Left-Bottom: parameter a in function of parameter using classical Gompertz model for the fusion of the 2 free growth cohorts. Right- Bottom: fits using reduced Gompertz model for the fusion of the 2 free growth cohorts.</figcaption></figure>

**Electroporation** Figure 6 shows MSE associated to each considered initial time for EF500, EF1000-B and EF1000-A. Using this data, checking the stability of the estimated parameters and the quality of fits, we fix the value at $t_{init}$ = 90 h (resp. 90 h and 43 h) for cohort EF500 (resp. EF1000-B and EF1000-A). The corresponding fits are given in Figure 7 and the corresponding estimated parameters in Table 3. Figure 8-Left shows the Gompertz evolution for the three cohorts using the mean of estimated parameters between 0 h and 250 h. This figure can be directly compared to the data given in Figure 3-Right.

<figure><img src="image_7.png" /><figcaption>Figure 6: Mean squared error (MSE) depending on $t_{init}$ using reduced Gompertz model. Left: EF500. Middle: EF1000-B. Right: EF1000-A.</figcaption></figure>
<table><caption>Table 3: Estimated parameters of reversed reduced Gompertz model after electrical shock. Left: EF500. Middle: EF1000B. Right: EF1000A.</caption><thead><tr><th>Parameters</th><th>EF500</th><th>EF1000-B</th><th>EF1000-A</th></tr></thead><tbody><tr><td>t<sub>init</sub></td><td>90h</td><td>90h</td><td>43h</td></tr><tr><td>V<sub>last</sub></td><td>0.24 (mean) - 0.039 (std)</td><td>0.24 (mean) - 0.017 (std)</td><td>0.19 (mean) - 0.051 (std)</td></tr><tr><td>k</td><td>3.17 (mean) - 0 (std)</td><td>3.32 (mean) - 0 (std)</td><td>10.7 (mean) - 0 (std)</td></tr><tr><td>b</td><td>0.013 (mean) - 0.0021 (std)</td><td>0.0086 (mean) - 0.0047 (std)</td><td>0.011 (mean) - 0.0014 (std)</td></tr><tr><td>MSE</td><td>4.6 × 10<sup>-5</sup></td><td>1.7 × 10<sup>-5</sup></td><td>2.9 × 10<sup>-5</sup></td></tr><tr><td>err<sub>1</sub></td><td>0.0050</td><td>0.0032</td><td>0.0035</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>Figure 7: Volume evolution after electrical shock. Left: fits using reduced Gompertz model for EF500 ($t_{init} = 90$ h). Middle: fits using reduced Gompertz model for EF1000-B ($t_{init} = 90$ h). Right: fits using reduced Gompertz model for EF1000-A ($t_{init} = 43$ h).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 8: Simulations using the mean of estimated parameters for the 3 cohorts. Top: volume models (Gompertz after electroporation process). Bottom: 3D models.</figcaption></figure>
## 4.2 3D spheroid fits

For the sake of simplicity, the tilde notations will be removed in this section.

**Free growth** The parameters $a = kb$ and $b$ have been fixed at the means of the estimated parameters obtained using the 2 cohorts of control, see Table 2. To reproduce heterogeneity similar to the data available in [15] (Figure 6, CAPAN-2 control), we fix the values $d = 0.3$, $\tau_{PtoQ}^{begin} = 0.04$, $\tau_{PtoQ}^{end} = 0$ and $s = 0.05$. We consider an initial constant value of $P$ of 0.8. None of these 3 parameters ($d, \tau_{PtoQ}^{begin}, \tau_{PtoQ}^{end}$ and $s$) or of the initial condition $P(0, \cdot)$ can be estimated using only the volume of the spheroid as the equation on $R$ (or on $V$) does not depend on them. Figure 9 shows the evolution of $r \mapsto P(\cdot, r)$ at different time snapshots. Figure 10 (first column) shows the spheroid evolution with distribution of proliferative cells. These simulations are obtained by using the mean values of estimated values.

<figure><img src="image_10.png" /><figcaption>Figure 9: Free growth - time snapshots of $r \mapsto P(\cdot, r)$.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 10: Spheroid evolutions with distribution of proliferative cells. Black (resp white) = 100 (resp. 0) % of proliferative cells. Column 1: free growth, column 2: EF500, column 3: EF1000B and column 4: EF1000A and column 5: EF2000.</figcaption></figure>

**Electroporation process** Table 4 gives the estimation of the 6 parameters for the 4 cohorts. As expected, values of parameters $V_0, k$ and $b$ are close to the values estimated with control
data. The estimated parameters of cohort EF1000-B are very similar to the parameters of cohort EF500. The values of $\lambda$ increase with the intensity of the pulse. Concerning the value of $p$, it is interesting to note that there is no (resp. small) decrease of the volume for EF500 (resp. EF1000-B $V \cdot cm^{-1}$) while for the group 1000-A $V \cdot cm^{-1}$ (resp. 2000 $V \cdot cm^{-1}$), a decreasing of 80% (resp. almost 100%) is observed. Concerning $m$, the value increases with the value of the electrical field. Concerning the MSE and the $l^1$-norm, the values are close to the values obtained with the free growth cohort, see Table 2 for example. Fits are given in Figure 11. Figure 10 shows the spheroid evolutions with distribution of proliferative cells for the 4 cohorts (column 2 to column 5). Figure 8-Right shows the Gompertz evolution for the four cohorts using the mean of estimated parameters between 0 h and 250 h. This figure can be directly compared to the data given in Figure 3-Right.

<table><caption>Table 4: Estimated parameters of electroporation system (8)-(9) using volume data for the 4 cohorts: EF500, EF1000B, EF1000A and EF2000.</caption><thead><tr><th>Parameters</th><th>EF500</th><th>EF1000-B</th></tr></thead><tbody><tr><td>V<sub>0</sub></td><td>0.063 (mean) - 2.14 &times; 10<sup>-6</sup> (std)</td><td>0.062 (mean) - 0.00011 (std)</td></tr><tr><td>b</td><td>0.0071 (mean) - 0.0014 (std)</td><td>0.0086 (mean) - 0.00055 (std)</td></tr><tr><td>k</td><td>2.12 (mean) - 0 (std)</td><td>1.83 (mean) - 0 (std)</td></tr><tr><td>&lambda;</td><td>0.63 (mean) - 0.015 (std)</td><td>0.66 (mean) - 0.061 (std)</td></tr><tr><td>p</td><td>0.99 (mean) - 9.89 &times; 10<sup>-5</sup> (std)</td><td>0.86 (mean) - 0.15 (std)</td></tr><tr><td>m</td><td>1.36 (mean) - 0.039 (std)</td><td>1.98 (mean) - 0.32 (std)</td></tr><tr><td>MSE</td><td>0.00013</td><td>6.17 &times; 10<sup>-5</sup></td></tr><tr><td>err<sub>1</sub></td><td>0.0089</td><td>0.0060</td></tr></tbody><thead><tr><th>Parameters</th><th>EF1000-A</th><th>EF2000</th></tr></thead><tbody><tr><td>V<sub>0</sub></td><td>0.055 (mean) - 1.10 &times; 10<sup>-5</sup> (std)</td><td>0.057 (mean) - 0.0014 (std)</td></tr><tr><td>b</td><td>0.011 (mean) - 6.04 &times; 10<sup>-5</sup> (std)</td><td>0.0078 (mean) - 0.0019 (std)</td></tr><tr><td>k</td><td>2.30 (mean) - 0 (std)</td><td>2.46 (mean) - 0 (std)</td></tr><tr><td>&lambda;</td><td>0.90 (mean) - 0.0029 (std)</td><td>0.99 (mean) - 5.36 &times; 10<sup>-6</sup> (std)</td></tr><tr><td>p</td><td>0.20 (mean) - 0.059 (std)</td><td>0.0012 (mean) - 0.00019 (std)</td></tr><tr><td>m</td><td>4.89 (mean) - 0.67 (std)</td><td>5.81 (mean) - 6.43 (std)</td></tr><tr><td>MSE</td><td>4.38 &times; 10<sup>-5</sup></td><td>4.83 &times; 10<sup>-7</sup></td></tr><tr><td>err<sub>1</sub></td><td>0.0044</td><td>0.00024</td></tr></tbody></table>
<figure><img src="image_12.png" /><figcaption>Figure 11: Fits of electroporation system (8)–(9) using volume data for the 4 cohorts: EF500 (top-left), EF1000B (top-right), EF1000A (bottom-left) and EF2000 (bottom-right).</figcaption></figure>

# 5 Discussion and Conclusion

We first focused on the control cohort and we validated the reduced Gompertz model a more stable version of the very classical Gompertz model (1) using biological data corresponding to multicellular spheroids in free growth. We also showed that the presence of the electrodes did not impact the growth by comparing the control cohort with the data in the electroporation cohort, where electric pulses were not applied. It is possible to obtain reasonable fits using the reduced Gompertz model, see Figure 7. Even if this volume approach shows encouraging preliminary results, this approach did not allow to properly quantify the impact of the electroporation and did not integrate the heterogeneity of the tumor, i.e. the fact that the tumor was composed of proliferative and quiescent cells. These two limitations motivate the introduction of a PDE system enabled to represent the tumor evolution upon submission to electrical shocks. Three densities of tumor cells were considered: the proliferative cells, the quiescent cells and the cells whose functioning was altered by the effect of electric field. This model was confronted to biological data corresponding to multicellular spheroids submitted to electric pulses with different intensities. The quality of the fits on the spheroid volumes see Figure 11 and the errors see Table 4 validated our model. These results allowed determining the percentage of destroyed cells and the percentage of cells whose functioning was altered by the effect of electroporation for each value of the electrical field. One could see that for an electrical field close to $1000 \text{ V} \cdot \text{cm}^{-1}$ (partially irreversible electroporation and partially reversible electroporation), the increasing of the electrical field seemed to boost the resumption of spheroids growth, see Figure 8-Right.

Our results provide a good insight into the impact of electroporation on tumor growth. In clinics, when irreversible electroporation is used to ablate tumors, the treatment may induce cell death by irreversible electroporation. This effect is observed in a specific zone, namely,
the zone where the electric field strength is the highest at the optimal place in respect to the electrodes. The efficacy of completely irreversible electroporation was well observed in spheroids group exposed to 2000 V·cm⁻¹.

Nevertheless, adjacent to the epicenter of the highest electric field, where the tumor is ablated, we have zones that underwent milder electric field strengths exposure, and we have the occurrence of partially irreversible electroporation. Within this *shadow* zone, not all cells are killed, and the ones that survive will resume their growth. This effect was clearly present in the multicellular spheroids exposed 1000 V·cm⁻¹. This field amplitude generated thus in two radically different behaviors: either death or growth boost. We hypothesize that this variable response rate is due to the non-homogeneity of the electric field. If the spheroid (or tumor) was not efficiently attained, we have a relapse. If the positioning of the spheroid (or tumor) is in the center of the electrodes, the treatment is efficient. This is a recurrent problem also in clinics, where clinicians observe a different response rate either efficient treatment of a patient with a decreased tumor growth or tumor relapse [20]. The advantage of our model, is that soon after the electroporation, we can detect the trend of tumor growth. In practice, this model can thus allow to assess which tumors were inefficiently treated, and that the patient is on his or her way to relapse. This would indicate to the treating doctor that the treatment was not efficient and that an additional cycle of electroporation is necessary to re-ablate the tumor. Concerning the growth boost, it is particularly intriguing and certainly requires further inquiries, from the biological, physical and mathematical perspective.

## Acknowledgements

This work has been partly performed in the scope of the projects NUMEP (PC201615) and MECI n°21CM119-00 granted by Plan Cancer Inserm.

## Conflict of interest

The Authors declare that there is no conflict of interest.

## References

[1] Monolix version 2019R1. (Lixoft SAS, Antony, France), 2019.

[2] 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.

[3] R. Cannon, S. Ellis, D. Hayes, G. Narayanan, and R. C. Martin. Safety and early efficacy of irreversible electroporation for hepatic tumors in proximity to vital structures. *Journal of surgical oncology*, 107(5):544-549, 2013.

[4] A. E. Casey. The experimental alteration of malignancy with an homologous mammalian tumor material: I. results with intratesticular inoculation. *The American Journal of Cancer*, 21(4):760-775, 1934.

[5] L. Chopinet, L. Wasungu, and M.-P. Rols. First explanations for differences in electro-transfection efficiency in vitro and in vivo using spheroid model. *International journal of pharmaceutics*, 423(1):7-15, 2012.
[6] 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.

[7] S. K. Frandsen, L. Gibot, M. Madi, J. Gehl, and M.-P. Rols. Calcium electroporation: evidence for differential effects in normal and malignant cell lines, evaluated in a 3d spheroid model. *PLoS One*, 10(12):e0144028, 2015.

[8] L. Gibot, L. Wasungu, J. Teissié, and M.-P. Rols. Antitumor drug delivery in multicellular spheroids by electropermeabilization. *Journal of controlled release*, 167(2):138-147, 2013.

[9] B. Gompertz. XXIV. On the nature of the function expressive of the law of human mortality, and on a new mode of determining the value of life contingencies. In a letter to Francis Baily, Esq. FRS &c. *Philosophical transactions of the Royal Society of London*, pages 513-583, 1825.

[10] J. Kolosnjaj-Tabi, L. Gibot, I. Fourquaux, M. Golzio, and M.-P. Rols. Electric field-responsive nanoparticles and electric fields: physical, chemical, biological mechanisms and therapeutic prospects. *Advanced drug delivery reviews*, 138:56-67, 2019.

[11] E. Kuhn and M. Lavielle. Maximum likelihood estimation in nonlinear mixed effects models. *Computational Statistics & Data Analysis*, 49(4):1020-1038, 2005.

[12] A. K. Laird. Dynamics of tumour growth. *British journal of cancer*, 18(3):490, 1964.

[13] A. K. Laird. Dynamics of tumour growth: comparison of growth rates and extrapolation of growth curve to one cell. *British Journal of Cancer*, 19(2):278, 1965.

[14] M. Marty, G. Sersa, J. R. Garbay, J. Gehl, C. G. Collins, M. Snoj, V. Billard, P. F. Geertsen, J. O. Larkin, D. Miklavcic, I. Pavlovic, S. M. Paulin-Kosir, M. Cemazar, N. Morsli, D. M. Soden, Z. Rudolf, C. Robert, G. C. O'Sullivan, and L. M. Mir. Electrochemotherapy – an easy, highly effective and safe treatment of cutaneous and subcutaneous metastases: Results of esope (european standard operating procedures of electrochemotherapy) study. *European Journal of Cancer Supplements*, 4(11):3-13, 2006. Electrochemotherapy.

[15] 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.

[16] L. Miller, J. Leor, and B. Rubinsky. Cancer cells ablation with irreversible electroporation. *Technology in Cancer Research & Treatment*, 4(6):699-705, 2005. PMID: 16292891.

[17] L. Norton. A gompertzian model of human breast cancer growth. *Cancer research*, 48(24 Part 1):7067-7071, 1988.

[18] S. Pelofy, H. Bousquet, L. Gibot, M.-P. Rols, and M. Golzio. Transfer of small interfering rna by electropermeabilization in tumor spheroids. *Bioelectrochemistry*, page 107848, 2021.

[19] B. Ribba, N. H. Holford, P. Magni, I. Trocóniz, I. Gueorguieva, P. Girard, C. Sarr, M. Elishmereni, C. Kloft, and L. E. Friberg. A review of mixed-effects models of tumor growth and effects of anticancer drug treatment used in population analysis. *CPT: pharmacometrics & systems pharmacology*, 3(5):1-10, 2014.
[20] H. J. Scheffer, K. Nielsen, M. C. de Jong, A. A. van Tilborg, J. M. Vieveen, A. R. Bouwman, S. Meijer, C. van Kuijk, P. M. van den Tol, and M. R. Meijerink. Irreversible electroporation for nonthermal tumor ablation in the clinical setting: a systematic review of safety and efficacy. *Journal of Vascular and Interventional Radiology*, 25(7):997-1011, 2014.

[21] C. Vaghi, A. Rodallec, R. Fanciullino, J. Ciccolini, J. P. Mochel, M. Mastri, C. Poignard, J. M. Ebos, and S. Benzekry. Population modeling of tumor growth curves and the reduced gompertz model improve prediction of the age of experimental tumors. *PLoS computational biology*, 16(2):e1007178, 2020.

# Appendix

## Validation of the SAEM algorithm using few measurements

To check whether the SAEM algorithm works correctly with a low number of multicellular spheroids and a low number of measurements, Figure 12 shows the error $||V_0 - V_0^{obj}, k - k^{obj}, b - b^{obj}||_2$ using various cohorts with $y$ individuals, drawn at random from the control group and $x$ consecutive times. The values ${V_0, k, b}^{obj}$ correspond to the values estimated with the maximum number of measurements and multicellular spheroids (top, right of Figure 12).

<figure><img src="image_13.png" /><figcaption>Figure 12: Errors $||V_0 - V_0^{obj}, k - k^{obj}, b - b^{obj}||_2$ according to the number of consecutive measurements (x-axis) and to the number of multicellular spheroids (y-axis).</figcaption></figure>

## Radial equation of multicellular spheroids evolution in case of free growth

When considering free growth, we have $F(t, x) = 0, \forall t, \forall x$ leading to study the following system

$$ \begin{aligned} \partial_t P + \nabla \cdot (\vec{v}P) &= \tau_G(P+Q) - \tau_{PtoQ}P, \\ \partial_t Q + \nabla \cdot (\vec{v}Q) &= \tau_{PtoQ}P, \\ \nabla \cdot \vec{v} &= \tau_G(P+Q), \\ \tau_G(t) &= ae^{-bt}, \end{aligned} \quad (12) $$

coupled with the dynamics of $\tau_{PtoQ}$.

The first equation of System 12 can be rewritten as

$$ \partial_t P + \vec{v} \cdot \nabla P = \tau_G(P+Q)(1-P) - \tau_{PtoQ}P. \quad (13) $$
Simplified equation for the density of proliferating cells - We will search the simplified equation followed by P

$$ \begin{aligned} \partial_t P(t, x) &= \frac{d}{dt} \tilde{P}(t, r(t, x)) \\ &= \partial_t \tilde{P}(t, r(t, x)) - r(t, x) \frac{R'(t)}{R(t)} \partial_r \tilde{P}(t, r(t, x)) \end{aligned} $$

Concerning the gradient term $\nabla P$, we have

$$ \nabla P(t, x) = \nabla \tilde{P}(t, r(t, x)) = \nabla r(t, x) \partial_r \tilde{P}(t, r(t, x)). $$

As $\nabla r(t, x) = \frac{1}{R(t)} \frac{x}{\|x\|}$, we obtain

$$ \nabla P(t, x) = \frac{1}{R(t)} \frac{x}{\|x\|} \partial_r \tilde{P}(t, r(t, x)) $$

and using the fact that the velocity is radial and Equation 13, we obtain

$$ \partial_t \tilde{P}(t, r(t, x)) + \frac{\tilde{v} - rR'(t)}{R(t)} \partial_r \tilde{P}(t, r(t, x)) = \tau_G(\tilde{P} + \tilde{Q})(1 - \tilde{P}) - \tilde{\tau}_{PtoQ} \tilde{P}. \quad (14) $$

Using Eq. (7), this equation can be rewritten

$$ \partial_t \tilde{P}(t, r(t, x)) + \left( \frac{\tilde{v}}{R(t)} - \frac{r}{3} \tau_G(t) \right) \partial_r \tilde{P}(t, r(t, x)) = \tau_G(\tilde{P} + \tilde{Q})(1 - \tilde{P}) - \tilde{\tau}_{PtoQ} \tilde{P}. $$

Equation for the radial velocity - The equation $\nabla \cdot \tilde{v} = \tau_G(P+Q)$ becomes in radial coordinates

$$ \frac{1}{R(t)} \frac{1}{r^2} \partial_r (r^2 \tilde{v}(t, r)) = \tau_G(\tilde{P} + \tilde{Q}). $$

By integration of the previous equation between 0 and r (the invariance by rotation implies that $\tilde{v}(t, 0) = 0$), we obtain

$$ \frac{\tilde{v}(t, r)}{R(t)} = \frac{1}{r^2} \int_0^r \tilde{r}^2 \tau_G(\tilde{P} + \tilde{Q}) d\tilde{r}. $$

and this gives as $\tilde{P} + \tilde{Q} = 1$ inside the tumor

$$ \frac{\tilde{v}(t, r)}{R(t)} = \frac{r}{3} \tau_G(t). $$

This implies that $\tilde{P}$ follows for each $r \in [0, R(t)]$, the following EDO

$$ \tilde{P}' = \tau_G(1 - \tilde{P}) - \tilde{\tau}_{PtoQ} \tilde{P}. $$