Without cells

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

With cells

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

Figure 3.1 – Top: Bode plots of the impedances $Z_{B,M}^{(c)}$. The error bars represent the 95% of the impedance measured at a given frequency. Bottom: Scheme of the experiment and bode plots of the $Z_{C,BP,M}^{(5,12)}$ and $Z_{C,AP,M}^{(5,12,t)}$ with measurements before the pulse (in green) and those done after the pulse during five minutes in blue – increasingly clear at the passage of the time, as explained by the scheme. Left: Corresponding time evolving curves for the frequency 1 kHz. Right: Zoom on the beginning.
**Quick presentation of the data** The full setup used to obtain the bioimpedance measurements is presented in [29]. We introduce the following datasets:

$$F_{kHz} = \{1, 3, 5, 7, 11, 15, 19, 25, 33, 41, 51, 63, 77, 95, 117, 143, 173, 209, 253, 307, 371, 447, 539, 649, 781, 939\} \text{kHz (frequencies)},$$

$$C_B = \{0, 25, 50, 75, 100\} \text{ mM (concentration without cells)},$$

$$C_C = \{0, 5, 10, 20, 30\} \text{ mM (concentration with cells)},$$

$$V = \{12, 30\} \text{ V (considered voltages)},$$

$$T = \{t_n = n \times 210 \text{ ms}, n = 0, 1, \dots, 1500\} \text{ (time after pulse).}$$

The following notations are also used: B for buffer, $\mathcal{M}$ for measures, BP for before pulse and AP for after pulse. The measured impedances given in Figure 3.1 can be summarized as:

(1) *Impedance without cells* in the buffer denoted $Z_{B,\mathcal{M}}^{(c)} : f \mapsto Z_{B,\mathcal{M}}^{(c)}(f)$, for all $(c, f) \in C_B \times F_{kHz}$, see Figure 3.1 (top),

(2) *Impedance with cells* in the buffer:

(a) 100 ms (only the mean is considered) in pre-electroporation context denoted $Z_{C,BP,\mathcal{M}}^{(c,v)} : f \mapsto Z_{C,BP,\mathcal{M}}^{(c,v)}(f)$, for all $(c, v, f) \in C_B \times V \times F_{kHz}$, see Figure 3.1 (bottom, green curves),

(b) ~5 min *after one biphasic microsecond electroporation pulse of v \in V* denoted $Z_{C,AP,\mathcal{M}}^{(c,v,t)} : f \mapsto Z_{C,AP,\mathcal{M}}^{(c,v,t)}(f)$, for all $(c, v, t, f) \in C_B \times V \times T \times F_{kHz}$, see Figure 3.1 (bottom, blue curves).

For almost all cases, several experiments were performed. The measured impedances $Z_{B,\mathcal{M}}^{(c)}$, $Z_{C,BP,\mathcal{M}}^{(c,v)}$ and $Z_{C,AP,\mathcal{M}}^{(c,v,t)}$ differ from the true sample impedances – respectively denoted by $Z_{B,T}^{(c)}$, $Z_{C,BP,T}^{(c,v)}$ and $Z_{C,AP,T}^{(c,v,t)}$ – because of distortion due to the experimental setup. The strategy to deal with this issue is presented in [29].

**Equivalent circuit model** Our strategy to study the impact of electroporation on the electrical properties of the cells consists in analyzing impedance data using an equivalent electrical model of the experimental setup. It is composed of two parallel branches: one for the buffer compartment and one for the cell compartment. An important – and hardly measurable – phenomenon lies in the ion exchanges between the cells and the buffer. This phenomena is particularly important at low buffer conductivity. Therefore the buffer impedances obtained by the calibration are modified – and thus need to be estimated – in the experiments with cells. Denoting by $\sigma_{0,b}$ the buffer conductivity without KCl, the buffer conductivity as a function of the KCl concentration is deduced from the Kohlrausch law:

$$\sigma_b([KCl]_b) = \sigma_{0,b} + [KCl]_b (\Lambda_0 - K \sqrt{[KCl]_b}).$$

To account for the cells-buffer ion exchange, we set $\sigma_{0,b} = \sigma_0 + \Delta\sigma_0$, where $\Delta\sigma_0$ is the change of conductivity due to ions exchange between cells and the buffer that has to be fit. More precisely, in Figure 3.2, one can remark the data variability (dot points) between two same experiments, as for example for $c = 10$ or 20 mM with $v = 12$ V (top, left) or for $c = 5$ mM with $v = 30$ V (bottom, left). We assume that this variability is due to different states of the cells and exchanges between the cell monolayer and the buffer. The parameter $\Delta\sigma_0$ has been added to the buffer conductivity in order to account for this phenomenon.
<figure><img src="image_31.png" /><figcaption>Figure 3.2 – Bode plot of $Z_{C,BP,T}^{(c,12)}$ (top-left) and $Z_{C,AP,T}^{(c,12,t)}$ at different times for $c \in C$.</figcaption></figure>

The buffer resistance is then modeled as

$$R_b^r = \frac{R_b}{\phi} = \frac{k}{\sigma_b([KCl]_b)\phi}, \quad \text{with } k = 40 \text{ m}^{-1},$$

where $1 - \phi$ is the volume fraction of the cells.

The influence of the cells on the measured impedance is mainly due to the cytoplasm and the membrane of the cells similarly to [FS98]. The cytoplasm compartment is assumed only conductive while the membranes compartment is modeled by a constant phase element (CPE) of parameters ($C_m, \alpha_m$) in parallel with a resistance denoted by $R_m$ following [Lvo12, Chapter 3]. Before the pulse, $R_m$ is fixed at the value 2.34 MΩ. The impedance of our equivalent circuit with cells $Z_{C,th}$ (before and after pulse) is defined by

$$\frac{1}{Z_{C,th}(f, \theta)} = \frac{\sigma_b([KCl]_b)\phi}{k} + \frac{1}{\frac{R_{cyt}}{(1-\phi)} + \frac{R_m}{(1 + R_m C_m (2i\pi f)^{\alpha_m})}}, \quad (3.1)$$

where $\theta$ is a vector which concatenates the 6 parameters: $\theta = (\phi, \Delta\sigma_0, [KCl]_b, R_{cyt}, \alpha_m, C_m)$.

Before pulse, the parameters are estimated under constraints. The constraints and the estimated values are given in [29]. After pulse, the parameters are fixed at the estimated values before pulse except $R_m$ and $[KCl]_b$ which are fit by the minimization algorithm without any model.

Figure 3.2 shows the calibrated data $Z_{C,BP,T}^{(c,12)}$ (dot points and top-left) fitted with the electrical circuit before pulse i.e. with $R_m$ fixed at 2.34 MΩ (straight line) and the calibrated data $Z_{C,AP,T}^{(c,12,t)}$ (dot points, top-right and bottom) fitted with the electrical circuit after pulse (straight line) for $t = 30, 90$ and $300$ s.
As expected, the generated true impedances with cells before the pulse delivery $Z_{C,BP,T}^{(c,v)}$ display a flat modulus for low frequencies. This is due to the shielding effect of the membrane at low frequency. These results validate simultaneously our calibration procedure and the electrical circuit combined with Kohlrausch's conductivity to model the buffer. The curves with $v = 30$ V are given in [29]. One can see again that the impedances are accurately determined by the model. The errors on the phase are smaller than 2 degrees (for 12 V and 30 V) and the relative errors on the amplitudes are smaller that 1% for 12 V and for 30 V, which is very satisfactory for such a complex experimental setup.

<figure><img src="image_32.png" /><figcaption>Figure 3.3 – Evolutions of the membrane resistance $R_m$ within the first 2 seconds after the permeabilizing pulse (12V and 30V resp.).</figcaption></figure>

The previous calibration enables us to follow the influence of electroporation on the changes in the parameters corresponding to the cell membrane of the equivalent electrical circuit proposed. Figure 3.3 (resp. Figure 3.4) presents the evolution of $R_m$ (resp. [KCl] concentration) during 2 seconds (resp. 300 seconds) after electroporation. Let us emphasize that this postpulse calibration of $R_m$ and [KCl]$_b$ is performed without any *a priori* on the evolution of these 2 parameters.

As for the membrane resistance, the dynamics show a very fast and strong decrease due to the pulse, followed by a small increase (see Figure 3.3) within the first milliseconds after the pulse. These initial dynamics on the membrane could be interpreted as the opening and closing of the pores of the membrane due to the pulse. It is interesting to note that the instantaneous drop in membrane resistance at 30 V is much higher than at 12 V (see Figure 3.3), which again is consistent with the classical view of electropermeabilization (the pulses with higher magnitude produce more defects than the pulses with lower magnitude). As for the evolution of the membrane resistance within 5 min (see [29] for the figure), the calibration shows that the membrane compartment does not recover its original resistance after 5 min to the value before the pulse, which is consistent with several experimental observations showing that electroporation takes several minutes [PKW+07, MNS+22]. However, the membrane resistance we estimate continues to decrease after several tens of seconds. We think that this unanticipated behaviour is due to ion exchange between the buffer and the cells. Indeed, the cells and the buffer exchange more ions than just the KCl ions, so it is more likely that $\sigma_{0,b}$ increases with time. In [29] we explore a way to answer this problem and conclude that time point measurements of buffer conductivity after the pulse would allow us to better track the long-term behaviour of the membrane after the pulse.

In [29], we show that the evolution of [KCl] concentrations can be described by the following function
<figure><img src="image_33.png" /><figcaption>Figure 3.4 – $[KCl]_b$ concentrations and fits obtained when fitting them with the function given in Equation (3.2).</figcaption></figure>

$$[KCl]_b^{\text{fit}} = [KCl]_b^{\text{begin}} + ([KCl]_b^{\text{end}} - [KCl]_b^{\text{begin}})(1 - e^{-t/\tau_{[KCl]_b}}). \quad (3.2)$$

The calibrations are plotted in Figure 3.4. The fits are very closed to the data, validating the model. From the values given in [29], one sees that for a given experiment, the release of $[KCl]$ –that is $[KCl]_b^{\text{end}} – [KCl]_b^{\text{begin}}$– is higher after the pulse at 30 V than at 12 V, and the time constant of the release $\tau_{[KCl]_b}$ is shorter at 30 V than at 12 V, meaning the dynamics of the release is faster at 30 V than at 12 V. Another interesting observation is the fact that the time constants $\tau_{[KCl]_b}$ are similar between all the experiments for a given voltage, which is consistent with the fact that the permeabilisation depends on the magnitude of the pulse. We eventually note that for higher initial concentration in the buffer, the release is smaller due to a lower gradient between the cell and the buffer $[KCl]$ concentrations.

**Conclusion** In summary, thanks to a robust setup calibration procedure, the data sets without and with cells were evaluated sequentially to appropriately calibrate the electrical circuit corresponding to the setup. Interestingly, we were able to distinguish between the release of $[KCl]$ and the dynamics of the membrane resistance within the first seconds after the pulse. This work provides, for the first time, numerical quantification of membrane resealing after an electroporation pulse using bioimpedance measurements. Our calibration suggests that additional measurements of buffer conductance should be made to track membrane dynamics during the 5 min following the pulse. We emphasize that the electrical equivalent circuit approach is a very interesting step to decipher membrane resealing immediately after the pulse. However, a 3D simulation of the entire experimental setup, which is currently under investigation, could provide more insight. This could involve writing a model that includes an intracellular region corresponding to the cell layer, which means that the layer is considered as a very large cell. To do this, modeling electroporation at the cell scale and, in particular, membrane behavior is a necessary step. This topic will be explored in the next section. This means that we will focus on the local behavior of the membrane.
### 3.1.2 Membrane modeling

Even though the phenomenon of electroporation has been discovered in the late 1960s, the ways the membranes become permeable by the effect of an intense electric field is still not well understood. In the late 1990s, the teams of Chizmadzhev and Weaver [CCSZ95, WC96] proposed a description of the emergence of water pores in a sea of lipids, under very constraint geometries (the pores are cylindrical). Then Krassowska and Neu proposed a Smoluchovskii equation for the population of cylindrical pores $n(r, t)$ of radius $r > 0$ at the time $t$, and they derived an asymptotic analysis to link the transmembrane potential to the total density of pores $N(t) = \int_{\mathbb{R}^+} n(r, t)dr$, which satisfies an ordinary differential equation [NK99, DK98]. From this times, only slight modifications of the models have been proposed by Weaver [WSE+12, SSGW14]. Kavian et al. proposed then a phenomenological version of the Krassowska and Debruin model, by limiting the number of parameter [KLPW14]. Leguèbe et al. proposed in [LSMP14] an extension by introducing the surface reaction-diffusion of oxidised lipids. This model – coupled with well-chosen boundary conditions for the electrodes used to measure, see Section 3.2.2 – could be used in a harmonic context to model the 3D configuration of the setup presented in the previous section. However, this model seems closer to the observations but it is phenomenological and lacks of any physical basis.

The aim of this section is to present a new model of membrane electroporation that combines the water content of the membrane (initially 0 for non porated membrane) and the transmembrane voltage. Interestingly, thanks to a well defined free-energy of the membrane, we somehow generalise the approach of Chizmadzhev and Weaver, getting rid of the geometrical cylindrical assumption. Our approach is physically relevant and we recover a surface diffusion equation of the lipid phase that can be seen as a justification of the model of Leguèbe et al.

**Quick presentation of the model** We model the state of the membrane by a function $\phi: \Gamma \to \mathbb{R}$ measuring the amount of lipid in the membrane $\Gamma$ at any position $x$ at time $t$. The state $\phi = 0$ represents the ideal pure lipid phase and $\phi = 1$ the ideal pure water phase. The free-energy of the membrane subjected to a transmembrane voltage $u^{(m)}$ is the functional energy so-called Gibbs energy given as [KP14]

$$ \mathcal{E}(\phi, u^{(m)}) = \frac{\kappa}{2} \int_{\Gamma} |\nabla\phi|^2 + \int_{\Gamma} W_m(\phi) - \frac{1}{2} \int_{\Gamma} C_m(\phi)(u^{(m)})^2, $$

where the first term corresponds to the membrane diffusion, the second to the potential energy of the membrane and the third to the influence of the transmembrane voltage. More precisely, $\kappa > 0$ is linked with the interfacial tension [FMZB20] and $W_m$ is the double-well potential energy that describes the stable states of the membrane. The term $|\nabla\phi|^2$ regularises $\phi$ and helps control the thickness of the water-lipid interface, as explained by Bray [Bra02]. The term $\frac{1}{2}C_m(\phi)(u^{(m)})^2$ is the electrical energy. The effect of the transmembrane voltage acting on the membrane is modeled as the energy associated to the change in specific capacitance of the membrane [AAC+79]. In fact, Weaver and Chizmadzhev [WC96] consider a similar term taking into account the difference of dielectric coefficients for water and lipids. Roughly speaking, this term affects the potential energy to favor the phase $\phi = 1$ leading to the entry of water within the membrane (see Figure 3.5).

As the water volume fraction is not conserved (defects can be created and disappear from the membrane), we consider the non-conserved dynamics associated to this energy functional (also called model-A by Bray [Bra02]). The evolution of $\phi$ is determined by the $L^2$-gradient flow associated with the energy functional above. It corresponds to the Euler-Lagrange equation for the energy $\phi \mapsto \mathcal{E}(\phi, u^{(m)})$. In other words, for a given kinetics coefficient $\alpha > 0$ – also called phase field mobility – the order parameter $\phi$
<figure><img src="image_34.png" /><figcaption>Figure 3.5 – Electrostatic term $W_m(\phi) - \frac{1}{2}C_m(\phi)(u^{(m)})^2$ versus $\phi$ when the transmembrane potential is assumed to be constant throughout the membrane. Each curve corresponds to a different value of $u^{(m)}$ between 0 and 1 V.</figcaption></figure>

satisfies the following Allen-Cahn equation on the membrane surface $\Gamma$ [AC79, Bra02, CH58]:

$$ \partial_t \phi - D_0 \Delta \phi = -\alpha \frac{\partial W_m}{\partial \phi}(\phi) + \frac{\alpha}{2} \frac{\partial C_m}{\partial \phi}(\phi) v^2, \quad \forall t > 0, \qquad (3.3) $$

$$ \phi(0, \cdot) = \phi^0(\cdot), \qquad (3.4) $$

where $D_0 := \alpha\kappa$ is a diffusion coefficient in $m^2.s^{-1}$, and $\alpha$ is in $m^2.J^{-1}.s^{-1}$.

Throughout the work, we consider for $W_m$ the simple smooth double-well potential given by

$$ W_m(\phi) := 16a_1\phi^2(1-\phi)^2 + 8a_2(\phi + \frac{1}{2})(\phi - 1)^2, \quad \forall \phi \in [0, 1], \qquad (3.5) $$

where $4a_1 > 3|a_2|$, so that $W_m$ satisfies

$$ W_m(1) = 0, \quad W'_m(1) = W'_m(0) = 0, \quad W''_m(0) > 0, \quad \text{and} \quad W''_m(1) > 0. $$

The above energy potential describes the fact that the pure lipid ($\phi = 0$) and pure water ($\phi = 1$) phases are the only two stable phases with an energy barrier roughly equal to $W_m(1/2)$ between them. In the above expression $a_2$ mainly controls the value of the first local minimum ($W(0) = 4a_2$) and $a_1$ the size of the local maximum (at $\phi_0 := 1/2 - 3a_2/8a_1$) which is the height of the energy barrier to pass from 0 to 1.

**Remark 5.** Concerning the theoretical aspects, a well-posedness result for the model of the membrane under some known electric field influence has been established in [25]. Furthermore, there is an extensive literature studying the qualitative mathematical properties of (3.3) when $u^{(m)} = 0$. One can cite, for example [AC79, ESS92, AHM08].
The equation on $\phi$ is coupled to the microscopic bidomain model already introduced in Chapter 1 (System (1.6)) in a cardiac context. In this work, $\Gamma$ is a closed 2D-surface of $\mathbb{R}^3$ without boundary. We denote by $\Omega$ the bounded domain of $\mathbb{R}^3$ in which $\Gamma$ is embedded, and $\Omega^{(e)}$ and $\Omega^{(e)}$ are the 2 connected components subsets of $\Omega$ separated by $\Gamma$. We denote as previously $u^{(i)}$ (resp. $u^{(e)}$) the intra-cellular (resp. extra-cellular) potential and $\sigma^{(i)}$ (resp. $\sigma^{(e)}$) its associated diffusion coefficient. The transmembrane potential $u^{(m)}$ is then as previously defined as the jump between $u^{(i)}$ and $u^{(e)}$ along the membrane $\Gamma$ i.e. $u^{(m)} = [u^{(i)} - u^{(e)}]_{\Gamma}$. We will consider in this electroporation context the following bidomain model

$$ \left\{ \begin{array}{ll} \vec{\nabla}_{\vec{x}} \cdot (\sigma^{(\alpha)} \vec{\nabla}_{\vec{x}} u^{(\alpha)}) = 0, & \text{in } \Omega^{(\alpha)} \\ \sigma^{(i)} \vec{\nabla}_{\vec{x}} u^{(i)} \cdot \underline{n}^{(i)} = \sigma^{(e)} \vec{\nabla}_{\vec{x}} u^{(e)} \cdot \underline{n}^{(i)}, & \text{on } \Gamma \\ -\sigma^{(i)} \vec{\nabla}_{\vec{x}} u^{(i)} \cdot \underline{n}^{(i)} = C_m(\phi) \partial_t u^{(m)} + S^{(m)}(\phi) u^{(m)}, & \text{on } \Gamma. \end{array} \right. \quad (3.6) $$

One can see that the differences with System (1.6) (when the mechanical deformations are neglected) concern the last equation which has been modified. Indeed the last equation in System (1.6) was

$$ -\sigma^{(i)} \vec{\nabla}_{\vec{x}} u^{(i)} \cdot \underline{n}^{(i)} = C_m \partial_t u^{(m)} + I^{\text{ion}}(u^{(m)}, w), $$

meaning that the main differences are:

* The capacitance per unit area of the membrane $C_m$ depends on $\phi$ (as one could have remarked in Equation (3.3)).

* Instead of the ionic term, we consider the conductance $S^{(m)}$ which measures the opposition to the flow of electric current. It models the electroporation of the cell meaning the opening of the pores. The ionic term is not considered as we do not consider cardiac cells in this context. However, a second term could be added to model the permeabilization of the membrane, see the conclusion paragraph of this subsection.

Furthermore, this system is then coupled to Dirichlet boundary condition on the electric potential

$$ u^{(m)}(t, \cdot) = g(t, \cdot), \quad \partial\Omega, $$

where $g$ corresponds to the applied voltage.

The analysis of the coupled model (3.6-3.3) has been done in [25]. It extends the results obtained in [KLPW14] to time-varying and space dependent capacitance.

The first step consists in eliminating the unknown electric potentials following an idea developed in [BCP09] for the bidomain equation at the macroscopic level or in [KLPW14] for an electroporation model at the microscopic level and already quickly introduced in Chapter 1 (see Equation 1.13). We show that there exist two regular operators $\Lambda$ and $G$ defined as a combination of Dirichlet-to-Neumman type of operators such that solutions of System (3.6) verify

$$ C_m(\phi) \partial_t u^{(m)} + (S^{(m)}(\phi) + \Lambda) = Gg. \quad (3.7) $$

**Spherical vs flat membranes** Always in [25], we characterize the operator $\Lambda$ in two important configurations: a spherical membrane $\mathcal{S}$ and a flat periodic membrane defined $(\mathbb{R}/(L\mathbb{Z}))^2$, see Figure 3.6 for the definitions and the notations. These two cases are of importance, because cells in suspension are mostly
<figure><img src="image_35.png" /><figcaption>Figure 3.6 – Two situations of membrane. Left: Diagram of spherical cell model. The cell membrane $\Gamma$ is the sphere of radius $R_0$ and the boundary condition at $\partial\Omega$ describes the effect of the electric field applied. Right: Bi-periodic square flat membrane of size $L$. The two surfaces $\partial\Omega^{\pm}$ parallel to $\Gamma$ and equidistant from it, represent the electrodes. This very simplified geometry can be seen as a volume extraction close to the boundary $\Gamma$ of the spherical case (see blue draw).</figcaption></figure>

round, while molecular dynamic simulations deals mostly with flat membrane. It is thus important to compare these two settings, in particular to understand the main differences in terms of order of magnitude of the time constants. The operators in this two particular cases can be explicitly diagonalised in terms of the eigenfunctions of the Laplace-Beltrami operator on the sphere. The results are given in [25] and interestingly it allows to determine the charging times of these two cases. In order to consider the periodic flat membrane as a zoomed-in flat patch of the spherical cell model, it is thus necessary to match them accordingly. This means that to compare both settings, we need to consider the following condition

$$H = \frac{R_0 \sigma^{(i)} + 2\sigma^{(e)}}{2 \sigma^{(i)} + \sigma^{(e)}}.$$

It is worth noting that this condition is far from being satisfied by the setup of the molecular dynamics simulations of membrane electroporation, for which the simulation boxes are a few nanometers thick, while cell radii are about several of micrometers. Therefore it is highly probable that the time constants of the molecular dynamic simulations are too small. It is clear that in the case of pulsed electric fields, this condition cannot be avoided by changing the duration or magnitude of pulses in each setting (spherical cell or flat membrane). In fact, as much as the charge time is affected by the geometry, so is the discharge time, therefore even between pulses different behaviors arise in both settings.

Aside from the height $H$, the other important geometric factor of the flat membrane is its length $L$. Its influence on the transmembrane potential is not so easy to quantify and compare to the spherical case. In [25], we show that small patches tend to flatten the transmembrane potential artificially even if we keep the same proportion of water in the membrane. This heuristic calculation suggests that large patches of flat membrane are needed so that flat torus can be seen as a zoom of the spherical setting.
**Membrane conductance and capacitance models** The membrane conductance $S_m$ directly impacts on the amplitude of the transmembrane voltage $u^{(m)}$ whereas $C_m$ directly affects the order parameter $\phi$ as it is directly responsible (along with the transmembrane potential) for pore formation, which is why we must also consider the membrane dynamics in this modeling. We consider for $S^{(m)}$ a sigmoid function which interpolates between the conductivity of the membrane $\sigma_l$ (lipid conductivity) and the conductivity of a pore $\sigma_w$ with $k_0 = 100$ the steepness of the transition and $\phi_{th} = 0.5$ the transition threshold. For $C_m$, we consider the modeling given in [Loo65] whose the parameter values are given in [25].

We emphasize that these choices will have little effect on the pore evolution once the membrane is discharged, so we can expect the same pore sealing behavior.

**Numerical simulations** We perform numerical simulations in the case of a flat toroidal membrane. To determine the values of the parameters of our coupled model, we either use the literature or compare our model to others or examine limiting cases, see [25] for the values.

An accurate splitting scheme combined with Fast Fourier Transforms for efficient computations of the model is presented in [25]. It is worth noting that the solving the problem in a geometry close to the cell geometry is still challenging. Indeed, a pore edge is somehow comparable in size to the width of the cell membrane $h \sim 10^{-9}$ m and the characteristic length of a cell membrane is comparable with the cell radius $R_0 \sim 10^{-5}$ m. Then to capture the evolution of pores we need a discretized membrane with a minimum amount of nodes of around $(\frac{10^{-5}}{10^{-9}})^2 = 10^8$. Additionally, calculating Dirichlet-to-Neumann operators is particularly slow even in the case where we know how to diagonalize it. Therefore such simulations require specific numerical strategies that will be addressed in a near future.

We simulate a flat biperiodic membrane of size $L = 200$ nm under the influence of a 12 ns square pulse of a uniform electric field of intensity $|E| = 3.2 \cdot 10^6$ V.m$^{-1}$. After the membrane discharges its transmembrane potential, we continue to simulate it for another $10 \ \mu s$. We pick as the initial condition $\phi_{\diamond} : (x, y) \mapsto \mu(x, y)$, where $|\mu| < 1.5 \cdot 10^{-2}$ is a Gaussian random noise (which has been smoothed for numerical stability purposes).

During the initial charge time the membrane evolves more or less uniformly in space until about $t = 12$ ns. At this point the electric field is turned off, but the membrane remains charged. It then discharges like a capacitor, slowly, while the order parameter continues to evolve under the influence of the remaining transmembrane potential. At around $t = 27$ ns, the conductive effect of the membrane comes into play. At this point, the membrane is still charged, but now discharges mainly through the sites of high conductivity resulting from the affected order parameter ($\phi > 1/2$). These sites are better seen at the same time in the transmembrane potential distribution (bottom). However, because the membrane charge is low, the membrane discharges within nanoseconds after this phenomenon occurs. This sudden discharge of the transmembrane potential leads to the assumption that electroporation has occurred.

However, a definite membrane configuration does not occur until about $1 \ \mu s$ later. In this case, the membrane is almost completely discharged by $t = 33$ ns. From this point on, the transmembrane voltage no longer plays a significant role, and only the Allen-Cahn equation controls the membrane dynamics. Figure 3.7 shows the evolution of the transmembrane potential shortly after the initial electric pulse. The lifetime of the pores depends on their size. In our simulations, most pores disappeared before $10 \ \mu s$. Of course, the evolution of these pores is consistent with the Allen-Cahn equation, since they minimize the mean curvature. Therefore, we see that the pores become rounder as they shrink.
<figure><img src="image_36.png" /><figcaption>Figure 3.7 – The first row of simulations corresponds to the membrane order parameter $\phi$ while the second row gives the associated transmembrane potential $u^{(m)}$.</figcaption></figure>

<figure><img src="image_37.png" /><figcaption>Figure 3.8 – The evolution of the membrane $\phi$ for 10 µs after the end of the pulse delivery.</figcaption></figure>
**Conclusion** We present a new model of membrane electroporation that combines membrane water content and transmembrane voltage. Interestingly, thanks to a well-defined membrane free energy, we generalise the pioneering approach of Chizmadzhev, Weaver, and Krassowska by getting rid of the geometric cylindrical assumption on which most of the current electroporation models are based. Our approach is physically relevant and we recover a surface diffusion equation of the lipid phase proposed by Leguèbe et al. in an earlier phenomenological model. We also perform a fine analysis of the non-local operators involved in two simple configurations (a spherical membrane and a flat periodic membrane), which allows us to compare the time constants of the phenomenon in spherical and flat membranes. An accurate splitting scheme combined with Fast Fourier Transforms is developed for efficient computation of the model. Our numerical results allow us to establish a link between molecular dynamics simulations of membrane permeabilization and experimental observations on vesicles and cells.

From the modeling point of view, it is incomplete at it concerns only cell electroporation and not cell permeabilization. Indeed, it has been shown by L. M. Mir's group that the intense electric field generates an oxidation of the lipids, that could explain the long duration of the permeabilised state of the membrane. This means that the pores are closed but the membrane is still permeable. To integrate this phenomenon, a possible extension of our model is to consider two order parameters: ones which follows the (non conserved) Allen-Cahn equation for the water content and a second ones which follows Cahn-Hilliard equation on the lipid phases (oxidised or not).

Furthermore, the changes in the ionic and adenosine triphosphate (ATP) concentrations of the permeabilization modeling, and in the cellular volume need to be addressed if we would like to make this model compatible with electroporation of cardiac cells. More details will be given in my research project given in Chapter 4.

## 3.2 Tissue modeling

In parallel with the work – presented on the previous section – on membrane modeling, we examine modeling at the tissue scale. It is the subject of this new section. To obtain electroporation modeling at the tissue scale, one could consider homogenization of the microscopic electroporation bidomain model. The formal derivation of this model based on 2-scale convergence – see Section 1.2.2 of Chapter 1 – has not been done. This homogenization will allow to determine a model at the tissue scale compatible with the model at the cell scale. This point is on the list of my research project, see Chapter 4 for details.

However, other electroporation models are available in the literature. The objective of this section is to first present the classical models at the tissue scale which already exist and are used in the literature. We also introduce a macroscopic bidomain system inspired from our knowledge of the phenomenon but not validated through the homogenization of a microscopic model. More precisely, in [15], we compare them and we investigate their calibrations for electrostatic descriptions of a tissue from bioelectrical measurements on rabbit livers performed by Sel et al. [SCB+05]. This work is presented in Section 3.2.1.

However, either the considered model, the choice of boundary conditions of the – very conductive and thin – passive electrodes is a difficult question because, unlike the active electrodes on which either the potential or the current is imposed – corresponding to Dirichlet or Neumann boundary conditions – for the inactive electrodes the value of the equipotential is an unknown of the PDE. This leads to a non-local constraint of the flux along the conductor interface, the so-called floating potential problem. We develop a strategy to address this problem presented in Section 3.2.2 and published in [18]. One could note that this problem appears in the bioimpedancemetry experiment presented in Section 3.1.1. Indeed, the setup includes six spiral micro-electrodes, see [29] whose two correspond to active electrodes to deliver
the electroporating pulses and four to measure the voltage.

<figure><img src="image_38.png" /><figcaption>Figure 3.9 – The typical geometrical configuration consists of the tissue domain $\Omega$ in gray, deprived of the two needles $\mathcal{E}^{\pm}$, in white. The outer boundary $\Gamma_{\text{out}}$ is represented in bold.</figcaption></figure>

### 3.2.1 Comparison and calibration of different electroporation models

To compare the different models available in the literature, the following simple configuration is considered, see Figure 3.9. The domain $\tilde{\Omega}$ is a connected, open smooth subset of $\mathbb{R}^d$, $d \in \{2, 3\}$, and consists of two non intersecting subsets: $\Omega$, which corresponds to the tissue under consideration (rabbit liver for instance in the experiment of Sel et al. [SCB+05]) and the needles, represented by two parallel cylinders $\mathcal{E}^{\pm}$. The needles are set at the isopotential $\pm g$ respectively. The outer boundary is denoted by $\Gamma_{\text{out}} := \partial\Omega \setminus \mathcal{E}^{\pm}$.

#### 3.2.1.1 Different electroporation models

**The standard static model** The most used model to describe tissue electroporation consists of the phenomenological electrostatic problem. The tissue is described as a conductive medium, whose conductivity tensor $\sigma$ depends on the amplitude of the electric field $-\nabla u$. The model reads then

$$ \left\{ \begin{array}{ll} -\nabla \cdot (\sigma(\|\nabla u\|)\nabla u) = 0, & \Omega, \\ \partial_n u|_{\Gamma_{\text{out}}} = 0, & u|_{\mathcal{E}^{\pm}} = g^{\pm}, \end{array} \right. \qquad (3.8) $$

The tissue conductivity consists of a 4 parameters sigmoid function. Typically

$$ \forall \lambda \ge 0, \sigma(\lambda) = \sigma_0 + \frac{\sigma_1 - \sigma_0}{2} (1 + \operatorname{erf}(k_{\text{ep}}(\lambda - E_{\text{th}}))), \qquad (3.9) $$

where $\sigma_0$ is the conductivity of the non electroporated tissue, $\sigma_1$ is the tissue conductivity of the fully porated tissue, $E_{\text{th}}$ is the threshold amplitude for electroporation, and $k_{\text{ep}}$ is the slope of the nonlinearity. Here, erf is the Gauss error function.

**Remark 6.** It should be noted that this choice of $\sigma$ is largely phenomenological, and as shown in [IVM10], the available experimental data does not seem sufficient to characterize the dependence of the conductivity on the electric field.
**The biphasic dynamical model** In [VSM<sup>+</sup>18], Voyer et al. propose a biphasic dynamical model based on the description of an individual cell and surrounding matrix as an electric circuit. The ODEs at the cell level are formally generalized to PDEs at the tissue level. It describes the extra-cellular potential $u^{(e)}$ and the electric field inside cells $\mathbf{J}^{(i)}$. The parameters are the extra-cellular and intra-cellular electric conductivities, respectively $\sigma^{(e)}$ and $\sigma^{(i)}$. The conductivity of the cell membrane $\sigma^{(m)}$ depends on time in a way which mimics the effects of poration (appearance of holes on the membrane) and permeabilisation, (degradation of membrane molecules). Both phenomena increase the conductivity as discussed in the conclusion of the previous section in the context of microscopic scale. The resulting system reads

$$ 
\left\{ 
\begin{array}{lll}
\nabla \cdot (\sigma^{(e)} \nabla u^{(e)} + \mathbf{J}^{(i)}) = 0, & \Omega, \\
\varepsilon_m \partial_t \mathbf{J}^{(i)} + (\sigma^{(m)}(t, \|\mathbf{E}^{(m)}\|) + \sigma^{(i)}) \mathbf{J}^{(i)} = \sigma^{(i)} \sigma^{(m)} \nabla u^{(e)}, & \Omega, \\
\sigma^{(m)}(t, \|\mathbf{E}^{(m)}\|) = \sigma_0^{(m)} + \sigma_1^{(m)} X_1(t, \|\mathbf{E}^{(m)}\|) + \sigma_2^{(m)} X_2(t, \|\mathbf{E}^{(m)}\|), & \Omega, \\
\partial_n u^{(e)}|_{\Gamma_{\text{out}}} = 0, \quad u^{(e)}|_{\varepsilon^{\pm}} = g^{\pm}, \quad \mathbf{J}^{(i)}|_{\Gamma_{\text{out}}} = 0,
\end{array}
\right.
\quad (3.10) $$

where $\mathbf{E}^{(m)} = \nabla u^{(e)} - (\sigma^{(i)})^{-1}\mathbf{J}^{(i)}$. The functions $X_1$ and $X_2$ are the respective degrees of poration and degrees of permeabilisation. They satisfy

$$ \dot{X}_1 = \frac{\beta_1(\|\mathbf{E}^{(m)}\|) - X_1}{\tau_1} \quad \text{and} \quad \dot{X}_2 = \frac{\beta_2(X_1) - X_2}{\tau_2}, \quad (3.11) $$

where $\beta_1$ and $\beta_2$ are 2-parameters sigmoid functions of the form $\beta_l(x) = \frac{1}{2}(1 + \text{erf}(k_l(x - x_l^{\text{trans}})))$, where $k_l$ is the stiffness of the sigmoid and $x_l^{\text{trans}}$ the transition threshold.

A simplified static version of this model allows to compare it with the standard model. Indeed, the static version of this model reads

$$ 
\left\{ 
\begin{array}{lll}
\nabla \cdot (\sigma^{(e)} \nabla u^{(e)} + \mathbf{J}^{(i)}) = 0, & \Omega \\
(\sigma^{(m)}(\|\mathbf{E}^{(m)}\|) + \sigma^{(i)}) \mathbf{J}^{(i)} = \sigma^{(i)} \sigma^{(m)} (\|\mathbf{E}^{(m)}\|) \nabla u^{(e)}, & \Omega \\
\partial_n u^{(e)}|_{\Gamma_{\text{out}}} = 0, \quad u^{(e)}|_{\varepsilon^{\pm}} = g^{\pm}, \quad \mathbf{J}^{(i)}|_{\Gamma_{\text{out}}} = 0,
\end{array}
\right.
\quad (3.12) $$

where $\sigma^{(m)}$ is a sigmoid function similar to (3.9). Written in $u^{(e)}$ and $\mathbf{E}^{(m)}$, the equations become

$$ 
\left\{ 
\begin{array}{lll}
\nabla \cdot \left( (\sigma^{(e)} + \sigma^{(i)}) \nabla u^{(e)} - \sigma^{(i)} \mathbf{E}^{(m)} \right) = 0, & \Omega, \\
(\sigma^{(m)}(\|\mathbf{E}^{(m)}\|) + \sigma^{(i)}) \mathbf{E}^{(m)} = \sigma^{(i)} \nabla u^{(e)}, & \Omega, \\
\partial_n u^{(e)}|_{\Gamma_{\text{out}}} = 0, \quad u^{(e)}|_{\varepsilon^{\pm}} = g^{\pm}, \quad \mathbf{E}^{(m)}|_{\Gamma_{\text{out}}} = 0.
\end{array}
\right.
\quad (3.13) $$

System (3.13) makes it possible to define the equivalent tissue conductivity. Actually, the potential $u^{(e)}$ satisfies

$$ \nabla \cdot \left( \left( \sigma^{(e)} + \frac{\sigma^{(i)} \sigma^{(m)} (\|\mathbf{E}^{(m)}\|)}{\sigma^{(i)} + \sigma^{(m)} (\|\mathbf{E}^{(m)}\|)} \right) \nabla u^{(e)} \right) = 0. $$

The equivalent tissue conductivity can thus be defined as

$$ \sigma_{\text{eq}}(\|\mathbf{E}^{(m)}\|) = \sigma^{(e)} + \frac{\sigma^{(i)} \sigma^{(m)} (\|\mathbf{E}^{(m)}\|)}{\sigma^{(i)} + \sigma^{(m)} (\|\mathbf{E}^{(m)}\|)}. $$

A simplified version of the problem consists in assuming that $\sigma^{(m)}$ depends on $\|\nabla u^{(e)}\|$ instead of $\|\mathbf{E}^{(m)}\|$.
Then the simplified problem reads

$$ \left\{ \begin{array}{l} \nabla \cdot \left( \left( \sigma^{(e)} + \frac{\sigma^{(i)} \sigma^{(m)} (\|\nabla u^{(e)}\|)}{\sigma^{(i)} + \sigma^{(m)} (\|\nabla u^{(e)}\|)} \right) \nabla u^{(e)} \right) = 0, \quad \Omega, \\ \partial_n u^{(e)}|_{\Gamma_{\text{out}}} = 0, \quad u^{(e)}|_{\varepsilon^{\pm}} = g^{\pm}, \end{array} \right. \quad (3.14) $$

which is similar to the standard model (3.8) since the function $\sigma_{eq}$ is a sigmoid function

$$ \forall \lambda \ge 0, \quad \sigma_{eq}(\lambda) = \sigma^{(e)} + \frac{\sigma^{(i)} \sigma^{(m)}(\lambda)}{\sigma^{(i)} + \sigma^{(m)}(\lambda)}. $$

**The static bidomain model** As explained in Section 1.2.2 of Chapter 1, the bidomain model has been proven to be the homogenization limit of the cell scale electric potential. We consider the following static bidomain written in terms of the extra-potential $u^{(e)}$ and the transmembrane potential $u^{(m)}$

$$ \left\{ \begin{array}{l} \nabla \cdot \left( (\sigma^{(e)} + \sigma^{(i)}) \nabla u^{(e)} - \sigma^{(i)} \nabla u^{(m)} \right) = 0, \qquad \Omega \\ A_m S^{(m)} (\|\nabla u^{(e)}\|) u^{(m)} - \nabla \cdot (\sigma^{(i)} \nabla u^{(m)}) + \nabla \cdot (\sigma^{(i)} \nabla u^{(e)}) = 0, \qquad \Omega \\ \sigma^{(e)} \nabla u^{(e)} \cdot \vec{n}|_{\Gamma_{\text{out}}} = 0, \quad u^{(e)}|_{\varepsilon^{\pm}} = g^{\pm}, \quad \sigma^{(i)} \nabla u^{(m)} \cdot \vec{n}|_{\partial \Omega} = \sigma^{(i)} \nabla u^{(e)} \cdot \vec{n}|_{\partial \Omega}, \end{array} \right. \quad (3.15) $$

where as always $A_m$ (m⁻¹) is the ratio of membrane area by unit volume and $S_m$ (S.m⁻²), the conductance of the cell membranes. To account for the electroporation phenomenon – corresponding to the increase of the membrane conductance – we assume that $S^{(m)}$ follows a sigmoid function of $-\nabla u^{(e)}$

$$ S^{(m)}(\|\nabla u^{(e)}\|) = S_{0}^{(m)} + \frac{1}{2}(S_{1}^{(m)} - S_{0}^{(m)}) \left[ 1 + \operatorname{erf}(k_{\text{ep}}(\|\nabla u^{(e)}\| - E_{\text{th}})) \right]. $$

In cardiac electrophysiology (see Chapter 1), the nonlinearity usually depends on the homogenised transmembrane voltage $u^{(m)} = u^{(i)} - u^{(e)}$ through the ionic term $I^{\text{ion}}$. However, the nonlinearity in $u^{(m)} = u^{(i)} - u^{(e)}$ is inconsistent with the electroporation phenomenon. Actually, assuming that $g^- = -g^+$ is constant, the symmetry along the axis ($Oy$) implies that $u^{(e)}, u^{(i)}$ and the transmembrane voltage vanishes along ($Oy$), which is inconsistent with the fact that electroporation occurs between the needles. For this reason, a nonlinearity depending on $\|\nabla u^{(e)}\|$ is chosen, which is relevant with the nonlinearity of the static nonlinear model. This choice has not been validated and it is one main perspective that I will develop in my research project (see Chapter 4) but we emphasize that this model seems more physiological than the previous phenomenological models, see Remark 6.

**Numerical simulations** The numerical discretizations that have been used to solve the different models are presented in [15]. One can remark that the strong non-linearity of the systems require the use of iteration schemes as for example the modified fixed point proposed in [BBK⁺¹⁵]. In the article, we also develop other strategies allowing to reduce the computational times.

### 3.2.1.2 Comparison on an experimental setup

We compare the previous model on a realistic situation coming from the experimental setup of Sel et al.. It consists of an ex-vivo cubic piece of rabbit liver in which two needles are inserted (Figure 3.10-Left).

Square pulses of different amplitudes are applied on the needles and the corresponding intensities are recorded. The data set consists of the measurements of the electric intensity that flows through one
<figure><img src="image_39.png" /><figcaption>Figure 3.10 – Left: Experimental setup of Sel et al. The piece of rabbit liver is located below the plexiglass. Right: Intensity measurements done at different voltages and with different size of electrodes.</figcaption></figure>

needle, say $\mathcal{E}^+$ for instance. Taking advantage of the symmetries of the liver piece, the computational domain is restricted to a quadrant with homogeneous Neumann conditions at the right, top, and bottom borders, homogeneous Dirichlet conditions at the left border, and the non-homogeneous Dirichlet conditions at the needle. The expression of the numerical electric intensity depends on the model. Denoting by $\ell_{\mathcal{E}^+}$ the length of the electrode in the perpendicular plane to the 2D simulation plane, one has

$$
\begin{align*}
I_{\mathcal{E}}^{1} &= \ell_{\mathcal{E}^{+}} \int_{\mathcal{E}^{+}} \sigma (\|\nabla u\|) \nabla u \cdot \vec{n} \, ds, && \text{where } u \text{ is the solution to (3.8)}, \\
I_{\mathcal{E}}^{2} &= \ell_{\mathcal{E}^{+}} \int_{\mathcal{E}^{+}} \sigma_{\text{eq}}(\|\mathbf{E}^{(m)}\|) \nabla u^{(e)} \cdot \vec{n} \, ds, && \text{where } (u^{(e)}, \mathbf{E}^{(m)}) \text{ is the solution to (3.13), or to (3.10)} \\
I_{\mathcal{E}}^{2,s} &= \ell_{\mathcal{E}^{+}} \int_{\mathcal{E}^{+}} \sigma_{\text{eq}}(\|\nabla u^{(e)}\|) \nabla u^{(e)} \cdot \vec{n} \, ds, && \text{where } u^{(e)} \text{ is the solution to (3.14)}, \\
I_{\mathcal{E}}^{3} &= \ell_{\mathcal{E}^{+}} \int_{\mathcal{E}^{+}} \sigma^{(e)} \nabla u^{(e)} \cdot \vec{n} \, ds, && \text{where } (u^{(e)}, u^{(m)}) \text{ is the solution to (3.15)}.
\end{align*}
$$

The above definition of the electrical intensities of the models involve the gradient of the solution along the needle. This is numerically unstable since the electric field is the most intense nearby the electrodes. To smoothen this behavior, it is possible to introduce a specific function $w$ such that the surface integral is replaced by a volume integral. More precisely, let $w$ be defined by

$$ - \Delta w = 0, \quad \text{in } \Omega, \qquad (3.16a) $$

$$ w|_{\mathcal{E}^{\pm}} = \pm 1, \quad \nabla w \cdot \vec{n}|_{\Gamma_{\text{out}}} = 0. \qquad (3.16b) $$
Then, thanks to appropriate integration by parts one has

$$I_{\varepsilon}^{1} = \frac{\ell_{\varepsilon^{+}}}{2} \int_{\Omega} \sigma (\|\nabla u\|) \nabla u \cdot \nabla w \, dx, \quad \text{where } u \text{ is the solution to (3.8)},$$

$$I_{\varepsilon}^{2} = \frac{\ell_{\varepsilon^{+}}}{2} \int_{\Omega} \sigma_{\text{eq}}(\|\mathbf{E}^{(m)}\|) \nabla u^{(e)} \cdot \nabla w \, dx, \quad \text{where } (u^{(e)}, \mathbf{E}^{(m)}) \text{ is the solution to (3.13)},$$

$$I_{\varepsilon}^{2,s} = \frac{\ell_{\varepsilon^{+}}}{2} \int_{\Omega} \sigma_{\text{eq}}(\|\nabla u^{(e)}\|) \nabla u^{(e)} \cdot \nabla w \, dx, \quad \text{where } u^{(e)} \text{ is the solution to (3.14), or to (3.10)}$$

$$I_{\varepsilon}^{3} = \frac{\ell_{\varepsilon^{+}}}{2} \int_{\Omega} \left( (\sigma^{(e)} + \sigma^{(i)}) \nabla u^{(e)} - \sigma^{(i)} \nabla u^{(m)} \right) \cdot \nabla w \, dx, \quad \text{where } (u^{(e)}, u^{(m)}) \text{ is the solution to (3.15).}$$

Figure 3.11 shows the characteristic solutions of all models with an electrode diameter of 0.7 mm and a voltage of 800 V (see [15] for the values of the parameters). The results for the standard static model (3.8) (see top row), the dynamic biphasic model (3.10) at $t = 8$ ms (see middle row), and the bidomain model (3.15) (see bottom row) are shown. Although the unknowns of the two systems are not identical, we can make a comparison by comparing the potentials (see left column), the electric flux (see middle column), and the electroporated areas (see right column). For the potential (see left column), we could consider $u(\cdot)$, the solution of (3.8), $u^{(e)}(8 \text{ ms}, \cdot)$, the solution of (3.10), and $u^{(e)}(\cdot)$, the solution of (3.15). We can see that the potentials are very similar. Regarding the electric flux, we decide to compare $\|\nabla u(\cdot)\|$, $\mathbf{E}^{(m)}$ and $\|\nabla u^{(e)}(\cdot)\|$. The results are still comparable (see middle column). Finally, one could try to define the electroporated region. It can be considered as the region where $\sigma(\|\nabla u(\cdot)\|)$ (or $\sigma^{(m)}(\|\mathbf{E}^{(m)}(8 \text{ ms}, \cdot)\|)$) reaches a high value for the standard static (resp. dynamic biphasic) model. More precisely, one could define thresholds to select the irreversible or reversible electroporated domains. For the bidomain model, the situation is quite different, since the extracellular conductivity does not depend on the electric flux. We then consider the conductance. Of course, as mentioned earlier, this highlights the need for better comparison of these tissue models, for example by asymptotic analysis. This point will be discussed in my research project. Nevertheless, we can conclude that the solutions of the different models are realistic and can be qualitatively compared. All models seem to be able to represent the electroporation phenomena at the tissue scale.

We will now focus on an important question: Can we calibrate the models using the available data? For this purpose and for simplicity, we will focus on the standard static model. The parameters are then the parameters of the sigmoid function: $E_{th}$, $\sigma_0$, $\sigma_1$, and $k_{ep}$. The data given in [SCB+05] are intensity measurements for 5 voltages and 3 electrode sizes. To estimate the parameters, we use a stochastic algorithm based on the unscented transform, see [15] for details. We first validate the estimation strategy on synthetic data in a framework close to the real data. We reach the following conclusions:

* only $E_{th}$, $\sigma_0$, and $\sigma_1$ can be estimated from the data, since there is no sensitivity associated with $k_{ep}$,

* the sensitivity of $\sigma_0$ is weaker and more iterations are needed (compared to $E_{th}$ and $\sigma_1$),

* without noise, the results are slightly better when three electrodes are considered,

* with noise, the more the data noise is important, the more it is difficult to well estimate the parameters and it becomes crucial to consider measurements performed with different electrode sizes.

As for the real data (see Figure 3.10-Right), the estimation of the parameters allows to perform simulations of the experiment and, in particular, to determine the electroporated area. Figure 3.12 illustrates the results using $\sigma(\|\nabla u\|)$ for all electrode sizes and for all voltages. These results are very encouraging as they show that we are able to estimate tissue conductivities that allow direct inference of reversible or irreversible electroporated regions based on intensity measures. This work can be considered as a first step in determining the treated zone of a patient's organ.
<figure><img src="image_40.png" /><figcaption>Figure 3.11 – Top: results for the standard static model (3.8). Middle: Results for the dynamic biphasic model (3.10) at 8 ms. Bottom: Results for the static bidomain model (3.15). Even though the unknowns of the systems are not the same, we can make a comparison by comparing the potentials (see left column), the electric flux (see middle column), and the electroporated areas (see right column).</figcaption></figure>
<figure><img src="image_41.png" /><figcaption>Figure 3.12 – $\sigma(\|\nabla u\|)$ in S/m for applied voltages 200 V and 1000 V and the three geometries for the parameters estimated on all geometries simultaneously.</figcaption></figure>

However, before we move on to clinical applications, there are some mathematical challenges. In particular, multiple electrodes are required in the clinical context. In the context of tumor ablation, the authors in [SCN<sup>+</sup>17] explain that three to six electrodes are used for inoperable hepatocellular carcinoma. They explain that this number of electrodes correlates with tumor size. In the context of cardiac ablation, multiple electrodes are included in a catheter, see, for example, the Faraware or Faraflex catheter developed by Farapulse Inc. [RDN<sup>+</sup>21] Because IRE is a multibipolar ablative technology, the electrodes are used in pairs. This means that during a pulse, two electrodes are active and the others are passive. The choice of boundary conditions for these very conductive and thin passive electrodes is a difficult question. While Dirichlet (resp. Neumann) boundary conditions model the imposed potential (resp. the current), the value of the equipotential on the inactive electrodes is an unknown of the PDE. This leads to a nonlocal constraint on the flux along the conductor interface, the so-called floating potential problem. One might also note that this problem also applies to the 3D simulation of the entire experimental setup currently under investigation, which is presented in Section 3.1.1. The next section is devoted to the presentation of our work on the floating potential problem in the context of electroporation.

### 3.2.2 Floating potential

As explained earlier, when multiple needles are inserted, the influence of the inactive electrodes on the distribution of the electric field must be carefully considered in order to accurately determine the ablation region. The focus of this section, which summarizes our article [18] – is to present an effective and rigorous way to compute the static electric field in the case of highly conductive thin inclusions.

Electroquasistatic theory states that the surface of a highly conductive conductor is an equipotential surface whose value is implicitly determined by the constraints of the problem. This is the so-called floating potential problem, which has been studied for several decades. In [ABOS14], Amann et al. have shown that the penalization method, which consists of imposing a high conductivity in the inclusion pro-
vides a less accurate electric potential than a well-designed numerical method for the floating potential problem. This result may seem strange, since the penalization is a kind of model of the real problem, while the floating potential is a perfect conductor approximation. Therefore, the question naturally arises as to how the floating potential approximates the real electric potential, whether there is a relationship between the size and the conductivity of the high conductive material that prevents the use of this approximation, and whether it is possible to increase the accuracy with an asymptotic analysis.

We first examine the observations of Amann et al. on a simplified case for which an explicit solution is available. We consider the case of a dielectric (low) conductive material Ω which is the annulus of radii $r_0 \in (0, 1)$ and 1 and with conductivity equal to 1, surrounded by a high conductive sheet $O_ε$ of thickness $ε$, and whose conductivity – after nondimensionalization – is of order $1/ε^ℓ$, where $ℓ = 1$ or 2 and $ε$ is a small parameter. The electroquasistatic potential $u_ε$ satisfies the following elliptic problem coupled with transmission conditions (TC) and boundary conditions (BC)

$$ egin{cases} 
\frac{1}{r}\partial_r(r\partial_r u_\varepsilon) + \frac{1}{r^2}\partial_\theta^2 u_\varepsilon = 0, & (\{r_0 < r < 1\} \cup \{1 < r < 1 + \varepsilon\}) \times \mathbb{R}/(2\pi\mathbb{Z}), \\ 
u_{\varepsilon|r=1-} - u_{\varepsilon|r=1+} = 0, & (\text{TC}), \\ 
\partial_r u_{\varepsilon|r=1-} - \frac{1}{\varepsilon^\ell}\partial_r u_{\varepsilon|r=1+} = 0, & (\text{TC}), \\ 
\partial_r u_{\varepsilon|r=1+\varepsilon} = 0, & (\text{BC}), \\ 
u_{\varepsilon|r=r_0} = 1 + e^{i\theta}, & (\text{BC}). 
\end{cases} 
\quad (3.17) $$

The corresponding floating potential (FP) problem consists in finding $(u, \alpha) \in H^1(\Omega) \times \mathbb{R}$ such that

$$ egin{cases} 
\frac{1}{r}\partial_r(r\partial_r u) + \frac{1}{r^2}\partial_\theta^2 u = 0, & \{r_0 < r < 1\} \times \mathbb{R}/(2\pi\mathbb{Z}), \\ 
u|_{r=1} = \alpha, & (\text{FP}), \\ 
\int_0^{2\pi} \partial_r u(1, \theta) d\theta = 0, & (\text{FP}), \\ 
u|_{r=r_0} = 1 + e^{i\theta}, & (\text{BC}). 
\end{cases} 
\quad (3.18) $$

The expression of the exact solution is given in [18]. This solution is then compared with the numerical resolution by standard second order finite difference scheme. This enables us to compare simultaneously how the solution of System (3.17) is approached by the solution of System (3.18), and how accurate is a standard second order numerical scheme for System (3.17). Numerical results are shown in Figure 3.13. Two main observations arise from these simplistic simulations.

First, the floating potential for $ℓ = 1$ does not approach the solution of System (3.17), while for $ℓ = 2$ it approximates with an accuracy of the order of $O(ε)$, which means that the floating potential cannot be used to approximate the electric potential when the ratio $R_{length}$ of the thickness of the conductor divided by the characteristic length of the dielectric is of the same order of magnitude as the ratio $R_{cond}$ of the conductivity of the dielectric divided by the conductivity of the high conductive sheet. Second, one can see that the numerical solution of System (3.17) is inaccurate when $ε$ becomes too small compared to the mesh grid. Indeed, the condition number of the matrix for the discretization grows considerably as $ε$ approaches zero. For small discretization steps, the instability increases and the numerical solution has lower accuracy, which provides an explanation for the statement of Amann et al. that the penalization method is less accurate than the floating potential.

As shown on this simplistic example, the direct resolution of the electroquasistatic problem in a domain with highly conductive inclusion leads to ill-conditioned matrix and floating potentials are preferred to avoid the computational cost. The numerical resolution of such floating potential problems has been studied for several decades. One can cite for instance the paper by Dular et al. [DHM+97], where the authors proposed a finite element method, which consists in enriching the finite elements space with spe-
<figure><img src="image_42.png" /><figcaption>Figure 3.13 – Left: Comparison of the exact solution $u_{\epsilon}^{\text{exact}}$ to (3.17) with the floating potential problem for $\ell = 1$ (blue line) and 2 (red line), and comparison of the explicit and the numerical solutions (3.17) (green line) as $\epsilon$ tends to zero. Right: Numerical convergence with the steps $(d\theta, dr) \in \{(314, 200), (628, 400), (1256, 800), (2512, 1600)\}$ of the 2nd order scheme to solve System (3.17) for $\epsilon = 1$ (blue line) and $\epsilon = 0.01$ (red line).</figcaption></figure>

cific functions defined on the nodes of the interfaces $\Gamma_k$. Amann et al. proposed in [ABOS14] a boundary element method to tackle the problem using single boundary layer integral formulation of the solution. Note that recently, a hybrid Galerkin method has been proposed by Sala et al. for a similar problem in the context of ocular hemodynamic, the electric potential being replaced by the Darcy pressure [SPGS18].

The aim of our work on floating potential is threefold. First, we give the proof of the well-posedness of the floating potential problem in the case of $N$ multiple highly conductive inclusions. Second, we propose a new numerical strategy to tackle the floating potential problem. Unlike the previous works cited above, our numerical strategy does not require any new specific numerical method. More precisely, it consists in characterizing the solution to the floating potential problem as a linear combination of $N+1$ explicit Dirichlet problems thanks to the definition of well designed Steklov-Poincaré operators. Finally, we propose an asymptotic analysis of the electroquasistatic potential in the case of a highly conductive thin inclusion, in the asymptotic regime where the ratio of the conductivities $R_{cond}$ is of order $\epsilon^2$, while the ratio $R_{length}$ is of order $\epsilon$. In particular, we prove the convergence of the asymptotic approximation at any order as $\epsilon$ goes to 0. We quickly give the results in what follows.

### 3.2.2.1 Problem setting and analysis

Let us state precisely the problem. Let $\mathcal{O}$ be a domain of $\mathbb{R}^d$, $d = 2, 3$ and let $(\mathcal{O}_k)_{k=1}^N$ be $N$ highly conductive inclusions embedded in $\mathcal{O}$. We denote by $\Gamma_{\text{out}}$ the outer boundary of $\mathcal{O}$, and by $\Gamma_k$ the boundary of $\mathcal{O}_k$ for $k = 1, \dots, N$. Define $\Omega = \mathcal{O} \setminus \cup \mathcal{O}_k$. Let $\sigma \in L^\infty(\Omega)$ be the conductivity map of $\Omega$ which satisfies for a given constant $a > 0$: $a \le \|\sigma\|_{L^\infty(\Omega)} \le 1/a$. Given $(g_k)_{k=1}^N \in \mathbb{R}$ and $f \in H^{-1}(\Omega)$, the floating potential problem<fnref n="1" /> consists in finding the $N+1$-uple $(u, \alpha_1, \dots, \alpha_N) \in H^1(\Omega) \times \mathbb{R}^N$ such that

$$ \left\{ \begin{array}{ll} -\nabla \cdot (\sigma \nabla u) & = f, & \text{in } \Omega, \\ u|_{\Gamma_{\text{out}}} & = 0, & \text{(BC)}, \\ u|_{\Gamma_k} = \alpha_k, \quad \int_{\Gamma_k} \sigma \partial_n u \, ds & = g_k, & \Gamma_k, \forall k = 1, \dots, N. \end{array} \right. \qquad (3.19) $$

<footnote>(*) value not available due to computational cost.</footnote>

<footnote n="1">¹Note that if the inclusion $\mathcal{O}_k$ is isolated, then $g_k$ is nothing but 0.</footnote>
Even though the well-posedness of System (3.19) has been addressed by Amann et al. in [ABOS14] for one inclusion, we present a variant proof for N inclusions that will lead to our simple numerical strategy. For $i = 1, \dots, N$, we consider the following Steklov-Poincaré operators defined as

$$ \Lambda_{\text{out}}^{(i)} : H^{-1}(\Omega) \to H^{-1/2}(\Gamma_i) \\ f \quad \mapsto \quad \sigma \partial_n v|_{\Gamma_i} \quad \text{s. t.} \quad \left\{ \begin{array}{l} -\nabla \cdot (\sigma \nabla v) = f \quad \text{in } \Omega, \\ v|_{\Gamma_{\text{out}}} = 0, \quad v|_{\Gamma_{\ell}} = 0, \quad \text{for } \ell = 1, \dots, N. \end{array} \right. $$

For $k = 1, \dots, N$, we define $\Lambda_k^{(i)}$ by

$$ \Lambda_k^{(i)} : H^{1/2}(\Gamma_k) \to H^{-1/2}(\Gamma_i) \\ \gamma \quad \mapsto \quad \sigma \partial_n v|_{\Gamma_i} \quad \text{s. t.} \quad \left\{ \begin{array}{l} -\nabla \cdot (\sigma \nabla v) = 0 \quad \text{in } \Omega, \\ v|_{\Gamma_i} = \gamma, \\ v|_{\Gamma_{\text{out}}} = 0, \quad v|_{\Gamma_{\ell}} = 0, \quad \text{for } \ell \neq i. \end{array} \right. $$

If it exists, the solution $(u, \alpha_1, \dots, \alpha_N)$ to System (3.19) satisfies $\sigma \partial_n u|_{\Gamma_i} = \sum_{\ell=1}^N \alpha_\ell \Lambda_\ell^{(i)}(1) + \Lambda_{\text{out}}^{(i)}(f)$, for $i = 1, \dots, N$ and the nonlocal constraints of System (3.19) read

$$ g_i = \sum_{\ell=1}^{N} \alpha_{\ell} \int_{\Gamma_i} \Lambda_{\ell}^{(i)}(1) ds + \int_{\Gamma_i} \Lambda_{\text{out}}^{(i)}(f) ds, \quad \text{for } i = 1, \dots, N. \quad (3.20) $$

By defining $\mathcal{M} = (\mathcal{M}_{ij})_{i,j=1,\dots,N}$ and $\mathcal{B} = (\mathcal{B}_i)_{i=1,\dots,N}$, Equality (3.20) can be rewritten as

$$ \mathcal{M} \begin{pmatrix} \alpha_1 \\ \vdots \\ \alpha_N \end{pmatrix} = \mathcal{B}, \quad \text{where } \mathcal{M}_{ij} = \int_{\Gamma_i} \Lambda_j^i(1) ds \quad \text{and} \quad \mathcal{B}_i = g_i - \int_{\Gamma_i} \Lambda_{\text{out}}^i(f) ds. \quad (3.21) $$

The proof – given in [18] – of the well-posedness of System (3.19) is reduced to proving the invertibility of $\mathcal{M}$.

**Proposition 7.** Let $\Omega$ be a domain of $\mathbb{R}^d$, $d = 2, 3$. Let us endow the space $(L^2(\Omega))^d$, $d = 2, 3$ with the scalar product $\langle \cdot, \cdot \rangle_{(L^2(\Omega))^d}$ defined by

$$ \langle F, G \rangle_{(L^2(\Omega))^d} = \int_{\Omega} \sigma F \cdot G \, dx, \quad \forall (F, G) \in (L^2(\Omega))^d. $$

The matrix $\mathcal{M}$ is a Gram matrix of the linearly independent vectors $\nabla v_1, \dots, \nabla v_N$ of $(L^2(\Omega))^d$, where the functions $(v_\ell)_{\ell=1}^N$ are defined by

$$ \left\{ \begin{array}{l} -\nabla \cdot (\sigma \nabla v_{\ell}) = 0 \quad \text{in } \Omega, \\ v_{\ell}|_{\Gamma_{\ell}} = 1, \\ v_{\ell}|_{\Gamma_{\text{out}}} = 0, \quad v_{\ell}|_{\Gamma_k} = 0, \quad \text{for } k \neq \ell. \end{array} \right. \quad (3.22) $$

Therefore $\mathcal{M}$ is invertible and there exists a unique $N+1$-uple $(u, \alpha_1, \dots, \alpha_N) \in H^1(\Omega) \times \mathbb{R}^N$ solution of System (3.19).
### 3.2.2.2 Numerical strategy to solve the floating potential problem

Proposition 7 leads to a simple characterization of the solution of System (3.19), and thus a simple numerical strategy. We first compute $v_\ell$ solution to (3.22) for $\ell = 1, \dots, N$ and $v_{\text{out}}$, solution to

$$ \begin{cases} -\nabla \cdot (\sigma \nabla v_{\text{out}}) = f & \text{in } \Omega, \\ v_{\text{out}}|_{\partial \Omega} = 0. \end{cases} \tag{3.23} $$

Second, we compute ($M_{ij}$) and ($B_{ij}$) given by (3.21) or equivalently

$$ M_{ij} = \int_{\Omega} \sigma \nabla v_i \cdot \nabla v_j \, dx, \quad B_i = g_i - \int_{\Omega} \sigma \nabla v_{\text{out}} \cdot \nabla v_i \, dx, $$

and deduce $(\alpha_1, \dots, \alpha_N)$ by solving the linear system (3.21). Finally, the solution $u$ to the floating potential problem (3.19) is obtained by the following linear combination

$$ u = v_{\text{out}} + \sum_{\ell=1}^{N} \alpha_{\ell} v_{\ell}, $$

meaning that to compute $u$, one just has to solve $N+1$ Dirichlet problems, which can be parallelized.

### 3.2.2.3 Asymptotic analysis

Always in [18], we have proposed an asymptotic analysis to approach accurately the solution to the electroquasistatic potential in a smooth domain with a highly conductive inclusion. We have shown that the so-called floating potential approaches the electroquasistatic potential with a first order accuracy, and we have given the expansion at any order. For the sake of simplicity, we have only considered the case where the relative thickness of the inclusion is of order $\epsilon$ and the ratio of the conductivities (the conductivity of the conductive inclusion divided by the conductivity of the domain) is of order $\frac{1}{\epsilon^2}$. Following the classical strategy:

(1) derive uniform estimates of the solution,

(2) make formal expansion (after defining local coordinates and Laplace operator),

(3) and prove the convergence,

we show that the so-called floating potential approaches the electroquasistatic potential.

### 3.2.3 Conclusion

In this section, we briefly present our work on modeling electroporation at the tissue level. We first give a comparison of the classical models available in the literature and compare them with a bidomain model that now needs to be justified mathematically. Second, we propose a mathematical strategy for dealing with inactive electrodes which is compatible with all these models. The work of this section is closely related to the work of the first section on cell modeling, while the last section, instead of modeling the electroporation phenomenon, focuses on the effects of electroporation on tumor cells.
## 3.3 Growth of tumor spheroids subjected to pulsed electric field

The aim of this section is to investigate the effects of irreversible and reversible electroporation on tumor growth. It is closely related to the work presented in Chapter 2. The aim of this work was to study the growth specificities of cancer cells in spheroids exposed to a pulsed electric field. Data were obtained by our collaborators Jelena Kolosnjaj, Muriel Golzio and Marie-Pierre Rols (IPBS, Toulouse University). Our strategy is based on an advection model (compatible with System 2.1 of Chapter 2) and well-fitted estimation strategies.

While irreversible electroporation is efficient in most cases, with a complete ablation rate of approximately 75% of treated tumor nodules after an irreversible electroporation procedure [CEH<sup>+</sup> 13], it is known that in some patients ablation is unsuccessful and the event can lead to an accelerated tumor growth. To improve our understanding of the effects of pulsed electric fields on tumor response, we use multicellular tumor spheroids as *in vitro* micro-tumor model. Multicellular tumor spheroids are 3D structures composed of cancer cells that can accurately reproduce the behavior of small solid tumors in their avascular early stages. Multicellular spheroids are becoming an important tool in cancer research as they represent an intermediate stage between 2D monolayer cell models and *in vivo* solid tumors. Spheroids resemble solid microtumors in many ways, such as heterogeneous architecture, internal gradients of signaling factors, nutrients, and oxygenation. In addition, cells within multicellular spheroids interact with each other via cell-cell interactions, and their growth kinetics resemble those observed in tumors *in vivo*. These similarities offer great potential for studying the biological properties of tumors, so they are often used for screening drugs and determining therapeutic efficacy. Such models have already demonstrated their relevance for the study of electroporation phenomena [GWTR13, FGM<sup>+</sup> 15]. The aim of this work was to characterize the regrowth of multicellular spheroids after pulse delivery as a function of pulse magnitude.

Thanks to mathematical models and a calibration strategy, we fit the model to the experimental data to evaluate the global behavior of the spheroid population. Two approaches were considered. The first is to study only the resumption of growth after exposure to the pulsed electric field using the Gompertz model (see Section 2.1.2.2 of Chapter 2). The second model is a partial differential equations as introduced in Section 2.1 of Chapter 2. Three populations of tumor cells are 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 requires additional parameters but allows us to represent the heterogeneity of the tumor cell response within the multicellular spheroids.

We begin with a presentation of the data in Section 3.3.1. In Section 3.3.2 we present a volume approach that shows encouraging preliminary results, and in Section 3.3.3 we present a 3D model that accounts for cell population heterogeneity. These first three sections correspond to the work published in [21]. Finally, in Section 3.3.4, we illustrate a population-based estimator for PDE systems that we developed ourselves (see [28]) on this system.

### 3.3.1 Biological experiments

Multicellular HCT-116-GFP spheroids were exposed to different electric field intensities – more specifically, 80 unipolar pulses of 100 µs and a frequency of 1 Hz at 0 V.cm⁻¹ (referred to as EF0) or 500 V.cm⁻¹ (referred to as EF500 group) or 1000 V.cm⁻¹ (referred to as EF1000 group) or 2000 V.cm⁻¹ (referred to as EF2000 group) – and their volumes were monitored by fluorescence and bright-field microscopy. The field magnitudes reported are the *nominal magnitudes*, i.e., the voltage between the electrodes divided by their distance. It should be noted that this is not the electric field seen by the cells. There are two control
<figure><img src="image_43.png" /><figcaption>Figure 3.14 – Left: Micrographs of multicellular spheroids followed by fluorescence microscopy of treated spheroids, after the application of electric pulses (80 unipolar pulses of 100 µs and a frequency of 1 Hz). Right: Volume evolutions of multicellular spheroids. Top: Control experiment (24 cases). Bottom: Electroporation experiment (59 cases: 13 for EF0, 16 for EF500, 14 for EF1000 and 16 for EF2000)</figcaption></figure>

groups: free growth or 0 V.cm⁻¹. In free growth condition, the spheroids are kept in growth medium. For 0 V.cm⁻¹ condition, the spheroids are transferred to the low conduction buffer for 10 minutes, following which the spheroids are rinsed in phosphate saline buffer and re-transferred to growth medium. The experiments were designed and performed by Jelena Kolosnjaj, Muriel Golzio, and Marie-Pierre Rols of the IPBS Institute at the University of Toulouse and are described in detail in [21]. Very large volume values were measured because of the osmotic shock in the cells forming the spheroids after the electric shock. Figure 3.14-Left shows for a spheroid of each of the 3 categories: EF500, EF1000, and EF2000 tumor development. Figure 3.14-Right shows the extracted volume of spheroids: free growth at the top and electroporation group at the bottom. In total 83 multicellular spheroids were tracked over 250 hours.

In [21] we first eliminate outliers and explain why the EF1000 group should be divided into two groups: in group A (6 spheroids) there are the multicellular spheroids very damaged by electric shock, whose volume of green fluorescing cells reached almost 0 mm³ after 48 hours, and in group B (5 spheroids) there are the spheroids whose dynamics are very similar to those observed in cohort EF500.

### 3.3.2 Modeling of volume evolution

**Control groups** We consider the Gompertz model already introduced in Chapter 2 for the volume evolution for the two control groups

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

where – as already said – $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
<table><caption>Table 3.1 – Estimated parameters of reversed reduced Gompertz model after electrical shock. Left: EF500. Middle: EF1000B. Right: EF1000A.</caption><thead><tr><th></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 &times; 10<sup>-5</sup></td><td>1.7 &times; 10<sup>-5</sup></td><td>2.9 &times; 10<sup>-5</sup></td></tr></tbody></table>

the initial growth rate. As with the meningioma cohort, we estimate parameters using a mixed-effects approach. Full details are provided in [21]. The population estimation strategy for the Gompertz model is applied to the 22 multicellular spheroids of the free growth cohort and the 13 multicellular spheroids EF0 of the electroporation experiment to test whether the cohorts are similar. The aim is to confirm the assumption that the electrodes have no effect on growth. As a first result, we show that the estimated parameter values are very close, implying that the electrodes do not interfere with growth (see [21]). Moreover, as already stated in Section 2.3.3 of Chapter 2, there is a strong correlation between a and b, and in [VRF+20] even a reduced Gompertz model is introduced, which consists in estimating b as before and to use $k = \frac{a}{b}$ as a constant parameter in the population. As a second result, we show in [21] that this reduced Gompertz model is valid in this context. The estimated parameters are $k = 2.18$ (mean) - 0.0 (std) and $b = 0.010$ (mean) - 0.0017 (std).

**Electroporation** 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.14-Right. 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.

Using the mean squared error (MSE) associated with each initial time considered for EF500, EF1000-B, and EF1000-A, and by checking the stability of the estimated parameters and the quality of the fits, we set the value to $t_{init} = 90h$ (resp. 90h and 43h) for cohort EF500 (resp. EF1000-B and EF1000-A). The corresponding fits are shown in Figure 3.15 and the corresponding estimated parameters in Table 3.1. As a third result, we can note that the regrowth still follows a Gompertz law and that it also seems to be an accelerated regrowth. Namely, while the estimated mean values of b are comparable to the estimated value of the control cohorts ($b = 0.010$), the estimated values of k are larger than the estimated value of the control cohorts ($k = 2.18$).
<figure><img src="image_44.png" /><figcaption>Figure 3.15 – Volume evolution after electrical shock. Left: fits using reduced Gompertz model for EF500 ($t_{init}$ = 90h). Middle: fits using reduced Gompertz model for EF1000-B ($t_{init}$ = 90h). Right: fits using reduced Gompertz model for EF1000-A ($t_{init}$ = 43h).</figcaption></figure>

Even if this volume approach shows encouraging preliminary results, this approach does not allow to properly quantify the impact of the electroporation and does not integrate the heterogeneity of the tumor, i.e. the fact that the tumor is 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.

### 3.3.3 Modeling the 3D evolution of spheroid

**Modeling** 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.

We consider an advanced mathematical model based on advection model which is compatible with System 2.1 of Chapter 2. 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 (3.24) $$

and the density of $Q$ follows

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

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 (3.26) $$

Using the saturation hypothesis $P + Q + F = 1$ in $\Omega(t)$, we have $\nabla \cdot \vec{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)} (\vec{v} \cdot \vec{n}) dX = \int_{\Omega(t)} \nabla \cdot \vec{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 [17] – 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, \quad (3.27) $$

where $\tau_G(t) = ae^{-bt}$ corresponds to the growth rate. 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. 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} \quad (3.28) $$

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** Inspiring from the work of [MFL+18], 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 (contrary to System 2.2). 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)}, $$
where R the radius of the spheroid. We 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}\vec{e}_r$. One can easily show that the volume $V$ (resp. the radius $R$) follows

$$ 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 (see the proof given in Appendix of [21]) and the fact that $V = \frac{4}{3\pi}R^3$, one can show that

$$ \left\{ \begin{aligned} R'(t) &= R(t)\tau_G(t)I(t, 1), \\ \partial_t \tilde{P}(t, r) + \tau_G(r^{-2}I(t, r) - rI(t, 1))\partial_r \tilde{P}(t, r) &= \tau_G(1 - \tilde{F})(1 - \tilde{P}) - \tilde{\tau}_{PtoQ}\tilde{P}, \\ \partial_t \tilde{F}(t, r) + \tau_G(r^{-2}I(t, r) - rI(t, 1))\partial_r \tilde{F}(t, r) &= -\tau_G(1 - \tilde{F})\tilde{F}, \\ \tilde{Q} &= 1 - (\tilde{P} + \tilde{F}). \end{aligned} \right. \qquad (3.29) $$

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.

**Parameters estimation** Six parameters need to be estimated: $V_0, k, b, p, m$ and $\lambda$ using volume evolution. We always want to use a population estimation strategy. However, the software Monolix [Ant19] cannot be used here to estimate the parameters because our model is a PDE system. We consider two different strategies to make the estimation. First in [21], we consider a Matlab library called nlmefitsa in which the SAEM algorithm is implemented. We will give the details obtained with the strategy right after. Second, we propose a population-based estimation for PDE systems in [28] and illustrate it using the PDE system and the data given here. We will develop this work in Section 3.3.4.

Very large volume values were measured because of the osmotic shock in the cells forming the spheroids after the electric shock. Because our model does not account for this phenomenon, these values are not used to estimate the parameters. Different strategies are possible: either not consider them or fix them at the minimum volume. To avoid identifiability issues, $b$ and $k$ can be estimated with the control cohort, and their mean values are used to estimate the other parameters in the 4 cases: EF500, EF1000A, EF1000B, and EF2000. Table 3.2 (first rows) gives the results obtained in this context.
<table><caption>Table 3.2 – Mean (mean) and standard deviation (std) of estimated parameters using real data. The first rows give the results obtained with the NLME algorithm using the Matlab function *nlmefitsa* and the second rows the results obtained using the population-based estimation presented in Section 3.3.4.</caption><thead><tr><th></th><th>Case</th><th>EF500 (12 spheroids)</th><th>EF1000-B (5 spheroids)</th></tr></thead><tbody><tr><td rowspan="2">λ</td><td>Function <i>nlmefitsa</i></td><td>0.22 (mean) - 0.092 (std)</td><td>0.37 (mean) - 0.012 (std)</td></tr><tr><td>Algo [28]</td><td>0.37 (mean) - 0.035 (std)</td><td>0.36 (mean) - 0.045 (std)</td></tr><tr><td rowspan="2">p</td><td>Function <i>nlmefitsa</i></td><td>0.89 (mean) - 0.014 (std)</td><td>0.82 (mean) - 0.086 (std)</td></tr><tr><td>Algo [28]</td><td>0.90 (mean) - 0.0028 (std)</td><td>0.86 (mean) - 0.012 (std)</td></tr><tr><td rowspan="2">m</td><td>Function <i>nlmefitsa</i></td><td>1.11 (mean) - 0.024 (std)</td><td>1.40 (mean) - 0.16 (std)</td></tr><tr><td>Algo [28]</td><td>1.14 (mean) - 0.18 (std)</td><td>1.21 (mean) - 0.16 (std)</td></tr></tbody></table><table><thead><tr><th></th><th>Case</th><th>EF1000-A (6 spheroids)</th><th>EF-2000 (14 spheroids)</th></tr></thead><tbody><tr><td rowspan="2">λ</td><td>Function <i>nlmefitsa</i></td><td>0.77 (mean) - 0.10 (std)</td><td>0.97 (mean) - 0.065 (std)</td></tr><tr><td>Algo [28]</td><td>0.69 (mean) - 0.04 (std)</td><td>0.92 (mean) - 0.0015 (std)</td></tr><tr><td rowspan="2">p</td><td>Function <i>nlmefitsa</i></td><td>0.17 (mean) - 0.010 (std)</td><td>0.0016 (mean) - 9×10<sup>-5</sup> (std)</td></tr><tr><td>Algo [28]</td><td>0.11 (mean) - 0.015 (std)</td><td>0.00097 (mean) - 2.1×10<sup>-6</sup> (std)</td></tr><tr><td rowspan="2">m</td><td>Function <i>nlmefitsa</i></td><td>5.19 (mean) - 1.02 (std)</td><td>x</td></tr><tr><td>Algo [28]</td><td>5.3 (mean) - 0.77 (std)</td><td>x</td></tr></tbody></table>

**Results** 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 small decrease of the volume for EF500 and EF1000-B while for the group EF1000-A (resp. EF2000), a decreasing of 75% (resp. almost 100%) is observed. Concerning $m$, the value increases with the value of the electrical field. Note that the post-pulse data for EF2000 do not really allow us to estimate the $m$ parameter, since the data for all spheroids are close to 0. Therefore, this parameter is not estimated. Figure 3.16 shows the spheroid evolutions with distribution of proliferative cells for the 4 cohorts (column 2 to column 5) compared to the control cohort (column 1).

**Discussion** The use of a PDE system makes it possible to determine the percentage of cells destroyed and the percentage of cells whose functioning was altered by the effect of electroporation for each value of the electrical field. It was seen that at an electrical field close to $1000 \text{ V.cm}^{-1}$ (partially irreversible electroporation and partially reversible electroporation), the increase of the electrical field seemed to boost the resumption of spheroid growth. Our results provide good insight into the effects of electroporation on tumor growth. When irreversible electroporation is used in clinics to ablate tumors, the treatment may induce cell death by irreversible electroporation. This effect is observed in a specific area, namely the area where the electric field strength is highest – at the optimal location with respect to the electrodes. The effectiveness of fully irreversible electroporation was well observed in the spheroid group exposed to $2000 \text{ V.cm}^{-1}$. Nevertheless, next 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 surviving cells resume growth. This effect was clearly present in the multicellular spheroids exposed to $1000 \text{ V.cm}^{-1}$. Thus, this field amplitude produced two radically different behaviors: either death or growth boost. We hypothesize that this variable response rate is due to the inhomogeneity of the electric field. If the spheroid (or tumor) was not reached efficiently, we have a relapse. If the positioning of the spheroid
<figure><img src="image_45.png" /><figcaption>Figure 3.16 – 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, column 4: EF1000A and column 5: EF2000.</figcaption></figure>

(or tumor) is in the center of the electrodes, the treatment is efficient. This is a recurrent problem even in clinics, where physicians observe a different response rate – either efficient treatment of a patient with a decreased tumor or a tumor relapse [SNdJ+14]. The advantage of our model is that we can determine the trend of tumor growth shortly after electroporation. So, in practice, this model can be used to determine which tumors have been treated inefficiently and that the patient is on the path to relapse. This would alert the treating physician that the treatment was not efficient and that an additional electroporation cycle is needed to remove the tumor again. As for the growth boost, it is particularly intriguing and certainly requires further investigation from a biological, physical, and mathematical perspective.

### 3.3.4 Population-based estimation for PDE systems

The estimation of the parameters given in Table 3.2 (first rows) requires several hours for each of the cohorts even if the cohorts are small (< 15 spheroids) and the number of parameters is small (= 3) and it is a 1D PDE system. It is possible that the Matlab function used for estimation could be optimized, but it remains a challenge for a mixed-effects strategy to increase the size of the underlying dynamics – for
example, when using PDE systems – with an acceptable complexity cost to the algorithm. To our knowledge, there is only one work that aims to develop a strategy for using population estimation strategies for PDE systems. In [GLV14], the PDE system is solved by interpolations of solutions precomputed on a carefully constructed mesh of the parameter domain.

In [28], we propose a formulation of a population-based estimator for PDE systems to address this challenge. The errors in the initial conditions are controlled by a Luenberger observer that corrects for the state, and the parameters are estimated using a population reduced-order Kalman-based filter restricted to the parametric space. This sophisticated method for solving this difficult problem is based on two important methodological strategies: (1) the population-based Kalman filter introduced in [23], (2) the joint state-parameter estimation introduced in [MCLT08]. More precisely, the population-based Kalman filter [23] captures all population information in a unified maximum likelihood estimation procedure assuming Gaussian disturbances. However even the proposed reduced version of this population-based Kalman filter – defined by combining clusters subpopulations with common observational backgrounds – results in the famous curse of dimensionality [Bel56] that makes it numerically intractable for partial differential equations (PDE). Assuming that the dimension of the parameters is smaller than the dimension of the system state, which is then responsible for the curse of dimensionality, one can restrict the population-based Kalman filter to parametric space and couple it with a Luenberger observer following the strategy presented in [MC11b]. This strategy assumes that a suitable Luenberger observer has been previously designed to control the trajectory.

The first section of [28] is devoted to the formalization of the estimation problem and the introduction of the population-based estimator and we refer to the paper for this part. The second section applies the complete strategy on spheroid electroporation. The problem, model, and data are the same as in Section 3.3.3. We first introduce and analyze a Luenberger observer, see Section 3.3.4.1 and we second perform numerical experiments on synthetic data to validate the approach, see Section 3.3.4.2. The choice of the 1D PDE model stems from the fact that the computational cost of solving the PDE system is still reasonable and allows to perform many replicates – 7500 simulations in total.

#### 3.3.4.1 A Luenberger observer for System (3.29)

We design the following Luenberger-type observer

$$ egin{cases} 
\hat{R}'(t) = \hat{R}(t)\tau_G(t)\hat{I}(t, 1) - \gamma_{\text{obs}}(t)(\hat{R} - R), \\ 
\partial_t \hat{P}(t, r) + \tau_G(r^{-2}\hat{I}(t, r) - r\hat{I}(t, 1))\partial_r \hat{P}(t, r) = \tau_G(1 - \hat{F})(1 - \hat{P}) - \hat{\tau}_{PtoQ}\hat{P}, \\ 
\partial_t \hat{F}(t, r) + \tau_G(r^{-2}\hat{I}(t, r) - r\hat{I}(t, 1))\partial_r \hat{F}(t, r) = -\tau_G(1 - \hat{F})\hat{F}, \\ 
\hat{Q} = 1 - (\hat{P} + \hat{F}). 
\end{cases} 
\quad (3.30) $$

with

$$ \hat{I}(t, r) = \int_{0}^{r} (1 - \hat{F}) \underline{r}^{2} d\underline{r}, \quad \text{and} \quad \hat{\tau}_{PtoQ}(t, r) = \tilde{\tau}_{PtoQ}^{\text{begin}} - \frac{\tilde{\tau}_{PtoQ}^{\text{begin}} - \tilde{\tau}_{PtoQ}^{\text{end}}}{1 + e^{\frac{(R(t)(1-r)-d)}{s}}}. $$

The function $\gamma_{\text{obs}}$ is the gain function. We assume that the function goes from $[0, T]$ to $]0, \infty[$ and is bounded i.e. there exists $\gamma_{\text{obs}}^{\min}, \gamma_{\text{obs}}^{\max} \in ]0, \infty[$ such that $\gamma_{\text{obs}}^{\min} \le \gamma_{\text{obs}}(t) \le \gamma_{\text{obs}}^{\max}$, for all $t \in [0, T]$. Initial conditions are given by

$$ \hat{R}(0) = \hat{R}_0, \quad \hat{P}(0, r) = \hat{P}_0 \in [0, 1], \quad \hat{Q}(0, r) = 1 - \hat{P}_0, \quad \hat{F}(0, r) = 0, \quad \forall r. $$
In this part, we only consider uncertainties reduced to the initial condition. The values of the parameters are fixed to the true values. We admit here that there exists unique solutions to System 3.29 and System 3.30 and these solutions are smooth (this can be proved by a Banach contraction mapping theorem). The proof of the following proposition is given in [28].

**Proposition 8.** We have

$$ (1) \quad \forall t \ge 0, \forall r \in [0, 1], U(t, r) \in [0, 1], \text{ for } U = \{\tilde{P}, \tilde{Q}, \tilde{F}, \hat{P}, \hat{Q}, \hat{F}\} \text{ and } (\tilde{P} + \tilde{Q})(t, r), (\hat{P} + \hat{Q})(t, r) \in [0, 1], $$

$$ (2) \quad \forall t \ge 0, \tilde{R}(t) \in \left[ \min(\tilde{R}(0), \tilde{R}(t_{as})), \tilde{R}(0)e^{\frac{ma_{bs}}{3b}} \right]. $$

In order to analyze the convergence of the observer we will study the estimation error defined by $\tilde{R} = R - \hat{R}$ and $\tilde{P} = P - \hat{P}$. One can prove the two following propositions (see [28] for the proofs and for numerical illustrations).

**Proposition 9.** If $\gamma_{\text{obs}}^{\min} > \frac{ma_{bs}}{3}$, the radius $t \mapsto \tilde{R}(t)$ converges exponentially to 0 when $t$ goes to $+\infty$.

**Proposition 10.** If $\gamma_{\text{obs}}^{\min} > \frac{ma_{bs}}{3} + T_O^{\text{end}}$ and $\tilde{P}(0, \cdot) \in H^s([0, r])$, the norm $t \mapsto \|\tilde{P}(t, \cdot)\|_{L^2([0,1])}^2$ converges exponentially to 0 when $t$ goes to $+\infty$.

### 3.3.4.2 State and parameter estimation

We would like to test the ability of our algorithm to estimate parameters. First, a synthetic data benchmark is considered to validate our implementation. To do so, we run different simulation series with different numbers of spheroids denoted by $N_S$ and with different measurement errors – i.e. with different values of the standard deviation of the additive Gaussian noise always denoted by $\sigma_{err}$ – for 6 scenarios: strong priors with $\gamma_{\text{obs}} = 0$ (str. pr.), strong priors with $\gamma_{\text{obs}} > 0$ (str. pr. - so), weak priors with $\gamma_{\text{obs}} = 0$ (w. pr.), weak priors with $\gamma_{\text{obs}} > 0$ (w. pr. - so), weak priors and false initial conditions with $\gamma_{\text{obs}} = 0$ (w. pr., f. IC), weak priors and false initial conditions with $\gamma_{\text{obs}} > 0$ (w. pr., f. IC - so). To limit the energy consumption of this validation, we will focus on estimating only two parameters: $b$ and $a_{bs}$ considering the model in the free growth case, i.e. without electroporation. More precisely, instead of estimating $a_{bs}$, we will estimate $k = \frac{a_{bs}}{b}$ which is assumed to be fixed in the population when using the reduced Gompertz model introduced in [VRF+20]. To validate our algorithm, we consider several statistical validation criteria whose the definitions are given in [28]: the relative bias (rBIAS), the relative mean squared error (rMSE), the mean estimated standard deviation (STD), the empirical standard deviation of estimates, defined as the standard deviation of estimates at convergence (ESTD) and the coverage (COV). For information on how to generate data, please refer to [28]. We generate $N_R = 100$ simulation replicates for $N_S = 1, 5, 10, 20$ and $\sigma_{err} = 0, 0.04, 0.08$ for the 6 scenarios. We also generate $N_R = 100$ simulation replicates for $N_S = 40$ but only for the scenario w. pr., f. IC - so.

To illustrate how the algorithm works, Figure 3.17 shows the time evolution of the states and parameters of two spheroids belonging to a cohort of 10 spheroids and $\sigma_{err} = 0.08$. The generated data and the two selected spheroids are indicated (top, left). The second line corresponds to Spheroid 2 and the third line to Spheroid 6. In both cases, the first column corresponds to radius, the second to proliferative and quiescent densities, and the third to parameter $b$. The parameter $k$ – estimated fixed in the population – is given in the upper right corner. The target states and target values of the parameters are shown in blue. The blue points correspond to the noisy observed data. The states and parameter values of the observer are shown in yellow. The dashed lines correspond to the 95% confidence interval when available (i.e. for
<figure><img src="image_46.png" /><figcaption>Figure 3.17 – Time evolution of states and parameters for two spheroids with $N_S = 10$ and $\sigma_{err} = 0.08$. The target states and target values of the parameters are shown in blue. The blue points in the graphs of the radius correspond to the noisy observed data. The states and parameter values with population estimation are shown in yellow. The dashed lines correspond to the 95% confidence interval. First line, left: generated data of a cohort of the 10 spheroids. The two selected spheroids are in dark blue. First line, right: evolution of parameter *k* (assumed fixed in the cohort). Second (resp. third) line: Spheroid 2 (resp. 6) with time evolution state of radius on the left, proliferative and quiescent cells at a time step in the middle and parameter *b* on the right.</figcaption></figure>

the parameters only). One can see that the yellow lines converge to blue lines over time, showing the ability of our algorithm to estimate well.

Figure 3.18-Left shows the rMSE, the rBIAS, and the COV concatenated for both parameters *b* and *k* and the STD, the ESTD, and the BMIXED for parameter *b* against the number of spheroids considered $N_S$ for a standard deviation error $\sigma_{err}$ of 0.08 for the 6 scenarios. First, one can see that, as expected, the results are better when strong priors and true initial conditions are considered. Second, using the state observer in conjunction with the Kalman observer for the parameters leads to better results when weak priors are considered and especially when the initial conditions are false. Since Scenario w. pr., f. IC - so is the most interesting, we will focus on it. Figure 3.18-Right shows the rMSE, the rBIAS, and the COV concatenated for both parameters *b* and *k* and the STD, the ESTD and the BMIXED for parameter *b* against the considered errors. As expected, increasing the number of spheroids gives better results and increasing the errors worsens them, although it can be seen that they are still acceptable.
<figure><img src="image_47.png" /><figcaption>Figure 3.18 – Performances evaluated over 100 replicates for the 6 scenarios when $\sigma_{err} = 0.08$ (left) and for Scenario w. pr., f. IC - so (right). First lines: rMSE (log-scale), rBIAS and COV concatenated for both parameters b and k as a function of either the errors (left) or the number of spheroids (right) considered. Second line: STD, ESTD and BMIXED for parameter b (k is fixed in the population) as a function of either the errors (left) or the number of spheroids (right) considered.</figcaption></figure>

<table><caption>Table 3.3 – Computation times given in minutes when estimating 2 parameters (b and k) in Scenario w. pr., f. IC - so. Comparison with the Matlab function *nlmefitsa* average of computation times for 10 replicates).</caption><thead><tr><th>Algo</th><th>Δ<sub>x</sub></th><th>N<sub>S</sub> = 1</th><th>N<sub>S</sub> = 5</th><th>N<sub>S</sub> = 10</th><th>N<sub>S</sub> = 20</th><th>N<sub>S</sub> = 40</th></tr></thead><tbody><tr><td>Algorithm [28]</td><td>0.05</td><td>0.03</td><td>0.3</td><td>1</td><td>4</td><td>12</td></tr><tr><td>Algorithm [28]</td><td>0.01</td><td>x</td><td>x</td><td>3.5</td><td>x</td><td>x</td></tr><tr><td>Algorithm [28]</td><td>0.001</td><td>x</td><td>x</td><td>60</td><td>x</td><td>x</td></tr><tr><td>Function nlmefitsa</td><td>0.05</td><td>x</td><td>80</td><td>x</td><td>x</td><td>x</td></tr></tbody></table>
These estimates are supplemented by computational times obtained using a 2.3 GHz quad-core Intel Core i7 computer and an implementation in MATLAB-2020-b, see Table 3.3. This gives an order of magnitude numerical complexity of the algorithm, albeit with a not fully optimized implementation. In particular, we give the computation times when no parallelization of the code is performed, although it is very easy to parallelize the loop on sigma points. Our strategy must solve $N_σ (= N_S × N_θ + 1$ using simplex sigma-points) times a time step of the PDE system at each time step. Therefore – even if the estimation part was reduced to the parameters with the Kalman filter, which means, that the dimension of the covariance matrix goes to $N_S × (N_θ + N_{Δx}) = N_S × (2 + 2 × \frac{1}{Δ_x} + 1) = 43N_S$ to $N_S N_θ = 2N_S$ – the computation times are always strongly related to the size of $N_{Δx}$. To illustrate this, we also shrink the spatial step $Δ_x$ for $N_S = 10$. We can see that the computation times for a non-optimized algorithm are still reasonable and this shows the possibility of applying it to 2D or 3D PDE systems.

Finally, for comparison, we use the Matlab function *nlmefitsa*, which is an implementation of the nonlinear mixed effects model (NLME) with the stochastic EM algorithm, but only for $N_S = 5$ because of the high computational cost. We report the mean of the computation times for 10 replicates. The rMSE – equal to 0.0067 (mean for only 10 replicates) – is very close to the rMSE – equal to 0.0074 (mean for 100 replicates) – obtained with our algorithm in exactly the same context.

Our approach shows to be successful on synthetic data. It is clear that using the state observer in conjunction with the Kalman observer for the parameters leads to better results when weak priors are considered and especially when the initial conditions are false. This implies that we recommend the use of a state observer especially when the confidence in the initial conditions is low. We are aware that building a state observer to handle errors in the initial conditions can be a very difficult step and is the main drawback of our strategy. We also perform estimation on real data – see second rows of Table 3.2 – and our results compare very well with the estimation procedures used in the literature as standard for estimating parameters of nonlinear mixed-effects models represented here through *nlmefistsa* function.

Having demonstrated the concept of our strategy and these performances in terms of accuracy, manageability and computation times, the main perspective of this work is to develop a code that can be easily used by the community, see more details in my research project given in Chapter 4.

## Conclusion

This chapter presents my contributions to electroporation modeling, an innovative therapy. More specifically, I present the work done to develop mathematical approaches to track membrane dynamics, the work done to compare and improve tissue-scale electroporation models, and then how we are beginning to investigate the effects of irreversible and reversible electroporation on tumor growth. Although our work to date has enabled great progress in the understanding of electroporation phenomenon, much work remains to be done, and I will present my future work in this area in my research project given in the next chapter.
<figure><img src="image_48.png" /></figure>
# Chapter 4

## Research project

In this very short chapter, some perspectives and future work are given. I plan, at least in the short and medium term, to continue working on mathematical models for biomedical problems.

### 4.1 Tumor growth modeling

I intend to continue work around tumor growth modeling, focusing on three main projects.

**Free growth of solid tumors** The first project is on free growth of solid tumors, which was studied as part of V. Montalibet's PhD (funded by the University of Bordeaux as part of the ANR project AI by UBx).

First, we will continue our work with our collaborators at the CHU of Bordeaux (J. Engelhardt and H. Loiseau, neurosurgeons) on meningiomas. The open clinical questions are: Can we segment the tumor fully automatically? Can we use the first medical image to determine clustering obtained with mathematical models? Can we perform 3D simulations within minutes? To answer at least the first two questions, we will use machine and deep learning strategies. We have recently gained this expertise in the MONC team by expanding our knowledge in these areas. To answer clinical questions, it is often essential to combine different mathematical branches.

Second, a new collaboration with IHU of Strasbourg (A. Venkatasamy, radiologist) has been started, dealing with acoustic neuroma, also called vestibular schwannoma. This is a noncancerous and usually slow-growing tumor that arises on the main nerve leading from the inner ear to the brain. The clinical questions are similar to meningiomas: Can we segment the tumor fully automatically? Can we model the volume dynamics? Can we classify schwannomas according to their aggressiveness? Can we predict the volume and shape of the tumor at a later time point with only two imaging examinations? This new project will test whether our strategies can be easily transferred to other contexts.

**Dynamics of circulating tumor cells in the blood vessels** As part of C. Etchegaray's collaboration with the Goetz Laboratory of Tumor Biomechanics at the University of Strasbourg, we are interested in the dynamics of circulating tumor cells (CTC) in blood vessels. This is a crucial phenomenon in the spread of metastases, which are responsible for about 90% of deaths in cancer patients. Indeed, some CTC manage to escape hydrodynamic forces and adhere to vessel walls before extravasating and further invading tissues and organs. Therefore, a deeper understanding of CTC dynamics will help to find new
therapeutic targets against invasion in the future. Our goal is to better understand the competition be-
tween hydrodynamic forces and cell adhesion forces and their effects on CTC circulation. To this end, we
aim to calibrate a mathematical model of cell dynamics using experimental data from *in vitro* (microflu-
dics) and *in vivo* (zebrafish embryos). This work will help provide insight into the biomechanical drivers
of CTC invasion. It is funded by Inria through the post-doc of G. Ciavolella.

**Effect of the drug propranolol on tumors** Propranolol belongs to a group of drugs called beta-blockers. Although it is mainly used to treat heart problems, anxiety, and to prevent migraines, it is also used in oncology. More specifically, propranolol is currently the first-line treatment for severe infantile hemangiomas, the most common childhood tumor affecting 1-5% of newborns [SdlRB<sup>+09</sup>]. The molecular mechanism of action of propranolol and the cellular and molecular targets of the drug are still being explored [MOKC<sup>+21</sup>]. Recently I started a new collaboration with a biologist from the University of Bordeaux (F. Moisan). The goal is to use 3D spheroids to better understand the antitumor mechanism of action of this beta-blocker and to identify new predictive markers of response. This challenge is very important because this new therapeutic strategy could also be useful for the treatment of other vascular tumors, especially angiosarcomas, which are rare malignant tumors with poor prognosis.

## 4.2 Electroporation modeling

As far as electroporation modeling is concerned, three main directions will be investigated in the com-
ing years. Most of this work will be done in a manner consistent with cardiac electroporation. One of the
major treatments for cardiac arrhythmias is catheter ablation, in which small areas of cardiac tissue are
destroyed to isolate or eliminate the cause of the rapid and irregular heartbeats. Most catheter ablation
therapies are performed thermally through the application of radiofrequency electromagnetic field (RFA),
and pulsed electric field ablation (PFA), which takes advantage of irreversible electroporation, is still in
its infancy. Preclinical evaluations of PFA in atrial fibrillation [KKI<sup>+19a</sup> and ventricular ablation in large
animal studies [COJ<sup>+20b</sup>] show successful results with possible transmural lesions, sparing vulnerable
adjacent structures. In terms of clinical performance, the recurrence of atrial fibrillation is on the order
of 15% [RDN<sup>+21</sup> with PFA and 30% with RFA [WN06]. Despite these promising results, there are still
some difficulties in applying PFA in routine clinical practice because of the technical complexity of this
novel approach. As in oncology, mathematical models and numerical strategies could be developed to
improve the understanding of PFA and to develop numerical criteria for treatment assessment thanks
to clinical data. This area of research offers me the opportunity to combine two of my main areas of
application: cardiac electrophysiology and electroporation modeling.

**Membrane modeling compatible with cardiac cells** The cell-level modeling presented in Section 3.1.2 of Chapter 3 is not compatible with cardiac cells because changes in ion and adenosine triphosphate (ATP) concentrations and cell volume were not considered. Therefore, the PDE system needs to become considerably more complex to incorporate other aspects of the electroporation phenomenon and adapt it to cardiac cells. In particular, the increase in cell membrane permeability when short high-voltage pulses are applied probably leads to cardiac cell death by Ca2+ uptake. To quantify what happens in cardiac cells under PFA – apoptosis or necrosis – the model should describe the *in vitro* process of internalization of extracellular species. This new model will additionally incorporate the formation of conduction defects, ion and molecular efflux and influx, and cell swelling. Since electroneutrality can be temporarily broken near the membrane due to ion exchange, the electroosmosis model should be
preferred. This means coming back to the Nernst-Planck equation corresponding to System 1.4 given in Chapter 1, and coupling it to the Stokes equation

$$-\eta \Delta_{\underline{x}} \underline{v}_s + \vec{\nabla}_{\vec{x}} P_s = q \sum_{k=1}^{N^{(\alpha)}} z_k c_k^{(\alpha)} \vec{\nabla}_{\vec{x}} u^{(\alpha)} \text{ with } \vec{\nabla}_{\vec{x}} \cdot \underline{v}_s = 0.$$

In this context, the velocity $\underline{v}_s$ does not correspond to the mechanical deformation, but to the velocity of the fluid. The viscosity is denoted by $\eta$. The transmission through the cell membrane has not yet been derived for this complete model, especially in the context of electroporation, but our goal is to couple it with the membrane equation of state $\phi$ and define a realistic electroporation current.

The goal is to include the currents of passive ion channels, Na/K-ATPase and calcium pumps, and Na/H channels, which are essential for the equilibrium of the cell. One may refer to the electrophysiological laws as in [Arm03, FH04, PSC+11]. In addition, specific nonlinearities must be introduced in the jump conditions to account for the various degrees of membrane permeabilization. We assume that we can derive such jump conditions thanks to the knowledge of electroporation acquired through classical electroporation modeling [KLPW14, LSMP14] and with a rigorous asymptotic analysis in the same vein as [PP13]. Our approach then generalizes the standard electroquasistatic models [LSMP14] by describing the ion and molecular fluxes and integrating the morphological properties of the cardiac cells through fitted geometries.

This work is funded by the ANR project MIRE4VTach, which I coordinate. A PhD student will be hired, and biological experiments to validate our modeling will be performed by our collaborators at IHU Lyric in Bordeaux (G. Caluori, cardiac arrhythmia researcher and P. Jaïs, cardiologist).

### A macroscopic bidomain model for electroporation

As explained in Chapter 3, System 3.15 was not formally derived by a homogenization procedure. Our goal is to use 2-scale convergence to derive a dynamic bidomain equation for electroporation that is compatible with cardiac cells i.e. accounting for ionic current. Preliminary results show that an effective model whose unknown is the sum of the first terms of the asymptotic development of the potentials could be a solution to this difficult problem, which is due to the scale difference between ion and electroporation currents. Of course, this work is valid even in the absence of ionic current, which means that the model is also suitable for electroporation of other organs and should make it possible either to justify System (3.15) or to propose a mathematically justified alternative.

Once this model is available, we would like to perform numerical simulations for ventricular tachycardia. Ventricular tachycardia (VTach) – common in both men and women – generates disorder of rapid heartbeat due to abnormal electrical signals in the ventricles. We will study the effects of cardiac geometry (wall thickness and fibrosis) and fibers orientation using a few patients. Regarding the catheter, one of the biggest challenges is the distance between the tissue and the catheter. Using CT or MRI images, we will consider the most realistic configuration and add appropriate variability to see the impact on the ablated area. The signal data recorded during the ablation procedure by our clinical collaborators will be used to derive the electroporation protocol (i.e. find out which electrodes are used for each pulse), to validate the model, and to estimate its parameters such as the conductivities. This last part is a challenge in terms of data assimilation.

This work is funded by the Fédération Française de Cardiologie and Inria through the PhD thesis of S. Bihoreau. Access to real data from IHU Lyric will be provided in the context of the ANR project MIRE4VTach. Interestingly, we should obtain depolarization maps, and we could use the estimation
strategy presented in Section 1.4.3 of Chapter 1 to complete them.

**Bidomain model of a cardiac domain containing a region ablated by PFA** In this third part, we aim to derive an electrophysiological model of a cardiac domain containing a region ablated by PFA. To this end, we aim to combine the bidomain model in the healthy parts of the heart with an adapted bidomain system in the electroporated region. Our main goal is to propose a numerical explanation for the higher rate of fibrillation recurrence after RFA compared with PFA. For PFA, the model results from an asymptotic analysis in which the small parameter is the thickness of the electroporated region and its low intracellular conductivity tensor. Preliminary results were obtained and published in the following conference proceeding [27].

## 4.3 Data assimilation aspects

I do not focus on this work in this manuscript, but an important part of my work is to develop efficient strategies for estimating the parameters of ODE or PDE systems. More specifically, together with P. Moireau (DR Inria) and M. Prague (CR Inria), we have developed a population-based Kalman filter in [23]. An extension to PDE systems using a coupling with a Luenberger observer to handle state errors was proposed in [28]. Important work has been done to validate the approach, and we have demonstrated its potential in two real-world and difficult examples: (1) to evaluate associations between governmental nonpharmaceutical interventions – such as lockdowns, curfews, or school closures – introduced to control the pandemic COVID-19 and transmission rates without prior knowledge of the dynamics of the transmission rates [22] (see Figure 2 in the introduction chapter for illustration) and (2) for very efficient estimation of the parameters of a 1D PDE system [28] (see Section 3.3.4.2 of Chapter 2). Since the proof of concept shows very interesting results, we would now like with my two collaborators to disseminate the method. This means that we need to develop a code and an associated interface that is easy to use, and very accurate documentation. We want to have code that is capable of solving ODE systems, and that is also compatible with any existing ODE or PDE solver. This is still under consideration, but a medium term solution is to seek funding to hire an engineer for this project.
My publications (peer-review journals and proceedings)

[1] D. Chapelle, A. Collin, and J.-F. Gerbeau. A surface-based electrophysiology model relying on asymptotic analysis and motivated by cardiac atria modeling. *M3AS*, 23(14):2749–2776, 2013. [Link](#).

[2] A. Collin, J.-F. Gerbeau, M. Hocini, M. Haïssaguerre, and D. Chapelle. Surface-based electrophysiology modeling and assessment of physiological simulations in atria. In *International Conference on Functional Imaging and Modeling of the Heart*, volume 7945, pages 352–359, 2013. [Link](#).

[3] D. Chapelle and A. Collin. Strong convergence results for the asymptotic behavior of the 3D-shell model. *Journal of Elasticity*, 115(2):173–192, 2014. [Link](#).

[4] A. Collin, D. Chapelle, and P. Moireau. A Luenberger observer for reaction–diffusion models with front position data. *Journal of Computational Physics*, 300:288–307, 2015. [Link](#). See Section 1.4.1 of Chapter 1 for details.

[5] A. Collin, D. Chapelle, and P. Moireau. Sequential state estimation for electrophysiology models with front level-set data using topological gradient derivations. In *International Conference on Functional Imaging and Modeling of the Heart*, pages 402–411. Springer, 2015. [Link](#). See Section 1.4.2 of Chapter 1 for details.

[6] E. Schenone, A. Collin, and J.-F. Gerbeau. Numerical simulation of electrocardiograms for full cardiac cycles in healthy and pathological conditions. *International Journal for Numerical Methods in Biomedical Engineering*, 2015. [Link](#).

[7] A. Collin, G. Sangalli, and T. Takacs. Analysis-suitable G1 multi-patch parametrizations for C1 isogeometric spaces. *Computer Aided Geometric Design*, 47:93–113, 2016. [Link](#).

[8] A. Collin and S. Imperiale. Mathematical analysis and 2-scale convergence of a heterogeneous microscopic bidomain model. *Mathematical Models and Methods in Applied Sciences*, 28(05):979–1035, 2018. [Link](#). See Section 1.2 of Chapter 1 for details.

[9] M. C. Rochoux, A. Collin, C. Zhang, A. Trouvé, D. Lucor, and P. Moireau. Front shape similarity measure for shape-oriented sensitivity analysis and data assimilation for eikonal equation. *ESAIM: Proceedings and Surveys*, 63:258–279, 2018. [Link](#).

[10] N. Tarabelloni, E. Schenone, A. Collin, F. Ieva, A. M. Paganoni, and J.-F. Gerbeau. Statistical assessment and calibration of numerical ECG models. *JP Journal of Biostatistics*, 15(2):151–173, 2018. [Link](#).
[11] A. Collin, S. Imperiale, P. Moireau, J.-F. Gerbeau, and D. Chapelle. Apprehending the effects of mechanical deformations in cardiac electrophysiology: A homogenization approach. *Mathematical Models and Methods in Applied Sciences*, 29(13):2377-2417, 2019. Link. See Sections 1.1 and 1.3 of Chapter 1 for details.

[12] A. Gérard, A. Collin, G. Bureau, P. Moireau, and Y. Coudière. Model assessment through data assimilation of realistic data in cardiac electrophysiology. In *International Conference on Functional Imaging and Modeling of the Heart*, pages 121-130. Springer, 2019. Link. See Section 1.4.3 of Chapter 1 for details.

[13] C. Zhang, A. Collin, P. Moireau, A. Trouvé, and M. C. Rochoux. Front shape similarity measure for data-driven simulations of wildland fire spread based on state estimation: Application to the rxcadre field-scale experiment. *Proceedings of the Combustion Institute*, 37(3):4201-4209, 2019. Link.

[14] C. Zhang, A. Collin, P. Moireau, A. Trouvé, and M. C. Rochoux. State-parameter estimation approach for data-driven wildland fire spread modeling: Application to the 2012 rxcadre s5 field-scale experiment. *Fire Safety Journal*, 105:286-299, 2019. Link.

[15] G. Jankowiak, C. Taing, C. Poignard, and A. Collin. Comparison and calibration of different electro- poration models. Application to rabbit livers experiments. *ESAIM: Proceedings and Surveys*, 67:242- 260, 2020. Link. See Section 3.2.1 of Chapter 3 for details.

[16] D. Voyer, S. Corridore, A. Collin, R. Scorretti, and C. Poignard. Numerical modeling of floating potentials in electrokinetic problems using an asymptotic method. *IEEE Transactions on Magnetics*, 56(4):1-4, 2020. Link.

[17] 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. Link. See Sections 2.1.2 and 2.2 of Chapter 2 for details.

[18] A. Collin, S. Corridore, and C. Poignard. Floating potential boundary condition in smooth domains in an electroporation context. In *Methods of Mathematical Oncology: Fusion of Mathematics and Biology, Osaka, Japan, October 26-28, 2020*, pages 91-106. Springer, 2021. Link. See Section 3.2.2 of Chapter 3 for details.

[19] A. Collin, V. Groza, L. Missenard, F. Chomy, T. Colin, J. Palussière, and O. Saut. A model-strengthened imaging biomarker for survival prediction in EGFR-mutated non-small-cell lung carcinoma patients treated with tyrosine kinase inhibitors. *Bulletin of Mathematical Biology*, 83(6):68, 2021. Link.

[20] A. Collin, T. Kritter, C. Poignard, and O. Saut. Joint state-parameter estimation for tumor growth model. *SIAM Journal on Applied Mathematics*, 81(2):355-377, 2021. Link. See Section 2.3 of Chapter 2 for details.

[21] 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. Link. See Section 3.3 of Chapter 3 for details.
[22] A. Collin, B. P. Hejblum, C. Vignals, L. Lehot, R. Thiébaut, P. Moireau, and M. Prague. Using population based Kalman estimator to model COVID-19 epidemic in France: estimating the effects of non-pharmaceutical interventions on the dynamics of epidemic. *International Journal of Biostatistics*, 2022. Link. See Figure 2 of the introduction chapter for the main result in terms of application.

[23] A. Collin, M. Prague, and P. Moireau. Estimation for dynamical systems using a population-based Kalman filter–Applications in computational biology. *MathematicS In Action*, 2022. Link.

[24] J. Engelhardt, V. Montalibet, H. Loiseau, O. Saut, and A. Collin. Evaluation of four tumour growth models to describe the natural history of meningiomas. *EBioMedicine*, 2023. Link. See Section 2.2 of Chapter 2 for details.

[25] P. Jaramillo-Aguayo, A. Collin, and C. Poignard. Phase-field model of bilipid membrane electroporation. *Journal of Mathematical Biology*, 87(1):18, 2023. Link. See Section 3.1.2 of Chapter 3 for details.

[26] P. Loubet, A. Vincent, A. Collin, C. Dejous, A. Ghiotto, and C. Jego. Life cycle assessment of ICT in higher education: a comparison between desktop and single-board computers. *The International Journal of Life Cycle Assessment*, pages 1–19, 2023. Link.

[27] S. Nati Poltri, G. Caluori, P. Jaïs, A. Collin, and C. Poignard. Electrocardiology modeling after catheter ablations for atrial fibrillation. In *International Conference on Functional Imaging and Modeling of the Heart*. Springer, 2023. Link.

[28] A. Collin. Population-based estimation for PDE systems – Applications in spheroids electroporation. *ESAIM: Control, Optimisation and Calculus of Variations*, 2023 (in press). Link. See Section 3.3.4 of Chapter 3 for details.

[29] A. Collin, T. García-Sánchez, S. Corridore, L. M. Mir, and C. Poignard. Deciphering immediate post-pulse membrane resealing from 4-electrode impedance measurements by numerical modeling. *Bioelectricity*, 2023 (in press). Link. See Section 3.1.1 of Chapter 3 for details.
<figure><img src="image_49.png" /></figure>
Bibliography

[AAC+79] I. G. Abidor, V. B. Arakelyan, L. V. Chernomordik, Y. A. Chizmadzhev, V. F. Pastushenko, and M. P. Tarasevich. Electric breakdown of bilayer lipid membranes. *Journal of Electro-analytical Chemistry and Interfacial Electrochemistry*, 104:37–52, jan 1979.

[ABE16] S. Alonso, M. Bär, and B. Echebarria. Nonlinear physics of electrical wave propagation in the heart: a review. *Reports on Progress in Physics*, 79(9):096601, 2016.

[ABOS14] D. Amann, A. Blaszczyk, G. Of, and O. Steinbach. Simulation of floating potentials in industrial applications by boundary element methods. *Journal of Mathematics in Industry*, 4(1):13, 2014.

[AC79] S. M. Allen and J. W. Cahn. A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening. *Acta Metallurgica*, 27(6):1085–1095, jun 1979.

[ACFS00] L. Ambrosio, P. Colli Franzone, and G. Savaré. On the asymptotic behaviour of anisotropic energies arising in the cardiac biodomain model. *Interfaces and Free Boundaries*, 2(3):213–266, 2000.

[ADH95] G. Allaire, A. Damlamian, and U. Hornung. Two-scale convergence on periodic surfaces and applications. *In Proceedings of the International Conference on Mathematical Modelling of Flow through Porous Media*, 1995.

[AHM08] M. Alfaro, D. Hilhorst, and H. Matano. The singular limit of the Allen–Cahn equation and the FitzHugh–Nagumo system. *Journal of Differential Equations*, 245(2):505–565, jul 2008.

[AHZ+13] I. Adeniran, J. C. Hancox, H. Zhang, et al. Effect of cardiac ventricular mechanical contraction on the characteristics of the ECG: A simulation study. *Journal of Biomedical Science and Engineering*, 6(12):47, 2013.

[All92] G. Allaire. Homogenization and two-scale convergence. *SIAM Journal of Mathematical Analysis*, 23(6):1482–1518, 1992.

[All12] G. Allaire. *Shape optimization by the homogenization method*, volume 146. Springer Science & Business Media, 2012.

[AM10] F. M. Andre and L. M. Mir. Nucleic acids electrotransfer in vivo: mechanisms and practical aspects. *Current gene therapy*, 10(4):267–280, 2010.

[Ant19] Antony, France: Lixoft SAS. Monolix version 2019R1, 2019.
[AP96] R. Aliev and A. Panfilov. A simple two-variable model of cardiac excitation. *Chaos, Solitons Fractals*, 3(7):293-301, 1996.

[Arm03] C. M. Armstrong. The Na/K pump, Cl ion, and osmotic stabilization of cells. *Proceedings of the National Academy of Sciences*, 100(10):6257-6262, 2003.

[BAB$^{+12}$] S. Benzekry, N. André, A. Benabdallah, J. Ciccolini, C. Faivre, F. Hubert, and D. Barbolosi. Modeling the impact of anticancer agents on metastatic spreading. *Mathematical Modelling of Natural Phenomena*, 7(1):306-336, 2012.

[BBB$^{+15}$] E. Baratchart, S. Benzekry, A. Bikfalvi, T. Colin, L. S. Cooley, R. Pineau, E. J. Ribot, O. Saut, and W. Souleyreau. Computational modelling of metastasis development in renal cell carcinoma. *PLoS computational biology*, 11(11):e1004626, 2015.

[BBK$^{+15}$] M. Breton, F. Buret, L. Krähenbühl, M. Leguèbe, L. M. Mir, R. Perrussel, C. Poignard, R. Scorretti, and D. Voyer. Non-linear steady-state electrical current modeling for the electropermeabilization of biological tissue. *IEEE Transactions on Magnetics*, 51(3):1-4, 2015.

[BBRRH$^{+19}$] F. Brugada-Bellsolà, P. T. Rodríguez, A. Rodríguez-Hernández, R. Garcia-Armengol, M. Tardáguila, A. Gonz'alez-Crespo, C. J. Domínguez, and J. M. Rimbau. Growth prediction in asymptomatic meningiomas: the utility of the aimss score. *Acta Neurochirurgica*, 161(11):2233-2240, nov 2019.

[BCF$^{+10}$] M. Boulakia, S. Cazeau, M.A. Fernández, J.-F. Gerbeau, and N. Zemzemi. Mathematical modeling of electrocardiograms: a numerical study. *Annals of Biomedical Engineering*, 38(3):1071-1097, 2010.

[BCG$^{+10}$] D. Bresch, T. Colin, E. Grenier, B. Ribba, and O. Saut. Computational modeling of solid tumor growth: the avascular stage. *SIAM Journal on Scientific Computing*, 32(4):2321-2344, 2010.

[BCP94] Y. Bourgault, Y. Coudière, and C. Pierre. A collocation-galerkin finite element model of cardiac action potential propagation. *IEEE Transactions on Biomedical Engineering*, 41(8):743-757, 1994.

[BCP09] Y. Bourgault, Y. Coudière, and C. Pierre. Existence and uniqueness of the solution for the bidomain model used in cardiac electrophysiology. *Nonlinear Analysis-Real World Applications*, 10(1):458-482, 2009.

[BDS93] J.-Z. Bao, C. C. Davis, and R. E. Schmukler. Impedance spectroscopy of human erythrocytes: system calibration, and nonlinear modeling. *IEEE Transactions on biomedical engineering*, 40(4):364-378, 1993.

[BDS$^{+12}$] M. Breton, L. Delemotte, A. Silve, L. M. Mir, and M. Tarek. Transport of siRNA through lipid membranes driven by nanosecond electric pulses: an experimental and computational study. *Journal of the American Chemical Society*, 134(34):13938-13941, 2012.

[BEB$^{+18}$] I. Baldi, J. Engelhardt, C. Bonnet, L. Bauchet, E. Berteaud, A. Grüber, and H. Loiseau. Epidemiology of meningiomas. *Neurochirurgie*, 64(1):5-14, 2018.

[Bel56] R. Bellman. Dynamic programming and Lagrange multipliers. *Proceedings of the National Academy of Sciences of the United States of America*, 42(10):767, 1956.
[Ben71] A. Bensoussan. *Filtrage optimal des systèmes linéaires*. Dunod, 1971.

[BGA<sup>+</sup>11] I. Baldi, A. Gruber, A. Alioum, E. Berteaud, P. Lebailly, A. Huchet, T. Tourdias, G. Kantor, JP. Maire, A. Vital, et al. Descriptive epidemiology of cns tumors in france: results from the gironde registry for the period 2000–2007. *Neuro-oncology*, 13(12):1370–1378, 2011.

[BHR04] M. Burger, B. Hackl, and W. Ring. Incorporating topological derivatives into level set methods. *Journal of Computational Physics*, 194(1):344–362, 2004.

[BLB<sup>+</sup>14] S. Benzekry, C. Lamont, A. Beheshti, A. Tracz, J. ML. Ebos, L. Hlatky, and P. Hahnfeldt. Classical mathematical models for description and prediction of experimental tumor growth. *PLoS computational biology*, 10(8):e1003800, 2014.

[BLP78] A. Bensoussan, J.-L. Lions, and G. Papanicolaou. *Asymptotic Analysis for Periodic Structures*. Elsevier, 1978.

[BOCF08] A. Bueno-Orovio, E. M. Cherry, and F. H. Fenton. Minimal model for human ventricular action potentials in tissue. *Journal of Theoretical Biology*, 253:544–560, 2008.

[BP11] M. Bishop and G. Plank. Bidomain ecg simulations using an augmented monodomain model for the cardiac source. *Biomedical Engineering, IEEE Transactions on*, 58(8):2297–2307, 2011.

[Bra02] A. J. Bray. Theory of phase-ordering kinetics. *Advances in Physics*, 51(2):481–587, mar 2002.

[Bri93] M. Briane. Three models of non periodic fibrous materials obtained by homogenization. *RAIRO-Modélisation mathématique et analyse numérique*, 27(6):759–775, 1993.

[BRS<sup>+</sup>09] F. Billy, B. Ribba, O. Saut, H. Morre-Trouilhet, T. Colin, D. Bresch, J.-P. Boissel, E. Grenier, and J.-P. Flandrois. A pharmacologically based multiscale mathematical model of angio-genesis and its use in investigating the efficacy of a new cancer treatment strategy. *Journal of theoretical biology*, 260(4):545–562, 2009.

[BV97] J. B. Braunstein and N. A. Vick. Meningiomas: the decision not to operate. *Neurology*, 48(5):1459–1462, 1997.

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

[CCJ<sup>+</sup>15] T. Colin, F. Cornelis, J. Jouganous, J. Palussière, and O. Saut. Patient-specific simulation of tumor growth, response to the treatment, and relapse of a lung metastasis: a clinical case. *Journal of Computational Surgery*, 2(1):1, 2015.

[CCSZ95] Y. A. Chizmadzhev, F. S. Cohen, A. Shcherbakov, and J. Zimmerberg. Membrane mechanics can account for fusion pore dilation in stages. *Biophysical Journal*, 69(6):2489–2500, dec 1995.

[CEH<sup>+</sup>13] R. Cannon, S. Ellis, D. Hayes, G. Narayanan, and R. CG. 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.
[CFPS14] P. Colli Franzone, L.F. Pavarino, and S. Scacchi. *Mathematical Cardiac Electrophysiology*, volume XIV of MS&A. Springer, 2014.

[CFPS16] P. Colli Franzone, L.F. Pavarino, and S. Scacchi. Bioelectrical effects of mechanical feedbacks in a strongly coupled cardiac electro-mechanical model. *Mathematical Models and Methods in Applied Sciences*, 26(01):27-57, 2016.

[CFPS17] P. Colli Franzone, L. F. Pavarino, and S. Scacchi. Effects of mechanical feedback on the stability of cardiac scroll waves: A bidomain electro-mechanical simulation study. *Chaos: An Interdisciplinary Journal of Nonlinear Science*, 27(9):1-15, 2017.

[CFPS18] P. Colli Franzone, L. F. Pavarino, and S. Scacchi. A numerical study of scalable cardiac electro-mechanical solvers on hpc architectures. *Frontiers in Physiology*, 9:1-16, 2018.

[CFPT05] P. Colli Franzone, L.F. Pavarino, and B. Taccardi. Simulating patterns of excitation, repolarization and action potential duration with cardiac Bidomain and Monodomain models. *Mathematical Biosciences*, 197(1):35-66, 2005.

[CFS02] P. Colli Franzone and G. Savaré. Degenerate evolution systems modeling the cardiac electric field at micro and macroscopic level. *Progress in Nonlinear Differential Equations and Their Applications*, 50:49-78, 2002.

[CH58] J. W. Cahn and J. E. Hilliard. Free energy of a nonuniform system. I. Interfacial free energy. *The Journal of Chemical Physics*, 28(2):258-267, feb 1958.

[CMLT09] D. Chapelle, P. Moireau, and P. Le Tallec. Robust filtering for joint state parameter estimation for distributed mechanical systems. *Discrete and Continuous Dynamical Systems-Series A*, 23(1-2):65-84, 2009.

[COJ$^{+}$20a] G Caluori, E. Odehnalova, T. Jadczyk, M. Pesl, I. Pavlova, L. Valikova, S. Holzinger, V. Novotna, V. Rotrekl, A. Hampl, et al. AC pulsed field ablation is feasible and safe in atrial and ventricular settings: a proof-of-concept chronic animal study. *Frontiers in bioengineering and biotechnology*, 8:552357, 2020.

[COJ$^{+}$20b] G. Caluori, E. Odehnalova, T. Jadczyk, M. Pesl, I. Pavlova, L. Valikova, S. Holzinger, V. Novotna, V. Rotrekl, A. Hampl, et al. AC pulsed field ablation is feasible and safe in atrial and ventricular settings: a proof-of-concept chronic animal study. *Frontiers in bioengineering and biotechnology*, 8:552357, 2020.

[Cou06] W. T. Couldwell. Asymptomatic meningiomas. *Journal of neurosurgery*, 105(4):536-537, 2006.

[CRN98] M. Courtemanche, R.J. Ramirez, and S. Nattel. Ionic mechanisms underlying human atrial action potential properties: insights from a mathematical model. *American Journal of Physiology*, (275):H301-H321, 1998.

[CV91] T.F. Chan and L.A. Vese. Active contours without edges. *IEEE Transactions on Image Processing*, 10(2):266-277, 1991.

[DHM$^{+}$97] P. Dular, F. Henrotte, B. Meys, A. Genon, and W. Legros. Une méthode naturelle de traitement des potentiels flottants associée à la méthode des éléments finis. *Journal de Physique III*, 7(11):2201-2209, 1997.
[DK98] K. A. DeBruin and W. Krassowska. Electroporation and shock-induced transmembrane potential in a cardiac fiber during defibrillation strength shocks. *Annals of Biomedical Engineering*, 26(4):584–596, July 1998.

[DMR05] R. V. Davalos, L. M. Mir, and B. Rubinsky. Tissue ablation with irreversible electroporation. *Annals of biomedical engineering*, 33(2):223, 2005.

[Don21] E. Donada. La moindre efficacité du vaccin contre la transmission remet-elle en cause le principe du pass sanitaire ? *Libération*, 2021.

[ESS92] L. C. Evans, H. M. Soner, and P. E. Souganidis. Phase transitions and generalized motion by mean curvature. *Communications on Pure and Applied Mathematics*, 45(9):1097–1123, oct 1992.

[Eve09] G. Evensen. *Data Assimilation: The Ensemble Kalman Filter*. Springer, 2009.

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

[FH04] J. A. Fraser and C. L.-H. Huang. A quantitative analysis of cell volume and resting potential determination and regulation in excitable cells. *The Journal of Physiology*, 559(2):459–478, 2004.

[Fit61] R. Fitzhugh. Impulses and physiological states in theoretical models of nerve membrane. *Biophysical Journal*, 1(6):445–466, 1961.

[FMZB20] D. Fraggedakis, M. Mirzadeh, T. Zhou, and M. Z. Bazant. Dielectric breakdown by electric-field induced phase separation. *Journal of The Electrochemical Society*, 167(11):113504, 2020.

[FPP+16] M. Favino, S. Pozzi, S. Pezzuto, F.W. Prinzen, A. Auricchio, and R. Krause. Impact of mechanical deformation on pseudo-ECG: a simulation study. *EP Europace*, 18(4):77–84, 2016.

[FS98] E. C. Fear and M. A. Stuchly. Modeling assemblies of biological cells exposed to electric fields. *IEEE Transactions on Biomedical Engineering*, 45(10):1259–1271, 1998.

[Geh08] J. Gehl. Electroporation for drug and gene delivery in the clinic: doctors go electric. *Electroporation Protocols*, pages 351–359, 2008.

[Gér19] A. Gérard. Modèles numériques personnalisés de la fibrillation auriculaire. PhD thesis, Université de Bordeaux, 2019.

[GLV14] E. Grenier, V. Louvet, and P. Vigneaux. Parameter estimation in non-linear mixed effects models with SAEM algorithm: extension from ODE to PDE. *ESAIM: Mathematical Modelling and Numerical Analysis*, 48(5):1303–1329, 2014.

[Gom25] 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*, (115):513–583, 1825.
[GSAL<sup>+15</sup>] T. García-Sánchez, A. Azan, I. Leray, J. Rosell-Ferrer, R. Bragos, and L. M. Mir. Interpulse multifrequency electrical impedance measurements during electroporation of adherent differentiated myotubes. *Bioelectrochemistry*, 105:123-135, 2015.

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

[GXCI<sup>+19</sup>] A. Gupta, Z. Xu, O. Cohen-Inbar, M. H. Snyder, L. K. Hobbs, C. Li, Q.-T. Nguyen, and J. P. Sheehan. Treatment of asymptomatic meningioma with gamma knife radiosurgery: Long-term follow-up with volumetric assessment and clinical outcome. *Neurosurgery*, 2019.

[HA76] J. E. Hoke and R. A. Anthes. The initialization of numerical models by a dynamic-initialization technique. *Monthly Weather Review*, 104(12):1551-1556, 1976.

[HAS07] H. L. P. Harpold, E. C. Alvord, and K. R. Swanson. The Evolution of Mathematical Modeling of Glioma Proliferation and Invasion. *J Neuropathol Exp Neurol*, 66(1):1-9, January 2007.

[Hec21] C. Hecketsweiler. Covid-19: face à la quatrième vague, les effets trop tardifs de l'accélération de la vaccination, selon l'Institut Pasteur. *Le Monde*, 2021.

[HH52] A.L. Hodgkin and A.F. Huxley. A quantitative description of membrane current and its application to conduction and excitation in nerve. *Journal of Physiology*, 117(4):500-544, 1952.

[HHI<sup>+09</sup>] T. Hashiba, N. Hashimoto, S. Izumoto, T. Suzuki, N. Kagawa, M. Maruno, A. Kato, and T. Yoshimine. Serial volumetric assessment of the natural history and growth pattern of incidentally discovered meningiomas. *Journal of neurosurgery*, 110(4):675-684, 2009.

[HKO07] L. He, C.Y. Kao, and S. Osher. Incorporating topological derivatives into shape derivatives based level set methods. *Journal of Computational Physics*, 225(1):891-909, 2007.

[HL09] M. Hintermüller and A. Laurain. Multiphase image segmentation and modulation recovery based on shape and topological sensitivity. *Journal of Mathematical Imaging and Vision*, 35(1):1-22, 2009.

[HM05] S. N. Healy and A. D. McCulloch. An ionic model of stretch-activated and stretch-modulated currents in rabbit ventricular myocytes. *Europace*, 7(s2):S128-S134, 2005.

[HP10] P. E. Hand and C. S. Peskin. Homogenization of an electrophysiological model for a strand of cardiac myocytes with gap-junctional and electric-field coupling. *Bulletin of Mathematical Biology*, 72:1408-1424, 2010.

[Hui98] G. Huiskamp. Simulation of depolarization in a membrane-equations-based model of the anisotropic ventricle. *Biomedical Engineering, IEEE Transactions on*, 45(7):847-855, 1998.

[IMM<sup>+19a</sup>] A. I. Islim, M. Mohan, R. DC. Moon, N. Srikandarajah, S. J. Mills, A. R. Brodbelt, and M. D. Jenkinson. Correction to: Incidental intracranial meningiomas: a systematic review and meta-analysis of prognostic factors and outcomes. *Journal of neuro-oncology*, 144(2):427-429, 2019.

[IMM<sup>+19b</sup>] A. I. Islim, M. Mohan, R. DC. Moon, N. Srikandarajah, S. J. Mills, A. R. Brodbelt, and M. D. Jenkinson. Incidental intracranial meningiomas: a systematic review and meta-analysis of prognostic factors and outcomes. *Journal of neuro-oncology*, 142(2):211-221, 2019.
[ITC<sup>+13</sup>] P. B. Ishai, M. S. Talary, A. Caduff, E. Levy, and Y. Feldman. Electrode polarization in dielectric measurements: a review. *Measurement Science and Technology*, 24(10):102001, 2013.

[IVM10] A. Ivorra, J Villemejane, and L. M. Mir. Electrical modeling of the influence of medium conductivity on electroporation. *Physical chemistry chemical physics : PCCP*, 12:10055-64, 2010.

[JKK<sup>+11</sup>] K.-W. Jo, C.-H. Kim, D.-S. Kong, H.-J. Seol, D.-H. Nam, K. Park, J.-H. Kim, and J.-I. Lee. Treatment modalities and outcomes for asymptomatic meningiomas. *Acta neurochirurgica*, 153(1):62-67, 2011.

[JOV<sup>+11</sup>] D. R. Johnson, J. E. Olson, R. A. Vierkant, J. E. Hammack, A. H. Wang, A. R. Folsom, B. A. Virnig, and J. R. Cerhan. Risk factors for meningioma in postmenopausal women: results from the Iowa Women's Health Study. *Neuro-oncology*, 13(9):1011-1019, 2011.

[JU97] S.J. Julier and J.K. 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.

[JU04] S. J. Julier and J. K. Uhlmann. Unscented filtering and nonlinear estimation. *Proceedings of the IEEE*, 92(3):401-422, 2004.

[KB61] R. Kalman and R. Bucy. New results in linear filtering and prediction theory. *Trans. ASME J. Basic. Eng.*, 83:95-108, 1961.

[KCM<sup>+10</sup>] E. Konukoglu, O. Clatz, B. H. Menze, B. Stieltjes, M.-A. Weber, E. Mandonnet, H. Delingette, and N. Ayache. Image guided personalization of reaction-diffusion type tumor growth models using modified anisotropic eikonal equations. *Medical Imaging, IEEE Transactions*, 29(1):77-95, 2010.

[Kee91] J. P. Keener. An eikonal-curvature equation for action potential propagation in myocardium. *Journal of Mathematical Biology*, 29:629-651, 1991.

[Kir10] B. J. Kirby. *Micro- and Nanoscale Fluid Mechanics: Transport in Microfluidic Devices*. Cambridge University Press, 2010.

[KKC<sup>+18</sup>] K. H. Kim, S. J. Kang, J.-W. Choi, D.-S. Kong, H. J. Seol, D.-H. Nam, and J.-I. Lee. Clinical and radiological outcomes of proactive Gamma Knife surgery for asymptomatic meningiomas compared with the natural course without intervention. *Journal of neurosurgery*, 130(5):1740-1749, 2018.

[KKI<sup>+19a</sup>] J. Koruth, K. Kuroki, J. Iwasawa, Y. Enomoto, R. Viswanathan, R. Brose, E. D. Buck, M. Speltz, S. R. Dukkipati, and V. Y. Reddy. Preclinical evaluation of pulsed field ablation: electrophysiological and histological assessment of thoracic vein isolation. *Circulation: Arrhythmia and Electrophysiology*, 12(12):e007781, 2019.

[KKI<sup>+19b</sup>] J. Koruth, K. Kuroki, J. Iwasawa, Y. Enomoto, R. Viswanathan, R. Brose, ED. Buck, M. Speltz, SR. Dukkipati, and VY. Reddy. Preclinical evaluation of pulsed field ablation: electrophysiological and histological assessment of thoracic vein isolation. *Circulation: Arrhythmia and Electrophysiology*, 12(12):e007781, 2019.
[KL05] E. Kuhn and M. Lavielle. Maximum likelihood estimation in nonlinear mixed effects models. *Computational Statistics & Data Analysis*, 49(4):1020-1038, 2005.

[KLPW14] O. Kavian, M. Leguèbe, C. Poignard, and L. Weynans. "Classical" electropromeabilization modeling at the cell scale. *Journal of mathematical biology*, 68(1-2):235-265, 2014.

[KNP<sup>+</sup>04] W. Krampla, S. Newrkla, W. Pfisterer, S. Jungwirth, P. Fischer, T. Leitha, W. Hruby, and K. H. Tragl. Frequency and risk factors for meningioma in clinically healthy 75-year-old patients: Results of the transdanube ageing study (VITA). *Cancer: Interdisciplinary International Journal of the American Cancer Society*, 100(6):1208-1212, 2004.

[KNP08] R. H. Keldermann, M. P. Nash, and A. V. Panfilov. Modeling cardiac mechano-electrical feedback using reaction-diffusion-mechanics systems. *Physica D*, 238:1000-1007, 2008.

[KP14] D. Kondepudi and I. Prigogine. *Modern thermodynamics: from heat engines to dissipative structures*. John Wiley & Sons, 2014.

[Kri18] T. Kritter. *Utilisation de données cliniques pour la construction de modèles en oncologie*. PhD thesis, Bordeaux, 2018.

[KS04] J. Keener and J. Sneyd. *Mathematical Physiology*. Springer, 2004.

[KSW<sup>+</sup>07] D. Keller, G. Seemann, D. Weiss, D. Farina, J. Zehelein, and O. Dössel. Computer based modeling of the congenital long-QT 2 syndrome in the visible man torso: From genes to ECG. In *Engineering in Medicine and Biology Society, 2007. EMBS 2007. 29th Annual International Conference of the IEEE*, pages 1410-1413, 2007.

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

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

[Lav14] M. Lavielle. *Mixed effects models for the population approach: models, tasks, methods and tools*. CRC press, 2014.

[LBG<sup>+</sup>03] G.T. Lines, M.L. Buist, P. Grottum, A.J. Pullan, J. Sundnes, and A. Tveito. Mathematical models and numerical methods for the forward problem in cardiac electrophysiology. *Comput. Visual. Sci.*, 5(4):215-239, 2003.

[LDAW18] K. A. Layne, P. I. Dargan, J. RH. Archer, and D. M. Wood. Gadolinium deposition and the potential for toxicological sequelae-A literature review of issues surrounding gadolinium-based contrast agents. *British journal of clinical pharmacology*, 84(11):2522-2534, 2018.

[Lev44] K. Levenberg. A method for the solution of certain non-linear problems in least squares. *Quarterly of applied mathematics*, 2(2):164-168, 1944.

[LKP<sup>+</sup>17] E. J. Lee, J. H. Kim, E. S. Park, Y.-H. Kim, J. K. Lee, S. H. Hong, Y. H. Cho, and C. J. Kim. A novel weighted scoring system for estimating the risk of rapid growth in untreated intracranial meningiomas. *Journal of neurosurgery*, 127(5):971-980, 2017.

[LL08] S. Lakshmivarahan and J.M. Lewis. Nudging methods: A critical overview. In S.K. Park and L. Xu, editors, *Data Assimilation for Atmospheric, Oceanic, and Hydrologic Applications*, volume XVIII. Springer, 2008.
[Loo65] H. Looyenga. Dielectric constants of heterogeneous mixtures. *Physica*, 31(3):401-406, 1965.

[LPPK17] E. J. Lee, J. H. Park, E. S. Park, and J. H. Kim. "Wait-and-see" strategies for newly diagnosed intracranial meningiomas based on the risk of future observation failure. *World neurosurgery*, 107:604-611, nov 2017.

[LSMP14] M. Leguebe, A. Silve, L. M. Mir, and C. Poignard. Conducting and permeable states of cell membrane submitted to high voltage pulses: mathematical and numerical studies validated by the experiments. *Journal of theoretical biology*, 360:83-94, 2014.

[Lue63] D. G. Luenberger. Determining the State of a Linear with Observers of Low Dynamic Order. PhD thesis, Stanford University, 1963.

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

[Lvo12] V. F. Lvovich. *Impedance spectroscopy: applications to electrochemical and dielectric phenomena*. John Wiley & Sons, 2012.

[LW82] N. M. Laird and J. H. Ware. Random-effects models for longitudinal data. *Biometrics*, pages 963-974, 1982.

[Mar63] D. W. Marquardt. An algorithm for least-squares estimation of nonlinear parameters. *Journal of the society for Industrial and Applied Mathematics*, 11(2):431-441, 1963.

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

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

[MCLT08] P. Moireau, D. Chapelle, and P. Le Tallec. Joint state and parameter estimation for distributed mechanical systems. *Computer Methods in Applied Mechanics and Engineering*, 197:659-677, 2008.

[MCLT09] P. Moireau, D. Chapelle, and P. Le Tallec. Filtering for distributed mechanical systems using position measurements: Perspectives in medical imaging. *Inverse Problems*, 25(3):035010 (25pp), 2009.

[MCZ+17] M. H. Mohammad, E. Chavredakis, R. Zakaria, A. Brodbelt, and M. D. Jenkinson. A national survey of the management of patients with incidental meningioma in the United Kingdom. *British journal of neurosurgery*, 31(4):459-463, 2017.

[MDFG12] V. Martin, A. Drochon, O. Fokapu, and J-F. Gerbeau. Magnetohemodynamics in the aorta and electrocardiograms. *Physics in Medicine and Biology*, 57:3177-3195, 2012.

[MFL+18] T. Michel, J. Fehrenbach, V. Lobois, 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.
[Mir01] L. M. Mir. Therapeutic perspectives of in vivo cell electropermeabilization. *Bioelectrochemistry*, 53(1):1-10, 2001.

[MLR05] L. Miller, J. Leor, and B. Rubinsky. Cancer cells ablation with irreversible electroporation. *Technology in cancer research & treatment*, 4(6):699-705, 2005.

[MNS+22] C. Merla, M. Nardoni, M. Scherman, S. Petralito, L. Caramazza, F. Apollonio, M. Liberti, P. Paolicelli, B. Attal-Tretout, and L. M. Mir. Changes in hydration of liposome membranes exposed to nanosecond electric pulses detected by wide-field Coherent anti-Stokes Raman microspectroscopy. *Bioelectrochemistry*, 147:108218, 2022.

[MOBP91] L. M. Mir, S. Orlowski, J. Belehradek, and C. Paoletti. Electrochemotherapy potentiation of antitumour effect of bleomycin by local electric pulses. *European Journal of Cancer and Clinical Oncology*, 27(1):68-72, 1991.

[MOKC+21] 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.

[Mor09] Y. Mori. From three-dimensional electrophysiology to the cable model: an asymptotic study. *arXiv preprint arXiv:0901.3914*, 2009.

[MP95] J. Malmivuo and R. Plonsey. *Bioelectromagnetism - Principles and Applications of Bioelectric and Biomagnetic Fields*. Oxford University Press, 1995.

[MS83] J. J. Moré and D. C. Sorensen. Computing a trust region step. *SIAM Journal on scientific and statistical computing*, 4(3):553-572, 1983.

[MS03] C. C. Mitchell and D. G. Schaeffer. A two-current model for the dynamics of cardiac membrane. *Bulletin Math. Bio.*, 65:767-793, 2003.

[MSB+12] D. Miklavčič, G. Serša, E. Brecelj, J. Gehl, D. Soden, G. Bianchi, P. Ruggieri, C. R. Rossi, LG. Campana, and T. Jarm. Electrochemotherapy: technological advancements for efficient electroporation-based treatment of internal tumors. *Medical & biological engineering & computing*, 50(12):1213-1225, 2012.

[NAY62] J. S. Nagumo, S. Arimoto, and S. Yoshizawa. An active pulse transmission line stimulating nerve axon. *Proceedings of the IEEE*, 50:2061-2071, 1962.

[NFF+98] A. Nygren, C. Fiset, L. Firek, J.W. Clark, D.S. Lindblad, R.B. Clark, and W.R. Giles. Mathematical model of an adult human atrial cell: the role of K+ currents in repolarization. *Circulation Research*, 82(1):63-81, 1998.

[NFN+05] S. Nakasu, T. Fukami, M. Nakajima, K. Watanabe, M. Ichikawa, and M. Matsuda. Growth pattern changes of meningiomas: long-term analysis. *Neurosurgery*, 56(5):946-955, 2005.

[Ngu89] G. Nguetseng. A general convergence result for a functional related to the theory of homogenization. *SIAM Journal on Mathematical Analysis*, 20(3):608-623, 1989.

[NK93] J. C. Neu and W. Krassowska. Homogenization of syncytial tissues. *Critical Reviews in Biomedical Engineering*, 21(2):137-199, 1993.
[NK99] J. C. Neu and W. Krassowska. Asymptotic model of electroporation. *Physical Review E*, 59(3):3471-3482, mar 1999.

[NNF<sup>+</sup>11] S. Nakasu, Y. Nakasu, T. Fukami, J. Jito, and K. Nozaki. Growth curve analysis of asymptomatic and symptomatic meningiomas. *Journal of neuro-oncology*, 102(2):303-310, 2011.

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

[NP04] M. P. Nash and A. V. Panfilov. Electromechanical model of excitable tissue to study reentrant cardiac arrhythmias. *Progress in Biophysics and Molecular Biology*, 85(2-3):501-522, 2004.

[NRM<sup>+</sup>03] M. Nakamura, F. Roser, J. Michel, C. Jacobs, and M. Samii. The natural history of incidental meningiomas. *Neurosurgery*, 53(1):62-71, 2003.

[OGF<sup>+</sup>13] Q. T. Ostrom, H. Gittleman, P. Farah, A. Ondracek, Y. Chen, Y. Wolinsky, N. E. Stroup, C. Kruchko, and J. S. Barnholtz-Sloan. CBTRUS statistical report: Primary brain and central nervous system tumors diagnosed in the united states in 2006-2010. *Neuro-oncology*, 15(suppl 2):ii1-ii56, 2013.

[OKSL11] S. Oya, S.-H. Kim, B. Sade, and J. H. Lee. The natural history of intracranial meningiomas. *Journal of neurosurgery*, 114(5):1250-1256, 2011.

[OVVR11] T. O'Hara, L. Virág, A. Varró, and Y. Rudy. Simulation of the undiseased human cardiac ventricular action potential: Model formulation and experimental validation. *PLoS Comput Biol*, 7(5):e1002061, 2011.

[PDG03] M. Potse, B. Dubé, and R.M. Gulrajani. ECG simulations with realistic human membrane, heart, and torso models. In *Engineering in Medicine and Biology Society, 2003. Proceedings of the 25th Annual International Conference of the IEEE*, volume 1, pages 70-73, 2003.

[PDV09] M. Potse, B. Dubé, and A. Vinet. Cardiac anisotropy in boundary-element models for the electrocardiogram. *Medical & biological engineering & computing*, 47(7):719-729, 2009.

[PKN05] A. V. Panfilov, R. H. Keldermann, and M. P. Nash. Self-organized pacemakers in a coupled reaction-diffusion-mechanics system. *Physical review letters*, 95(25):258104, 2005.

[PKW<sup>+</sup>07] A. G. Pakhomov, J. F. Kolb, J. A. White, R. P. Joshi, S. Xiao, and K. H. Schoenbach. Long-lasting plasma membrane permeabilization in mammalian cells by nanosecond pulsed electric field (nsPEF). *Bioelectromagnetics*, 28(8):655-663, 2007.

[PP13] R. Perrussel and C. Poignard. Asymptotic expansion of steady-state potential in a high contrast medium with a thin resistive layer. *Applied Mathematics and Computation*, 221:48-65, 2013.

[PSC<sup>+</sup>11] C. Poignard, A. Silve, F. Campion, L. M. Mir, O. Saut, and L. Schwartz. Ion fluxes, trans-membrane potential, and osmotic stabilization: a new dynamic electrophysiological model for eukaryotic cells. *European Biophysics Journal*, 40(3):235-246, 2011.

[PSCF05] M. Pennacchio, G. Savaré, and P. Colli Franzone. Multiscale modeling for the bioelectric activity of the heart. *SIAM Journal on Mathematical Analysis*, 37(4):1333-1370, 2005.
[Pta17] M. Ptashnyk. Multiscale modelling and analysis of signalling processes in tissues with non-periodic distribution of cells. *Vietnam Journal of Mathematics*, 45(1):295-316, 2017.

[QLRRB17] A. Quarteroni, T. Lassila, S. Rossi, and R. Ruiz-Baier. Integrated heart-coupling multiscale and multiphysics models for the simulation of the cardiac function. *Computer Methods in Applied Mechanics and Engineering*, 314:345-407, 2017. Special Issue on Biological Systems Dedicated to William S. Klug.

[RC11] G. Richardson and J. Chapman. Derivation of the bidomain equations for a beating heart with a general microstructure. *SIAM Journal on Applied Mathematics*, 71(3):657-675, 2011.

[RDN<sup>+</sup>21] V. Y. Reddy, S. R. Dukkipati, P. Neuzil, A. Anic, J. Petru, M. Funasako, H. Cochet, K. Minami, T. Breskovic, I. Sikiric, et al. Pulsed field ablation of paroxysmal atrial fibrillation: 1-year outcomes of IMPULSE, PEFCAT, and PEFCAT II. *Clinical Electrophysiology*, 7(5):614-627, 2021.

[RHM<sup>+</sup>14] B. Ribba, Nick H. Holford, P. Magni, I. Trocóniz, I. Gueorguieva, P. Girard, C. Sarr, M. Elishmereni, C. Kloft, and Lena 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.

[Ric09] G. Richardson. A multiscale approach to modelling electrochemical processes occurring across the cell membrane with application to transmission of action potentials. *Mathematical Medicine and Biology*, 26(3):201-224, 2009.

[RKP<sup>+</sup>12] B. Ribba, G. Kaloshi, M. Peyre, D. Ricard, V. Calvez, M. Tod, B. Čajavec-Bernard, A. Idbaih, D. Psimaras, L. Dainese, et al. A tumor growth inhibition model for low-grade glioma treated with chemotherapy or radiotherapy. *Clinical Cancer Research*, 18(18):5071-5080, 2012.

[Ros21] F. Rosier. COVID-19 : malgré cinq semaines de baisse continue en métropole, les épidémiologistes restent prudents. *Le Monde*, 2021.

[RRBP18] R. Romani, G. Ryan, C. Benner, and J. Pollock. Non-operative meningiomas: long-term follow-up of 136 patients. *Acta neurochirurgica*, 160(8):1547-1553, 2018.

[RSC<sup>+</sup>06] B. Ribba, O. Saut, T. Colin, D. Bresch, E. Grenier, and J.-P. Boissel. A multiscale mathematical model of avascular tumor growth to investigate the therapeutic benefit of anti-invasive agents. *Journal of theoretical biology*, 243(4):532-541, 2006.

[Sac04] F.B. Sachse. *Computational Cardiology: Modeling of Anatomy, Electrophysiology and Mechanics*. Springer-Verlag, 2004.

[SAM00] K. R. Swanson, E. C. Alvord, and J. D. Murray. A quantitative model for differential motility of gliomas in grey and white matter. *Cell Prolif.*, 33(5):317-329, 2000.

[SB93] D. R. Stauffer and J.-W. Bao. Optimal determination of nudging coefficients using the adjoint equations. *Tellus A*, 45(5):358-369, 1993.

[SBMAJ03] K. R. Swanson, C. Bridge, J.D. Murray, and E. C. Alvord Jr. Virtual and real brain tumors: using mathematical modeling to quantify glioma growth and invasion. *Journal of the neurological sciences*, 216(1):1-10, 2003.
[SCB+05] D. Sel, D. Cukjati, D. Batiuskaite, T. Slivnik, L. M. Mir, and D. Miklavcic. Sequential finite element model of tissue electropermeabilization. *IEEE Transactions on Biomedical Engineering*, 52(5):816-827, 2005.

[Sch66] H. P. Schwan. Alternating current electrode polarization. *Biophysik*, 3(2):181-201, 1966.

[SCN+17] O. Sutter, J. Calvo, G. N'Kontchou, J.-C. Nault, R. Ourabia, P. Nahon, N. Ganne-Carrié, V. Bourcier, N. Zentar, F. Bouhafs, et al. Safety and efficacy of irreversible electroporation for the treatment of hepatocellular carcinoma not amenable to thermal ablation techniques: a retrospective single-center case series. *Radiology*, 284(3):877-886, 2017.

[SdlRB+09] V. Sans, E. D. de la Roque, J. Berge, N. Grenier, F. Boralevi, J. Mazereeuw-Hautier, D. Lipsker, E. Dupuis, K. Ezzedine, P. Vergnes, et al. Propranolol for severe infantile hemangiomas: follow-up report. *Pediatrics*, 124(3):e423-e431, 2009.

[SF68] H. P. Schwan and C. D. Ferris. Four-electrode null techniques for impedance measurement with high resolution. *Review of scientific instruments*, 39(4):481-485, 1968.

[SLCFS14] O. Saut, J.-B. Lagaert, T. Colin, and H. M. Fathallah-Shaykh. A multilayer grow-or-go model for GBM: effects of invasive cells and anti-angiogenesis on growth. *Bulletin of mathematical biology*, 76(9):2306-2333, 2014.

[SNdJ+14] H. J. Scheffer, K. Nielsen, M. C. de Jong, A. van Tilborg, J. M. Vieveen, A. Bouwman, S. Meijer, C. van Kuijk, P. 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.

[SPGS18] L. Sala, C. Prud'homme, G. Guidoboni, and M. Szopos. Ocular mathematical virtual simulator: A hemodynamical and biomechanical study towards clinical applications. *Journal of Coupled Systems and Multiscale Dynamics*, 6(3), 2018.

[SRA+10] M. E. Sughrue, M. J. Rutkowski, D. Aranda, I. J. Barani, M. W. McDermott, and A. T. Parsa. Treatment decision making based on the published natural history and growth rate of small meningiomas: a review and meta-analysis. *Journal of neurosurgery*, 113(5):1036-1042, 2010.

[SS90] D. R. Stauffer and N. L. Seaman. Use of four-dimensional data assimilation in a limited-area mesoscale model. Part I: Experiments with synoptic-scale data. *Monthly Weather Review*, 118(6):1250-1277, 1990.

[SSGW14] K. C. Smith, R. S. Son, T. R. Gowrishankar, and J. C. Weaver. Emergence of a large pore subpopulation during electroporating pulses. *Bioelectrochemistry*, 100:3-10, dec 2014.

[TDP+04] M.-C. Trudel, B. Dubé, M. Potse, R.M. Gulrajani, and L.J. Leon. Simulation of QRST integral maps with a membrane-based computer heart model employing parallel processing. *IEEE Transactions on Biomedical Engineering*, 51(8):1319-1329, 2004.

[TGR05] J. Teissie, M. Golzio, and MP. Rols. Mechanisms of cell membrane electropermeabilization: a minireview of our present (lack of?) knowledge. *Biochimica et Biophysica Acta (BBA)-General Subjects*, 1724(3):270-280, 2005.

[TM05] R. Temam and A. Miranville. *Mathematical Modeling in Continuum Mechanics*. Cambridge University Press, second edition, 2005.
[tTNNP04] K. H. W. J. ten Tusscher, D. Noble, P.J. Noble, and A. V. Panfilov. A model for human ventricular tissue. *American Journal of Physiology-Heart and Circulatory Physiology*, 286(4):H1573-H1589, 2004.

[tTP06] K. H. W. J. ten Tusscher and A. V. Panfilov. Alternans and spiral breakup in a human ventricular tissue model. *American Journal of Physiology-Heart and Circulatory Physiology*, 291(3):H1088-H1100, 2006.

[Tun78] L. Tung. *A bi-domain model for describing ischemic myocardial d-c potentials*. PhD thesis, Massachusetts Institute of Technology. Dept. of Electrical Engineering and Computer Science., 1978.

[Ver97] G. Verbeke. Linear mixed models for longitudinal data. In *Linear mixed models in practice*, pages 63-153. Springer, 1997.

[VIT<sup>+</sup>07] M. W. Vernooij, M. A. Ikram, H. L. Tanghe, A. JPE Vincent, A. Hofman, G. P. Krestin, W. J. Niessen, M. MB Breteler, and A. van der Lugt. Incidental findings on brain MRI in the general population. *New England Journal of Medicine*, 357(18):1821-1828, 2007.

[VLDP03] P.A. Vidard, F.-X. Le Dimet, and A. Piacentini. Determination of optimal nudging coefficients. *Tellus A*, 55(1):1-15, 2003.

[VRF<sup>+</sup>20] C. Vaghi, A. Rodallec, R. Fanciullino, J. Ciccolini, JP. Mochel, M. Mastri, C. Poignard, JML. 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.

[VSM<sup>+</sup>18] D. Voyer, A. Silve, L. M. Mir, R. Scorretti, and C. Poignard. Dynamical modeling of tissue electroporation. *Bioelectrochemistry*, 119:98-110, 2018.

[War75] F. Wartak. *Electrocardiogram interpretation*. Medical Education Systems, 1975.

[WC96] J. C. Weaver and Y. A. Chizmadzhev. Theory of electroporation: A review. *Bioelectrochemistry and Bioenergetics*, 41(2):135-160, dec 1996.

[WGRS12] S. T. Wall, J. M. Guccione, M. B. Ratcliffe, and J. S. Sundnes. Electromechanical feedback with reduced cellular connectivity alters electrical activity in an infarct injured left ventricle: a finite element model study. *American Journal of Physiology-Heart and Circulatory Physiology*, 302(1):H206-H214, 2012.

[WN06] F. H. M. Wittkampf and H. Nakagawa. RF catheter ablation: lessons on lesions. *Pacing and Clinical Electrophysiology*, 29(11):1285-1297, 2006.

[WOH<sup>+</sup>95] D. Wei, O. Okazaki, K. Harumi, E. Harasawa, and H. Hosaka. Comparative simulation of excitation and body surface electrocardiogram with isotropic and anisotropic computer heart models. *Biomedical Engineering, IEEE Transactions on*, 42(4):343-357, 1995.

[WSE<sup>+</sup>12] J. C. Weaver, K. C. Smith, A. T. Esser, R. S. Son, and T. R. Gowrishankar. A brief overview of electroporation pulse strength-duration space: A region where additional intracellular effects are expected. *Bioelectrochemistry*, 87:236-243, oct 2012.

[ZAD<sup>+</sup>08] L. A. Zeidman, W. J. Ankenbrandt, H. Du, N. Paleologos, and N. A. Vick. Growth rate of non-operated meningiomas. *Journal of neurology*, 255(6):891-895, 2008.
[ZLJ<sup>+</sup>15] L. Zeng, P. Liang, J. Jiao, J. Chen, and T. Lei. Will an asymptomatic meningioma grow or not grow? A meta-analysis. *Journal of Neurological Surgery Part A: Central European Neurosurgery*, 76(05):341-347, 2015.

[ZNLD92] X. Zou, I.M. Navon, and F.-X. Le Dimet. An optimal nudging data assimilation scheme using parameter estimation. *Quarterly Journal of the Royal Meteorological Society*, 118(508):1163-1186, 1992.