teaching and research institutions in France or abroad, or from public or private research centers.

Marginal regeneration in a horizontal film. Instability growth law in the non linear regime.

Alice Gros, Adrien Bussonnière, Sanjiban Nath and Isabelle Cantat Univ Rennes, CNRS, IPR (Institut de Physique de Rennes) -UMR 6251, F-35000 Rennes (Dated: February 25, 2021) Marginal regeneration is the process at the origin of the thinner patches produced at the bottom of vertical soap films and slowly rising to the top. It controls the film drainage, and thus its life time, and ultimately the life span of the whole foam sample. Despite its importance, the dynamical properties of this original hydrodynamical instability remains to elucidate. Using a deformable frame, we identify a way to trigger the instability in a horizontal film; as the gravity-induced rising motion of the patches is suppressed in this case, we are able to study its long time dynamics for the first time, in a well controlled situation. The predictions made by Mysels et al. are successfully tested on our data and are extended by an original coarsening model we develop. This model quantitatively reproduces the patch growth in the horizontal geometry and, when adapted to the vertical situation, is consistent with the known properties of the instability.

I. INTRODUCTION

When a liquid film has just been pulled out from a foaming solution, its thickness is usually much greater than its equilibrium value, which causes a drainage. The film is bounded by menisci which connect the film either to a solid frame, to the liquid reservoir or, in a foam, to the neighboring films. In this last case, the meniscus is called a Plateau border. These menisci are at a pressure lower than the film and thus drain it. The drainage dynamics, induced by this capillary suction or by the gravity, is slowed down by the interface elasticity (Gibbs elasticity) [START_REF] De Gennes | Young soap films[END_REF][START_REF] Kim | Marangoni elasticity of flowing soap films[END_REF] or interface viscosity, which reduces or suppresses interface stretching. The drainage governs the film life time, and is thus crucial for most applications involving foam or bubbles [START_REF] Cantat | Foams. Structure and Dynamics[END_REF].

The relevant systems to address the drainage dynamics are the film elements, defined as elementary volumes of film spanning from one interface to the other, of volume h dS. If the Poiseuille flows between the interfaces are negligible, which is usually the case outside a small domain close to the menisci, these elements are closed systems, of conserved volume. For large Gibbs elasticity, which will be assumed in this paper, the area dS is conserved too, leading to the conservation of h. One key degree of freedom remains nevertheless present in the film: the deformation and recirculation motion of the film elements in the plane of the film. This sliding puzzle dynamics is only damped by the interface shear viscosity and may consequently be driven by tiny forces [START_REF] Chomaz | The dynamics of a viscous soap film with soluble surfactant[END_REF][START_REF] Winkler | Exponentially fast thinning of nanoscale films by turbulent mixing[END_REF].

In their seminal studies, Mysels, Shinoda and Frankel [START_REF] Mysels | Soap films: Study of their thinning and a bibliography[END_REF] proposed a subtle mechanism to explain the drainage, dominated by exchanges of film elements between the film and the meniscus, the so-called marginal regeneration. This theory predicts that along a meniscus, one finds alternatively some thick film elements plunging into the meniscus, and, at the same time, thinner film elements being extracted from the meniscus : the invariance of the flow along the meniscus direction is broken by a hydrodynamical instability. The film area conservation imposes a vanishing net flux of interface between the film and the meniscus, film elements leaving the film being replaced by film elements entering the film. However, as the elements leaving the film are thicker than the ones entering the film, the volume flux toward the meniscus is positive, which results in an efficient drainage [START_REF] Lhuissier | Bursting bubble aerosols[END_REF][START_REF] Frostad | Dynamic fluid-film interferometry as a predictor of bulk foam properties[END_REF][START_REF] Poulain | Ageing and burst of surface bubbles[END_REF].

The determination of the instability wave length and of the extracted film thickness is a long standing question, and the respective roles of gravity, surface tension gradient and initial film thickness in the control of these parameters are still controversial [START_REF] Lhuissier | Bursting bubble aerosols[END_REF][START_REF] Bruinsma | Theory of hydrodynamic convection in soap films[END_REF][START_REF] Vincent | Marangoni flow driven instabilities and marginal regeneration[END_REF]. In this paper, we provide quantitative observations of a horizontal foam film drainage and confirm that marginal regeneration stands even without gravity effect. We observe a surprisingly regular pattern along the meniscus, consisting in semi-circular patches of thin film invading the initially present thick film, separated by small channels through which the thick film flows toward the meniscus (see Fig. 1). We identify the condition allowing to close Mysels et al. set of equations and we develop an original model which succeeds in predicting the patch thickness, shape and growth rate.

II. EXPERIMENTAL RESULTS

The foam film is produced on a deformable frame similar to the one used in [START_REF] Bussonnière | Dynamical coupling between connected foam films: Interface transfer across the menisci[END_REF], supporting three flat rectangular films and a free meniscus (see Fig. 1(a)). The film of interest is horizontal and limited by two immobile metallic edges (the lateral edges), by a free meniscus and by a mobile metallic edge of length W = 65 mm which can translate along the lateral edges. We vary the distance L between the free meniscus and the mobile edge to prepare the film and the dynamics is then observed with the motor at rest. The free meniscus and the meniscus on the mobile edge play the same role and are not distinguished in the article.

Using a dedicated preparation protocol we produce a horizontal film, invariant along the meniscus (xdirection), having a vanishing velocity and a uniform thickness in a large margin close to a meniscus (the green/yellow domains in Fig. 1(b)). The dark blue profile, shown in Fig. 1(c), is an example of the film profile P 0 just before the instability begins. It is obtained as follows. First, a film is produced by plunging the frame into a foaming solution ; the central part of P 0 of thickness h center ≈ 0.5 µm (hereafter the central film) is produced at this step. Secondly, the film is deformed by pulling the mobile edge of the frame ; the thick parts on both sides of P 0 (hereafter the thick films), of average thickness close to 1.2 µm in the example of Fig. 1(c), are extracted from the menisci, the film thickness being controlled by the motor pulling velocity. Lastly, the film deformation is stopped by a pushing step. The whole film thickness increases during the pushing step and reaches its final profile P 0 at a time close to the instability onset. The average thick film thickness obtained at this stage is denoted h in . The details of the procedure are given in Appendix A. In particular, we detail there the role of the pushing step, we provide the used parameter values, and describe the averaging procedure used to define h in .

The foaming solution is made of deionized water, sodium dodecyl sulfate (SDS, cmc = 2.36 g/L) at 5.6 g/L, dodecanol at 0.05 g/L, glycerol at 15 % in volume and fluorescein sodium salt at 0.8 g/L. The viscosity is η = 1.52 mPa.s, and the equilibrium surface tension is γ = 32 mN/m. After the motor has stopped, and the film has stabilized at large scale, we measure the evolution of the film thickness close to the menisci from the emission intensity of the fluorescein dissolved in the foaming solution. The calibration of the relationship between the fluorescence intensity and the thickness has been performed by measuring the film thickness along a line using a spectral camera (see [START_REF] Bussonnière | Dynamical coupling between connected foam films: Interface transfer across the menisci[END_REF]).

At the end of the film preparation, a pinch is visible close to the menisci and it grows with time, in qualitative agreement with the asymptotic growth laws predicted in [START_REF] Aradian | marginal pinching" in soap films[END_REF] (the fast destabilization of the pinch and our space resolution do not allow us to make quantitative comparisons). This pinch is initially invariant in the x-direction, as shown in Fig. 1 (d,A). At some point, chosen as reference time t = 0, the pinch destabilizes (Fig. 1(d,B)) and turns into circular spots of very thin film. These initial spots, of typical diameter 300 µm, seem to define nucleation points for the instability observed afterwards (see the movie presented in Appendix D). After this pinch destabilization, very regular semi-circular patches of thin film, of typical thickness 1 µm, appear along the meniscus. Theses patches are extracted from the meniscus at a velocity v out . This velocity represents the film element velocity at the frontier between the patch and the meniscus. It will be shown to be homogeneous along the meniscus and close to the patch radius time derivative. The radius of the patches increases with time and they progressively invade the film, as shown by the evolution of average thickness profiles in Fig. 1(c). When two patches become too close they coalesce and a larger patch is produced, with the total patch area being conserved (see Fig. 1(d,D)). The patch radii R are measured by image processing of the top views. The mean patch radius R, averaged over all patches, is plotted as a function of time in Fig. 2(a), for two thicknesses h in of the thick film surrounding the patches (see Appendix B for details on the image processing and on the averaging procedure of the patch radii and extraction velocities). The wavelength of the marginal regeneration is 2 R + b ≈ 2 R, with b the average width of the channel between two patches, and thus evolves with time during this coarsening process. The dynamics is faster for thicker films, and sublinear with respect to time. As d R/dt ∼ v out (see eq. 10 in section III), this indicates that the extraction velocity slows down with time.

As shown in Fig. 1(c), the frontier between the central film and the thick film stays at the same position, and the central film thickness does not vary. This indicates an absence of compression/dilatation in this central region of the film during the patch growth, and thus a conservation of the local interface area. Importantly, the stretching or compression deformations are homogeneous in a foam film [START_REF] Durand | Relaxation time of the topological t1 process in a two-dimensional foam[END_REF]. The patch dynamics, occurring on the lateral sides of the same film, thus verifies the interface area conservation too: the assumption of incompressible interface is verified during the patch growth.

Consequently, after its extraction from the meniscus, the thickness of a film element remains the same during the whole dynamics, and the patch thickness distribution deduced from the color scale in Fig. 1 contains the memory of patch formation history: the darkest spots, which are the thinnest part of the patches, are produced by the initial pinch destabilization; around these small spots, we find concentric domains of equal thicknesses. These domains correspond to the film extracted at successive times, from the thickest domains close to the dark spot, to the thinnest ones far from it, and later reorganized by the sliding puzzle-like recirculation motion in the patches. The darkest spots provide passive tracers and clearly show this recirculation in Fig. 2(c): the central part of the patch moves away from the meniscus at a typical velocity of the order of 0.2 mm/s, while its side parts move toward the meniscus at a higher velocity, of the order of 2 mm/s. Note that patches can contain several dark spots as patches coalesce over time. Despite this spatial reorganization, the thickness of the film elements extracted from the meniscus at the time t, denoted h out (t), can be measured from the comparison of the thickness distribution in the patches between two successive times (see the details of the image processing in Appendix B). As shown in Fig. 2(b), the extracted thickness averaged over all patches hout (t) decreases with time and is thinner for a thinner surrounding film h in .

When some film is extracted from the meniscus at the velocity v out , additional interface area is brought to the film and the interface area conservation in the film thus imposes that the same amount of area flows back into the meniscus. This is fulfilled by the thick film flowing back toward the meniscus through the small channels between the patches. The channel velocity at the meniscus is denoted v in . This velocity is more difficult to measure directly but is close to the velocity measured on the sides of the patches (see Fig. 2(c)), by velocity continuity. The channel width b is much smaller than the patch diameter, imposing that v in is significantly larger than v out .

The theoretical predictions plotted in Fig. 2 are established from the model developed in the next sections, which extends and closes the model proposed by Mysels et al. [START_REF] Mysels | Soap films: Study of their thinning and a bibliography[END_REF].

III. MODEL A. Lubrication equations

The flow in the domain located between a meniscus of radius r m and a thin film of uniform thickness h is a A c c e p t e d m a n u s c r i p t classical problem [START_REF] Mysels | Soap films: Study of their thinning and a bibliography[END_REF]. Let us consider (i) a steady situation, (ii) with incompressible interfaces and (iii) invariant in the x-direction. The static meniscus of uniform curvature is reached at y → -∞, the asymptotic flat film is reached at y → ∞ and the whole dynamics occurs in the so-called dynamical meniscus, spanning over a length r m from the static meniscus, in which the interface curvature varies. In this simplest case, the film profile ĥ(y) is solution of the equation

v 0 8 ĥ3 ∂ 3ĥ ∂y 3 + v(h -ĥ) = 0 , (1)
where v 0 = γ 3η is the visco-capillary velocity and v is the uniform interface velocity defined such as v < 0 when the interface flows towards the meniscus ("in"), and v > 0 otherwise ("out") (see Fig. 3).

The velocity v out ∼ d R/dt varies slowly with time (see Fig. 2) and the capillary number Ca = v out /v 0 is at most of the order of 10 -4 in our experiments. The assumption of incompressible interface is valid if the Gibbs modulus E verifies E/γ Ca 2/3 ∼ 2 .10 -3 , which is the case for our foaming solution [START_REF] Bussonnière | Dynamical coupling between connected foam films: Interface transfer across the menisci[END_REF][START_REF] Seiwert | Theoretical study of the generation of soap films: role of interfacial viscoelasticity[END_REF][START_REF] Champougny | Surfactant-induced rigidity of interfaces: a unified approach to free and dip-coated films[END_REF]. Our experiments therefore verify the theoretical assumptions (i) and (ii). The observed flow is however obviously not invariant along the meniscus and the problem should in principle be solved numerically in the (x, y) plane. In order to provide simple analytical predictions, we assume in the following that the velocity profile in the dynamical meniscus is binary, i.e. that the flow is invariant along x inside each patch on one hand (with v = v out > 0), and in the channels between the patches on the other hand, (with v = -v in < 0) (see Fig. 3). This simplification allows us to use the eq. 1 to predict the dynamics in each domain.

The solution for the pulling case (v > 0) verifies [START_REF] Mysels | Soap films: Study of their thinning and a bibliography[END_REF]

h out = 1.28 r m v out v 0 2/3 , (2)
δγ out = 1.85 γ v out v 0 2/3 , (3)
with δγ out ∼ 1 mN/m, hereafter called film tension, the surface tension difference between the film and the meniscus assumed to be at a uniform reference tension γ.

To the best of our knowledge, only the pulling case has been studied experimentally so far, leading to a quantitative agreement with the theoretical prediction of eq. 2 [START_REF] Saulnier | What is the mechanism of soap film entrainment?[END_REF][START_REF] Seiwert | Extension of a suspended soap film: A homogeneous dilatation followed by new film extraction[END_REF].

When v = 0 a steady solution is never reached and the evolution of the film profile with time in this case is given in [START_REF] Aradian | marginal pinching" in soap films[END_REF].

In the pushing case (v < 0), the solution depends on the imposed velocity v in and on a given asymptotic film thickness h in , which is, in contrast to the pulling case, a control parameter. The film tension is given by

δγ in = 2.5 γ v in v 0 2/3 H(α) (4)
where

α = h in 2r m v in v 0 -2/3 . ( 5
)
The function H does not have an analytical expression, but its numerical value is well fitted by H(α) = λα 1/3 -1.05, as proposed in [START_REF] Mysels | Soap films: Study of their thinning and a bibliography[END_REF] with λ = 1 (Fig 5 .5 of the reference). The film tension is negative for small values of α, which is achieved for thin films, large menisci or high velocities. In such cases, the tension in the film is lower than the tension in the meniscus, and is therefore a driving force for the motion toward the meniscus. An important transition occurs at α = 1.16/λ 3 : at greater values, the tension in the film is higher than the meniscus tension, i.e. δγ in > 0 as in the pulling case, and the driving force is the meniscus capillary suction only. The surface tension gradients, at the origin of the surface tension difference δγ, are localized in the dynamical meniscus spanning over ∼ r m (Ca out ) 1/3 ∼ 100µm. For greater values of y (further from the meniscus), the Laplace pressure becomes negligible. As the dynamical meniscus width is much smaller than the patch radii ( R), the pressure in the patches is thus the reference pressure and the tension obeys η s ∆ 2D v s + ∇ 2D γ = 0, with η s ∼ 10 -7 kg/s [START_REF] Drenckhan | Fluid dynamics of rivulet flow between plates[END_REF] the interface shear viscosity, and v s the interface velocity in the (x, y) plane. The observed recirculations in the film occur over a characteristic distance equal to the patch radius R, and with a velocity v in ∼ 10 -3 m/s (see Fig. 2(c)). The associated tension variations thus scale as η s v in /R ∼ 10 -7 N/m.

The surface tension in the film is also modified by the line tension T ∼ 10 -10 N induced by the film thickness variation across the frontier between the thin and thick A c c e p t e d m a n u s c r i p t parts of the film, in a transition region of width δ ∼ 100 µm [START_REF] Shabalina | Rayleigh-taylor-like instability in a foam film[END_REF]. This line tension imposes a surface tension jump across the frontier, which is the 2D equivalent of the 3D Laplace pressure jump: in a patch of radius R, the tension is slightly smaller than in the thick film on the other side of the frontier, and the tension jump is T /R ∼ 10 -7 N/m.

The semi-circular patch shapes we experimentally observed indicate that the line tension T dominates in the film, over the viscous interface stress associated to η s : the patches are very mobile and adjust immediately their shapes to insure the minimization of their perimeters. Anyway, as both surface tension variations associated to T and to η s are much smaller than δγ, the patch growth can be predicted by assuming that the tension is uniform in the entire film and that it only varies across the dynamical meniscus from γ in the static meniscus to γ + δγ in the film.

We deduce from eq. 3 that the extraction velocity v out , which depends only on δγ out , is homogeneous along the frontier with the meniscus. The relationship between δγ in and v in involves in contrast the additional parameter h in . In the following, we assume that the film thickness in the channels is homogeneous and equal to the thickness of the incoming film, which thus imposes a uniform value of v in . As discussed below, this assumption may need further theoretical and experimental improvements. Finally, using δγ in = δγ out in eq. 3 and 4, we obtain 1.85

v out v 0 2/3 = 2.5 v in v 0 2/3 λ α 1/3 -1.05 , (6)
with α given by eq. 5.

Using our experimental values r m = 320 µm, h in ∼ 1 µm, v in ∼ 2 mm/s, we get α = 0.36, leading to a negative value of H if λ = 1, as established in [START_REF] Mysels | Soap films: Study of their thinning and a bibliography[END_REF]. This corresponds to a negative value of δγ in , which can therefore not be equal to δγ out ∼ 10 -3 N/m obtained from eq. 3. The same inconsistency between the theory and the experimental observations has already been noticed in § 5.12 of [START_REF] Mysels | Soap films: Study of their thinning and a bibliography[END_REF]. We believe that the channel width is not large enough for the assumption of invariance along x (assumption (iii)) to be verified. The x-invariant model predicts thickness oscillations along the steady profile [START_REF] Mysels | Soap films: Study of their thinning and a bibliography[END_REF]. As we observe that the x-invariant pinch shown in Fig. 1 (d,A) is unstable, we believe that the pinch of this theoretical profile may be unstable too. The actual profile along the y direction in the middle of the thick film channels may thus differ from the x-invariant profile, the thinnest parts of this invariant profile being repelled on the lateral sides of the channel. Note that without pinch in the thickness profile, we expect the viscous dissipation to be lowered, and the film tension to be higher. The non-uniformity of the film thickness across the channels also implies, from eq. 4, the non-uniformity of v in . In this simple model, v in should thus be considered as an average value across the channels (identical in each channel).

These remarks led us to the introduction of the adhoc free parameter λ in the fitting expression of H (nonexistent in [START_REF] Mysels | Soap films: Study of their thinning and a bibliography[END_REF] and thus taken as 1 when introduced in eq. 4), arbitrarily keeping the functional dependency proposed in [START_REF] Mysels | Soap films: Study of their thinning and a bibliography[END_REF] for H. A good agreement between theory and experiments is obtained for a value of λ around 1.8, as presented in the following results.

B. Interface area conservation

The interface area conservation imposes that the net interface area flux between the film and the meniscus vanishes (see Fig. 3), leading to the relation

v out ξ -v in (1 -ξ) = 0 , (7)
with ξ the length fraction of meniscus in contact with the patches. Using eq. 6, we get

v in = v F(ξ) ; v out = v 1 -ξ ξ F(ξ) (8)
with v = v 0 λ 3 h in 2 r m 3/2
and

F(ξ) = 1.85 2.5 1 -ξ ξ 2/3 + 1.05 -9/2 .
The film tension δγ is eventually expressed as a function of ξ solely, using eqs. 3 and 8. As shown in Fig. 4, the family of solutions (v in , v out ), obtained for given ξ values in the range [0 -1], has a film tension in the range [0 δγ max ], with δγ max = 0.17 γ λ 3 h in /r m . For tension differences out of that range, the coexistence between the patches and the thick film is not possible, and an x-invariant solution should be observed. The patch thickness, given by eq. 2, can be expressed as

h out = 0.64λ 3 h in 1 -ξ ξ F(ξ) 2/3 < 0.11 λ 3 h in . (9)
For λ = 1.8, the upper bound of thickness ratio is h out /h in < 0.7 < 1. For any value of ξ, the marginal regeneration thus leads to a decrease of the film volume, as expected from the drainage behavior.

C. Growth model

The whole dynamics is eventually governed by the sole dynamical parameter ξ, which remains to be determined. Mysels et al. proposed to close the problem by assuming that the selected ξ value should minimize the dissipation, leading to the prediction α mysels = 3.5 in eq. 6 (with λ = 1). Here, we propose instead an argument of kinematic nature to predict ξ as a function of time, from which we deduce all other observable quantities. and v in /v (dashed line) given by eq. 8, as a function of ξ. The red curve is δγ/γ, given by eq. 3, with γ = γλ 3 h in /(2rm) and λ = 1.8. The shaded domain highlights the experimentally relevant values of ξ.

We first model the dynamic based on the averaged quantities and discuss the patch distribution in the next section. Note that patch size distribution is peaked, with a maximum close to its average, which justifies this approach. In this section we therefore assume a large number N of patches of homogeneous radius R along a meniscus of length W , separated by channels of width b (see Fig. 3). The parameter

ξ equals 2N R/W = 1 -N b/W . Assuming b 2 R as observed experi- mentally, W ≈ 2N R and ξ ≈ 1 -b/(2 R). The patch area is Ā = π R2 /2 and varies between coalescences as d Ā/dt = π R(d R/dt) grow = 2 Rv out .
The patch size increases until the channel width reaches a critical value b min , leading to the coalescence of two neighboring patches. The channel width observed on the image just before the coalescence is b min ≈ 0.2 mm which corresponds to a few pixels. We expect the coalescence to be triggered by an attractive force, of capillary nature [START_REF] Ho | Direct measurement of capillary attraction between floating disks[END_REF], arising when two interface deformations associated to the thick/thin film frontier overlap. The characteristic length b min should thus scale as the frontier thickness 2δ ∼ 0.2 mm so the orders of magnitude are consistent.

Let us consider the coalescence of a single pair of patches. It obviously breaks the patch size uniformity and lead to new average quantities, which are however considered as the new uniform size of all patches just after the coalescence. The coalescence leads to a new patch, of area A n = 2 Ā and radius R n = √ 2 R. The average radius jumps from R to R +δR c , with δR c = ( √ 2-1) R/N . The average distance between patches increases from b min just before the coalescence to b min + δb c just after, with δb c = 2(2 -√ 2) R/N . Assuming that the free space between patches is equally redistributed between all patches, the time δt c before the next coalescence is the time needed for b to reach b min again. So

δt c = δb c /(4v out /π) = π/(2v out )(2 - √ 2
) R/N . The average patch size increase is the sum of the growing and coalescence contributions, so d R/dt = d R/dt grow + δR c /δt c , leading eventually to the rule

d R dt = 2 π(2 - √ 2) v out ∼ 1.1v out . ( 10
)
The value of ξ, which determines v out , varies between 1-b min /(2 R) and 1-(b min +δb c )/(2 R). Experimentally, R varies in the range [0.5-3] mm. With W = 65 mm, the total number of patches is in the range [100 -20]. Consequently δb c is much smaller than b min for small patches and becomes of the same order for the largest ones. We thus get ξ ≈ 0.8 at the beginning of the experiment and ξ is in the range [0.94 -0.97] at the end (see the shaded domain in Fig. 4).

Simple laws are obtained at large R where b → b min with b min R. In this limit, we have ξ → 1 -b min /(2 R) and F(ξ) → 0.8 and by writing the eqs. 8 and 9 we get :

d R dt = 2 π(2 - √ 2) v out = 0.29 γ η b min 2 R λ 3 h in 2 r m 3/2 , (11)
R(t) = 0.29 γb min η λ 3 h in 2 r m 3/2 t + R(0) 2 1/2 , (12)
hout (t) = 0.55 h in b min 2 R(t) 2/3 λ 3 . (13)
The computed values of these scaling laws ( R(t) and hout (t)) are compared to the experimental mean radius and thickness of the patches in Fig. 2.

The theoretical prediction shown in Fig. 2(a) is obtained by fitting the parameters λ and R(0) (for both series independently) in the equation 12, assuming b min = 200 µm. The uncertainty on λ is deduced from the error bar on R(t). For hout (t), shown in Fig. 2(b), we determine λ from the time average hout (t) R(t) 2/3 /h in t . The error bar on λ can be obtained from the standard deviation of hout (t) R(t) 2/3 /h in at the different times, or from the error bar on hout (t), which are of the same order. The prediction shown in Fig. 2(b) is obtained using eq. 13, with the fitted λ value and the experimental value of R(t). The agreement between experiment and theory is good and the obtained λ values are all in the range 1.6 -1.8. This value depends on the value of b min which is poorly determined, and on the arbitrary choice of the function H. However, this adjustable parameter is close to one and does not significantly depend on h in , nor on the quantity considered ( R(t) and hout ). We thus believe that this simple model captures the main physical processes governing the dynamics.

Note that the thickness slightly decreases with y in the thick film, as visible in Fig. 1(c), potentially leading to the decrease of the film thickness in the channel h in as the patches invade the film. We carefully checked that this effect is negligible. The sublinear behavior of R(t) is thus due to the predicted dynamics and not to this thickness heterogeneity. Normalized distribution of the patch area for different times given in seconds in the legend, averaged over all the series having h in = 1.2 µm. The light gray line is the distribution obtained using the N0 patches measured at t = 0.2 s and making random coalescences between them until the patch number of the time t = 4.2 s is reached (we averaged the process over 100 realizations). The dark gray line is the exponential law corresponding to a random coalescence process for an infinite number of patches.

To investigate further than the average values, we now focus on the radius distribution, plotted in Fig. 5. At initial time the size resolution is very poor, which artificially produces a very peaked distribution. At later times however the distribution is well resolved and remains surprisingly stable in time. This significantly differs from the behavior expected from a random coalescence process, leading to an exponential size distribution, with a maximum at the initial patch size [START_REF] Aldous | Deterministic and stochastic models for coalescence (aggregation and coagulation): a review of the mean-field theory for probabilists[END_REF], as shown in Fig. 5 (see Appendix C for further details). This proves that the probability of coalescence depends on the patch size and that the smallest patches are more likely to coalesce. We identify below the physical process at the origin of this patch size stabilization.

The flow lines in the vicinity of two patches are represented in Fig. 3. The patches, labeled i and i + 1, are separated by a channel of width b i . The patch centers are respectively at the position x i and x i+1 and their radii are R i and R i+1 . The frontiers x = x i and x = x i+1 are local symmetry axes of the flow, so the x-component of the velocity vanishes along these lines. The interface flux balance (eq. 7) can therefore be written locally in the subsystem limited by these frontiers as

v out (R i + R i+1 ) -v in b i = 0 . (14)
The velocities v out and v in are determined by the uniform film tension, and thus by the global value of ξ: they do not depend on the pair of patches we consider. This means that the local quantity

ξ i = 1 -b i /(R i + R i+1 )
does not depend on i and is everywhere equal to ξ.

This analysis leads to the prediction that the channel width is proportional to the average radius of the neighboring patches: the smallest patches are limited by the narrowest channels and thus coalesce with a higher probability. This explains the narrow patch size distribution observed and justify our approach based on averaged quantities.

Finally, let us mention the stabilizing effect of ξ on the channel widths. In the range of ξ observed experimentally, the film tension decreases with ξ, as shown in Fig. 4. Consequently, if a fluctuation decreases the distance between two patches (thus imposing locally ξ i > ξ), then δγ i < δγ, and the lateral forces on the patches are unbalanced: the two patches are pulled apart, with counteracts the initial fluctuation. The additional free space δb c induced by a coalescence should be rapidly redistributed equally between all patches by this mechanism, thus supporting the assumption made in the model.

E. Extension of the model to the vertical case

In a vertical film, gravity induces an effective buoyancy force [START_REF] Adami | Surface tension profiles in vertical soap films[END_REF]: the thinnest film elements rise in the film, which is the most visible and popular signature of the marginal regeneration occurring at the bottom of a film. At small scale, however, gravity is negligible [START_REF] Lhuissier | Bursting bubble aerosols[END_REF] and our model should be relevant: we suggest that the patches may grow according to the dynamics described in this paper, until they reach a critical size R max at which they are torn away the meniscus by the gravity. At this critical size, the line tension should balance the buoyancy force, which is compatible with the recent observations made in [START_REF] Miguet | The role of marginal regeneration on bubble thinning[END_REF]. The capillary condition (b min ) is thus replaced by a gravitational condition (R max ) in vertical films.

The removal of a patch from the meniscus produces an empty space in which a new patch appears. Note that the apparition of such a new patch in a free space is neglected in our model because it is rare enough in horizontal films, but is clearly visible in some examples (see Fig. 1(C) for example). Importantly, the qualitative observation of such events suggests that the thickness of the new patch is the current extracted thickness h out (t) imposed by the dynamics in the rest of the film (see eq. 2), and not the thickness observed at the onset of the instability, governed by a different law. The steady dynamics of marginal regeneration under gravity would in that case obeys the laws established in this paper, the average patch radius being stabilized at R max .

An important remark is that the pattern shown in [START_REF] Lhuissier | Bursting bubble aerosols[END_REF] and [START_REF] Miguet | The role of marginal regeneration on bubble thinning[END_REF] have a distance between patches and a patch size of the same order of magnitude, so a value of ξ ∼ 0.5. For this value, v out does not vary much with ξ and stays close to its maximal value (see Fig. 4). The patch thickness should thus be close to h out = 0.7h in (see eq. 9), which is very close to the observations h out = 0.8h in reported in [START_REF] Vincent | Marginal regeneration in thin vertical liquid films[END_REF] and used in subsequent models [START_REF] Lhuissier | Bursting bubble aerosols[END_REF][START_REF] Miguet | The role of marginal regeneration on bubble thinning[END_REF]. Using an original set-up and a well controlled horizontal foam film, we provide the first direct and quantitative observations of the alternate velocities associated to marginal regeneration and we exhibit a coarsening process inducing a time dependency of the instability wave length. Film thicknesses and velocities are consistent with Mysels' model if a phenomenological correction in the viscous resistance of the film is introduced. Further 3D simulations would be required to predict the value of the added pre-factor λ, which is determined here by fitting the experimental data.

We believe that our theoretical work shed a new light on the marginal regeneration process. It provides predictions for the patch size and thickness in fairly good agreement with our experimental data. Moreover, first considerations for an adaptation to a vertical film provide a good agreement with data of the literature for the thickness ratio between the foam film and the marginal regeneration patches. This ratio is the key ingredient controlling the film drainage rate.

Finally, the precise physical process which initially triggers the destabilization of the invariant pinch remains to be elucidated, and the domains of stability and/or metastability of the x-invariant solutions of the eq. 1 (or its unsteady variant) still need to be identified.

APPENDIX A -INITIAL CONDITION IN THE HORIZONTAL FILM

Film preparation

A film is first produced with L 0 = 2.5 mm by plunging the frame into the foaming solution (the central film of thickness h center ). After a short drainage time, L is increased up to L max at a velocity V pull . During this step two films, called the thick films in the article, are extracted from the two menisci, with a thickness h in controlled by the extraction velocity [START_REF] Mysels | Soap films: Study of their thinning and a bibliography[END_REF][START_REF] Seiwert | Extension of a suspended soap film: A homogeneous dilatation followed by new film extraction[END_REF], which depends on V pull and of the central film elasticity. If the motor is stopped just after this pulling step, the extraction keeps going during few seconds, at a slower velocity. Indeed the central film, which was stretched during the pulling step, relaxes to a smaller area. To avoid this relaxation, and thus to get a better film thickness uniformity close to the meniscus, we reverse the motor motion when reaching L max and decrease L at the velocity V push down to L = L max -∆L. This pushing step stops the film extraction and systematically triggers the instability observed. It induces a transient film compression which modifies h center and h in .

We used L max = 18.5 mm for all series. The final thickness h in is 0.66 ± 0.08 µm for V pull = 2 mm/s and 1.2±0.1 µm for V pull = 3 mm/s. We used V push = 2 or 3 mm/s and ∆L in the range [START_REF] Mysels | Soap films: Study of their thinning and a bibliography[END_REF][START_REF] Lhuissier | Bursting bubble aerosols[END_REF][START_REF] Frostad | Dynamic fluid-film interferometry as a predictor of bulk foam properties[END_REF][START_REF] Poulain | Ageing and burst of surface bubbles[END_REF][START_REF] Bruinsma | Theory of hydrodynamic convection in soap films[END_REF] mm. This modifies the time at which the film velocity vanishes but not the final thickness. No correlation has been observed between the parameters V push and ∆L and the measured dynamics, in the investigated range.

Condition of vanishing velocity

The condition of vanishing velocity and vanishing extension/compression in the thick film is difficult to obtain and requires to tune the frame deformation used to prepare the film by trial and error. Both conditions are coupled. Indeed, if we denote L 0 the width (in the y direction) that the horizontal film would have at the equilibrium tension and ε its extension with respect to L 0 , we have L = εL 0 and

dL dt = L 0 dε dt + 2 v men , (15)
with v men the average film velocity along the menisci (with the convention v > 0 if oriented toward the meniscus). Once the motor has stopped, dL/dt = 0 by definition, and v men can be determined from a measure of dε/dt. The mass conservation in the central film imposes

dε dt = - h center 0 (h center ) 2 dh center dt , (16)
with h center 0 the central film thickness at long time (taken as arbitrary reference). The measure of the central film thickness thus provides information on the evolution of the extension (or compression) state of the film. Using eq. 15 it provides the average velocity close to the meniscus. In this paper, we only kept data series for which the variation of ε is less than 5% during the 8 s of patch evolution. This corresponds to an average velocity in the thick film at the meniscus below 0.1 mm/s during the measure.

Thickness profile in the thick film

The thick film characteristics are measured on the last image before the destabilization of the pinch (defined as t = 0). For each series, the profile is first averaged in the x direction to obtain h av1 (y) (see the first profile in Fig. 1(c)). On each side of the film, we define

h av2 = h av1 (y) [y1y2] , (17)
with y 1 and y 2 respectively at the distance 1 mm and 3 mm from the meniscus. The obtained value h av2 is finally averaged over all series sharing the same pulling velocity, and on both sides of each film, leading to

h in = h av2 series . (18)
The error bar given in the text is the standard deviation of h av2 on the different series.

A c c e p t e d m a n u s c r i p t

The mean radii and extracted film thicknesses are determined from the experimental data, to be compared to the output of our model. Axes are in pixels.

Patch radius

We analyze all the patches in the domain of width ∆x = 30 mm in the center of the film, in which the boundary effects are negligible and the lightening homogeneous. The position of the meniscus is automatically detected as well as the frontier between the thin and thick parts of the film. The channel positions are determined as the points of the frontier at minimal distance from the meniscus (see Fig. 6). The frontier shape between two successive channels are close to half ellipses, which radii R x and R y (in the directions x and y) are measured. In average, large patches are slightly elongated in the y direction, with R y being up to 20% larger than R x . This effect is not discussed in the paper and the patch shape is considered as being semi circular. We use arbitrarily R y in the analysis.

In order to filter artifactual small patches produced by the automatic contour detection, we suppressed from the averaging process the patches of radius smaller than 20% of the largest ones. We checked that the results do not depend on the threshold chosen. Note that, after the destabilization of the initial pinch, the appearance of new patches are rare events. However when a fluctuation produces a large channel during an especially long time, a new patch can be produced in the middle of the channel. These patches, suppressed by the filtering process, coalesce rapidly with their neighbor and do not participate to the global dynamics. An example of patch radius distribution is shown in Fig. 7.

Extracted film thickness

The image processing used to determine the thickness h out (t) of the film extracted from the meniscus between the times t and t -∆t is the following. We select all the pixels between the meniscus and the frontier separating the thin and the thick parts of the film (i.e. the area inside the patches) and make a histogram of film thicknesses in this domain at time t: this quantifies the film area A(h, t, dh) having a thickness in the range [h; h+dh] at time t.

The film extracted from the meniscus between the time t -∆t and t has the thickness distribution A extr (t, ∆t, h, dh) obtained as A extr (t, ∆t, h, dh) = A(h, t, dh) -A(h, t -∆t, dh) . [START_REF] Drenckhan | Fluid dynamics of rivulet flow between plates[END_REF] As the patches grow and as each film element in the patch keeps its area and thickness, this quantity is expected to be positive. The negative parts of the histogram shown in Fig. 8 are due to the experimental noise. The average thickness of the film extracted between the time t -∆t and t is given by hout (t) = h A extr (t, ∆t, h, dh)dh A extr (t, ∆t, h, dh)dh .

APPENDIX C-DISTRIBUTION INDUCED BY A RANDOM COALESCENCE PROCESS

Random coalescences lead generically to exponential distributions (see the review [START_REF] Aldous | Deterministic and stochastic models for coalescence (aggregation and coagulation): a review of the mean-field theory for probabilists[END_REF]) when the initial distribution is peaked and the initial and final numbers of object N 1 and N 2 verifies N 1 N 1 -N 2 1. As we have experimentally relatively small initial and final numbers of patches we performed numerical simulations with the experimental values N 1 = 300 (at t=0.2 s) and N 2 = 140 (at t=4.2 s) used in the Fig. 5 of the article (i. e. all the patches measured in the region of interest of all series having h in = 1.2µm).

A c c e p t e d m a n u s c r i p t

We generate a list of N 1 patches with a Gaussian area distribution having the same average and standard deviation than the experimental patch area distribution at t = 0.2 s. Then we sum the areas of randomly chosen neighboring patches, until only N 2 patches remain. The distribution shown in Fig. 5 of the main article is the final distribution obtained after averaging over 100 identical simulations. It is already close to the asymptotic exponential law, and clearly differs from the experimental distribution.

APPENDIX D -CAPTION OF THE MOVIE

The movie [26] shows the top view of the film, from t = 0 to t = 8.4 s. The ROI is cropped at the limit of the meniscus (at the top of the image) and its dimensions are 4 × 42 mm 2 . The thick film (at the bottom of the image) has a thickness h in = 1.2 µm.

FIG. 1 .

 1 FIG. 1. (a) Scheme of the set-up. Three foam films (in blue) are produced on a metallic frame and are connected to each other along a free meniscus. The red parts of the frame are immobile whereas the gray edge can move, which modifies the length L of the horizontal film. (b) Top view of the film, before the instability onset. Colors indicate the film thickness (same color scale as in (d)). (c) Evolution of film thickness profile as a function of y, averaged in the x-direction at every 0.8 s, from t = 0 (instability onset, dark blue profile P0), to t = 8.8 s (yellow). Same experiment as in (d). (d) Top view of the film, close to one meniscus, for the times (0, 0.2, 1, 2.6, 7.4) s, from (A) to (E). A meniscus appears in yellow (saturation of the thickness scale) at the bottom of the images and a homogeneous thick film at the top. A pinch (in blue) is visible at the meniscus boundary in (A) and it destabilizes in (B). Then patches of thin film grow along the meniscus (C-E). The local structure of the flow close to the meniscus is depicted in (E).

FIG. 2 .

 2 FIG. 2. (a)Average patch radius as a function of time, for h in = 0.66 ± 0.08 µm (blue, averaged over 4 series) and h in = 1.2 ± 0.1 µm (red, averaged over 12 series). The shaded areas are the standard deviations. The darkest lines are the power law eq. 12 with the adjustable parameter λ = 1.8 ± 0.05 (blue) and λ = 1.6 ± 0.05 (red). All series have been synchronized at the time t0.5 at which R = 0.5 mm (experimentally less noisy than t = 0). (b) Thickness of the film extracted from the meniscus in the period of time [t -∆t, t], as a function of time, with ∆t = 1.4 s. The circles are the experimental results, with the same convention as for (a). The darkest solid lines are the power law eq. 13 with λ = 1.8 ± 0.05 for both series. (c) We selected a patch without coalescence between t = 3 s and t = 8 s, and rescaled their images at successive times to keep the patch radius and position constant. This allowed to follow trajectories, with respect to the patch boundary position: the colored circles are the successive rescaled positions of a passive tracer, shown on the last image. We also computed (on the original, non rescaled, images) the tracer velocity, which amplitude is indicated by the color code. On this example, it varies from 0.2 mm/s (blue) to 2 mm/s (red).

FIG. 3 .

 3 FIG. 3. Scheme of the flow and notations used in the text. Each black semi-circle represents a patch indexed by i, of radius Ri, and bounded on one side by a dynamical meniscus of width , and on the other side by its frontier with the thick film, of width δ. The gray lines shown on the left are flow lines: the dashed lines are open lines beginning at the meniscus in the patch and ending at the meniscus in the channel; the solid lines are closed lines, associated to the recirculations visible in the patches. The uniform velocities at the meniscus in the patches (resp. in the channel) are symbolised by the red (resp. blue) arrows. The channel between the patches i and i + 1 is indexed by i and has a width bi.

3 FIG. 4 .

 34 FIG.[START_REF] Chomaz | The dynamics of a viscous soap film with soluble surfactant[END_REF]. The black curves are the velocities v out /v (solid line) and v in /v (dashed line) given by eq. 8, as a function of ξ. The red curve is δγ/γ, given by eq. 3, with γ = γλ 3 h in /(2rm) and λ = 1.8. The shaded domain highlights the experimentally relevant values of ξ.

  FIG.5. Normalized distribution of the patch area for different times given in seconds in the legend, averaged over all the series having h in = 1.2 µm. The light gray line is the distribution obtained using the N0 patches measured at t = 0.2 s and making random coalescences between them until the patch number of the time t = 4.2 s is reached (we averaged the process over 100 realizations). The dark gray line is the exponential law corresponding to a random coalescence process for an infinite number of patches.

FIG. 6 .

 6 FIG.6. On this raw image, the color indicates the thickness and the frontier automatically detected is shown in yellow. Axes are in pixels.

25 FIG. 7 . 8 FIG. 8 .

 25788 FIG. 7. (left) Size of all patches observed at a time t after the instability onset, for one series. (right) Patch number as a function of time for the same series.

ACKNOWLEDGMENTS

This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 725094). We thank M. Le Fur for technical support and A. Saint-Jalmes and E. Rio for fruitful discussions.


Simona Dobre

Thierry Bastogne

Muriel Barberi-Heyob

Denise Bechet

Jacques Didelon

Alain Richard
email: alain.richard@cran.uhp-nancy.fr

System identification of the intracellular photoreaction process induced by photodynamic therapy

Photodynamic therapy (PDT) is an alternative treatment for cancer that involves the administration of a photosensitizing agent, which is activated by light at a specific wavelength. This illumination causes a sequence of photoreactions, which -in the presence of molecular oxygen -is supposed to be responsible for the death of the tumor cells. The PDT efficiency stems from the optimal interaction between these three factors (light, drug and oxygen). In this paper, a new approach is proposed to estimate photophysical parameters which characterize the ability of a photosensitizing drug to produce singlet oxygen. This approach is based on system identification techniques. This model-based method would allow biologists to estimate all the photophysical parameters from spectro-fluorescence data generated by only one experiment. Secondly, contrary to usual techniques which are restricted to in vitro studies, this approach can be directly applied to in vivo data.

I. INTRODUCTION

Photodynamic therapy (PDT) [START_REF] Moser | Photodynamic Tumor Therapy: 2nd and 3rd Generation[END_REF], [START_REF] Goldman | Photodynamic Therapy[END_REF], [START_REF] Dougherty | Photodynamic Therapy[END_REF] is a treatment of displastic tissues such as cancers. Mainly, it involves the selective uptake and retention of the photosensitizing drug (photosensitizer, PS) in the tumor, followed by its illumination with light of appropriate wavelength. The PS activation is thought to produce, after multiple intermediate reactions, singlet oxygen at high doses (in the presence of molecular oxygen) and thereby to initiate apoptotic and necrotic death of tumor.

Currently, the efficiency of the photodynamic phase (production of singlet oxygen), as well as the comparison of different PS, are based on three photophysical characteristics:

• the absorption coefficient, σ S , of the PS in the spectral region of the excitation light. σ S characterizes the PS- light interaction; • the quantum yield of the PS triplet state, Φ T and the triplet state lifetime, τ T ; • the quantum yield of singlet oxygen production, Φ ∆ which characterizes the PS-oxygen interaction. These photophysical coefficients were till now separately determined from in vitro dedicated experiments ( [START_REF] Braun | Applications of singlet oxygen reactions: mechanistic and kinetic investigations[END_REF], [START_REF] Bonnett | Photobleaching of 5,10,15,20-tetrakis(m-hydroxyphenyl)porphyrin (m-THPP) and the corresponding chlorin (m-THPC) and bacteriochlorin(m-THPBC). A comparative study[END_REF]), but they cannot always predict the in vivo cytotoxic action.

In this paper, a model-based method, based on system identification techniques [START_REF] Walter | Identification of Parametric Models from experimental data[END_REF], is proposed as an alternative and complementary approach for the estimation of PDT photophysical characteristics. Such an approach could provide two main benefits: [START_REF] Moser | Photodynamic Tumor Therapy: 2nd and 3rd Generation[END_REF] the estimation of the photophysical coefficients directly from in vivo experiments and (2) a significant cut of the experimental cost by decreasing the number of trials. So far, the estimation of the photophysical coefficients by system identification techniques has never been assessed in practice.

Five main difficulties can be addressed for the application of such a solution to in vivo spectro-fluorescence data:

• the nonlinearity structure of the photopohysical model; • the large range of possible values (up to four decades) for the unknown parameters; • some initial state variables are unknown;

• the lack of information (only one measured variable over six state variables); • the limited degrees-of-freedom for the choice of the laser light stimulus (input variable), see paragraph II-C for more information about the experimental setup.

A dynamic model of the photoreaction phase, the relationships between the photophysical characteristics and the model parameters, and the experimental setup are introduced in section II. Section III is devoted to the experimetal modeling methodology while section IV deals with the selection of practically identifiable parameters. The estimation of the model parameters, the model falsification tests and the determination of the parameter uncertainty are presented in sections V, VI and VII respectively. A discussion of the results is provided before to conclude.

II. PROBLEM STATEMENT

As previously mentioned in the introduction, the interaction between the light (at a specific wavelength) and the PS already incorporated in the tumor, determines a sequence of photoreactions which are presumed to have as result the death of the tumor cells (through the production of singlet excited oxygen -a very reactive specie).

A. Model Description

Based on the kinetics equations describing the type-II reactions specific to PDT [START_REF] Moser | Photodynamic Tumor Therapy: 2nd and 3rd Generation[END_REF], [START_REF] Goldman | Photodynamic Therapy[END_REF], [START_REF] Bastogne | Modélisation systéme de la thérapie photodynamique[END_REF], a nonlinear state-space model of the intracellular photoreaction process may be introduced. This model is composed of six state variables and eleven unknown photophysical parameters, and is described as follows,

                                                 d[S 0 ] dt = u P + k f [S 1 ] -k Pb [S 0 ] 1 O 2 + +k p [T 1 ] + k T [T 1 ] 3 O 2 -k A u L [S 0 ] d[S 1 ] dt = k A u L [S 0 ] -k f + k ISC [S 1 ] - -k SM [S 1 ] [M] d[T 1 ] dt = k ISC [S 1 ] -k p [T 1 ] -k T [T 1 ] 3 O 2 - -k T M [T 1 ] [M] d[ 3 O 2 ] dt = u O -k T [T 1 ] 3 O 2 + k l 1 O 2 d[ 1 O 2 ] dt = k T [T 1 ] 3 O 2 -k l 1 O 2 - -k ox [M] 1 O 2 -k Pb [S 0 ] 1 O 2 d[M] dt = -k ox [M] 1 O 2 -k SM [S 1 ] [M] - -k T M [T 1 ] [M] y (t i ) = γ [S 0 ] (t i ) + v (t i ) v (t i ) ∼ N (0, s 2 ) (1)
where

x (t, p) = [S 0 ] , [S 1 ] , [T 1 ] , 3 O 2 , 1 O 2 , [M] T ∈
}. p = k ISC , k Pb , k l , k f , k p , k A , k ox , k T M , k SM , k T , γ T ∈
R 11 are the unknown parameters of the model, with

k f = k F + k IC , k p = k P + k T S and k l = k nr l + k r l .
Tab. I contains the model variables description whereas Tab. II presents all the parameters with their initial guess (see step 2 of the experimental modeling procedure) and the ones of the state variables.

B. Photophysical Parameters

The relationships between the model parameters and the photophysical coefficients are:

k A = σ S hν A ⇒ σ S = h • υ A • k A Φ T = k ISC k ISC +k f ; τ T = 1 k ISC +k f Φ ∆ = Φ T • φ et = k ISC k ISC +k f • k T [ 3 O 2 ] k T [ 3 O 2 ]+kp (2)
where h = 6.026 • 10 -34 J • s is Planck's constant, ν A is the frequency of the incident light, φ et is the efficiency of the energy transfer to molecular oxygen. Note that for a given concentration of molecular oxygen

[ 3 O 2 ], all the

Symbol Definition

Values

k A PS absorption rate at u L = 76mW /cm 2 19 k Pb Photobleaching rate 1.2 • 10 8 mM -1 s -1 k T Biomolecular rate constant for 10 5 the reaction of 3 O 2 with T 1 mM -1 s -1 k f Rate constant for the reaction S 1 → S 0 2 • 10 7 s -1 k p Rate constant for the reaction T 1 → S 0 1250s -1 k ISC Rate constant for the reaction S 1 → T 1 8 • 10 7 s -1 k l Rate constant for the reaction 1 O 2 → 3 O 2 10 6 s -1 k ox Bimolecular rate constant for 2.6 • 10 9 the reactions 1 O 2 with M mM -1 s -1 k SM Bimolecular rate constant for 10 -2 the reaction of S 1 with M mM -1 s -1 k T M Bimolecular rate constant for 1 mM -1 s -1 the reaction of T 1 with M γ measurement coefficient 43.15 a.u. [S 0 ] 0 Initial condition for [S 0 ] 8.5 • 10 -3 mM [S 1 ] 0 Initial condition for [S 1 ] 0 mM [T 1 ] 0 Initial condition for [T 1 ] 0 mM 3 0 2 0 Initial condition for 3 O 2 83 • 10 -3 mM 1 O 2 0 Initial condition for 1 O 2 0 mM [M] 0 Initial condition for [M] 830 • 10 -3 mM
photophysical coefficients can be deduced from p, the vector of the model parameters.

C. Problem statement

As previously pointed out, the in vivo determination of σ S , Φ T , τ T and Φ ∆ can be expressed as a problem of identifia- bility and estimation of p from the in vivo measurement of u L and y.

D. Experimental Setup

Fig. 1 shows the in vivo experimental setup. Female athymic Foxn1 nude mice were used for this study (Harlan, Gannat, France). The mice were used for tumor implantation when they were 7 -9 weeks old (20 -25g). The model of human malignant glioma was obtained using U87 cells, as described in [START_REF] Pinel | Erythropoietin-induced reduction of hypoxia before and during fractionated irradiation contributes to improvement of radioresponse in human glioma xenografts[END_REF]. Tumors were treated when they reached 5 ± 1mm (15 ± 5 days after tumor grafting). The photosensitizing molecule is TPC-Ahx-ATWLPPR [START_REF] Tirand | Metabolic profile of a peptide-conjugated chlorintype photosensitizer targeting neuropilin-1: an in vivo and in vitro study[END_REF] the administrated dose, u P = 1.75µmol • kg -1 , was injected to the mice via the tail vein. The mice were kept in the dark for 4 h, and anesthetized. Irradiation of tumors was carried out at λ = 652nm, using a dye laser, pumped with an argon laser. Light was delivered through an optical fiber and illuminated the tumor surface over a 1cm diameter beam spot. The incident fluence rate on the tumors was previously evaluated using an optical power meter. The fluence rate (irradiance) u L (t) is a square wave signal. Its magnitude reaches 95 mW •cm -2 and its period is equal to T = 1 mn with a 50% duty cycle. The fluorescence intensity y(t) is collected by a spectrofluorimeter only during the dark phases (laser off) and with a sampling period of time equal to T s = 7 s, i.e. three points for each dark half cycle (technical limit of the data acquisition system). The total duration of the experiment is limited to about 42 min.

III. EXPERIMENTAL MODELING METHODOLOGY

The experimental modeling procedure used in this application study has been presented in [START_REF] Bastogne | Selection and identification of physical parameters from passive observation. application to a winding process[END_REF]. The latter approach is composed of six steps, 1) Defining the model structure and assessing its theoretical (or structural) identifiability [START_REF] Walter | Identification of Parametric Models from experimental data[END_REF] deals with the possibility to give an unique value to each parameter of a mathematical model structure. The uniqueness of this solution is assessed in an idealized or theoretical framework where the process and the model have identical structures, the data is noise-free, and where the input signals and the measurement times can be chosen at will. In this application, the global structural identifiability of model ( 1) was checked by the similarity transformation method. 2) A prior analysis, often based on a literature analysis, is carried out to provide the initial guess of p. In order to improve the numerical properties of the Fisher Information Matrix (see section IV), a rescaling of the parameters is done: p i = p 0 i • Θ i , with p 0 initial guess for the parameters, Θ the normalized parameters to be estimated and i = 1, • • • , 11.

3) The third step is devoted to the implementation of the model into a simulation environment. In this case, the simulation model is represented by a block-diagram implemented into Simulink c 1 . 4) Structural identifiability is a qualitative concept. When the model is proved structurally identifiable, an admissible experiment has to be selected to collect the estimation data. The design of optimal experiments may be regarded as maximizing a measure of identifiability, thereby transforming the intitial qualitative problem into a quantitative one, also called practical identifiability [START_REF] Walter | Identifiability of Parametric Models[END_REF], [START_REF] Vanrolleghem | Practical identifiability of a biokinetic model of activated sludge respiration[END_REF]. In a fourth step, a practical identifiability analysis is locally carried out through the evaluation of the output sensitivity functions and the corresponding Fisher Information Matrix [START_REF] Dochain | Dynamical Modelling and Estimation in Wastewater Treatment Processes[END_REF], [START_REF] Bastogne | Selection and identification of physical parameters from passive observation. application to a winding process[END_REF].

The output sensitivity functions are computed by the software Diffedge c which enables the sensitivity analysis of block diagrams by computer algebra [START_REF] Masse | Differentiation, analyze sensibilite and identification of the hybrid models decrits under Simulink[END_REF]. 5) In the next step, the practically identifiable parameters Θ are selected according to numerical properties of an empirical Fisher information Matrix. 6) The sixth step deals with the parameter estimation. The estimation criterion is defined as follows:

J (Θ) = N ∑ j=1 e (t j , Θ) T • e (t j , Θ) (3)
where Θ is a subvector of Θ containing the selected parameters to be estimated and e(t, Θ) is the output error defined by:

e(t, Θ) = y S (t) -y(t, Θ) ∼ N (0, σ 2 ), (4)
with y S (t) the output system measurement. {e(t j , Θ)| j = 1, • • • , N} is assumed to be a realization sequence of independent and identically distributed gaussian variables. {t j | j = 1, • • • , N} is the sequence of the measurement time instants. 7) The last step consists in estimating the empirical probability density functions of the identifiable model parameters in order to characterize their uncertainty.

IV. SELECTION OF IDENTIFIABLE PARAMETERS

The selection of the practically identifiable parameters relies on an empirical Fisher Information Matrix defined by

FIM ∝ N ∑ j=1 ∂ y (t j , Θ) ∂ Θ T ∂ y (t j , Θ) ∂ Θ (5)
The number of practically identifiable parameters is given by the numerical rank of FIM. Given u L and the measurement time instants defined in section II-D, the numerical rank of FIM was equal to 11, meaning that all the parameters are practically identifiable.

V. PARAMETER ESTIMATION

The parameter estimation step was based on a least squares method applied to the output error. The minimization of the estimation criterion (3) with model (1) as a constraint was solved using an Active-Set algorithm implemented in Matlab. Parameter estimates are given in Tab. (III). Figure 2 shows the measured values vs the estimated responses.

VI. MODEL FALSIFICATION

The model falsification test used herein consists in checking the modeling assumption, stated in section II-A, about the gaussian distribution of the output residuals. Two statistical tests were implemented. The first one is a quantilie-quantile plot and the second test is based on the Pearsons chi-square test.

Let us firstly formulate the null hypothesis (H 0 ) and the alternative hypothesis (H 1 ) such that:

H 0 : e ∼ N (0, σ 2 ) H 1 : e ∼ N (0, σ 2 ) (6)
where e denotes a random variable of which an observed sequence is {e(t

i )|i = 1, • • • , N}. A. Q-Q test
A Quantile-Quantile plot is a graphical method for diagnosing differences between the probability distribution of a statistical population from which a random sample has been taken and a comparison distribution. Herein, The Q-Q test The data are plotted against a theoretical normal distribution in such a way that the points should form an approximate straight line. Departures from this straight line indicate departures from normality. The results are represented in fig. 3. The normal probability plot shows a reasonably linear pattern in the center of the data, therefore the null hypothesis H 0 is sustained.

B. Pearsons chi-square test

As shown in Fig. 4, the sequence of output residuals is discretized into 60 bins. The number of degrees of freedom is r = 18. The test statistic computed from the observed frequencies (histogram in Fig. 4) χ 2 = 30.94 is lower than the critical value χ 2 α,r = 32.34. Accordingly, the null hypothesis cannot be rejected for α = 0.02. In other words, the model based on the estimated parameters is not falsified.

VII. CRITICAL ANALYSIS OF THE RESULTS

A possible problem with the proposed approach is that the numerical search may terminate far from the true minimum

p 6 = k 0 A • Θ 6 Θ 1 [0.7661, 1.2636] p 1 = k 0 ISC • Θ 1 Θ 5 [0.8517, 1.3626] p 5 = k 0 p • Θ 5 Θ 7 [0.8227, 1.3471] p 7 = k 0 ox • Θ 7 Θ 4 [0.7653, 1.3682] p 4 = k 0 f • Θ 4 Θ 2 [0.7762, 1.4060] p 2 = k 0 Pb • Θ 2 Θ 3 [0.7806, 1.4539] p 3 = k 0 l • Θ 3 Θ 10 [0.7856, 1.4716] p 10 = k 0 T • Θ 10 Θ 9 [4.6648, 5.4799] p 9 = k 0 SM • Θ 9
for several reasons, such as inadequate choice of the optimization algorithm coefficients or initial parameter guess, or trapping in a local minima.

VIII. UNCERTAINTY ANALYSIS

Because of the low signal-to-noise ratio, determining the optimal values of the parameters Θ with respect to the chosen criterion [START_REF] Walter | Identification of Parametric Models from experimental data[END_REF], [START_REF] Marsili-Libelli | Confidence regions of estimated parameters for ecological systems[END_REF], [START_REF] Pauw | Optimal Experimental Design for Calibration of Bioprocess Models: A Validated Software Toolbox[END_REF] is not enough. Herein, a Monte-Carlo approach was carried out to determine the empirical probability density functions of the parameter estimates. The implemented procedure is the following:

• Based on the experimental data, compute a parameter estimate Θ; • Check the validity of the hypothesis on the noise e (see section VI); • Generate vectors of fictitious data by running the model with parameters Θ for N sim = 1000 realizations of i.i.d.

N (0, σ 2 ) variables v(t i ) with σ 2 = 1 N -1 N ∑ i=1 e (t i ) 2 (7)
• Estimate N sim series of estimated parameters Θ based on each of these vectors of fictitious data ỹ and using the same method as for the experimental data; • Given the N sim estimates of Θ i , compute the 25th and 975th empirical 1000-quantiles (permillages), noted q 0.025 (Θ i ) and q 0.975 (Θ i ). Then, the empirical 95% confidence interval for Θ i is given by [q 0.025 (Θ i ); q 0.975 (Θ i )].

The empirical 95% confidence interval for all the parameters estimates are given in Tab. IV. Fig. 5 illustrates all the histograms obtained with the Monte-Carlo method.

IX. CONCLUSION In this study, the estimation of the photophysical parameters included in photoreactions induced by PDT is examined. This objective has never been reached in an in vivo framework. An original model-based approach, expressed as a system identification problem, is proposed. The latter approach involves to study the practical identifiability of the photoreaction model and the estimation of its parameters. The main difficulties of this inverse problem are (i) the non linearity of the model structure and (ii) the lack of information in the experimental data (small data sets, small sampling rate, no measurement during the treatment phase (Laser ON), few degrees of freedom concerning the choice of the input signal and only one output variable for six state variables). The local practical identifiability test is based on the conditioning of the Fisher Information Matrix. The parameter estimation relies on a least-squares method applied to the output error. Despite the previously mentioned difficulties, the results obtained in practice from in vivo experiments have shown promising results.

X. PERSPECTIVES

Two main developments could enhance the accuracy of the estimates:

• designing new experiment protocols, i.e. seeking the optimal input signal u L and increasing significantly the sampling rate f S during the dark phases (laser OFF) of the test. A new research platform equipped with a laser diode is currently developed in order to both increase the signal-to-nose ratio and the degrees of freedom concerning the choice of the input signal; • using global sensitivity approaches instead of the local approach used herein.

Fig. 2 .

 2 Fig. 2. Measured vs Estimated Responses

Fig. 4 .

 4 Fig. 3. Q-Q plot of output residuals

Fig. 5 .

 5 Fig. 5. Results of the Uncertainty Analysis using Monte Carlo method for the estimated parameters

  and

	u L (t)	t ON	T	fiber	head	Spectro-
					fiber	fluorimeter
				t	lig ht be am	y(t)	u L (t)
		Laser source	u L (t)
		In vivo tumor (nude mouse)			t
						t i
				Fig. 1. in vitro experimental set up

Simulink builds upon Matlab c , The MathWorks, Inc

ACKNOWLEDGMENTS

We gratefully acknowledge L. Denis-Vidal (Université des Sciences et Technologie de Lille, France) for her scientific support in structural identifiability in this study.