Nada El Berdan

Ahmad Makki

Morgan Pierre

FINITE ELEMENT APPROXIMATION OF A SINGULAR ELLIPTIC PROBLEM WITH DISCONTINUOUS COEFFICIENTS

Keywords: finite element method, singular elliptic problem, monotone operator, adaptive mesh refinement. Mathematics Subject Classification: 65N30, 35J70, 47H05, 65N50

We propose a finite element discretization of a singular elliptic problem with discontinuous coefficients. We use a "regularize then discretize" approach. We show that our method converges in 1, 2 and 3 space dimensions. We also perform numerical simulations in two space dimensions with FreeFem++, using an adaptive mesh strategy to deal with the singularity. The simulations confirm the validity of our approach.

Introduction

We consider the following singular elliptic problem,

       -div(A(x)∇u) = c(x) u γ in Ω, u > 0 in Ω, u = 0 on ∂Ω. (1.1)
Here, Ω is a bounded domain of R d , γ > 0, the function c satisfies c ∈ L ∞ (Ω), c ≥ 0, c = 0, (

and A = (a ij ) 1≤i,j≤d is a bounded elliptic matrix, i.e.

a ij ∈ L ∞ (Ω), for all 1 ≤ i, j ≤ d, (1.3)
and there exists α > 0 such that

A(x)ξ • ξ ≥ α|ξ| 2 , (1.4)
for all ξ ∈ R d , for almost every x ∈ Ω. Problem (1.1) is a generalization of the singular elliptic problem

       -∆u = c(x) u γ in Ω, u > 0 in Ω, u = 0 on ∂Ω, (1.5)
which has been extensively studied, along with many of its variants such as (1.1) (see, e.g., [START_REF] Boccardo | Semilinear elliptic equations with singular nonlinearities[END_REF][START_REF] Crandall | On a Dirichlet problem with a singular nonlinearity[END_REF][START_REF] Díaz | Elliptic problems on the space of weighted with the distance to the boundary integrable functions revisited[END_REF][START_REF] Berdan | The uniform Hopf inequality for discontinuous coefficients and optimal regularity in BMO for singular problems[END_REF][START_REF] Godoy | Elliptic problems with singular nonlinearities of indefinite sign[END_REF][START_REF] Gui | Regularity of an elliptic problem with a singular nonlinearity[END_REF] and references therein). Questions related to existence, uniqueness and regularity of solutions have been investigated. Such problems arise in the modeling of boundary layers [START_REF] Callegari | Some singular, nonlinear differential equations arising in boundary layer theory[END_REF][START_REF] Charkaoui | Weak periodic solution for semilinear parabolic problem with singular nonlinearities and L1 data[END_REF][START_REF] Gie | Singular perturbations and boundary layers[END_REF].

If A, c and ∂Ω are regular enough, it is well-known that (1.1) has a unique classical solution u ∈ C 2 (Ω) ∩ C 0 (Ω), which cannot generally belong to C 2 (Ω). For instance, the function u(x) = 1 -x 2 (1.6) satisfies (1.5) on the interval (-1, 1) with c(x) = 1 and γ = 3, that is -u (x) = 1 u 3 in (-1, 1), u > 0 in (-1, 1) and u(-1) = u(1) = 0.

(1.7)

We note that the function defined by (1.6) has the regularity C 1/2 ([-1, 1]) and does not belong to H 1 0 (-1, 1). These two remarks have been generalized to problem (1.5). In particular, if Ω is of class C 2 , if c is Hölder continuous and bounded from below by a positive constant, a result of Gui and Lin [START_REF] Gui | Regularity of an elliptic problem with a singular nonlinearity[END_REF] shows that If 0 < γ < 1, then u ∈ C 1,ν (Ω), where ν = 1 -γ,

If γ = 1, then u ∈ C ν (Ω), for all ν ∈ (0, 1), (1.9)

If γ > 1, then u ∈ C ν (Ω), where ν = 2 1 + γ .

(1.10)

Under the same assumptions, the solution u to (1.5) belongs to H 1 0 (Ω) if and only if γ < 3 (see [START_REF] Lazer | On a singular nonlinear elliptic boundary-value problem[END_REF] and [1, Lemma 2.1]).

In [START_REF] Barrett | Finite element approximation of a semilinear elliptic problem with a singular nonlinearity[END_REF], Barrett and Schneider performed a very nice a priori estimate analysis for a P 1 finite element approximation of (1.5) in 1, 2 and 3 space dimensions. In particular, they obtained a priori estimates in H 1 norm for γ < 3 and L ∞ norm for all γ > 0. They also provided extensive numerical simulations in space dimension 1. They used uniform subdivisions for the computations, but they pointed out in their conclusion that adaptive meshes would be more appropriate. This is illustrated by the example (1.6) or by the result (1.10): since the C ν regularity occurs near the boundary of Ω, it seems natural to use a mesh which is refined near the boundary.

Our purpose here is to propose a finite element approximation of problem (1.1). This is done on a "regularize then discretize" basis. We show that our method converges for d = 1, 2 and 3 and we also perform numerical simulations in two space dimensions with FreeFem++. In order to deal with the singularity, we use an adaptive mesh strategy.

Our paper is organized as follows. We first introduce in Section 2 a regularized version of problem (1.1) by using a local regularization of the nonlinearity s → 1/s γ , as in [START_REF] Barrett | Finite element approximation of a semilinear elliptic problem with a singular nonlinearity[END_REF]. A monotonicity argument shows that the regularized problem has a unique solution u ε in H 1 0 (Ω) (Proposition 2.2). By adapting the approach in [START_REF] Boccardo | Semilinear elliptic equations with singular nonlinearities[END_REF][START_REF] Berdan | The uniform Hopf inequality for discontinuous coefficients and optimal regularity in BMO for singular problems[END_REF], we show in Section 3 that u ε converges to a solution u to problem (1.1) in appropriate sense as ε goes to 0. In particular, when γ ≤ 1, u ε → u weakly in H 1 0 (Ω) (Theorem 3.4) and when γ > 1, u ε → u strongly in L 2 (Ω) (Theorem 3.7). Then, we propose in Section 4 a discretization of the regularized problem. For a fixed ε, we prove that the solution to the discretized problem converges strongly in H 1 0 (Ω) to u ε as the mesh size goes to 0 (Theorem 4.5). The assumptions on the family of triangulations are quite general. In Section 5, several numerical simulations show the validity of our approach.

The regularized problem

Throughout the paper, we assume that γ > 0, c satisfies (1.2) and A satisfies (1.3)- (1.4). In Sections 2 and 3, we only assume that Ω is a bounded domain of R d where d is any positive integer. In Section 4, we make further assumptions on Ω and d.

For the regularized problem, we work in the Sobolev space H 1 0 (Ω) equipped with the Hilbertian norm

v 1 = Ω |∇v| 2 dx 1/2 , which is equivalent to the standard H 1 norm. The norm in L 2 (Ω) is denoted • 0 . The first eigenvalue of -∆ is λ 1 = inf v∈H 1 0 (Ω), v =0 v 2 1 v 2 0 > 0. (2.1)
The duality product between H 1 0 (Ω) and its topological dual

H -1 (Ω) is denoted •, • . Following [1], we define the following C 1 regularization f ε of f (s) = s -γ , for each ε ∈ (0, 1]: f ε (s) = f (s) if s > ε, f (ε) + f (ε)(s -ε) if s ≤ ε. (2.2)
We note that f ε is decreasing and positive on R.

Remark 2.1. Another regularization of f (s) could be fε (s) = f (ε + s). (2.3)
However, as pointed out in [START_REF] Barrett | Finite element approximation of a semilinear elliptic problem with a singular nonlinearity[END_REF], f ε performs numerically better. The reason is that fε is a global regularization which always effects the practical computation, whereas f ε is a local regularization which has a better computational effect, cf. Section 5.2.

The regularized version of problem (1.1) reads

-div(A(x)∇u ε ) = c(x)f ε (u ε ) in Ω, u ε = 0 on ∂Ω. (2.4)
We have Proposition 2.2. Problem (2.4) has a unique solution u ε ∈ H 1 0 (Ω). Moreover, u ε > 0 a.e. in Ω.

Proof. We use a monotonicity argument. We consider the nonlinear operator

B ε v = Lv -c(•)f ε (v)
, where Lv = -div(A∇v).

(2.5)

We will show that i) B ε is a nonlinear continuous operator from

H 1 0 (Ω) into its dual H -1 (Ω). ii) B ε is strictly monotone. iii) B ε is coercive, i.e. lim v 1 →+∞ B ε v, v v 1 → +∞.
It follows from the theory of maximal monotone operators (see, e.g. [START_REF] Zeidler | Nonlinear functional analysis and its applications[END_REF]Theorem 26.A]) that there exists a unique u ε ∈ H 1 0 (Ω) such that B ε u ε = 0. The strong maximum principle [START_REF] Gilbarg | Elliptic partial differential equations of second order[END_REF]Theorem 8.19] implies that u ε > 0 on Ω.

We start with i). The linear operator L which is defined by duality through

Lv, w = Ω A(x)∇v • ∇wdx is continuous from H 1 0 (Ω) into H -1 (Ω) since A ∈ L ∞ (Ω; R d×d ). By definition (2.2) of f ε , we see that f ε is bounded on [ε, ∞) and affine on (-∞, ε]. Thus, there exist positive constants C 1 , C 2 (which depend on ε) such that |f ε (s)| ≤ C 1 |s| + C 2 (s ∈ R).
(2.6)

As a consequence, the operator v → a(•)f ε (v) is continuous from L 2 (Ω) into L 2 (Ω)
(see, e.g., [17, 16.1 Lemme]). This shows assertion i).

Next, we see that for all v, ṽ ∈ H 1 0 (Ω) with v = ṽ, we have

B ε v -B ε ṽ, v -ṽ = Ω (A(x)∇(v -ṽ)) • ∇(v -ṽ) - Ω a(x)[f ε (v) -f ε (ṽ)](v -ṽ)dx ≥ α v -ṽ 2 1 > 0,
where in the second line, we used (1.4). This proves assertion ii).

For the coercivity of B ε , we note that

-sf ε (s) = γs 2 ε -(γ+1) -s(1 + γ)ε -γ if s ≤ ε, -s 1-γ if s > ε, (2.7)
and we consider three cases, namely γ > 1, γ = 1 and γ < 1.

If γ > 1, then we have

-sf ε (s) ≥ -C 3 , ∀s ∈ R, (2.8)
for some positive constant C 3 which depends on ε. From (1.4), we deduce that

B ε v, v = Ω A(x)∇v • ∇vdx - Ω c(x)f ε (v(x))v(x)dx ≥ α v 2 1 -c L ∞ (Ω) |Ω|C 3 , (2.9
)

for all v ∈ H 1 0 (Ω). If γ ≤ 1, the minimum value of -sf ε (s) on (-∞, ε] is attained at s ε = ε. Thus, -sf ε (s) ≥ -s ε f ε (s ε ) = -ε 1-γ ≥ -1, ∀s ∈ (-∞, ε],
(2.10) since ε ∈ (0, 1]. If γ = 1, the estimate (2.10) holds for all s ∈ R (see (2.7)) and we recover the estimate (2.9) but with a constant C 3 = 1 independent of ε.

If γ < 1, we deduce from Young's inequality when s > ε and from (2.10) when s ≤ ε that for all β > 0, there exists C β > 0 independent of ε such that

-sf ε (s) ≥ -βs 2 -C β (s ∈ R).
(2.11)

Using (1.4) and (2.1), we obtain

B ε v, v = Ω A(x)∇v • ∇vdx - Ω c(x)f ε (v(x))v(x)dx ≥ α 2 v 2 1 + αλ 1 2 v 2 0 -c L ∞ (Ω) |Ω|(β v 2 0 + C β ),
for all v ∈ H 1 0 (Ω). Thus, for β > 0 small enough, we have

B ε v, v ≥ α 2 v 2 1 -C 4 , (2.12)
for some constant C 4 independent of ε.

Summing up, we have proved that (2.12) holds for a constant C 4 independent of ε ∈ (0, 1] when γ ≤ 1, and for a constant which depends on ε when γ > 1. This shows assertion iii) and the proof is complete.

The limit ε → 0

In order to understand the behaviour of (u ε ) as ε goes to 0, we adapt the approach in [START_REF] Boccardo | Semilinear elliptic equations with singular nonlinearities[END_REF][START_REF] Berdan | The uniform Hopf inequality for discontinuous coefficients and optimal regularity in BMO for singular problems[END_REF]. We first have:

Proposition 3.1. If 0 < ε ≤ δ ≤ 1, then u ε ≥ u δ a.e. in Ω. ( 3.1)
For all ω ⊂⊂ Ω, there exists c ω independent of ε such that

u ε (x) ≥ c ω > 0 for a.e. x ∈ ω, for every ε ∈ (0, 1]. (3.2)
Proof. Let 0 < ε ≤ δ ≤ 1. We have, in the sense of distributions,

-div(A(x)∇u δ ) = c(x)f δ (u δ ), -div(A(x)∇u ε ) = c(x)f ε (u ε ), (3.3) so that -div(A(x)∇(u δ -u ε )) = c(x)[f δ (u δ ) -f ε (u ε )]. (3.4) We choose (u δ -u ε ) + (which belongs to H 1 0 (Ω)
) as a test function, and we note that

[f δ (u δ ) -f ε (u ε )] (u δ -u ε ) + ≤ 0. (3.5) Indeed, f δ is nonincreasing so that [f δ (u δ ) -f δ (u ε )] (u δ -u ε ) + ≤ 0, and δ → f δ (s) is nonincreasing for all s, so that [f δ (u ε ) -f ε (u ε )] (u δ -u ε ) + ≤ 0.
On summing these two inequalities, we obtain (3.5). Thus, the choice (u δ -u ε ) + as a test function in (3.4) yields

α (u δ -u ε ) + 2
1 ≤ 0, thanks to assumption (1.4). We conclude that (u δ -u ε ) + = 0, which shows (3.1).

Next, we note that for all ε > 0, u ε belongs to L ∞ (Ω). This follows from [20, Théorème 4.2] since the right-hand side of (3.3) belongs to L ∞ (Ω). In particular, for ε = 1, there exists a constant C such that

u 1 L ∞ (Ω) ≤ C,
and we have

-div(A(x)∇u 1 ) = c(x)f 1 (u 1 ) ≥ c(x)f 1 ( u 1 L ∞ (Ω) ) ≥ c(x)f 1 (C).
Since c(x)f 1 (C) is not identically zero, the strong maximum principle [START_REF] Gilbarg | Elliptic partial differential equations of second order[END_REF]Theorem 8.19] implies that u 1 > 0 a.e. in Ω and that (3.2) holds for ε = 1. On choosing δ = 1 in (3.1), we obtain that (3.2) holds for all ε ∈ (0, 1].

In order to pass to the limit ε → 0, we need a priori estimates. We consider two cases, namely γ ≤ 1, in which case uniqueness holds [9, Remark 3], and γ > 1.

3.1. The case γ ≤ 1. We have:

Lemma 3.2. If γ ≤ 1, the family (u ε ) is bounded in H 1 0 (Ω). Proof. We choose v = u ε in (2.12). Since B ε u ε = 0, this yields 0 ≥ α 2 u ε 2 1 -C 4 .
The constant C 4 is independent of ε since γ ≤ 1. This proves the claim.

We recall the following representation formula.

Lemma 3.3 (see [3]). Let v ∈ L 1 loc (Ω) ∩ H -1 (Ω) such that v ≥ 0 a.e. in Ω. Then for all ϕ ∈ H 1 0 (Ω), we have vϕ ∈ L 1 (Ω) and v, ϕ = Ω vϕdx. (3.6)
Proof. Let ϕ ∈ H 1 0 (Ω). If ϕ ≥ 0 a.e. in Ω, the result is a consequence from [START_REF] Brézis | A property of Sobolev spaces[END_REF]. In the general case, we write ϕ = ϕ + -ϕ -where ϕ + and ϕ -are nonnegative functions in H 1 0 (Ω).

Our convergence result reads as follows.

Theorem 3.4. Assume that γ ≤ 1. Then there exists a unique solution u ∈ H 1 0 (Ω) to (1.1) in the sense that

∀ ω ⊂⊂ Ω, ∃c ω > 0 such that u(x) ≥ c ω > 0, for a.e. x ∈ ω, (3.7) c u γ ∈ L ∞ loc (Ω) ∩ H -1 (Ω), (3.8
)

Ω A(x)∇u • ∇ϕdx = Ω c(x) u γ ϕdx, ∀ϕ ∈ H 1 0 (Ω). (3.9)
Moreover, u ε → u weakly in H 1 0 (Ω) and strongly in L 2 (Ω), as ε → 0. We note that the right-hand side of (3.9) is well-defined because we may apply Proof. We first prove uniqueness of u. If ũ ∈ H 1 0 (Ω) is another solution of (1.1) in the sense of (3.7)-(3.9), then u -ũ satisfies

Ω A(x)∇(u -ũ) • ∇ϕdx = Ω c(x) 1 u γ - 1 ũγ ϕdx,
for all ϕ ∈ H 1 0 (Ω). We choose ϕ = u -ũ as a test function, we use that f (s) = s -γ is decreasing and assumption (1.4). This yields

α u -ũ 2 1 ≤ Ω A(x)∇(u -ũ) • ∇(u -ũ)dx ≤ 0.
Thus, ũ = u as claimed.

Next, we consider the behaviour (u ε ) as ε tends to 0. By Lemma 3.2, (u ε ) is bounded in H 1 0 (Ω). Thus, up to a subsequence, (u ε ) converges weakly in H 1 0 (Ω) to a function u ∈ H 1 0 (Ω). We know by Rellich's theorem that L 2 (Ω) is compactly imbedded into H 1 0 (Ω), so that (up to a subsequence) u ε → u strongly in L 2 (Ω) and a.e. in Ω. Proposition 3.1 implies that (3.7) holds for u. In turn, this implies that u > 0 a.e. in Ω and

c u γ ∈ L ∞ loc (Ω). Now, we use that B ε u ε = 0, that is Ω A∇u ε • ∇ϕdx = Ω c(x)f ε (u ε )ϕdx, ∀ϕ ∈ H 1 0 (Ω). (3.10) Let ϕ ∈ C ∞ c (Ω) in (3.10). Thanks to (3.1), f ε (u ε ) is bounded in L ∞ loc (Ω) and f ε (u ε ) → f (u) a.e.
in Ω. Since ϕ has compact support in Ω, we may pass to the limit in (3.10). We obtain that (3.9) holds for all ϕ ∈ C ∞ c (Ω). Equation (3.10) also shows that

c(x)f ε (u ε ) belongs to a bounded set of H -1 (Ω) since Ω c(x)f ε (u ε )ϕdx ≤ C u ε 1 ϕ 1 ≤ C ϕ 1 , for all ϕ ∈ H 1 0 (Ω)
, where the constants C and C are independent of ε. This shows that c(x)f (u) ∈ H -1 (Ω). Finally, we use that C ∞ c (Ω) is dense in H 1 0 (Ω) and the representation formula of Brezis and Browder (Lemma 3.3), and we find that equation (3.9) holds for all ϕ ∈ H 1 0 (Ω). The uniqueness of the limit u implies that the whole family (u ε ) converges to u, weakly in H 1 0 (Ω) and strongly in L 2 (Ω). This concludes the proof.

Remark 3.5. The results in [START_REF] Boccardo | Semilinear elliptic equations with singular nonlinearities[END_REF][START_REF] Berdan | The uniform Hopf inequality for discontinuous coefficients and optimal regularity in BMO for singular problems[END_REF] show that u ) is bounded in

H 1 0 (Ω) and (u ε ) is bounded in H 1 loc (Ω).
As a consequence, we can state the following convergence result: Theorem 3.7. Assume that γ > 1 and let (ε n ) be a decreasing sequence of positive real numbers which tends to 0. Then the pointwise limit u of (u εn ) solves (1.1) in the sense that

∀ ω ⊂⊂ Ω, ∃c ω > 0 such that u(x) ≥ c ω > 0, for a.e. x ∈ ω, (3.11)
u ∈ H 1 loc (Ω), (3.12)

Ω A(x)∇u • ∇ϕdx = Ω c(x) u γ ϕdx, ∀ϕ ∈ C 1 c (Ω), (3.13)
u 1+γ 2
∈ H 1 0 (Ω) (this is the meaning of u = 0 on ∂Ω).

(3.14)

Moreover, u εn → u weakly in H 1 (ω), for all ω ⊂⊂ Ω, u 1+γ 2 εn → u 1+γ 2
weakly in H 1 0 (Ω) and u εn → u strongly in L 2 (Ω).

Proof. We note that by (3.1), the sequence (u εn ) n is nondecreasing so that its pointwise limit u(x) = lim n→+∞ u εn (x) is well defined a.e. in Ω. Property (3.11) is a consequence of (3.2). Let ϕ ∈ C 1 c (Ω) and let ω = {x ∈ Ω : ϕ(x) = 0}.

We note that ω ⊂⊂ Ω. By Lemma 3.6, (u εn ) is bounded in H 1 (ω), so that, up to a subsequence, u εn → u weakly in H 1 (ω), strongly in L 2 (ω) and a.e. in Ω. The uniqueness of the limit u implies that the whole sequence (u εn ) converges.

By definition of u εn , we have

Ω A(x)∇u εn • ∇ϕdx = Ω c(x)f εn (u εn )ϕdx. (3.15)
Since ω ⊂⊂ Ω, we may pass to the limit in (3.15) and we obtain that u satisfies (3.13).

Finally, we know by Lemma 3.6 that (u

1+γ 2
εn ) is bounded in H 1 0 (Ω). This shows that u satisfies (3.14) and that (u

1+γ 2
εn ) converges weakly in H 1 0 (Ω) to u 1+γ 2 . Assume now that d ≥ 3 and let s < 2 = 2d d-2 . Thanks to the Sobolev imbeddings [START_REF] Gilbarg | Elliptic partial differential equations of second order[END_REF], we know that

H 1 0 (Ω) is compactly imbedded in L s (Ω), so (u 1+γ 2 εn ) → u 1+γ 2
strongly in L s (Ω). There exists g ∈ L s (Ω) such that, up to a subsequence, (u

1+γ 2 εn ) → u 1+γ 2
a.e. in Ω and 0 ≤ u (Ω). The dominated convergence theorem shows that u εn → u in L q (Ω), with q = s(1+γ) 2 < d(1+γ) d-2 . Since d(1+γ) d-2 > 2 and s < 2 is arbitrary, we may choose q = 2. Note that by uniqueness of the limit a.e., the whole sequence converges.

If d = 2, a similar argument holds for all s < +∞ and we conclude similarly.

If d = 1, then H 1 0 (Ω) is compactly imbedded in C 0 (Ω), so (u 1+γ 2 εn ) → u 1+γ 2
strongly in C 0 (Ω). In particular, there exists M > 0 such that 0 ≤ u εn (x) ≤ M for all x ∈ Ω and for all n. Since the function t → t 2 1+γ is uniformly continuous on [0, M ], we conclude that u εn → u strongly in C 0 (Ω), and therefore in L 2 (Ω). This concludes the proof.

The finite element approximation

Throughout Section 4, we assume that ε ∈ (0, 1] is fixed and d = 1, 2 or 3. For the finite element approximation of the regularized problem (2.4), we assume that Ω is a bounded convex d-polyhedron, that is a bounded interval if d = 1, a convex polygon if d = 2 and a convex polyhedron if d = 3.

We consider a regular family {T h } h>0 of triangulations of Ω into d-simplices and such that each triangulation covers Ω exactly. For each triangulation T h , the P 1 finite element space is

V h = {v h ∈ C 0 (Ω), : v h is affine on every d-simplex of T h and v h = 0 on ∂Ω}. (4.1)
The space V h is a finite dimensional subspace of V = H 1 0 (Ω). We denote I h the interpolation operateur at the vertices of the triangulation. It is well-known [START_REF] Ern | Éléments finis: théorie, applications, mise en oeuvre[END_REF] that for every v ∈ H 2 (Ω), we have

I h v → v in V, as h → 0. (4.2) Since H 2 (Ω) is dense in V , this implies inf v h ∈V h u -v h 1 → 0, (4.3)
for all u ∈ V .

Remark

h ∈ V h such that Ω A(x)∇u h • ∇v h dx = Ω c(x)f ε (u h )v h dx, ∀v h ∈ V h . (4.4)
We have Proposition 4.2. Problem (4.4) has a unique solution.

Remark 4.3. The function u h = u ε,h defined in (4.4) depends also on ε. Since ε ∈ (0, 1] is fixed throughout Section 4, we omit the subscript ε.

Proof. We use a monotonicity argument, as in the continuous case. Let B h : V h → V h be the nonlinear operator defined by

Ω (B h v h )w h dx = Ω A(x)∇v h • ∇w h - Ω c(x)f ε (v h )w h dx. (4.5)
The same arguments as in the proof of Proposition 2.2 show that B h is continuous, strictly monotone and coercive. Thus, there exists a unique u h ∈ V h such that B h u h = 0.

Remark 4.4. We do not know if u h ≥ 0 in every case. If a discrete maximum principle holds then u h ≥ 0, but this requires additional assumptions on the triangulation and the matrix A [START_REF] Ciarlet | Maximum principle and uniform convergence for the finite element method[END_REF]. If Lv = -∆v and if the triangulation is weakly acute, then a discrete maximum principle holds and we have u h ≥ 0, see [START_REF] Barrett | Finite element approximation of a semilinear elliptic problem with a singular nonlinearity[END_REF][START_REF] Jiang | Effect of numerical integration for elliptic obstacle problems[END_REF][START_REF] Vanselow | About Delaunay triangulations and discrete maximum principles for the linear conforming FEM applied to the Poisson equation[END_REF]. If Lv = -∆v and d = 1, the discrete maximum principle holds without any additional assumption on the mesh [START_REF] Vejchodský | Discrete maximum principle for higher-order finite elements in 1[END_REF].

We recall that the parameter ε > 0 is fixed throughout this section. The following convergence holds: Theorem 4.5. We have u h → u ε strongly in V as h → 0.

Proof. We notice that

Ω (B h v h )w h dx = B ε v h , w h , ∀v h , w h ∈ V h , (4.6)
where B h and B ε are defined in (4.5) and (2.5), respectively. We apply this relation with v h = w h = u h . Using B h u h = 0 and (2.12), we obtain

0 ≥ α 2 u h 2 1 -C 4 .
This shows that the sequence (u h ) is bounded in V .

We may therefore assume that, up to a subsequence, u h converges weakly in V to some function ũ. Let ϕ ∈ C ∞ c (Ω). We choose v h = u h and w h = I h ϕ in (4.6). We find

0 = B ε u h , I h ϕ = Ω A(x)∇u h • ∇(I h ϕ)dx - Ω c(x)f ε (u h )I h (ϕ)dx. (4.7)
The quadratic term tends to Ω A(x)∇ũ•∇ϕdx thanks to (4.2). Since V is compactly imbedded into L 2 (Ω), we know that u h → ũ strongly in L 2 (Ω). Using the continuity of the nonlinear operator

v → c(•)f ε (v) from L 2 (Ω) into L 2 (Ω) (cf. (2.6
)), we may pass to the limit in (4.7) and we find

0 = B ε ũ, ϕ .
This is true for all ϕ ∈ C ∞ c (Ω). Thus, ũ = u ε and by uniqueness of the limit, the whole sequence (u h ) converges to ũ, weakly in V .

Next, we denote g(x) = c(x)f ε (u(x)) and g h (x) = c(x)f ε (u h (x)). We have already seen that g h → g in L 2 (Ω). The definition of u h and u ε reads

Ω A(x)∇u h • ∇v h dx = Ω g h v h dx, ∀v h ∈ V h , Ω A(x)∇u • ∇vdx = Ω gvdx, ∀v ∈ V.
We introduce the solution ũh ∈ V h to the linear problem

Ω A(x)∇ũ h • ∇v h dx = Ω gv h , ∀v h ∈ V h .
The assumptions on A ensure that ũh exists and is unique. Cea's lemma [START_REF] Ern | Éléments finis: théorie, applications, mise en oeuvre[END_REF] and (4.3) imply that ũh → u strongly in V . On the other hand, we have

α u h -ũh 2 1 ≤ Ω A(x)∇(u h -ũh ) • ∇(u h -ũh )dx = Ω (g h -g)(u h -ũh )dx ≤ g h -g 0 u h -ũh 0 , so u h -ũh 1 ≤ 1 αλ 1/2 1 g -g h 0 ,
which shows that u h -ũh 1 → 0. On writing u h = (u h -ũh ) + ũh , we obtain that u h → u strongly in V , as claimed.

Numerical simulations

We performed numerical simulations with the FreeFem++ software [START_REF] Hecht | New development in FreeFem++[END_REF] in two space dimension. We used a mesh adaptation algorithm available in FreeFem++ with the keyword adaptmesh. Our algorithm reads as follows, for a bounded domain Ω of R 2 with smooth boundary and a given ε > 0.

Algorithm 5.1. Input: an initial mesh T 0 h which gives a triangulation of Ω h (a polygonal approximation to Ω) and an initial guess ũ0 h of the solution u. For k from 0 to k max :

(1) build the space

V k h ⊂ H 1 0 (Ω h ) of continuous P 1 finite element associated to T k h , cf. (4.1); (2) Compute the solution u k h ∈ V k h to the variational problem Ω A(x)∇u k h • ∇v k h dx = Ω c(x)f ε (u k h )v k h dx, ∀v k h ∈ V k h ,
using a Newton algorithm with the initial value ũk

h ; (3) Build a new triangulation T k+1 h of Ω h adapted to the solution u k h ; (4) Compute the P 1 interpolate ũk+1 h of u k h on T k+1 h . Output: T kmax h and u kmax h .
We recall that f ε is the regularization of f (s) = s -γ given by (2.2) (except in Table 2 where we choose fε defined by (2.3)). For a polygonal domain Ω, we have Ω h = Ω but in the simulations below we have Ω h = Ω because the boundary ∂Ω is of class C 2 .

5.1.

A radial test case for A = Id. We first consider problem (1.5) on the ring Ω = {(x, y) ∈ R 2 : e -2 < x 2 + y 2 < e 2 }, with c(x, y) = 1 x 2 + y 2 and γ = 3.

(

Using polar coordinates, it is easy to check that the exact solution to this problem is

u(x, y) = 1 -log 2 ( x 2 + y 2 ) 1/2 . This solution belongs to C ∞ (Ω) ∩ C 1/2 (Ω) (cf. (1. 10
)), and it does not belong to

H 1 0 (Ω) since γ > 3 [1, Lemma 2.1].
The exterior (interior, respectively) circle of Ω with radius e (1/e, resp.) is approximated by a regular polygon having 170 (46, resp.) vertices; its edges have a length of approximately 0.1 (0.05, resp.). This yields a polygonal domain Ω h which approximates the ring Ω. Since Ω is not convex, Ω h ⊂ Ω and the approximation is nonconforming, that is H 1 0 (Ω h ) ⊂ H 1 0 (Ω). We performed Algorithm 5.1 with k max = 5 mesh adaptations, starting with an initial triangulation T 0 h of Ω h and the initial guess ũ0 h = 0. The FreeFem++ command in step (3) was Th=adaptmesh(Th,u,err=0.005,nbvx=9e4). This command means that T k+1 h is built using T k h and u k h with a P 1 interpolation error level of 0.005 and a maximum number of vertices equal to 9e4.

The final mesh T 5 h obtained for ε = 0.1 is shown in Figure 1. It has 11084 vertices and it has been obtained by successive mesh refinements by starting from an initial mesh T 0 h with 4181 vertices. The mesh T k h at iteration k = 0, In Table 1, we present the L 2 error u-u h L 2 (Ω h ) and the L ∞ error u-u h L ∞ (Ω h ) obtained for different value of ε (with u h = u 5 h ). The discretization parameters are the same as for the case ε = 0.1, which means that the initial triangulation T 0 h is the same, the initial guess is ũ0 h = 0 and the parameters in the adaptmesh command are the same. We checked numerically that the solution u h is nonnegative on Ω h , up to computer accuracy. For each ε, Table 1 gives the number of vertices in the final mesh T kmax h , as well as the minimum and maximum length of the edges of this triangulation.

On writing

u -u h = (u -u ε ) + (u ε -u ε,h ), (5.2)
the error u -u ε,h can be interpreted as an error u -u ε due to the regularization and an error u ε -u ε,h due to the discretization (we write u h = u ε,h to stress that u h depends on ε, cf. Remark 4.3).

In line 5 of Table 1, we have computed for each ε the ratio of the error in the L ∞ norm for 2ε and ε, namely

r ∞,ε = u -u 2ε,h L ∞ (Ω h ) u -u ε,h L ∞ (Ω h ) .
We see that for ε varying from 0.2 down to 0.025, the ratio r ∞,ε is very close to 2. This means that the error is of order O(ε) and the error due to the regularization is dominant in (5.2). For ε between 0.025 and 0.00625, the error due to the discretization becomes dominant.

In comparison, the L 2 error decreases more rapidly than the L ∞ error for ε between 0.4 and 0.025 (cf. line 3 in Table 1, which represents the same ratio as r ∞,ε , but in L 2 norm). This can be related to the choice of a local regularization f ε for f : the error u(x, y) -u h (x, y) is mostly located near the boundary ∂Ω (cf. Remark 2.1). ε 0.1 0.05 0.025 0.0125 6.25e-3 3.125e-3 1. We see that for ε from 0.1 down to 6.25e-3, the L ∞ error and the L 2 error are both dominated by the regularization error which is of order O(ε). The L 2 error is much larger than for the regularization f ε because we use a global regularization of f : the error u(x, y) -u h (x, y) is located everywhere on Ω and not just near the boundary. This is confirmed by the last line of Table 2, which shows for each ε the normalized ratio

u -u ε,h L 2 (Ω h ) |Ω| 1/2 u -u ε,h L ∞ (Ω h ) .
This ratio is close to 1 as long as the regularization error dominates.

5.3.

A radial test case with regular coefficients. We consider now problem (1.1) on the ring Ω = {(x, y) ∈ R 2 : 1 < x 2 + y 2 < 2 2 }, with A(x, y) = 1

x 2 + y 2 1 0 0 1 , c(x, y) = 1 4 x 2 + y 2 and γ = 3.

(5.

3)

The exact solution to this problem is

u(x, y) = v( x 2 + y 2 ), (5.4)
where

v(r) = [(r -1)(2 -r)] 1/2 . Indeed, it is straighforward to check that -v (r) = 1 4v 3 (r) for r ∈ (1, 2), and v(1) = v(2) = 0.
The ODE can be rewritten

- 1 r ∂ ∂r ra(r) ∂v(r) ∂r = 1 (4r)v 3 (r) for r ∈ (1, 2), with a(r) = 1 r . (5.5)
For a radial function u, this corresponds to the PDE (1.1) with the values (5.3). Since the coefficients are regular, the solution u to (1.1) is unique in the class C 2 (Ω)∩ C 0 (Ω) [START_REF] Crandall | On a Dirichlet problem with a singular nonlinearity[END_REF]. It is therefore radial. It is easy to check that u given by (5.4) belongs to For the numerical resolution of (5.3), we proceeded as previously. The interior and exterior circle of radius 1 and 2 were both approximated by a regular polygon with 251 vertices, yielding a polygonal approximation Ω h of Ω (with Ω h ⊂ Ω since Ω is not convex). The sides of the small polygon had a length of approximately 0.025, whereas the sides of the large polygon had a length of approximately 0.05.

C ∞ (Ω) ∩ C 1/2 (
We performed Algorithm 5.1 with k max = 5 iterations, starting with an initial triangulation of Ω h and the initial guess ũ0 h = 0. The parameters in the adaptmesh command were set as previously. Table 3 is the analogue of Table 1 for this test case.

Due to the regularity of the coefficients, the results are very similar to the first test case. In particular, for ε from 0.4 down to 0.05, the L ∞ error is O(ε) and the error due to the regularization is dominant. Then, for smaller values of ε, the error due to the discretization is dominant.

5.4.

A numerical simulation with a discontinuous matrix. In order to illustrate our method in the case of a discontinuous elliptic matrix, we computed the numerical solution obtained on the unit disc Ω in R 2 for the values A(x, y) = a(x, y) 1 0 0 1 , c(x, y) = 1 and γ = 1.1, (5.6) with a(x, y) = 0.2 if x 2 + y 2 < 0.2 2 , a 1 (x, y) if x 2 + y 2 ≥ 0.2 2 , (

where

a 1 (x, y) =         
1.05 -y if x ≥ 0 and y ≥ 0, 2.10 -2x if x ≥ 0 and y < 0, 1.05 + x if x < 0 and y ≥ 0, 2.10 + 2y if x < 0 and y < 0.

(5.8)

We note that a is discontinuous accross the circle of radius 0.2 and along the x and y axis when 0.2 < |x| < 1 or 0.2 < |y| < 1.

We use Algorithm 5.1 with ε = 0.05, k max = 5 and ũ0 h = 0. The unit circle is approximated by a regular polygon with 314 sides having a length of approximately 0.02. This yields a polygonal approximation Ω h of Ω and an initial triangulation of Ω h . Since Ω is convex, we have Ω h ⊂ Ω and we are in the situation described in Remark 4.1. The parameters in the adaptmesh command are the same as previously except for err whose value is set to 0.0025.

The output mesh T kmax h is shown in Figure 2. We observe that it is refined near the unit circle, due to the singularity in the PDE, but also along the discontinuities of the matrix A. The corresponding output solution u kmax h is shown in Figure 3.

Lemma 3 .

 3 3 to the nonnegative function v = c u γ .

1+γ 2 ∈

 2 H 1 0 (Ω) and u ∈ L ∞ (Ω), for all γ ∈ (0, 1]. 3.2. The case γ > 1. By arguing as in [2, Lemma 4.1], we have: Lemma 3.6. Assume that γ > 1. Then (u 1+γ 2 ε

1+γ 2 εn≤ 2 1+γ

 22 g, for all n. Thus, 0 ≤ u εn ≤ g

Figure 1 .

 1 Figure 1. Mesh for case (5.1) and ε = 0.1 (γ = 3)

Figure 2 .

 2 Figure 2. Output mesh for case (5.6) and ε = 0.05 (γ = 1.1)

Figure 3 .

 3 Figure 3. Output solution for case (5.6) and ε = 0.05 (γ = 1.1)

  4.1. Similar properties hold if we assume that Ω is a bounded convex open set with C 2 boundary. In this case, the d-polyhedron Ω h defined by the triangulation satisfies Ω h ⊂ Ω and the vertices of the boundary of Ω h lie on the boundary of Ω. We use thatV h ⊂ H 1 0 (Ω h ) ⊂ H 1 0 (Ω) by setting v h = 0 on Ω \ Ω h . Properties (4.2)-(4.3) are still true (see, e.g.,[START_REF] Pierre | Convergence of exponential attractors for a finite element approximation of the Allen-Cahn equation[END_REF] Section 4.1]), so the results in this section are also valid in such a case.

	The finite element approximation of problem (2.4) reads: find u

Table 1 .

 1 1, 2, 3, 4 and 5 contains 4181, 9356, 11758, 11366, 11097 and 11084 vertices respectively. L 2 error, L ∞ error and mesh characteristics for case(5.1)

	ε	0.4	0.2	0.1	0.05	0.025 0.0125 0.00625
	L 2 error	0.18	0.053	0.014 0.0036 0.0026 0.0029 0.0031
	L 2 ratio	•	3.4	3.8	3.9	1.4	0.9	0.9
	L ∞ error	0.090	0.045	0.023	0.012 0.0057 0.0035 0.0048
	L ∞ ratio	•	2.0	2.0	1.9	2.1	1.6	0.7
	Nb vertices	6790	8877	11084 13350 17370 21121	50717
	h min	2.7e-03 9.3e-04 3.3e-04 1.5e-04 6.7e-05 3.9e-05 2.7e-05
	h max	0.28	0.26	0.31	0.29	0.28	0.29	0.28

Table 2 .

 2 L 2 and L ∞ error for the regularization fε and case (5.1) 5.2. About the choice of the regularization. In Table2, we have computed the L 2 and L ∞ error with the regularization fε defined by (2.3) instead of f ε given by (2.2). The test case is (5.1) and the algorithm and the discretization parameters are exactly as previously. We have divided the values of ε by 4 in order to have a L ∞ error which is easily comparable to the values in Table

	1.5625e-3

Table 3 .

 3 L 2 error, L ∞ error and mesh characteristics for case(5.3)

	ε	0.4	0.2	0.1	0.05	0.025 0.0125 0.00625
	L 2 error	0.17	0.049	0.020 0.016 0.016	0.016	0.016
	L 2 ratio	•	3.5	2.5	1.3	1.0	1.0	1.0
	L ∞ error	0.098	0.046	0.023 0.012 0.0098 0.0099 0.0099
	L ∞ ratio	•	2.1	2.0	1.9	1.2	1.0	1.0
	Nb vertices	8652	12387 15633 19520 24442 36397	69629
	h min	1.0e-02 3.4e-03 1.2e-03 4.5e-4 1.8e-4 8.1e-05 4.7e-05
	h max	0.22	0.25	0.25	0.25	0.25	0.24	0.25

Ω) and that it does not belong to H 1 0 (Ω).

Lebanese University, Nabatieh, Lebanon.

Laboratoire de Mathématiques et Applications, Université de Poitiers, CNRS, F-86073 Poitiers, France1


Johannes Karges

Yih Shi Kuang

Ching Ong

Hui Chao

Gilles Gasser
email: gilles.gasser@chimieparistech.psl.eu

] Yih

Photon Phototherapeutic Effects of Ru(II) Polypyridine Complexes in the Hypoxic Centre of Large Multicellular Tumour Spheroids and Tumour-Bearing Mice

Keywords: Bioinorganic Chemistry, Medicinal Inorganic Chemistry, Metals in Medicine, Photodynamic Therapy, Photosensitizers

come   L'archive ouverte pluridisciplinaire

INTRODUCTION

Cancer has emerged as one of the deadliest diseases worldwide. Photodynamic Therapy (PDT) is a minimal invasive medical technique to treat this disease, often in combination with other methods (i.e., surgery, chemotherapy, radiation therapy or immunotherapy). In PDT, a preferably non-toxic photosensitizer (PS) is activated at a specific wavelength to generate reactive oxygen species (ROS). The majority of clinically approved PSs act by an energy transfer of the PS to molecular oxygen ( 3 O2) to generate singlet oxygen ( 1 O2). As ROS and 1 O2 are highly reactive, they can rapidly interact with their biological surroundings to trigger cell death. [1] Despite recent research advances, the clinical applications of PDT remains sometimes unsatisfactory. Many tumours and especially large tumours show, after a PDT treatment, a regression of their outer spheres while their core remain intact, causing, after a certain time, a cancer relapse. [2] This incomplete tumour eradication is due to several factors: 1) the PS is not able to penetrate a 3D cellular architecture; 2) the oxygen concentration in the centre of the tumour is low (i.e., hypoxia), which is hampering the production of the therapeutically necessary ROS and 1 O2; 3) light is not delivered to the tumour centre. [3] Consequently, there is a need for the evaluation of new PSs in large tumour models that can better simulate the conditions found in the tumours of patients.

Multicellular tumour spheroids (MCTSs) can mimic the pathological conditions found in clinically-relevant tumours. MCTSs with a diameter of 200 μm are frequently used to simulate 3D cell-cell as well as cell-matrix interactions. They can therefore be utilised to model the drug delivery of a compound. [4] This is important since many investigated drug candidates have failed the translation from promising results in a 2D monolayer model to a 3D or an in vivo model due to compromised drug delivery. [5] Furthermore, MCTSs are also able to simulate hypoxia and proliferation gradients to the centre. Since the vast majority of PDT agents act by an oxygen-dependent mechanism, the treatment of hypoxic tumors remains a major medical challenge. [6] Despite these limitations, the successful treatment of the hypoxic regions of melanoma MCTSs with a diameter of 280 μm was recently reported. [7] Studies have shown that the oxygen diffusion inside of MCTSs is limited to approximately 200 μm. [8] Therefore, larger MCTS are able to generate a hypoxic core, which can better simulate the pathological conditions found in clinically-relevant tumours. Despite these benefits, studies of the antitumor effects of compounds in large MCTSs remain very rare. [9] To overcome this limitation and promote a better understanding of the biological effects of new compounds, there is an urgent need for in-depth studies of the ability of new compounds in large MCTSs. Ru(II) polypyridine complexes are gaining increasing attention as PDT PSs due to their attractive chemical and photophysical properties (e.g., high water solubility, chemical stability and photostability). [10] Despite their remarkable characteristics, the majority of PSs are activated using either UV or blue light, limiting the light penetration inside the tissue and therefore their application to treat deep seated or large tumours. To overcome this limitation, there is a need for the development of Ru(II) polypyridine complexes with an absorption in the biological spectral window (600-900 nm). [11] This aim could be achieved by a red-shift of the 1-Photon (1P) absorption or the use of a 2-Photon (2P) process, in which the compound absorbs two photons of low energy/high wavelengths simultaneously. Next to a deeper tissue penetration, a 2P excitation correlates with a reduced photodamage as well as enhanced spatial resolution. However, the Ru(II) polypyridine complexes reported so far were found to have relatively poor 2P photophysical properties [7, 9a, 12] compared to porphyrins, porphyrin oligomers or expanded porphyrinoids, [13] limiting the application of this technique with these complexes. In this context, we have recently reported the rational design of the complexes 1) as effective PSs for 1P (540 nm) and 2P (800 nm) PDT. [14] Worthy of note, the use of (E,E')-4,4´-bisstyryl-2,2´-bipyridine type Ru(II) polypyridine complexes were previously reported as effective 2P absorbing chromophores. [15] The lead compound of our study, namely [Ru(bipy)2(L-OMe)] 2+ was found to be active in vivo. [14] Despite the impressive biological property of this complex, its lack of 1P absorption in the biological spectral window and poor 2P absorption is limiting the application for deep-seated or large tumors. With this in mind, in this work, the analogous complexes that bear different ancillary ligands [Ru(phen/bphen)2(L-NMe2/L-OMe)] 2+ (1-4, phen = 1,10-phenenthroline, bphen = 4,7diphenyl-1,10-phenanthroline) are reported as effective PSs. The lead compound of this new study was found to have a red-shifted 1P and strong 2P absorption, enabling a phototherapeutic effect upon 1P (595 nm) or 2P (800 nm) excitation. An in-depth analysis of their ability to treat the hypoxic core in extraordinary large MCTS is presented. Further this complex was found to eradicate an adenocarcinomic human alveolar basal epithelial tumour inside a mouse model.

[Ru(bipy)n(L-NMe2/L-OMe)3-n] 2+ (bipy = 2,2´-bipyridine, L-NMe2 = (E,E')-4,4´-bis[p-(N,N- dimethylamino)styryl]-2,2´-bipyridine, L-OMe = (E,E')-4,4´-bis[p-methoxystyryl]-2,2´- bipyridine, n = 0, 1, 2, Figure
R = NMe2 OMe R = NMe2 OMe R = NMe2 OMe N N N N N N Ru R R R 2+ R R R N N N N N N Ru R R 2+ R R N N N N N N Ru R 2+ R R = NMe2 OMe 1 2 N N N N N N Ru R 2+ R R = NMe2 OMe 3 4 N N N N N N Ru R 2+ R
Previously studied compounds:

New compounds:

RESULTS AND DISCUSSION

To the best of our knowledge, the synthesis of the complexes 1-4 (Figure 1) has not been reported yet. The ligands L-NMe2 and L-OMe were synthesized as previously reported by our group. [14] The precursors Ru(phen)2Cl2 and Ru(bphen)2Cl2 were prepared as described in the literature. [16] Finally, the chloride substituents of the respective precursors were substituted with L-NMe2 or L-OMe, yielding the desired complexes 1a-4a. All compounds were characterized using 1 H-, 13 C-NMR, HR-MS analysis and the purity of the compounds verified by HPLC and elemental analysis. Details on the synthesis and characterization of the complexes can be found in the SI (Scheme S1-S2, Figure S1-S12).

To evaluate the potential of the compounds, their photophysical properties (Table S1) were determined. All complexes were found to have a strong red-shifted 1P absorption as indicated by their high extinction coefficients and absorption tail towards longer wavelengths (Figure S13a). Strikingly, the complexes have an exceptionally high 2P absorption up to ~1600 GM (Figure S13b), which is an order of magnitude higher than the majority of previously reported Ru(II) polypyridine complexes (~40-250 GM). [7, 9a, 12] Importantly, the compounds have a 1P

absorption tail towards and 2P absorption within the biological spectral window (600-900 nm), potentially allowing them to be utilized for the treatment of deep-seated or large tumours. The maximum of the emission of the complexes (Figure S14) was determined to be at 698 nm for the L-NMe2-coordinated complexes (1a, 3a) and at 663 nm for the L-OMecoordinated complexes (2a, 4a), resulting in a large Stokes shift for all investigated compounds. The comparison of the excitation and absorption spectra of all compounds showed no significant differences. The L-NMe2 or L-OMe ligands play a significant role in the absorption and emission characteristics of the complexes since complexes 1 and 3 as well as 2a and 4a have similar absorption and emission profiles. The analysis of the excited states of the analogous bipy derivatives suggested that these transitions stem from metal-to-ligand charge transfer/ligand-to-metal charge transfer excitations. [14] The comparison between the complexes revealed that the L-OMe-coordinated analogous (2a: 1.9 %, 4a: 2.7 %) have a significantly higher luminescence quantum yield than the L-NMe2-coordinated compounds (1a: 0.4 %, 3a: 0.6 %). Worthy of note, other previously published Ru(II) polypyridine complexes with terminal dimethylamine groups were also found to be poorly luminescent. [11c, 23d] Theoretical calculations for the [Ru(bphen)2((E,E')-4,4'-bis(N,N'-dimethylaminovinyl)-2,2'-bipyridine)] 2+ complex indicated that the energetically lowest lying excitations are based on non-luminescent ligand-centred transitions. [23d] The bphen coordinated compounds (3a, 4a) showed a stronger luminescence than the phen coordinated compounds (1a, 2a). This observation is in agreement with the comparison between the corresponding parent complexes

[Ru(phen)3] 2+ and [Ru(bphen)3] 2+ . [17] All complexes were found to have excited state lifetimes in the nanosecond range (Figure S15-S18) in an aerated environment (i.e., 72-193 ns) and in a degassed environment (i.e., 334-981 ns). These values are in the same range than for other reported Ru(II) polypyridine complexes. [18] As the lifetimes drastically decrease in the presence of air, it indicates that the excited state can interact with a component of air. Capitalizing on this, the type of ROS generated upon light irradiation was then investigated by electron spin resonance spectroscopy using 2,2,6,6tetramethylpiperidine as a singlet oxygen ( 1 O2) scavenger and 5,5-dimethyl-1-pyrroline Noxide as a • OOH or • OH radical scavenger. In contrast to the characteristic 1 O2-induced triplet signal of 2,2,6,6-tetramethylpiperidinyloxyl in CH3CN and PBS for all complexes (Figure S19-S22), no significant signal was observed for the formation of • OOH or • OH radicals. The 1 O2 quantum yields upon light exposure were determined using two complementary methods:

1) directly by measuring the phosphorescence signal of 1 O2 upon excitation at 450 nm; 2)

indirectly by capturing 1 O2 with a reporter molecule and monitoring its change by absorbance spectroscopy upon excitation at 450 or 540 nm. The compounds (Table S2) were found to have 1 O2 quantum yields in CH3CN between 43-92 % and 2-14 % in an aqueous solution. The comparison between the compounds revealed that the 1 O2 production of 2a and 4a is drastically higher than of 1a and 3a, which is expected due to their superior photophysical properties. The 1 O2 generation of 2a and 4a is superior than their structurally related parent complexes[Ru(phen)3] 2+ and [Ru(bphen)3] 2+ [17] as well as their analogous bipy derivatives. [14] The stability of a compound is a crucial factor in PDT applications. Previous studied have shown some metal complexes are not stable under physiological conditions, leading to undesired side effects. [19] The compounds were incubated in human plasma at 37 °C and extracted after 48 h. No significant differences between the chromatograms (Figure S23-S26)

were observed, indicating that all complexes were stable in a biological environment.

Following this, the stability upon light exposure at 540 nm was investigated by monitoring their absorption spectra. Previous studies have shown that this could be an issue for some Ru(II) complexes. [20] It is important to note that clinically used PSs, represented here by protoporphyrin IX (PpIX), are generally associated with a fast photodegradation. [21] The compounds investigated in this study were found to be highly photostable (Figure S27-S30)

whereas PpIX completely decomposed (Figure S31) under identical experimental conditions.

To study the effect that these complexes have on cells, we first assessed the time they need to be taken up by human cervical carcinoma (HeLa) cells using inductively coupled plasma mass spectrometry (ICP-MS). Within 6 h, the asymptotic maximum concentration (Figure S32-S35) was reached. The comparison between the complexes (Figure S36) shows that the bphen-coordinated complexes (3a, 4a) are slightly better taken up by cells than the phen coordinated complexes (1a, 2a). This result is in agreement with the determination of the logP values of the complexes (Table S3), which demonstrate their strong lipophilic character.

Following this, the uptake mechanism of the compounds was investigated by blocking different pathways (Figure S37-S40) by preincubation with metabolic (2-deoxy-D-glucose and oligomycin), cationic transporter (tetraethylammonium chloride) and endocytotic (ammonium chloride or chloroquine) inhibitors as well as incubation at lower temperature (4 °C). [22] As the incubation with tetraethylammonium chloride had only a negligible effect on the uptake, the internalization by this pathway was ruled out. The change to lower temperatures as well as the incubation with metabolic inhibitors significantly decreased the uptake, indicating that the mechanism is energy dependent. As the preincubation with ammonium chloride or chloroquine strongly decreased the internalization of the complexes 1a-4a, it indicates that all compounds are taken up through an energy dependent endocytosis pathway. The subcellular localization was then examined using confocal laser scanning microscopy (Figure S41). The compounds were incubated in HeLa cells with commercial dyes for the major cellular organelles (i.e., nucleus, mitochondria, lysosomes, Golgi apparatus, endoplasmic reticulum) and their colocalization in the cells compared. As no significant overlap was observed, it suggests that the compounds do not majorly localize in these organelles. To further investigate the localization, the cell organelles were separately extracted and the amount of Ru determined by ICP-MS (Figure S42). All compounds were mainly found in the cytoplasm with small amounts of unselective accumulation. The structurally related bipy derivatives previously studied were also found to majorly localize in this organelle. [14] To study the (photo-)cytotoxic effects of 1a-4a, these compounds were incubated in the dark as well as upon light exposure at 480 nm (10 min, 3.1 J/cm 2 ) and 540 nm (40 min, 9.5 J/cm 2 ) in non-cancerous retinal pigment epithelium (RPE-1), HeLa, mouse colon carcinoma (CT-26)

and human glioblastoma astrocytoma (U373) cells. PpIX, the anticancer drug cisplatin and the parent complexes [Ru(phen)3] 2+ and [Ru(bphen)3] 2+ were used as controls. The light doses and irradiation times used in this work were optimized to the survival of the cells when treated only with the light source. All cell lines were tested identically to investigate the influence of the compounds on different types of cancer. The results (Table S4-S5) show that the phen-based compounds (1a, 2a) as well as [Ru(phen)3] 2+ are non-toxic in the dark (IC50 > 100 mM), which is a crucial requirement for PDT applications. On the contrary, the bphencoordinated complexes (3a, 4a) were found to be cytotoxic in the low micromolar range (IC50 = 5.2 -20.8 μM) in all tested cell lines. The observed cytotoxicity in the dark for the bphencoordinated Ru(II) complexes could also be observed for [Ru(bphen)3] 2+ and is in agreement with recent studies of structurally related bphen-coordinated compounds. [23] Upon light irradiation, all compounds were able to generate 1 O2, causing cell death in all investigated cell lines. The phen-coordinated compounds (1a, 2a) were found to have a phototoxic effect in the low micromolar range (IC50 = 0.9 -15.6 μM), while [Ru(phen)3] 2+ only had a negligible effect. The bphen-coordinated complexes (3a, 4a) as well as [Ru(bphen)3] 2+ have, next to a cytotoxic effect in the dark, a phototoxic effect, as demonstrated by their IC50 values in the high nanomolar to low micromolar range (IC50 = 0.5 -2.1 μM). The L-OMe-based compounds (2a, 4a) were more phototoxic than the L-NMe2-based compounds (1a, 3a). This observation is attributed to their superior photophysical properties (Table S1-S2), notably their better 1 O2 production. Overall, 2a was identified as the best compound of the series with no observed toxicity in the dark and a phototoxicity in the high nanomolar range in CT-26 cells (IC50,540nm = 0.9 ± 0.4 μM, PI540nm > 111). These highly promising results of 2a are an order of magnitude lower compared to the IC50 values of the PS PpIX, the parent complex [Ru(phen)3] 2+ and cisplatin tested under identical experimental conditions.

The phototoxic profile of the lead compound 2a was further explored using longer wavelengths. 2a was also active upon irradiation at 595 nm (60 min, 11.2 J/cm 2 ) in CT-26 cells with an IC50 value in the low micromolar range (IC50,595nm = 2.4 ± 0.3 μM, PI595nm > 41.7). In contrast, no phototoxicity was observed upon exposure to light at 620 nm (30 min, 3.3 J/cm 2 ). The cell death mechanism caused upon irradiation at 480 nm (10 min, 3.1 J/cm 2 ) was investigated by determining the cell viability upon preincubation with autophagy (3methyladenine), apoptosis (Z-VAD-FMK), paraptosis (cycloheximide) and necrosis (necrostatin-1) inhibitors (Figure S43). Since preincubation with autophagy and necrosis inhibitors did not significantly influence the cell survival, these pathways were ruled out for all compounds. The preincubation with a paraptosis inhibitor slightly increased cell survival, while the preincubation with an apoptosis inhibitor highly increased cell survival. This indicates that the cell death is mainly caused by the apoptosis pathway, with minor contribution from the paraptosis pathway. Some previously studied PS were triggering cell death by the same mechanism. [24] To further investigate the apoptosis mechanism, its dependency on caspases 3/7 was studied. Caspases 3/7 are well known executers of the intrinsic and extrinsic apoptosis mechanism. [25] The caspase activity was measured in HeLa cells upon irradiation at 480 nm (10 min, 3.1 J/cm 2 ). The results show highly increased caspase levels upon irradiation (Figure S44), indicating that the phototoxic effect of all compounds is caused by a caspase 3/7 pathway. Some previously studied Ru(II) polypyridine complexes were found to exert their phototoxic effect by the same mechanism. [26] Following this, the capacity of the compounds to generate 1 O2 upon irradiation inside the hypoxic centre of the MCTSs was investigated. For this purpose, HeLa MCTSs were incubated with 2 ´,7 ´-dichlorofluorescein diacetate (DCFH-DA), which is converted into the highly fluorescent 2 ´,7 ´-dichlorofluorescein in the presence of ROS. While no green fluorescence signal was observed in the dark, upon 2P irradiation (800 nm, 2 J/cm 2 , section interval of 5 μm), a strong signal in the whole MCTSs was detected (Figure 3). Strikingly, the signal was also observed in the centre of the MCTSs, suggesting the generation of 1 O2 despite the hypoxic conditions, which are typically observed in MCTS. To study the (photo-)cytotoxic effect in MCTSs, 1a-4a (20 μM), cisplatin (10 μM and 30 μM)

and the well-characterized PS tetraphenylporphyrin (H2TPP, 20 μM) were exposed to 1P (500 nm, 10 J/cm 2 ) and 2P (800 nm, 10 J/cm 2 , section interval of 5 μm) irradiation and their tumour growth monitored (Figure 4). As observed in 2D monolayer cells, 3a and 4a showed a weak tumour growth inhibition effect in the dark similar to cisplatin at 10 μM while 1a and 2a

did not significantly influence the tumour growth in the dark -the MTCSs grew in similar manner as the control. Upon 1P or 2P light treatment, all compounds 1a-4a were able to cause a phototoxic effect inside the MCTSs leading to a decrease of the volume of the MCTSs. In particular, the compounds 2a and 4a drastically decreased the tumour growth upon light exposure, demonstrating their ability to act efficiently as a PS inside a 3D cellular architecture. Under identical conditions, the MCTSs treated with H2TPP did not show any significant effects, highlighting the superior ability of these compounds. After this time, MCTSs were exposed to a 1P (500 nm, 10 J/cm 2 ) or a 2P irradiation (800 nm, 10 J/cm 2 , section interval of 5 μm). After 2 days, the cell survival was assessed by measurement of the fluorescence of calcein (λex = 495 nm, λem = 515 nm) and cell death by measurement of the fluorescence of propidium iodide (λex = 536 nm, λem = 617 nm).

To further study and quantify the (photo-)cytotoxic effect of the complexes, their IC50 values in MCTSs in the dark as well as upon 1P (500 nm, 10 J/cm 2 ) and 2P (800 nm, 10 J/cm 2 , section interval of 5 μm) excitation were determined by measuring their ATP concentration.

The obtained results (Table 1) are in agreement with the tumour growth inhibition assay,

showing that the bphen-coordinated complexes (3a, 4a) have a cytotoxic effect in the low micromolar range in the dark while the phen-coordinated complexes (1a, 2a) showed no dark toxicity even up to high micromolar range (IC50 > 300 μM). Upon exposure to 1P or 2P irradiation, all compounds were able to cause a phototoxic effect in the micromolar range (IC50,1P = 3.8 -32.6 μM, IC50,2P = 0.8 -27.8 μM). As the lead compound of this study, 2a demonstrated its remarkable ability as a PS with PI values of >40 upon 1P and >250 upon 2P irradiation. Further, 2a was also found to be phototoxic in HeLa MCTSs upon irradiation at 595 nm (60 min, 11.2 J/cm 2 ) in the low micromolar range (IC50,595nm = 16.7 ± 1.2 μM, PI595nm > 18). Overall, the lead compound 2 was found to have a two order of magnitudes higher phototoxicity in comparison to H2TPP, indicating that the compounds reported here can act using very low drug and light doses.

Table 1. IC50 values in the dark as well as upon 1P (500 nm, 10 J/cm 2 ) and 2P (800 nm, 10 J/cm 2 , section interval of 5 μm) irradiation of 1a-4a as well as of cisplatin and H2TPP in HeLa 3D MCTS. Average of three independent measurements. n.d. = not determinable.

For a further investigation, the biological properties of 2 inside a mouse model was investigated. As the counter ion of a metal complex can have a significant effect on the overall ability of the compound including its water solubility or toxicity [27] , 2a was converted into the chloride salt 2b using a counter ion exchange resin. The biodistribution of the compound inside of nude mice bearing an adenocarcinomic human alveolar basal epithelial (A549) tumour upon intravenous tail-injection was time-dependently (0.5, 1, 2, 4 h) studied. At each time point, the mice were sacrificed, all major organs collected (i.e., blood, spleen, intestine, gastric, liver, kidney, lung, heart, brain, tumour), ground and the amount of Ru inside each organ determined by ICP-MS (Figure S51). The compound was rapidly absorbed from the blood stream, accumulating mainly within the liver and kidney. The compound also accumulated inside the tumour with a concentration maximum 2 h after the injection.

Capitalizing on these results, in vivo PDT experiments were performed using a 1P (500 nm, 10 mW/cm with a 2P light source did not show a significant growth after the treatment, the tumours treated with a 1P light source did indicate a small tumour growth. We assume that this could be caused by the deeper light tissue penetration of 2P irradiation compared to 1P. In contrast, the tumours only treated by light or with the compound did not show any tumour inhibition effect (Figures 6a and6c). The animals treated with the compound behaved normally, without signs of pain, stress or discomfort and did not lose or gain weight (Figure 6b). After 17 days (representative picture: Figure S52), all mice were sacrificed and the tumour and organs separated. As excepted, using a histological examination with an H&E stain the tumour tissue showed noticeable pathological alterations (Figure S53) while all major organs did not show any pathological effect (Figure S54). Strikingly, after the treatment, the MCTSs were fully eradicated, including in their large hypoxic centres. While the lead compound of this study, 2 was found to be non-toxic in the dark, it is highly phototoxic in MCTSs in the very low micromolar range upon 1-Photon (595 nm) or 2-Photon (800 nm) irradiation. Worthy of note, in comparison to the previously published 2,2´-bipyridine analogous compounds, these complexes were found to be photoactive upon even longer wavelengths by 1-Photon irradiation. In addition, 2 was also able to eradicate an adenocarcinomic human alveolar basal epithelial tumour inside a mouse model upon clinically relevant 1-Photon (500 nm) or 2-Photon (800 nm) excitation, demonstrating its high potential as a PS. Overall, the compounds presented in this study can overcome the limitations of currently applied PSs. We strongly believe that these complexes, especially 2, has great potential for 1-and 2-Photon excited PDT and is a suitable candidate for further clinical investigations.

Figure 1 .

 1 Figure 1. Chemical structures of previously studied and new compounds investigated in this work. The followed letter a corresponds to the compounds as a hexafluorophosphate salt whereas the letter b indicated the compound as a chloride salt.

Figure 2 .

 2 Figure 2. 1P (λex = 458 nm, λem = 600 -750 nm) and 2P (λex = 800 nm, λem = 600 -750 nm) excited Z-stack images in HeLa MCTS after incubation of 2a (10 μM) for 12 h. a/c) Z-axis

Figure 3 .

 3 Figure 3. Confocal fluorescent images of HeLa MCTSs incubated with DCFH-DA and the compounds 1a-4a (2 μM) kept in the dark and after a 2P irradiation (800 nm, 2 J/cm 2 , section interval of 5 μm).

Figure 4 .

 4 Figure 4. Tumour growth inhibition assay. Change of the volume in HeLa MCTSs in correlation to treatment time. The MCTSs were treated with the compounds 1a-4a (20 μM), H2TPP (20 μM) and cisplatin (10 μM and 30 μM). The MCTSs were a) strictly kept in the dark b) exposed to a 1P irradiation (500 nm, 10 J/cm 2 ) c) exposed to a 2P irradiation (800 nm, 10 J/cm 2 , section interval of 5 μm) on day 1. Representative pictures of the MCTSs can be found in Figure S48-S50. The error bars correspond to the standard deviation of the three replicates.

Figure 5 .

 5 Figure 5. Representative image of the viability assay with HeLa MCTSs kept in the dark and exposed to light. MCTS were treated with the compounds 1-4 (20 μM) in the dark for 12 h.

Figure 6 .

 6 Figure 6. In vivo PDT study of 2b using 1P (500 nm, 10 mW/cm 2 , 15 min) or 2P (800 nm, 50 mW, 1 kHz, pulse width 35 fs, 5 s/mm) excitation on nude mice bearing an adenocarcinomic human alveolar basal epithelial cancer tumour. a) Tumour growth inhibition curves upon treatment. b) Average body weights of the tumour-bearing mice. c) Representative photographs of the tumour after different treatments on day 17. *p<0.05, **p<0.01.

  Encouragingly, the PDT treated tumours drastically shrank until they were nearly eradicated. It is important to mention that while the tumours treated

	light treatment on days 1, 4 and 7.
		dark	1P	PI	2P	PI
	1	>300	32.6 ± 2.5	>9.2	27.8 ± 3.1	>10.8
	2	>300	7.5 ± 0.2	>40.0	1.2 ± 0.3	>250.0
	3	27.8 ± 3.6	8.9 ± 0.7	3.1	3.1 ± 0.6	9.0
	4	29.3 ± 2.9	3.8 ± 0.4	7.7	0.8 ± 0.5	36.6
	H2TPP	>100	>100	n.d.	>100	n.d.
	cisplatin	18.6 ± 1.3	-	-	-	-

2

, 15 min) or 2P irradiation (800 nm, 50 mW, 1 kHz, pulse width 35 fs, 5 s/mm)

ACKNOWLEDGEMENTS

We thank Dr. Philippe Goldner for access to state-of-the-art laser apparatus. This work was financially supported by an ERC Consolidator Grant PhotoMedMet to G.G. (GA 681679), has received support under the program "Investissements d' Avenir" launched by the French Government and implemented by the ANR with the reference ANR-10-IDEX-0001-02 PSL (G.G.), the National Science Foundation of China (Nos. 21525105 and 21778079 for H.C.) and the 973 Program (No. 2015CB856301 for H.C.).


Aurélien Appriou
email: aurelien.appriou@inria.fr

Towards measuring states of epistemic curiosity through electroencephalographic signals 1 st

Keywords:

Understanding the neurophysiological mechanisms underlying curiosity and therefore being able to identify the curiosity level of a person, would provide useful information for researchers and designers in numerous fields such as neuroscience, psychology, and computer science. A first step to uncovering the neural correlates of curiosity is to collect neurophysiological signals during states of curiosity, in order to develop signal processing and machine learning (ML) tools to recognize the curious states from the non-curious ones. Thus, we ran an experiment in which we used electroencephalography (EEG) to measure the brain activity of participants as they were induced into states of curiosity, using trivia question and answer chains. We used two ML algorithms, i.e. Filter Bank Common Spatial Pattern (FBCSP) coupled with a Linear Discriminant Algorithm (LDA), as well as a Filter Bank Tangent Space Classifier (FBTSC), to classify the curious EEG signals from the non-curious ones. Global results indicate that both algorithms obtained better performances in the 3-to-5s time windows, suggesting an optimal time window length of 4 seconds (63.09% classification accuracy for the FBTSC, 60.93% classification accuracy for the FBCSP+LDA) to go towards curiosity states estimation based on EEG signals.

I. INTRODUCTION

Curiosity is an important mental state associated with spontaneous exploration, active learning, facilitated memorization and sustained engagement [START_REF] Oudeyer | Intrinsic motivation, curiosity, and learning: Theory and applications in educational technologies[END_REF]. Recent research in psychology [START_REF] Kidd | The psychology and neuroscience of curiosity[END_REF] and neuroscience [START_REF] Gottlieb | Towards a neuroscience of active sampling and curiosity[END_REF] has shown its pervasive role across multiple dimensions of human cognition and learning. The design of robotic and computer architectures can also be guided by our understanding and conceptualizations of curiosity [START_REF] Oudeyer | Intrinsic motivation systems for autonomous mental development[END_REF], [START_REF] Gordon | Can children catch curiosity from a social robot?[END_REF]. Human-computer interaction (HCI) researchers have investigated how systems can elicit curiosity in order to inform the design of persuasive, engaging, and playful interactions [START_REF] Tieben | Curiosity and interaction: making people curious through interactive systems[END_REF], or to motivate and incentivize crowd workers [START_REF] Law | Curiosity killed the cat, but makes crowdwork better[END_REF]. Finally, as curiosity-driven learning has been argued to be a crucial feature for efficient education [START_REF] Oudeyer | Intrinsic motivation, curiosity, and learning: Theory and applications in educational technologies[END_REF], [START_REF] Freeman | Active learning increases student performance in science, engineering, and mathematics[END_REF], learning technologies are being developed to promote curiosity and motivate curiosity-driven behaviours in students [START_REF] Ceha | Expression of curiosity in social robots: Design, perception, and effects on behaviour[END_REF].

Brain-Computer Interfaces (BCIs) enable users to interact with computers by using brain activity, as measured by elec-troencephalography (EEG) signals [START_REF] Clerc | Brain-Computer Interfaces 1: foundations and methods[END_REF]. There are two types of BCIs: active BCIs, where brain activity is converted into a command to an application, and passive BCIs, where the goal is to monitor the user's mental state in real time and adapt the interaction accordingly [START_REF] Zander | Towards passive brain-computer interfaces: applying brain-computer interface technology to human-machine systems in general[END_REF].

The premise of this paper is that passive BCIs that can monitor, through neurophysiological signals, the level of curiosity of users can be both a useful tool for understanding curiosity, and beneficial to designers of interactive systems, who wish to adapt the interaction paradigm or application content to users' levels of curiosity. In the context of active BCIs, which are known to be notoriously difficult for novices to use [START_REF] Jeunet | Towards a cognitive model of MI-BCI user training[END_REF], it would be beneficial to adapt the BCI training tasks to the mental states of these users, e.g., their curiosity level, and adapt the BCI tasks to users' interests in order to prevent boredom and improve learning. In this work, we conducted an experiment in which subjects were given a trivia questionbased task, designed to elicit different levels of curiosity. We used signal processing and machine learning tools to analyze the EEG signals collected through this protocol, and assessed how well we can estimate curiosity through EEG signals using Machine Learning (ML) classifiers. Results show that F1 score can reach 63% from only a few seconds of EEG, depending on the pre-processing and ML algorithms chosen.

The first to propose a definition of curiosity was William James, who said that "curiosity is an instinct that evolved to facilitate survival and adaptation through active exploration of the environment" [START_REF] Gruber | 16 Curiosity and Learning[END_REF]. Daniel Berlyne then introduced a multi-dimensional model of curiosity, characterizing this psychological state in two dimensions, i.e., perceptual/epistemic curiosity, and specific/diversive curiosity [START_REF] Berlyne | A theory of human curiosity[END_REF]. Following this theory, perceptual curiosity refers to "a drive which is aroused by novel stimuli and reduced by continued exposure to these stimuli", whereas epistemic curiosity refers to "a desire to acquire knowledge, and applies mainly to humans" [START_REF] Loewenstein | The psychology of curiosity: A review and reinterpretation[END_REF]. On the second dimension, specific curiosity is defined by the desire for a particular piece of information, whereas diversive curiosity is defined as a more general seeking for stimulation that is related to boredom. Although no scientific consensus has been reached concerning the definition of curiosity, certain types of stimuli are known to trigger it, i.e., those of surprising, novel, or intermediate complexity, as well as activities that are characterized by a knowledge gap or errors in prediction [START_REF] Oudeyer | Intrinsic motivation, curiosity, and learning: Theory and applications in educational technologies[END_REF]. Such triggers can lead to different momentary states of curiosity, including epistemic [START_REF] Mussel | Epistemic curiosity and related constructs: Lacking evidence of discriminant validity[END_REF], [START_REF] Litman | Encyclopedia of the Sciences of Learning[END_REF] and perceptual curiosities [START_REF] Gruber | 16 Curiosity and Learning[END_REF]. Epistemic curiosity has been particularly studied through psychological experiments in the last decades [START_REF] Brod | Lighting the wick in the candle of learning: generating a prediction stimulates curiosity[END_REF], and has more recently been the object of neuroscientific experiments.

Several studies have been run to better understand the neural mechanisms underlying states of curiosity. In [START_REF] Kang | The wick in the candle of learning: Epistemic curiosity activates reward circuitry and enhances memory[END_REF], the authors scanned participants with functional magnetic resonance imaging (fMRI) in order to study the brain areas activated during the triggering of curiosity states with trivia questions. They observed that the curiosity induced by the trivia questions was correlated to brain activity in the caudate region, an area shown to be associated with anticipated rewards [START_REF] Adcock | Reward-Motivated Learning: Mesolimbic Activation Precedes Memory Formation[END_REF]. Moreover, they also found a correlation between surprising new information and activation of brain areas linked to memory. In [START_REF] Gruber | States of Curiosity Modulate Hippocampus-Dependent Learning via the Dopaminergic Circuit[END_REF], the question was raised as to whether curiosity enhances longterm memory, similar to the way anticipated rewards do. They likewise conducted a fMRI study with trivia question tasks, and found a correlation between curiosity levels and variations of activation in the right hippocampus and bilateral nucleus accubens, both involved in long-term memory improvements [START_REF] Gruber | 16 Curiosity and Learning[END_REF]. Using frontal EEG asymmetry-a common tool for measuring engagement and motivation- [START_REF] Lima | Curiosity, frontal EEG asymmetry , and learning[END_REF] investigated the relation between curiosity and learning. Participants performed trivia question tasks, similar to those used in [START_REF] Kang | The wick in the candle of learning: Epistemic curiosity activates reward circuitry and enhances memory[END_REF], [START_REF] Gruber | States of Curiosity Modulate Hippocampus-Dependent Learning via the Dopaminergic Circuit[END_REF], while the EEG signals from the frontal cortex were recorded. Here, researchers found a correlation between frontal brain asymmetry (FBA) and memory recall; however they did not observe any correlation between FBA and self-reported curiosity.

In summary, these neuroscientific studies, mainly based on fMRI, support the existence of a correlation between epistemic curiosity and memory/learning, as well as correlations between states of curiosity and activation in specific brain regions. However, these prior works did not perform continuous state monitoring. In this work, we propose to measure states of epistemic curiosity using EEG. Compared to fMRI, which can be expensive and difficult to use outside of the lab, EEG provides a usable, portable and affordable tool for measuring the temporal activation of brain states associated with curiosity, making it suitable for applications such as BCIs.

II. MATERIALS & METHODS

We describe in detail the experimental setup, including the participants involved, the protocol, the signal pre-processing procedure, the EEG-based curiosity classifiers used, and the evaluations performed.

A. Participants

Twenty-seven participants (N=27) were recruited through ads posted on social media and the local university mailing lists (13F/14M; aged 28.7 ± 4.0). Levels of education varied between high school diploma and Ph.D. To be included in the experiment, people had to be at least 18 years old, speak French, and consent to the study1 . Non-inclusion criteria include bad vision, heart condition, neurological or psychological diseases, and emotion-related problems. Finally, as per typical EEG studies protocol, we asked participants not to drink coffee or tea within the 2 hours prior to the experiment.

B. Protocol

Each participant participated in a single session that took place at Inria Bordeaux Sud-Ouest, lasting approximately 2 hours. Participants were asked to fill-in a pre-session questionnaire, assessing personal characteristics (such as gender, age and education). Next, EEG and ElectroOculoGram (EOG) electrodes were placed on the participants. For recordings, we used a BrainProduct ActiCHamp amplifier (EEG, 61 active electrodes in a 10/20 system, and EOG, 3 active electrodes to measure ocular artifacts). To our knowledge, except for [START_REF] Lima | Curiosity, frontal EEG asymmetry , and learning[END_REF] where they attempted -but did not succeed -to measure curiosity with electrodes placed onto the scalp recovering the frontal cortex, there has been no study seeking to classify curiosity levels from EEG signals. Due to the exploratory nature of our study, we covered the entire scalp with a relatively high EEG spatial resolution (61 electrodes). EEG signals were recorded and visually inspected using OpenViBE [START_REF] Renard | Open-ViBE: An open-source software platform to design, test, and use braincomputer interfaces in real and virtual environments[END_REF]. Sensors for ElectroDermal Activity (EDA), breathing and heart rate (HR) were also installed; though these signals were not used in the present study. Following the setup, a 3 minute baseline, consisting of measuring EEG signals from participants at rest with opened eyes, was recorded. Finally, participants were asked to perform 4 runs of curiosity tasks (described below), around 10 minutes each, with 5 minute breaks between them. Each run consisted of a series of trivia questions and answers, inducing different levels of curiosity. Before each question, a fixation cross was presented on the screen for 3 seconds, in order to get the participant ready. Figure 1 illustrates the experiment flow.

C. Materials

Prior work using trivia questions to elicit curiosity presented participants with question/answer pairs that did not have any link to previously viewed or future questions [START_REF] Kang | The wick in the candle of learning: Epistemic curiosity activates reward circuitry and enhances memory[END_REF], [START_REF] Gruber | States of Curiosity Modulate Hippocampus-Dependent Learning via the Dopaminergic Circuit[END_REF], [START_REF] Lima | Curiosity, frontal EEG asymmetry , and learning[END_REF]. In our study, we introduce a novel protocol consisting of chains of trivia question/answer pairs-i.e., if participants were curious about the answer to a certain question, the following question would follow on the same topic. The assumption is that if participants were curious about a certain question, it was likely they would be curious about this topic in general. This new method allowed us to: 1) record a large enough amount of EEG signals from both curious and non-curious states, to then be able to train the curiosity classifiers with a balanced set of EEG examples, and 2) check the assumption that curiosity could be a mental state that increases over time, in the same way intrinsic motivation and self-directed learning increase when time spent in flow state increases [START_REF] Hektner | A Longitudinal Exploration of Flow and Intrinsic Motivation in Adolescents[END_REF].

The trivia questions used to elicit curiosity in this study came from an online trivia question dataset (https://www. randomtriviagenerator.com). The questions were grouped following a two-level categorization system. Questions from the website were already classified into classical trivia categories (referred to as first-order categories), such as Science, History, Geography, Arts, General Knowledge and Sports. We further classified the questions in each first-order category into groups of 4 to 20 questions based on 800 extracted keywords from the first-order category questions. For example, in the Geography category, we identified a sub-category based on the keyword "Nile". All questions from the Geography category with the word "Nile" were then grouped into a single chain of questions. Figures 1 and2 show another example of a question chain based on the keyword "World War One".

All chains containing less than 4 or more than 20 questions were not used in this study. A minimum of 4 was chosen in order to have enough questions to define a subcategory. Conversely, we chose 20 as the maximum to ensure the subcategory was not too vague. Overall, the dataset consisted of 2000 questions/answers divided between 250 subcategories.

D. Curiosity Task

The curiosity task was presented on a computer screen using images -black background with white text. Each run was set up as follows: a first question from a random chain from a random category was displayed to the participant -for 7 to 9 seconds -right after a 3-second fixation cross. The participant then had 2-3 seconds to choose to display the answer or not, as presented in Figure 1. The display time was determined based on the length of the question or answer in terms of number of characters. The participant chooses to display the answer by tapping on the keyboard space bar: this question Fig. 2. The trivia questions/answers system. In the example, a question about World War One is presented: if they choose to display the answer, they will stay on the "World War One" topic, but will continue on another topic -here scientific questions about ants -if they skip the answer.

was then flagged as "curious", the concerned subcategory was considered interesting for the participant, and the answer was directly displayed on the screen for a few seconds. Participants were asked to rank their level of curiosity for the question on a 1-7 scale using a number pad right after the answer had been displayed. Following the rating, a new question randomly selected from the same subcategory was displayed.

If the participant did not choose to display the answer (by tapping on the keyboard) before the decision time elapsed, the question was labeled as "non-curious", the subcategory was not considered to be interesting to the participant. The curiosity rating scale was immediately administered without revealing the answer, followed by a fixation cross and a new random question from the next subcategory. A run ended only if at least 15 trials with questions marked as curious and 15 trials with questions marked as non-curious had been displayed. Thus, for each participant, we obtained at least 30 trials per run, i.e., a total of 120 trials in 4 runs: 60 trials per state of curiosity (curiosity & no curiosity).

E. EEG signal processing & classification

Our system aims to discriminate curious from non-curious states using EEG signals. To do so, we employed machine learning approaches based on state-of-the-art algorithms developed for BCIs [START_REF] Lotte | A Review of Classification Algorithms for EEG-based Brain-Computer Interfaces: A 10-year Update[END_REF]- [START_REF] Appriou | Modern machine learning algorithms to classify cognitive and affective states from electroencephalography signals[END_REF].

1) EEG Pre-processing: We first pre-processed EEG signals into N-second windows, in order to create 5 different data sets, with 1-, 2-, 3-, 4-and 5-second windows. The EEG signals for each 10 minute run were divided into approximately 30 trials, i.e., one trial per trivia question displayed. More precisely, an EEG trial was defined as ending at the time the question disappeared from the screen, and starting N-seconds earlier, as represented by the blue arrows on Figure 3. Note that no artifact removal algorithm has been used in this study. 2) Classification Labels: Each trial was flagged as "answer" or "no-answer", based on the participant's choice to reveal the answer or not, respectively. Based on the flags and participants' self-reported curiosity ratings, we labeled each trial as "curious" if the trial was flagged as "answer" and the rating was higher than the mean of the participant's ratings; otherwise as "non-curious", if the trial was flagged as "no-answer" and the rating was lower than the mean of the participant's ratings. We obtained around 50 trials per curiosity level per participant (mean number of trials in class "non-curious" = 50.18±12.58; mean number of trials in class "curious" = 46.85±12.11).

3) Machine Learning Algorithms: We used two ML algorithms that are exploring multiple frequency bands: (1) a Filter Bank Common Spatial pattern (FBCSP) coupled with a Linear Discriminant Analysis (LDA) [START_REF] Ang | Filter bank common spatial pattern algorithm on BCI competition IV datasets 2a and 2b[END_REF], and (2) a Filter Bank Tangent Space Classifier (FBTSC) [START_REF] Appriou | Modern machine learning algorithms to classify cognitive and affective states from electroencephalography signals[END_REF]. Both algorithms proved effective for mental state classification from EEG [START_REF] Appriou | Modern machine learning algorithms to classify cognitive and affective states from electroencephalography signals[END_REF].

Prior work studying curiosity through EEG signals [START_REF] Lima | Curiosity, frontal EEG asymmetry , and learning[END_REF], only used two electrodes that were placed on the frontal cortex, and the recorded signals were band-pass filtered in the alpha frequency band (8-13Hz). They did not find any correlation between activity in the frontal cortex and curiosity. It was therefore of interest to extract information from multiple electrodes (here 61) and frequency bands (here ten).

The FBCSP+LDA algorithm, which won the fifth international BCI competition [START_REF] Ang | Filter bank common spatial pattern algorithm on BCI competition IV datasets 2a and 2b[END_REF], works as follows: first, the training phase consists of optimally identifying and extracting both spatial and spectral features. For the spectral dimension, i.e., the frequency bands, EEG signals are filtered into ten 4Hzwide frequency bands (in 1-4 Hz, 4-8 Hz, ..., 36-40 Hz) as in [START_REF] Ang | Filter bank common spatial pattern algorithm on BCI competition IV datasets 2a and 2b[END_REF]. For each band, the band-pass filtered EEG trials are used. Then spatial filters are built for each band using the Common Spatial Pattern (CSP) algorithm [START_REF] Blankertz | Optimizing spatial filters for robust EEG single-trial analysis[END_REF], which optimizes the EEG signal-to-noise ratio: the variance of spatially filtered signals is maximized for one class and minimized for the other class. In our study, 4 CSP filters (2 pairs) have been optimized for each frequency band, resulting in 40 features (4 CSP filters * 10 frequency bands). From those 40 features, 4 were selected using the maximum Relevance Minimum Redundancy (mRMR) feature selection algorithm [START_REF] Peng | Feature Selection Based On Mutual Information: Criteria of Max-Dependency,Max-Relevance, and Min-Redundancy[END_REF] to train an LDA classifier to discriminate curious from non-curious trials.

The second algorithm, FBTSC, [START_REF] Appriou | Modern machine learning algorithms to classify cognitive and affective states from electroencephalography signals[END_REF] represents EEG signals as covariance matrices and manipulates them with Riemannian geometry [START_REF] Yger | Riemannian approaches in Brain-Computer Interfaces: a review[END_REF]. Here, each trial is first band-pass filtered in the same ten 4Hz-wide frequency bands we used with the FBCSP (1-4Hz, 4-8Hz, ..., 36-40Hz). To design the classifier, the average spatial covariance matrix for each class (curious and non-curious) is computed for each frequency band, and all covariance matrices are then projected in tangent space at the point defined as the matrices mean. We then used the softmax function-based probabilistic output of a Logistic Regression (LR) that has been trained in the tangent space, to determine probabilities of belonging to each class. Since we have two classes and a bank of ten frequency bands, 10 pairs of probabilities were computed. From these pairs of probabilities, the four most discriminant were selected using mRMR on the training set, and then multiplied together in order to obtain 2 final probabilities, i.e., one for the "curious" class and one for the "non-curious" class. The class assigned to a test EEG trial was decided according to the highest probability.

4) Evaluation:

We assessed the performances of both algorithms using a within-participant study with five-fold stratified Cross-Validation. This means the data from each participant was divided into five parts: four parts were used for training the classifier and the fifth one for testing the resulting curiosity classifier for that participant. This process was repeated five times, with each part used exactly once as the testing set.

III. RESULTS

The F1-score, which is the weighted average of the precision and the recall, for each classifier and each time window length are reported in Figure 5. As a reference, the statistical chance levels using [START_REF] Combrisson | Exceeding chance level by chance: The caveat of theoretical chance levels in brain signal classification and statistical assessment of decoding accuracy[END_REF] was estimated at 51.59% (100 trials per participants on average, 27 participants).

The boxplot of the performances obtained by each algorithm with the different time window lengths are reported in Figure 6. We performed a 2-way ANOVA with repeated measures to evaluate the performance of the factor Time Window according to the factor Algorithm (FBCSP+LDA vs FBTSC). Note that we checked the data sphericity, as well as the normality, and used Greenhouse-Geisser (GG) correction in ANOVA if needed. The ANOVA revealed a main effect of Time Figure 4 shows the F1-score, precision and recall for each participant and algorithm, with the 4-seconds time window length, i.e., the time window with which both classification algorithms obtained the best performances. Still using the 4seconds time window length, we also studied the percentage of time that each frequency band was selected (by mRMR) by each algorithm, as reported in Figure 7. Those results show that the ML algorithms mainly used 4 frequency bands-i.e., delta (1-4Hz), theta (4-8Hz), alpha (8-12Hz) and low beta (12-16Hz)-to classify states of epistemic curiosity.

IV. DISCUSSION, CONCLUSION AND FUTURE WORK

In this paper, we conducted an experiment aimed at collecting EEG signals during states of curiosity that were triggered using chains of trivia questions, and using ML to distinguish curious vs non-curious states. We used two ML algorithms, i.e., the FBCSP+LDA and FBTSC, to classify EEG signals at five different time-window lengths. The best results were obtained with the FBTSC, reaching about 62% of F1score for the 3, 4 and 5s time-windows (respectively 62.22, 63.09 and 62.90%), significantly outperforming the chance level (51.59%). Results in those 3 time-window lengths also significantly outperformed the results in the 2s time window (55.05%), but not the ones in the 1s time window (57.36%). The FBCSP+LDA reached an F1-score of 59% for the 3, 4 and 5s time windows, with respectively 59.90, 60.09 and 59.81%, significantly outperforming both the chance level (51.59%) and the results in the 1s time window (53.39%), but not the ones in the 2s time window (57.34%). Overall, results indicate that both algorithms obtained better performances in the 3-to-5s time windows, suggesting a minimum time of 3 seconds to go towards curiosity state estimation based on EEG signals. Moreover, ML algorithms mostly used a range of frequency bands from delta to low beta in order to classify states of epistemic curiosity, suggesting variations of EEG activity in the low frequencies during states of epistemic curiosity.

As future work, we will explore ways to measure curiosity states through EEG without participant-dependent calibration. While trivia questions were used as a trigger of curiosity, new tools (e.g., social robots [START_REF] Ceha | Expression of curiosity in social robots: Design, perception, and effects on behaviour[END_REF], video games) or stimuli (e.g., videos of magic tricks) could be used in future experiments. We will also perform deeper neurophysiological analysis to identify the EEG sensors and sources mostly modulated by curiosity levels. For example, so far, results suggested that most of the information for discriminating curiosity levels are found in theta, alpha and low beta. These frequency bands are similar to the ones used to estimate levels of workload and levels of engagement [START_REF] Dehais | A Neuroergonomics Approach to Mental Workload, Engagement and Human Performance[END_REF]. This is interesting given that both states share some common characteristics with curiosity, such as implications in long-term memory improvements. Their similarities and differences would thus need to be further studied. Further analyses can also be done to compare our results against those obtained with fMRI [START_REF] Kang | The wick in the candle of learning: Epistemic curiosity activates reward circuitry and enhances memory[END_REF], [START_REF] Gruber | States of Curiosity Modulate Hippocampus-Dependent Learning via the Dopaminergic Circuit[END_REF], in order to gain a better understanding of the neurological markers underlying curiosity states. Finally, analyses of physiological signals, e.g., breathing, electrodermal activity or heart rate, could be added to enhance curiosity decoding [START_REF] Knapp | Physiological signals and their use in augmenting emotion recognition for human-machine interaction[END_REF].

Fig. 1 .

 1 Fig. 1. Experiment flow: 1) fixation cross 2) question presentation 3) choice to reveal the answer 4) answer presentation 5) curiosity rating.

Fig. 3 .

 3 Fig. 3. Diagram representing the way we epoched the signals into 1, 2, 3, 4 and 5-seconds time windows (TW).

Fig. 4 .

 4 Fig. 4. F1 score, precision and recall for the 4-seconds time window length, for each subject and for each algorithm.

Fig. 5 .

 5 Fig. 5. Average classification performances (F1-score) across participants obtained by each algorithm with the different time window lengths.

Fig. 6 .

 6 Fig. 6. F1-score for the different time window lengths, for each algorithm.

Fig. 7 .

 7 Fig. 7. Percentage of time that each frequency band was selected by each algorithm, with the 4-seconds time window length.

This study was approved by the ethics committee of Inria Bordeaux Sud-Ouest (COERLE approval number

.

ACKNOWLEDGMENT

This work received support from the European Research-Council (grant ERC-2016-STG-714567) and the Bordeaux-Waterloo Joint Grant 2019-2021.