Yves Coudière
email: yves.coudiere@inria.fr

Charlie Douanla-Lontsi
email: charlie.douanla-lontsi@inria.fr

Charles Pierre
email: charles.pierre@univ-pau.fr

Rush-Larsen time-stepping methods of high order for stiff problems in cardiac electrophysiology

Keywords: stiff equations, explicit high-order multistep methods, exponential integrators, stability and convergence, Dahlquist stability Subject classification: 65L04, 65L06, 65L20, 65L99

Rush-Larsen time-stepping methods of high order for stiff problems in cardiac electrophysiology

Yves Coudière, Charlie Douanla Lontsi, Charles

Introduction

This article concerns the problem of time integration of stiff reaction-diffusion equations, in particular when they are coupled to a system of ordinary differential equations (ODE). As developed below, for such problems, the matters of stability and accuracy are of first importance. As a systemic example of these questions, we will consider the monodomain model in cardiac electrophysiology [START_REF] Clements | Activation dynamics in anisotropic cardiac tissue via decoupling[END_REF][START_REF] Colli-Franzone | Monodomain simulations of excitation and recovery in cardiac blocks with intramural heterogeneity[END_REF][START_REF] Colli-Franzone | Simulating patterns of excitation, repolarization and action potential duration with cardiac bidomain and monodomain models[END_REF]. Given the heart domain Ω and the time interval [0, T ], it has the general form

∂v ∂t = Av + f 1 (v, ζ) + s(x, t), dζ dt = f 2 (v, ζ), (1)
where A is a diffusion operator. The unknown function v : Ω × [0, T ] → R is the transmembrane potential. The unknown function ζ : Ω × [0, T ] → R p+q gathers p + q variables describing the state of the cell membrane. It incorporates p gating variables and q ionic concentrations. The source term s(x, t) is an applied stimulation current. The reaction terms f 1 and f 2 model ionic currents across the cell membrane, and are called ionic models. Ionic models have originally been developed by Hodgkin and Huxley [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF] in 1952.

Highly detailed ionic models specific to cardiac cells have been designed since the 1960's, such as the Beeler and Reuter (BR) model [START_REF] Beeler | Reconstruction of the action potential of ventricular myocardial fibres[END_REF] or the ten Tusscher, Noble, Noble and Panfilov (TNNP) model [START_REF] Tusscher | A model for human ventricular tissue[END_REF]. A comprehensive review is available in [START_REF] Spiteri | Stiffness analysis of cardiac electrophysiological models[END_REF].

There are two major difficulties for numerical simulations in cardiac electrophysiology. First, the non-linear functions f 1 and f 2 in equation ( 1) induce expensive computations of the mappings (v, ζ) → f i (v, ζ). For example, the TNNP model [START_REF] Tusscher | A model for human ventricular tissue[END_REF] involves the computation of 50 scalar exponentials, that have to be performed for each mesh node to approximate solutions of the partial differential equation [START_REF] Beeler | Reconstruction of the action potential of ventricular myocardial fibres[END_REF]. They represent the predominant computational load during numerical simulations, and their total amount needs to be maintained as low as possible. Fully implicit time-stepping methods, which require a non-linear solver, are therefore avoided. Second, the equations (1) are stiff, but since implicit methods are not affordable, numerical instabilities are challenging to manage. More precisely, the stiffness is caused by the presence of different space and time scales. The solutions of equation (1) display sharp wavefronts. Typically, the scaling factor between the fast and slow variables ranges from 100 to 1000. This is commonly coped with by resorting to very fine space and time discretization grids, associated with high computational costs.

In this context, our strategy for solving problem (1) is to use high order methods, so as to have accurate simulations with coarser space and time discretization grids. A high order time-stepping method that fulfills the two following conditions is required: it must have strong stability properties, and has to be explicit for the reaction terms. To this aim, we will focus on the time integration of stiff ODE systems of the form dy dt = f (t, y), y(0) = y 0 ,

in which the nonlinear function f : [0, +∞[×R N → R N (e.g. N = p + q for the ionic models presented above) may be split as f (t, y) = a(t, y)y + b(t, y). This leads to a  formulation more suited to our needs, dy dt = a(t, y)y + b(t, y), y(0) = y 0 .

It involves the non-linear term b(t, y) and the operator y ∈ R N → a(t, y)y ∈ R N , which can be easily linearized as e.g. a(t, y n )y. This term a(t, y) will be inserted into the numerical scheme in order to stabilize the computations. It will be called the stabilizer in the sequel.

In practice, the term a(t, y) may be related to the Jacobian of the system ∂ y f (t, y).

However, no a priori definition of the stabilizer is made (such as a(t, y) = ∂ y f (t, y)), because we plan to analyze the formulation in (3) in general. This will allow us for instance to define the stabilizer as an approximation of the Jacobian, for technical reasons detailed below. This approach is relevant in cardiac electrophysiology, where the fastest variables are gating variables that are given by the p first equations of the ODE system dζ/ dt = f 2 (v, ζ) in [START_REF] Beeler | Reconstruction of the action potential of ventricular myocardial fibres[END_REF]. They have the general form

dζ i dt = ζ i,∞ (v) -ζ i τ i (v) ,
(see Section 4.1) that motivates the reformulation (3) with the diagonal stabilizer a = diag(-1/τ i ). Exponential integrators are well suited in this framework, we refer to [START_REF] Minchev | A review of exponential integrators for first order semilinear problems[END_REF][START_REF] Hochbruck | Exponential integrators[END_REF][START_REF] Hochbruck | A short course on exponential integrators[END_REF]] for general reviews. They have been widely studied for the semilinear equation ∂ t y = Ay + b(t, y), see e.g. [START_REF] Hochbruck | Exponential integrators for large systems of differential equations[END_REF][START_REF] Cox | Exponential time differencing for stiff systems[END_REF][START_REF] Hochbruck | Explicit exponential Runge-Kutta methods for semilinear parabolic problems[END_REF][START_REF] Hochbruck | Exponential Rosenbrock-type methods[END_REF][START_REF] Tokman | New adaptive exponential propagation iterative methods of Runge-Kutta type[END_REF][START_REF] Luan | Explicit exponential Runge-Kutta methods of high order for parabolic problems[END_REF]. Exponential integrators commonly define a time iteration based on the exact solution of an equation of the form ∂ t y = Ay + p(t) where p(t) is a polynomial. It is usually defined with the functions (ϕ k ) k≥0

ϕ 0 (z) = e z , ϕ j+1 (z) = ϕ j (z) -1/j! z , (4)
introduced by Nørsett [START_REF] Norsett | An A-stable modification of the Adams-Bashforth methods[END_REF]. In general, it requires to compute a matrix exponential applied to a vector, like e tA y.This is the supplementary cost associated with exponential integrators. A gain in stability is expected when A is the predominant stiff part of the equation.

The target equation ( 3) incorporates a non-constant linear part a(t, y), exponential integrators have been less studied in that case. Exponential integrators of Adams type for a non-constant linear part have been first considered by Lee and Preiser [START_REF] Lee | A class of non linear multistep A-stable numerical methods for solving stiff differential equations[END_REF] in 1978, and by Chu [START_REF] Chu | An automatic multistep method for solving stiff initial value problems[END_REF] in 1983. Recently, Ostermann et al., [START_REF] Hochbruck | Exponential multistep methods of Adams-type[END_REF] developed and analyzed the linearized exponential Adams method. In general, the original equation ( 2) is formulated after each time step as dy dt = J n y + c n (t, y), involving the Jacobian matrix J n = ∂ y f (t n , y n ), and the correction function c n (t, y) = f (t, y) -J n y. This has several drawbacks. It requires the computation of matrix exponential applied to a vector with a different matrix at each time step. Moreover, stabilization can be performed on the fast variables only, in case they are known in advance, e.g. because of modeling assumptions, or of our physical understanding of the problem. In this case, using the full Jacobian as the stabilizer will  cause unnecessary computational efforts. As an alternative, the stabilizer can be set to a part or an approximation of the Jacobian. This had already been proposed by Nørsett [START_REF] Norsett | An A-stable modification of the Adams-Bashforth methods[END_REF] in 1969 and has been analyzed in [START_REF] Tranquilli | Rosenbrock-Krylov methods for large systems of differential equations[END_REF], [START_REF] Rainwater | A new class of split exponential propagation iterative methods of Runge-Kutta type (sEPIRK) for semilinear systems of ODEs[END_REF], and [START_REF] Coudiére | Exponential Adams Bashforth integrators for stiff ODEs, application to cardiac electrophysiology[END_REF] for exponential Rosenbrock, exponential Runge-Kutta and exponential Adams methods, respectively. For exponential Adams methods, equation (3) is reformulated after each time step as dy dt = a n y + c n (t, y), with a n = a(t n , y n ), and c n (t, y) = f (t, y) -a n y. The resulting scheme with a time-step h > 0 is (see the details in [START_REF] Hochbruck | Exponential multistep methods of Adams-type[END_REF][START_REF] Coudiére | Exponential Adams Bashforth integrators for stiff ODEs, application to cardiac electrophysiology[END_REF])

y n+1 = y n + h (ϕ 1 (a n h) (a n y n + γ 1 ) + ϕ 2 (a n h)γ 2 + . . . ϕ k (a n h)γ k ) , (5)
where the numbers γ i are the coefficients of the Lagrange interpolation polynomial of c n (t, y) (in a classical k-step setting), and the functions ϕ j are given by (4). Independently, Perego and Veneziani [START_REF] Perego | An efficient generalization of the Rush-Larsen method for solving electro-physiology membrane equations[END_REF] presented in 2009 an innovative exponential integrator of order 2, of a different nature. They proposed a scheme of the form

y n+1 = y n + hϕ 1 (α n h) (α n y n + β n ) , (6)
involving two coefficients α n and β n to be computed at each time step. The resulting scheme has a very simple definition, and is in particular simpler than the exponential Adams integrators [START_REF] Colli-Franzone | Simulating patterns of excitation, repolarization and action potential duration with cardiac bidomain and monodomain models[END_REF]. The essential difference with the previous approaches is that α n = a(t n , y n ), but instead is fixed for the scheme to be consistent of order 2. Specifically, the coefficients α n and β n are given by α n = 3 2 a n -1 2 a n-1 and β n = 3 2 b n -1 2 b n-1 with a j = a(t j , y j ) and b j = b(t j , y j ). Perego and Veneziani presented their scheme as a "generalization of the Rush-Larsen method" in reference to the Rush-Larsen scheme [START_REF] Rush | A practical algorithm for solving dynamic membrane equations[END_REF] commonly used in electrophysiology.

This scheme resembles the Magnus integrator introduced by Hochbruck et al. in [START_REF] Hochbruck | On Magnus integrators for time-dependent Schrödinger equations[END_REF] for the time dependent Schrödinger equation iy = H(t)y, and extended by Gonzàlez et al. in [START_REF] González | A second-order Magnus-type integrator for nonautonomous parabolic problems[END_REF] to parabolic equations with time-dependent linear part y = a(t)y+b(t). The secondorder Magnus integrator also formulates as [START_REF] Coudiére | Exponential Adams Bashforth integrators for stiff ODEs, application to cardiac electrophysiology[END_REF], but with α n = a(t n+1/2 ) and b n = b(t n+1/2 ). The scheme of Perego and Veneziani generalizes the second-order Magnus integrator to the case where a = a(t, y) and b = b(t, y): it presents an approximation of the unknown terms a t n+1/2 , y(t n+1/2 ) and b t n+1/2 , y(t n+1/2 ) using a two-points interpolation.

In this paper we will study schemes under the form (6). We will show that they also exist at the orders 3 and 4, and will exhibit explicit definitions of the two coefficients α n and β n . The schemes will be referred to as as Rush-Larsen schemes of order k (shortly denoted by RL k ), in the continuation of the denomination used in [START_REF] Perego | An efficient generalization of the Rush-Larsen method for solving electro-physiology membrane equations[END_REF]. They will be shown to be stable under perturbation and convergent of order k. We also present the Dahlquist stability analysis for the RL k schemes. It is a practical tool that allows one to scale the time step h with respect to the variations of the function f (t, y) in problem [START_REF] Chu | An automatic multistep method for solving stiff initial value problems[END_REF], see e.g. [START_REF] Hairer | Solving Ordinary Differential Equations II[END_REF]. The splitting f (t, y) = a(t, y)y + b(t, y) may be arbitrary, but obviously the choice of the stabilizer term a(t, y) is critical for the stability of the method. When considering time-dependent stabilizers, the stability domain depends on this splitting. We compute stability domains numerically, and show that they are much larger if a(t, y) captures the variations of f (t, y), than in absence of stabilization (i.e., when a(t, y) = 0).



We finally evaluate the performances of the RL k methods as applied to the membrane equation in cardiac electrophysiology. They are compared to the exponential Adams integrators [START_REF] Colli-Franzone | Simulating patterns of excitation, repolarization and action potential duration with cardiac bidomain and monodomain models[END_REF]. The two methods have a very similar robustness with respect to stiffness, allowing stable computations with large time steps. For the considered test case, the RL 3 and RL 4 schemes are more accurate for large time steps.

The paper is organized as follows. The RL k schemes are derived in Section 1, and their numerical analysis is made in Sections 1 and 2. The Dahlquist stability analysis is completed in Section 3. The numerical results are presented in Section 4. The paper ends with a conclusion in Section 5.

In the sequel h denotes the time step, and t n = nh are the associated time instants, starting at t 0 = 0.

Definition of RL k schemes and consistency

Let us consider a solution y(t) of equation ( 3) on a time interval [0, T ]. It is recalled that the scheme ( 6) is consistent of order k if, given a time step h, a time instant kh ≤ t n ≤ T -h, and the numerical approximation y n+1 in ( 6) computed with y n-j = y(t n-j ) for j = 0, . . . , k -1, we have |y n+1 -y(t n + h)| ≤ Ch k+1 , for a constant C only depending on the data a, b, y 0 and T of the problem (3).

Lemma 1. Assume that the functions a(t, y) and b(t, y) are C k regular on [0, T ] × R N . Moreover, assume that a(t, y) is diagonal (a(t, y) = diag (a i (t, y))) or constant. Then the scheme in ( 6) is consistent of order k for k = 2, 3, 4 if

• for k = 2, we have

α n = a n + 1 2 a n h + O(h 2 )
, and

β n = b n + 1 2 b n h + O(h 2 );
• for k = 3, we have

α n = a n + 1 2 a n h + 1 6 a n h 2 + O(h 3 ),
and

β n = b n + 1 2 b n h + 1 12 (a n b n -a n b n )h 2 + O(h 3 );
• for k = 4, we have

α n = a n + 1 2 a n h + 1 6 a n h 2 + 1 24 a n h 3 + O(h 4 ),
and

β n = b n + 1 2 b n h + 1 12 (a n b n -a n b n )h 2 + 1 24 (b n + a n b n -a n b n ) h 3 + O(h 4 );  Yves Coudière,

Remark 1.

The assumption "a(t, y) is diagonal or constant" in Lemma 1 has the following origin. To analyze the consistency of the scheme, we will compute a Taylor expansion in h of the scheme in ( 6). This expansion is derived from Taylor expansions of α n and β n .

For the sake of simplicity, assume the simple form α n = α n,0 + hα n,1 . We need to expand ϕ 1 (α n h) as a power series in h, where the function ϕ 1 is analytic on C. However, in the matrix case, the equality,

ϕ 1 (M + N ) = ϕ 1 (M ) + ϕ 1 (M )N + . . . + ϕ (i) 1 (M )N i /i! + . .
. holds if M and N are commutative matrices. Therefore one cannot expand ϕ 1 (α n h) without the assumptions that α n,0 and α n,1 commute. This difficulty vanishes if a(t, y) is constant or a varying diagonal matrix.

y(t n + h) = y(t n ) + k j=1 s j j! h j + O(h k+1 ),
with s j = y (j) (t n ). Using that y = ay + b we get that

s 1 =a n y n + b n , s 2 =(a n + a 2 n )y n + a n b n + b n , s 3 =(a n + 3a n a n + a 3 n )y n + b n + a n b n + 2a n b n + a 2 n b n , s 4 =(a n + 4a n a n + 3a 2 n + 6a n a 2 n + a 4 n )y n + b n + b n a n + 3a n b n + 5a n a n b n + 3a n b n + a 3 n b n + a 2 n b n .
Series expansions in h for α n and for β n are introduced as

α n = α n,0 + α n,1 h + • • • + α n,k-1 h k-1 + O(h k ), β n = β n,0 + β n,1 h + • • • + β n,k-1 h k-1 + O(h k ).
If the matrix a(t, y) is diagonal or constant (see Remark 1), the Taylor expansion of the numerical solution y n+1 in (6) can be performed

y n+1 = y(t n ) + k j=1 r j j! h j + O(h k+1 ).



Rush-Larsen time-stepping methods of high order

A direct computation of the r j gives

r 1 =α n,0 y n + β n,0 , r 2 =(2α n,1 + α 2 n,0 )y n + 2β n,1 + α n,0 β n,0 , r 3 =(6α n,2 + α 3 n,0 + 6α n,0 α n,1 )y n + 3α n,1 β n,0 + 6β n,2 + α 2 n,0 β n,0 + 3α n,0 β n,1 , r 4 =(24α n,0 α n,2 + 24α n,3 + 12α n,1 α 2 n,0 + 12α 2 n,1 + α 4 n,0 )y n + 12α n,2 β n,0 + 24β n,3 + 12α n,0 β n,2 + 12α n,1 β n,1 + 4α 2 n,0 β n,1 + 8α n,0 α n,1 β n,0 + α 3 n,0 β n,0 ,
where y n denotes y(t n ). The condition to be consistent of order k is:

r i = s i for 1 ≤ i ≤ k.
The consistency conditions in Lemma 1 are obtained by solving recursively these relations.

We then can state our main result, which includes the definition of the RL k schemes.

Theorem 1. Assume (as in Lemma 1) that the functions a(t, y) and b(t, y) are C k regular on [0, T ] × R N , and that a(t, y) is diagonal or constant. Then, the three schemes defined for k = 2, 3, 4 by equation ( 6) and the coefficients,

• for k = 2,

α n = 3 2 a n - 1 2 a n-1 , β n = 3 2 b n - 1 2 b n-1 ,
• for k = 3,

α n = 1 12 (23a n -16a n-1 + 5a n-2 ), β n = 1 12 (23b n -16b n-1 + 5b n-2 ) + h 12 (a n b n-1 -a n-1 b n ),
• for k = 4,

α n = 1 24 (55a n -59a n-1 + 37a n-2 -9a n-3 ), β n = 1 24 (55b n -59b n-1 + 37b n-2 -9b n-3 ) + h 12 (a n (3b n-1 -b n-2 ) -(3a n-1 -a n-2 )b n ),
where a j = a(t j , y j ) and b j = b(t j , y j ), are consistent of order k.

The three methods stated above are called Rush-Larsen methods of order k, and denoted by RL k . They are explicit and k-step methods.



Remark 2. If the matrix a is a constant, a(t, y) = A, then we have α n = A for all three methods. In this case, the expressions of the coefficients β n in Theorem 1 for k = 3, 4 simplify as follows:

RL 3 case :

β n = 1 12 (23b n -16b n-1 + 5b n-2 ) - h 12 A(b n -b n-1 ).
RL 4 case :

β n = 1 24 (55b n -59b n-1 + 37b n-2 -9b n-3 ) - h 12 A(2b n -3b n-1 + b n-2 ).
Proof. It is a direct consequence of backwards differentiation formulas, that we first recall.

The derivatives of a real function f at the time instant t n can be approximated as follows (with common notations):

• first derivative,

f n = f n -f n-1 h + O(h) = 1 2h (3f n -4f n-1 + f n-2 ) + O(h 2 ) = 1 6h (11f n -18f n-1 + 9f n-2 -2f n-3 ) + O(h 3 );
• second derivative,

f n = 1 h 2 (f n -2f n-1 + f n-2 ) + O(h) = 1 h 2 (2f n -5f n-1 + 4f n-2 -f n-3 ) + O(h 2 );
• third derivative,

f n = 1 h 3 (f n -3f n-1 + 3f n-2 -f n-3 ) + O(h).
With these formulas, the consistency condition at order 3 on the coefficient α n becomes

α n = a n + 1 2 a n h + 1 6 a n h 2 + O(h 3 ) = a n + 1 4 (3a n -4a n-1 + a n-2 ) + 1 6 (a n -2a n-1 + a n-2 ) + O(h 3 ) = 1 12 (23a n -16a n-1 + 5a n-2 ) + O(h 3 ).
We retrieve the definition of α n for the RL 3 scheme. The same proof holds for β n , and extends to order 4.



We refer to [10, Ch. III-8] for the definitions of convergence and of stability under perturbation. For the analysis of time-stepping methods, it is commonly assumed that f in equation ( 2) is uniformly Lipschitz with respect to its second variable y. This hypothesis will be replaced by assumptions based on the formulation [START_REF] Clements | Activation dynamics in anisotropic cardiac tissue via decoupling[END_REF]. Precisely it will be assumed that a(t, y) is bounded, a(t, y), b(t, y) are uniformly Lipschitz in y.

The Lipschitz constants of a and b are denoted by L a and L b , respectively. The upper bound on |a(t, y)| is denoted by M a .

Proposition 1. If the assumption (7) holds, then, the RL k schemes are stable under perturbation, for k = 2, 3, 4. In addition, also for k = 2, 3, 4, if the consistency assumptions of Theorem 1 are satisfied (a(t, y) and b(t, y) are C k regular and a(t, y) is diagonal or constant), then the RL k scheme is convergent of order k.

Stability under perturbation together with consistency implies (nonstiff) convergence, see e.g. [START_REF] Hairer | Solving Ordinary Differential Equations I[END_REF], or [START_REF] Coudiére | Exponential Adams Bashforth integrators for stiff ODEs, application to cardiac electrophysiology[END_REF], where the current setting has been detailed. Therefore the proof of the convergence statement in Proposition 1 is immediate, and will not be recalled here.

The following definitions are necessary to prove Proposition 1. Equation ( 2

s t,h : Y = (y 1 , . . . , y k ) ∈ E k -→ s t,h (Y ) ∈ E, with s t,h (Y ) = y k + hϕ 1 (α t,h (Y )h) (α t,h (Y )y k + β t,h (Y )) ,
in such a way that the scheme in (6) reads y n+1 = s tn,h (y n-k+1 , . . . , y n ). The functions α t,h and β t,h map the vector Y of the k previous values to the values α n and β n given in Theorem 1. For instance, the function α t,h (Y ) for k = 3 (the RL 3 scheme) reads

α t,h (Y ) = 1 12 (23a(t, y 3 ) -16a(t -h, y 2 ) + 5a(t -2h, y 1 )).
A first technique to prove the stability under perturbation consists in showing that the function s t,h is globally Lipschitz in Y . To this aim, the derivative ∂ Y s t,h has to be analyzed. As developed in Remark 1, it implies restrictions on the function a(t, y): it has to be either diagonal or constant. A second technique consists in proving the following two stability conditions:

|s t,h (Y ) -s t,h (Z)| ≤ |Y -Z| ∞ (1 + Ch(|Y | ∞ + 1)) , ( 8
) |s t,h (Y )| ≤ |Y | ∞ (1 + Ch) + Ch, (9)
 for all Y and Z in E k , and where the constant C depends only on the data a, b, y 0 in equation ( 3), and on the final time T . These are sufficient conditions for the stability under perturbation, as proved in [6, Section 2]. We will use the conditions ( 8) and ( 9) here, because they are more general, and give rise to less computations. The core of the proof is the following property of the RL k scheme. For Y = (y 1 , . . . , y k ) ∈ E k , we have

s t,h (Y ) = z(t + h) for z = α t,h (Y )z + β t,h (Y ), z(t) = y k . (10)
It will be used together with the following Gronwall inequality (see [START_REF] Dragomir | Some Gronwall type inequalities and applications[END_REF]Lemma 196,p.150]). Suppose that z(t) is a C 1 function, and that there exists

M 1 > 0 and M 2 > 0 such that |z (t)| ≤ M 1 |z(t)| + M 2 for all t ∈ [t 0 , t 0 + h]. Then ∀t ∈ [t 0 , t 0 + h], |z(t)| ≤ e M 1 (t-t 0 ) (|z(t 0 )| + M 2 (t -t 0 )) . ( 11
)
Proposition 1. In this proof, we always assume that 0 ≤ h, t ≤ T , and denote by C i a constant that depends only on the data a, b and T of problem [START_REF] Clements | Activation dynamics in anisotropic cardiac tissue via decoupling[END_REF]. With the assumptions in [START_REF] Cox | Exponential time differencing for stiff systems[END_REF], and the definitions of α n (k = 2, 3, 4) in Theorem 1, the function α t,h is uniformly Lipschitz with a Lipschitz constant equal to L α . Moreover we have the uniform bound α t,h ≤ M α . Since the function b(t, y) is uniformly Lipschitz with respect to y, and since 0 ≤ t ≤ T , we have

|b(t, y)| ≤ |b(t, 0)| + |b(t, y) -b(t, 0)| ≤ K b (1 + |y|), (12)
with K b = max(L b , sup 0≤t≤T |b(t, 0)|). For the RL 3 scheme, we have

|β t,h (Y )| ∞ ≤ 11 3 K b (1 + |Y | ∞ ) + h 12 M a 2K b (1 + |Y | ∞ ) ≤ C 1 (1 + |Y | ∞ ).
The same inequality holds for the RL 2 , and RL 4 schemes. Afterwards, we can apply these bounds to the differential equation in ( 10)

|z | = |α t,h (Y )z + β t,h (Y )| ≤ M α |z| + C 1 (1 + |Y | ∞ ). The initial state is |z(t)| = |y k | ≤ |Y | ∞ .
Finally, the Gronwall inequality (11) yields, for

t ≤ τ ≤ t + h, |z(τ )| ≤ e Mαh (|Y | ∞ + hC 1 (1 + |Y | ∞ )) ≤ e Mαh (|Y | ∞ (1 + C 1 h) + C 1 h) ≤ |Y | ∞ (1 + C 2 h) + C 2 h, (13)
by bounding the exponential with an affine function for 0 ≤ h ≤ T . This gives the stability condition (9) for τ = t + h. For the RL 2 scheme, the function β t,h is uniformly Lipschitz. For the RL 3 scheme, for Y = (y 1 , y 2 , y 3 ) and Z = (z 1 , z 2 , z 3 ) in E 3 , we have

|β t,h (Y ) -β t,h (Z)| ∞ ≤ 11 3 L b |Y -Z| ∞ + h 12 |a(t, y 3 )b(t -h, y 2 ) -a(t, z 3 )b(t -h, z 2 )| + |a(t -h, y 2 )b(t, y 3 ) -a(t -h,
|F (Y ) -F (Z)| = |a(ξ, y 3 ) (b(τ, y 2 ) -b(τ, z 2 )) + (a(ξ, y 3 ) -a(ξ, z 3 )) b(τ, z 2 )| ≤ M a L b |Y -Z| ∞ + L a |Y -Z| ∞ |b(τ, z 2 )|.
With the inequality [START_REF] Hochbruck | A short course on exponential integrators[END_REF], this yields, for 0 ≤ τ, ξ ≤ T , and

Y , Z in E k , |F (Y ) -F (Z)| ≤ C 3 |Y -Z| ∞ (1 + |Z| ∞ ).
As a result, we have

|β t,h (Y ) -β t,h (Z)| ∞ ≤ C 4 |Y -Z| ∞ (1 + |Z| ∞ ) .
The same inequality holds for the RL 4 scheme. Finally we consider Y 1 and Y 2 in E k , and the notation α i = α t,h (Y i ), and

β i = β t,h (Y i ).
The property [START_REF] Hairer | Solving Ordinary Differential Equations I[END_REF]

shows that s t,h (Y 1 ) -s t,h (Y 2 ) = (z 1 -z 2 )(t + h), where z i is the solution to z i = α i z i + β i , with z i (t) = Y i,k .
On the first hand, with the inequality (13), we have

|z 2 (τ )| ≤ C 5 (1 + |Y 2 | ∞ ) for t ≤ τ ≤ t + h. On the second hand, on [t, t + h], we have |(z 1 -z 2 ) | ≤ |α 1 ||z 1 -z 2 | + |α 1 -α 2 ||z 2 | + |β 1 -β 2 | ≤ M α |z 1 -z 2 | + L α |Y 1 -Y 2 | ∞ C 5 (1 + |Y 2 | ∞ ) + C 4 |Y 1 -Y 2 | ∞ (1 + |Y 2 | ∞ ) ≤ M α |z 1 -z 2 | + C 6 |Y 1 -Y 2 | ∞ (1 + |Y 2 | ∞ ).
The initial condition yields

|(z 1 -z 2 )(t)| = |Y 1,k -Y 2,k | ≤ |Y 1 -Y 2 | ∞ .
As a consequence, the Gronwall inequality [START_REF] Hairer | Solving Ordinary Differential Equations II[END_REF] applied to these bounds shows that

|(z 1 -z 2 )(t + h)| ≤ e Mαh (|Y 1 -Y 2 | ∞ + hC 6 |Y 1 -Y 2 | ∞ (1 + |Y 2 | ∞ )) ≤ e Mαh |Y 1 -Y 2 | ∞ (1 + C 6 h(1 + |Y 2 | ∞ )) .
This last inequality implies the stability condition [START_REF] Dragomir | Some Gronwall type inequalities and applications[END_REF], again by bounding the exponential with an affine function for 0 ≤ h ≤ T .

Dahlquist stability

For the general ideas and definitions concerning the Dahlquist stability we refer to [START_REF] Hairer | Solving Ordinary Differential Equations II[END_REF].

The background for the Dahlquist stability of exponential integrators with a general varying stabilizer a(t, y) has been developed in [START_REF] Coudiére | Exponential Adams Bashforth integrators for stiff ODEs, application to cardiac electrophysiology[END_REF], following the ideas of Perego and Veneziani [START_REF] Perego | An efficient generalization of the Rush-Larsen method for solving electro-physiology membrane equations[END_REF]. The equation ( 2) is considered with the Dahlquist test function f (t, y) = λy, which is split into f (t, y) = a(t, y)y + b(t, y), in order to match the framework of equation ( 3), with a(t, y) = θλ, b(t, y) = λ(1 -θ)y.

For θ = 1, the methods are exact and thus A-stable. For θ 1, the exact linear part of f (t, y) in equation ( 2) is well approximated by a(t, y). The stability domain depends on θ, it is denoted by D θ . Given a value of θ, the region D θ is defined by the modulus of a stability function, with the same definition as for multistep methods, see e.g. [START_REF] Hairer | Solving Ordinary Differential Equations II[END_REF]. This stability function has been numerically computed, pointwise on a grid inside the complex plane C, for each of the three RL k schemes, k = 2, 3, 4. Order 2 Rush-Larsen

The stability domain for the RL 2 scheme has been analyzed in [START_REF] Perego | An efficient generalization of the Rush-Larsen method for solving electro-physiology membrane equations[END_REF]. The situation for this scheme is interesting, and we reproduced the results on Figure 1. We note the observations below.

• If 0 ≤ θ < 2/3, the stability domain D θ is bounded. Its size increases with θ, starting from the stability domain without stabilization for θ = 0, that corresponds to the Adams-Bashforth scheme of order 2.

• If θ = 2/3, the method is A(0) stable: R -⊂ D θ . The domain boundary is asymptotically parallel to the real axis, so that the method is not A(α) stable.

• If θ > 2/3, the stability domain is located around the y-axis: the method is A(α) stable. The angle α increases with θ, it goes to π/2 as θ → 1 -.

Rush-Larsen methods of orders 3 and 4

The situation is different for the Rush-Larsen methods of orders 3 and 4. The stability domains D θ are depicted on Figures 2 and3, for various values of θ, and for the orders 3 and 4, respectively. Excepted for the case θ = 1, the stability domain is always bounded: the scheme is not A(0)-stable. However, the stability domain for θ 1 is much larger than without stabilization (corresponding to the Adams-Bashforth schemes of orders 3 or 4). For the RL 3 scheme, the stability domain when θ = 0.85 is 25 times wider on the left than D θ|θ=0 , and when θ = 1.05 it is 400 times wider. For the RL 4 case, D θ|θ=1.05 is almost 300 times wider on the left than D θ|θ=0 . In the particular case θ = 0 (no stabilization, corresponding to the Adams-Bashforth scheme of order 3), the stability domain crosses the x-axis at x -0.3 (dark blue arrow).

Numerical results

In this section we present numerical experiments that illustrate the performances of the RL k methods. They will be compared to the exponential integrators of Adams type of order k defined by equation ( 5), shortly denoted by EAB k . The EAB k schemes have been numerically studied in [START_REF] Coudiére | Exponential Adams Bashforth integrators for stiff ODEs, application to cardiac electrophysiology[END_REF], for the resolution of the membrane equation in electrophysiology, as compared to several classical methods. In that context, they have been shown to be as stable as implicit methods with a much smaller cost. We present the same numerical tests here, so as to extend the comparison to the schemes benchmarked in [START_REF] Coudiére | Exponential Adams Bashforth integrators for stiff ODEs, application to cardiac electrophysiology[END_REF]. The cellular action potential for cardiac cells is described on Figure 4. This phenomenon displays a stiff behavior characterized by the presence of heterogeneous time scales. The electrical activity of cardiac cells is modeled with an ODE system called membrane equation. It has the form

The membrane equation

dv dt = -I ion (v, w, c) + I st (t), dw i dt = w ∞,i (v) -w i τ i (v) , dc dt = g(v, w, c), (14)
where w = (w 1 , . . . , w p ) ∈ R p is a vector of gating variables, c ∈ R q is a vector of ionic concentrations, and v ∈ R is the transmembrane potential, we refer to [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF][START_REF] Beeler | Reconstruction of the action potential of ventricular myocardial fibres[END_REF][START_REF] Luo | A dynamic model of the cardiac ventricular action potential[END_REF][START_REF] Tusscher | A model for human ventricular tissue[END_REF] for details. The gating variables describe the opening state (between 0 and 1) of various protein structures on the cell membrane, which control ionic transfers between the intra and extra-cellular media. Each gating variable w i evolves towards the state w ∞,i (v) at rate τ i (v). Specific ionic currents (sodium, potassium, ...) across the cellular membrane are computed with the help of the variables v, w and c. The sum of these currents defines the total ionic current I ion (v, w, c) across the membrane. The function I st (t) is a source term, it represents a stimulation current. The membrane equation corresponds to the ODE system in the monodomain model ( 1) with ζ = (w, c).

We will consider two such models: the BR model [START_REF] Beeler | Reconstruction of the action potential of ventricular myocardial fibres[END_REF] and the TNNP model [START_REF] Tusscher | A model for human ventricular tissue[END_REF].The BR model [START_REF] Beeler | Reconstruction of the action potential of ventricular myocardial fibres[END_REF] describes the membrane action potential of mammalian ventricular myocardial cells. It involves 6 gating variables (p = 6, they are denoted by m, h , j , d , f , x i ) and one ionic concentration (q = 1): the intra-cellular calcium [Ca i ]. The Nernst potential E Ca for the calcium ions then is time-dependent and a (slow inward) calcium current I s is modeled as I s = g s df (v -E Ca ) depending on the gating variables d, f , the transmembrane potential v and a constant g s . A fast inward sodium current I N a (depicted on Figure 4) that depends on the three gating variables m, h, j and on v is similarly described. Two outward currents are modeled: I x 1 that depends on x 1 and v and I K 1 that only depends on v. The total ionic current in ( 14) is the sum of these four currents

I N a + I s + I K 1 + I x 1 = -I ion (w, c, v).
The TNNP model is specifically designed for human ventricular myocytes. It is more sophisticated than the BR model but conserves the same general structure. It involves 12 gating variables and 4 ionic concentrations (p = 12 and q = 4). The total ionic current I ion is the sum of 15 specific ionic currents.

The membrane equation ( 14) can be reformulated in the form of (3) with

a(t, y) =   -1/τ (v) 0 0 0 0 0 0 0 0   , b(t, y) =   w ∞ (v)/τ (v) g(y) -I ion (y) + I st (t)   , (15)
for y = (w, c, v) ∈ R N (N = p + q + 1) and where -1/τ (v) is the p × p diagonal matrix with diagonal entries (-1/τ i (v)) i=1,...,p . The resulting matrix a(t, y) is diagonal.

Convergence

No analytical solution is available for the chosen application. A reference solution y ref for a reference time step h ref is computed with the Runge-Kutta scheme of order 4 to analyze the convergence properties of the RL k schemes. Numerical solutions y are compared to y ref for coarsest time steps h = 2 m h ref .

A numerical solution y consists in successive values y n at the time instants t n = nh. On every interval (t 3n , t 3n+3 ) the polynomial y of degree at most 3 so that y(t 3n+i ) = y 3n+i , i = 0, . . . , 3 is constructed. On (0, T ) y is continuous and piecewise polynomial of degree 3 its values at the reference time instants nh ref are computed. This provides a projection P (y) of the numerical solution y onto the reference grid. Then P (y) can be compared with the reference solution y ref .

The numerical error is defined by

e(h) = max |v ref -P (v)| max |v ref | , (16)
where the potential v is the last and stiffest component of y in equation [START_REF] Hochbruck | Exponential integrators for large systems of differential equations[END_REF]. The convergence graphs for the BR model are plotted on Figure 5. Each scheme displays the expected asymptotic behavior of Proposition 1: e(h) = O(h k ) as h → 0.

Stability

Spiteri et al. in [START_REF] Spiteri | Stiffness analysis of cardiac electrophysiological models[END_REF] have evaluated the stiffness of the BR and TNNP models along one cellular action potential (as depicted on Figure 4). The largest negative real part of the eigenvalues of the Jacobian matrix during the action potential is of -1170 and -82 for the TNNP and BR models, respectively. The TNNP model thus is 15 times stiffer than the BR model (15 1170/82).

Robustness to stiffness for the RL k schemes is evaluated by comparing the critical time steps for these two models. The critical time step ∆t 0 is defined as the largest time step such that the numerical simulations run without overflow for h < ∆t 0 . The results are presented in Table 1 An excellent robustness to stiffness can be observed. the RL k schemes are not A(α) stable, and the critical time step is expected to be divided by 15 in case of an increase of stiffness of magnitude 15. It is here divided by 2.7, 2.0 and 1.3 for k = 2, 3 and 4, respectively. A comparison with the EAB k schemes shows that the two schemes have similar robustness to stiffness. Loss of stability is induced by the non-stabilized part, whose eigenvalues are less modified by the change of model.

Accuracy

The RL k schemes are compared to the EAB k schemes in terms of accuracy. This is done using the relative error e(h) in equation ( 16), for the BR and TNNP models (we recall than the TNNP model is stiffer by a factor of 15). The results are collected in Tables 2 and3. For the RL 2 and the EAB 2 schemes, the accuracies are very close, the EAB 2 scheme being slightly more accurate for the BR model. For the orders 3 and 4, the RL k schemes are more accurate at large time steps. For smaller time steps, accuracies are almost the same. The RL k and EAB k have the same accuracy in the asymptotic convergence region.

Conclusion

In this paper, we have introduced two new ODE solvers, that we have called Rush-Larsen schemes of orders 3 and 4. They are explicit multistep exponential integrators. Their definition is simple inducing an easy implementation. We exposed the analysis of convergence and of stability under perturbation for these two schemes. We also analyzed their Dahlquist stability: they are not A(0) stable, but exhibit very large stability domains for sufficiently accurate stabilization. The numerical behavior of the schemes is analyzed for a complex and realistic stiff application. The RL k schemes are as stable as exponential integrators of Adams type, allowing simulations at large time step. On the presented example, the RL k schemes are more accurate for k = 3 and 4 than the exponential integrators of Adams type, when considering larger time steps. They are also shown to be robust to stiffness in terms of both stability and accuracy.

  Proof. Consider equation (3) on the closed time interval [0, T ], and its solution, the function y. Since the functions a and b are C k regular on [0, T ] × R N , the solution y(t) is C k+1 regular on [0, T ]. Its derivatives up to order k + 1 are bounded by constants only depending on the data of problem (3), and on T . The Taylor expansion of y at time instant t n is

  ) is considered on E = R N with the max norm | • |. A final time T > 0 is considered. The space of N × N matrices is equipped with the operator norm • associated with | • |. The space E k is equipped with the max norm |Y | ∞ = max 1≤i≤k |y i | with Y = (y 1 , . . . , y k ). The RL k scheme is defined by the mapping

  z 2 )b(t, z 3 )|  Rush-Larsen time-stepping methods of high order Let us bound the Lipschitz constant of a function of the type F (Y ) = a(ξ, y 2 )b(τ, y 3 ), for 0 ≤ τ, ξ ≤ T :

Figure 1 :

 1 Figure 1: Stability domain D θ for the RL 2 scheme for various values of θ. The stability domain for the particular case θ = 0 (no stabilization) is in gray, corresponding to the Adams-Bashforth scheme of order 2.

1 Figure 2 :

 12 Figure2: Stability domain D θ for the RL 3 scheme. In the particular case θ = 0 (no stabilization, corresponding to the Adams-Bashforth scheme of order 3), the stability domain crosses the x-axis at x -0.54 (dark blue arrow).

Figure 3 :

 3 Figure3: Stability domain D θ for the RL 4 scheme. In the particular case θ = 0 (no stabilization, corresponding to the Adams-Bashforth scheme of order 3), the stability domain crosses the x-axis at x -0.3 (dark blue arrow).

Figure 4 :

 4 Figure4: TNNP model illustration. Left, cellular action potential: starting from a (negative) rest value, the transmembrane potential v(t) has a stiff depolarization followed by a plateau and a repolarization to the rest value. Right, depolarization is induced by an ionic sodium current I N a , with obvious large stiffness.

Figure 5 :

 5 Figure 5: Relative error e(h) (definition (16)) as a function of the time step h for the RL k schemes, for k = 2 to 4 and in Log/Log scale.

  53 10 -3 5.96 10 -3 2.14 10 -2 7.34 10 -3 8.34 10 -3 0.0125 5.75 10 -3 8.05 10 -4 3.21 10 -4 5.11 10 -3 7.62 10 -4 3.70 10 -4

  Charlie Douanla-Lontsi and Charles Pierre where a n , a n , a n and b n , b n , b n denote the successive derivatives at time t

n of the functions t → a(t, y(t)) and t → b(t, y(t)).

Table 1 :

 1 . Critical time steps ∆t 0 for the RL k and EAB k schemes method RL 2

		RL 3	RL 4	EAB 2 EAB 3	EAB 4
	BR	0.323 0.200 0.149 0.424 0.203	0.123
	TNNP 0.120 0.148 0.111 0.233 0.108 7.56 10 -2

Table 2 :

 2 Relative error e(h) (eq. (16)) for the BR model. -2 5.86 10 -2 9.26 10 -2 9.17 10 -2 0.119 0.05 3.35 10 -2 6.34 10 -3 4.58 10 -3 2.31 10 -2 1.09 10 -2 8.96 10 -3 0.025 8.88 10 -3 7.57 10 -4 2.61 10 -4 5.39 10 -3 1.17 10 -3 4.33 10 -4

	h	RL 2	RL 3	RL 4	EAB 2	EAB 3	EAB 4
	0.2	0.251	0.147	-	0.284	0.516	-
	0.1	0.107	4.07 10

Table 3 :

 3 Relative error e(h) (eq. (16)) for the TNNP model. -2 4.54 10 -2 4.61 10 -2 9.01 10 -2 5.59 10 -2 8.93 10 -2 0.025 2.21 10 -2 6.

	h	RL 2	RL 3	RL 4	EAB 2	EAB 3	EAB 4
	0.1	0.177	0.305	0.421	0.351	0.530	-
	0.05	7.39 10

Acknowledgments. This study received financial support from the French government as part of the "Investissement d'Avenir" program managed by the "Agence Nationale de la Recherche" (ANR), grant reference ANR-10-IAHU-04. It also received fundings of the ANR project HR-CEM, grant reference 13-MONU-0004-04.


Laura Garcia De Dios

Marianne Boix

Sabine Simeoni-Sauvage

Iréa Touche

Roxelane Cakir

Ludovic Montastruc

Jose-Miguel Sanchez-Perez

 

 

 

Oatao

Laura Garcia De Dios3

Maii Anne Boix

Sabine Sauvage\ Fréa Touche3

Roxelane Cakir\

Ludovic Montastrnc

Jose Miguel Sanchez Perezb

Multiobjective optimization of eco-industrial parks: evaluation of environmental impacts at the watershed scale

Keywords: eco-industrial parks, multi-objective optimization, watershed scale, Material Flow Analysis

come   L'archive ouverte pluridisciplinaire

Introduction

Water is a fragile, precious and litnited resource and is not only itnportant for domestic or sanitaiy pw-poses, but every human economic activity relies on the existence of this element. In addition, leavit1g the human scope, every marine and ten-estrial ecosystem depends on water somehow. Tiùs fact underlines the huge impotiance of this resource. Nevertheless, many recent reports have demonstrated that water is in danger. Indeed, there is a clear evidence that in many rivers and aquifers there has been a reduction in its water quantity and quality. Water scarcity is a problem in many countries and many ecosystems have been modified until a point of no return. One of the most prominent problem is the pollution caused by agriculture and industry. This situation can cause significant damages, not only to ecosystems, but also to humans.

In conclusion, most of the problems evoked have a human origin and as a result, governments and society in general must act to reverse this situation. Efficiency in water use must be increased, technology must be developed and in this situation, research has an important role, because there are still many questions to answer and strategies to establish. All this concern about resources in general, not only water, led to the birth of a great deal of new concepts such as sustainable development or industrial ecology. Industrial ecology consists of understanding human industrial activities as a natural ecosystem (Frosh and Gallopolous, 1989). Indeed, a natural ecosystem represents equilibrium, sustainability and the use of every resource in the most efficient way as possible, reducing consumption and waste. One of the most important ideas in implementing Industrial Ecology is to establish eco-industrial parks. An eco-industrial park can be understood as a group of facilities in the same geographic region that try to optimise the use of energy, water and raw material flows, looking for economic benefits for all of them and for a better environmental performance [START_REF] Lowe | Creating by-product resource exchanges: strategies for eco-industrial parks[END_REF]. Nevertheless, among all the eco-industrial parks that have been built, most of them are working in suboptimal conditions because they have been built somehow, spontaneously. In other words, gradually some improvements were introduced, but a complete industrial park reconfiguration was never established. The objective of this research is to introduce the conception of robust mathematical methods to design eco-industrial parks while taking into account the effects of such a policy at a bigger scale: the watershed scale. In order to convince the industries owners to be a part of these changes and to obtain real profits from them through a multiobjective optimization.

Objectives and methodology

During this research, several EIPs are going to be established within a region, from existing facilities or existing industrial areas. The reconfiguration of these facilities into eco-industrial parks is going to be limited to their water networks. The idea is to investigate deeply the impact of this politic at a bigger scale, the scale of the whole The river flow along the watershed is modeled through a Material Flow Analysis (MFA), so its quality eau be studied before and after the reconfigura.fions of the ecoindustrial pru-ks. To strut with, the river is divided into several pruts (Figure 2). The main pruts are going to be called reaches ru1d the secondruy sections wi.11 be the tributru•ies. In general, this model is going to consist of some mass balru1ces per reach and tJ.ibuta1y, tlùs allows to solve several consecutive small problems rallier tllan a big non-lineru• one, more difficult to solve. This problem is fon.nulated as a Non-Linear Problem (NLP) that has been solved witliin GAMS enviromnent witli DICOPT solver. Its most imp01tant vrui.able is the pollutru1t concentration at the final disposa!. Witli tlie changes on this vrui.able before and after the establishment of the eco-industrial park their effects on tlie environment will be evaluated.

Tribuby

Moiti-objective optimization of eco-industrial parks

As it has been showed in the previous figure, witliin a single watershed there can be several industiial pai-ks or clusters of comp31Ùes. Moreover, these industrial complexes ru•e constituted by individual plants modelled as a group of process sources (i) and process sinks (j) to be reconfigured as EIP. Besides, a regeneration to clean the water that has been used is added. This unit can be called an interception unit as well. The interception unit consists of several interceptors (int). The superstructure defines the whole set of available solutions in a network so that all the possible links between units can be seen. For instance, the water used by a process unit leaves the sources to be used again in a process sink if some conditions, in relation to quality, are achieved. This water can also be sent to the regeneration unit. Furthermore, the flow leaving the sources can be sent directly to the environment without any previous treatment if the concentration constraints are respected. The process sinks can receive water from the sources, the interception unit or from fresh water sources directly. Finally, the interception unit is constituted by several interceptors, as it has been mentioned before. Each interceptor receives used water from process sources, it cleans the water and it can send it, either to the environment or to other process sinks so the water can be used again. The companies that constitute the industrial parks are going to be working on their own at the beginning. With this methodology, all the plants characterized by their sources, sinks and interceptors are going to start working together, by sharing their water networks. For example, the water of one source of one specific plant can feed any sink of any plant; it can go to any interceptor of the industrial complex, etc. This way, the industrial plants will integrate an eco-industrial park. This optimization step will give a new optimal topological configuration of the water network as a solution by minimizing the total cost and environmental impact (fresh water consumption). The problem is formulated as an MILP and solved within GAMS environment with CPLEX solver. The main equations are mass balances for each sink, source and interceptors, the total cost and environmental functions and more constraints that complete the model [START_REF] Boix | Industrial water management by multiobjective optimization: from individual to collective solution through eco-industrial parks[END_REF].

Case study: the Garonne watershed in France

Description of the case study

In this study, four possible eco-industrial parks are considered to be placed in the Garonne watershed. To model the Garonne watershed firstly, it is necessary to divide the watershed into different sub-watersheds, as the MFA model requires. Secondly, some information about the watershed was needed for each partition that had been done, for example the precipitations rate, industrial, residential and agricultural discharges, losses, usages, initial flows, etc. To reach this ptupose, a model is developed in ArcSWAT [START_REF] Arnold | Soil and Water Assessment Tool (SWAT), Introductory Manual[END_REF] for the watershed of Garonne. The watershed has been divided into 22 sub-basins and the final disposai has been placed in Tonneins. Also, it has been decided to place the possible eco-industrial parks in Toulouse, Montauban, Cahors and Agen (Figure 3). Firstly, with all this infonnation, the MFA model for the Garonne watershed with nonnal industriaI parks (without any collaboration) was developed and the concentration of pollutant at the final disposai is obtained: 2.03 ppm. Tuen, the multi-objective optirnization of all the fom industrial parks is catried out and the infonnation obtained in relation to the pollutant concentration and the wastewater is introduced a.gain into the MF A model to evaluate the changes in the pollutant concentration.

Reconfiguration of indust1ial parks into EIP by mufti-objective optimization

In this study, fotu• possible eco-industria.l pat-ks are considered to be placed in the Garonne watershed. For each of them a multi-objective optimisation process has been canied out. The optimisation variables were fresh water, regenerated water, munber of connections, total cost and the pollutant concentration expel to the envirolllllent. To solve this problem, for each eco-industiial pat•k a Pai•eto front has been drawn and, from this set of optimal solutions, the optimal TOPSIS solution has been chosen. In order to illusti-ate tllis step, Figure 4 represents one of the fotu• EIP Pareto fronts that have been built; the figure is related to EIP 11°1 and represents the fresh water consmnption (enviromnental objective) with the regenerated water (economic objective) for severnl ntunber of cotmexions.

3. Environmental impacts of EIP policy at the watershed scale

Severa! strategies were studied: first, 01ùy one of the fotu• EIP is working while the others stay as classica.l industiia.l parks, then some combinations of them and finally a.Il of them are working at fue same time. All these configtuc1tions were then compa.red to the pollutant concenti-ation at the final disposai, without a.ny eco-industiial park in the watershed. Figttre 5 smmnarizes the results obtained dming tllis study where the pollutant concenti•ation ate the fmaI disposai is repo1ted for the different strategies explored. Sorne imp01tant results must be tu1derlined:

A global diminution of the pollutant concentration at the estuaty is observed for every strategy ofEIP development

Conclusions

As a general conclusion, it must be highlighted that the implementation of industrial ecology strategies is always interesting. The development of eco-industrial parks reduces pollution and studying the eco-industrial parks effects on a bigger perspective or scale, such as the scale of a whole watershed. This study helps us to understand their influence or beneficial impacts in a better way and the strategy developed aims at implementing this strategy with a higher performance. The optimal location for these eco-industrial parks in the watershed has a major influence in their final environmental impact. Nevertheless, this task is not completely finished. Some improvements can be added to both models, the EIP's model and the MFA model. Starting with the first one, this model can be enhanced by making the code multi-pollutant, the cost equations can be improved, other ways of introducing the ecological effects into the code or more ecological parameters could be studied. Moreover, the eco-industrial park concept can be extended by sharing not only water but also more resources such as energy or other materials.

  watershed. The objectives of this work are summarized below: • To model a watershed. A Material Flow Analysis (MFA) to track and measure the quality of the river water is developed based on Lopez-Diaz et al. (2015). • To develop a multi-objective optimization model for the water networks reconfiguration. • To study the impacts of these changes on the watershed.• To look for data to apply, specifically, this model to the Adour-Garonne watershed in France. According to Figure1, the two models interact with firstly the run of the MFA model to obtain the parameters of the watershed without the eco-industrial park's effects. Then, in the optimization model (MILP), the data related to water network of the industries are changed to make the industries work as eco-industrial parks. Finally, the results of the EIPs reconfiguration are loaded again into the MFA model to obtain the effects of these modifications in the watershed, following the methodology developed by Lopez-Diaz et al .

Figure 2 :

 2 Figure 2: Decomposition of the river basin in reaches and tributaries (from Lapez-Diaz, 2015)

Figure 3 :

 3 Figure 3: Garonne watershed

Figure 4 :

 4 Figure 4: Pareto Front of the EIP 1

Figure 5 :

 5 Figure 5: Pollutant concentration at the final disposal under several strategies of EIP development

  -Each EIP does not have the same effect on the overall watershed, for example, EIP 3 leads to a greater reduction than EIP 4. This is mainly because EIPs have different sizes and are located at different places in the watershed what causes different environmental impacts.


Asmaa Amehraye
email: asmaa.amehraye@telecom-bretagne.eu.

Dominique Pastor
email: dominique.pastor@telecom-bretagne.eu.

Ahmed Tamtaoui
email: tamtaoui@inpt.ac.ma.

Driss Aboutajdine

From Maskee to Audible Noise in Perceptual Speech Enhancement

Keywords: Perceptual speech filtering, maskee to audible noise, distorsion, musical noise

A new analysis of perceptual speech enhancement is presented. It focuses on the fact that if only noise above the masking threshold is filtered, then noise below the masking threshold, but above the absolute threshold of hearing, can become audible after the masker filtering. This particular drawback of some perceptual filters, hereafter called the maskee-to-audible-noise (MAN) phenomenon, favours the emergence of isolated tonals that increase musical noise. Two filtering techniques that avoid or correct the MAN phenomenon are proposed to effectively suppress background noise without introducing much distortion. Experimental results, including objective and subjective measurements, show that these techniques improve the enhanced speech quality and the gain they bring emphasizes the importance of the MAN phenomenon.

I. INTRODUCTION

M ANY perceptual speech enhancement algorithms have been proposed in the last few decades. They are based on psychoacoustic models to take advantage of the masking phenomenon inherent to the human auditory system. For instance, in [START_REF] Virag | Single channel speech enhancement based on masking properties of the human auditory system[END_REF], the psychoacoustic model is used to control the parameters of the spectral subtraction in order to find the best trade-off between noise reduction and speech distortion; to make musical noise inaudible, the linear estimator proposed in [START_REF] Hu | Incorporating a psychoacoustic model in frequency domain speech enhancement[END_REF] incorporates the masking properties of the human auditory system. In general, the objective of perceptual speech enhancement is to improve the perceptual quality of the enhanced speech signal. Since human ears cannot perceive noise with level below the noise masking threshold, perceptual methods basically aim at reducing audible noise only. By so proceeding, these methods reduce speech distortion.

Although perceptual methods perform well in comparison with classical subtractive type algorithms, most of them still return some audible and annoying musical residual noise. The reasons are manifold. To begin with, biases are introduced by estimating the noise spectrum and the masking threshold. However, some experiments show that even when the noise spectrum and the masking threshold are known, musical noise is still present after denoising. The cornerstone of the letter is the following claim: the attenuation of speech components after perceptual speech enhancement lowers the masking threshold level and, therefore, may reveal noise components initially masked and not processed. This is hereafter called the maskee-to-audible-noise (MAN) phenomenon.

The first contribution of this letter is to highlight the relevance of the MAN phenomenon. The second contribution is the presentation of an elementary approach that takes into account this phenomenon to perform perceptual speech enhancement. This elementary approach involves weighting a standard perceptual filter.

Two weighting functions are considered for application to the same perceptual filter, chosen for its efficiency to process audible noise [START_REF] Lin | Speech denoising using perceptual modification of wiener filtering[END_REF], [START_REF] Amehraye | Perceptual improvement of wiener filtering[END_REF]. This leads to two weighted perceptual filters (WPF). The first one, WPF1, has been introduced in [START_REF] Amehraye | Perceptual improvement of wiener filtering[END_REF]. It is recalled here for comparison to the second one, WPF2, proposed below. In contrast to the former, the latter corrects the MAN phenomenon only at the specific frequencies that are candidates to the MAN phenomenon. By so proceeding, we aim to avoid introducing undesirable distortion.

It is worth noticing that WPF1 and WPF2 basically derive from the standard Wiener filter for two reasons. First, the Wiener filter is easy to implement and second, it can reasonably be expected that if we succeed in reducing the perception of residual noise resulting from some Wiener-like filtering, the quality of the denoised speech will be improved and yield a fairly satisfactory comfort of listening.

The organisation of the letter is as follows. Section II discusses the MAN phenomenon in perceptual speech enhancement. Section III presents the perceptual filter to which the weighting functions are applied and completes the presentation of WPF1 by the results of subjective tests. These experimental results motivate the design of WPF2, introduced in section IV. The two filters WPF1 and WPF2 are then compared by means of subjective and objective tests. Section V concludes this letter.

II. THE MAN PHENOMENON IN PERCEPTUAL SPEECH ENHANCEMENT

A. Frequency masking

The masking phenomenon derives from the limited frequency selectivity of the human auditory system. In this letter, we consider only the so-called frequency masking. This masking occurs when some powerful signal distorts the absolute threshold of hearing and, thus, makes inaudible weak signals that would be perceptible otherwise. How effective the masker is at increasing the masking threshold of hearing depends on the frequency of the maskee and on the frequency of the masker. The maximum masking effect occurs when the masker and the maskee are at the same frequency; the masking effect diminishes when the frequency of the maskee moves away from that of the masker.

B. MAN phenomenon

The main idea of perceptual speech enhancement is to incorporate the masking properties of the human auditory system to reduce audible noise only and, thus, avoid much distortion. Noise components that are not audible because of some maskers in the original noisy signal are still present after denoising. They can become audible if they are initially above the absolute threshold of hearing and their maskers are filtered. This is what we call the MAN phenomenon. It can affect the performance of perceptual filtering that processes audible noise only.

f 0 f 0 f 1 f 1 f 2 f 2

Fig. 1. Description of the MAN phenomenon

Figure 1 is an illustration of the MAN phenomenon located in the dashed area. At frequency f 0 a masker is present. The resulting masking threshold (dashed black curve) is below a first noise component at frequency f 1 and above another noise component, the maskee, located at frequency f 2 . We do not perceive the maskee at frequency f 2 . After filtering the additive noise around the audible frequency f 1 , the masker signal is attenuated, which normally lowers the masking threshold, so that noise at frequency f 2 becomes audible. This phenomenon may occur at each frequency where the energy of the noise maskee lies between the initial masking threshold and the absolute threshold of hearing ATH. These conditions can be easily satisfied physically and carrying out an experiment to illustrate the matter is simple. Consequently, the MAN phenomenon favours the emergence of isolated audible tones that contribute to musical noise. In fact, musical noise consists of rapidly changing random tones that are noticeable in the background of speech. They sound metallic and the denoised speech might be even more unpleasant than the original noisy signal.

III. PERCEPTUAL WEIGHTING BY WIENER FILTERING

A. Principle

We begin by presenting the perceptual filter on which WPF1 and WPF2 are based. Then, we recall the expression of WPF1. The notations introduced henceforth are kept throughout the rest of the paper with always the same meaning.

Let us assume that each frame of noisy signal contains the same number M of samples. Given frame k, Y k (ν) = S k (ν) + N k (ν) denotes the Discrete Fourier Transform (DFT) Y k of the noisy speech signal at frequency ν = 0, 1, . . . , M -1 where S k (resp. N k ) stands for the DFT of the speech signal (resp. additive and independent noise). The weighting functions proposed in this letter are applied to

G k (ν) = | S k (ν)| 2 | S k (ν)| 2 + max (γ k (ν) -T k (ν), 0) . ( 1
)
where

S k (ν) = Y k (ν)W k (ν) = Y k (ν)(ξ k (ν)/(1 + ξ k (ν))), W k is the Wiener filter, ξ k (ν) is the so-called decision-directed estimate of the a priori SNR [4], T k (ν)
is the masking threshold and γ k (ν) is the noise power spectrum estimate.

Filter G k performs a Wiener filtering of only the amount of noise that exceeds the masking threshold. In [START_REF] Lin | Speech denoising using perceptual modification of wiener filtering[END_REF], this method applies to the sub-band components returned by an auditory filter bank. Here, the filter G k results from the adaptation of this method to the usual case where the time-frequency analysis is performed by the standard DFT. The choice of this filter is motivated by the objective test results given in [START_REF] Amehraye | Perceptual improvement of wiener filtering[END_REF]. According to these results, G k outperforms the perceptual filters proposed in [START_REF] Beaugeant | New optimal filtering approaches for hands-free telecommunication terminals[END_REF], [START_REF] Hu | Incorporating a psychoacoustic model in frequency domain speech enhancement[END_REF] and the standard Wiener filter. Filter G k , as well as those in [START_REF] Virag | Single channel speech enhancement based on masking properties of the human auditory system[END_REF], [START_REF] Beaugeant | New optimal filtering approaches for hands-free telecommunication terminals[END_REF], [START_REF] Hu | Incorporating a psychoacoustic model in frequency domain speech enhancement[END_REF], are typical examples of perceptual filters which do not take into account the MAN phenomenon. In fact, they process audible noise only.

The weighting function considered in this section is the standard Wiener filter W k . The resulting WPF, namely WPF1, is thus specified by

G WPF1 k (ν) = W k (ν)G k (ν). (2)
The resulting filtering accentuates the denoising for the frequencies ν where noise is perceptually annoying, that is, when the noise power spectral density or spectrum γ k (ν) is above the masking threshold. By attenuating every frequency with W k , we avoid that the MAN phenomenon occurs.

In the next section, we experimentally assess the performance of WPF1 in comparison to the standard Wiener filter W k , G k (see Eq. ( 1)) and

H k (ν) = W k (ν) if γ k (ν) > T k (ν) 1 otherwise, (3)
which is the perceptually motivated filter of [START_REF] Beaugeant | New optimal filtering approaches for hands-free telecommunication terminals[END_REF]. This filter performs the Wiener denoising only for audible noise frequency components. A smoothing correlogram is applied to each perceptual filter so as to avoid discontinuities of the filter frequency response.

B. Subjective tests

Our purpose is now to assess the enhanced speech quality achieved by WPF1. If the MAN phenomenon is significant, we should notice some improvement by using WPF1 in comparison with filters that do not take into account the MAN phenomenon . Therefore, after the evaluation by means of objective tests presented in [START_REF] Amehraye | Perceptual improvement of wiener filtering[END_REF], we carried out subjective listening tests to compare WPF1, G k , H k and the Wiener filter.

Noise signals from the Noisex database (babble and car noise) were added with two SNRs (5 dB, 10 dB) to 10 sentences randomly chosen from the TIdigits database downsampled to 8 KHz. The experimental protocol was the following one. Short-time windows (32 ms) of noisy speech, with 50% overlap, were transformed into the frequency domain using the short-time Fast Fourier Transform (FFT). The auditory masking threshold was estimated by means of the Johnston model [START_REF] Johnston | Transform coding of audio signals using perceptual noise criteria[END_REF] applied to the Wiener estimate S k . The noise spectrum was assumed to be known in order to assess the filtering without taking the risk to introduce any bias due to noise spectrum estimation. The enhanced speech signal in the time domain was obtained using the overlap-and-add approach after transformation back into the time domain via the Short-Time Inverse FFT. We used the recently standardised methodology [START_REF] Itu-T | Subjective test methodology for evaluating speech communication systems that include noise suppression algorithm[END_REF] for the subjective evaluation. In short, this methodology requires the listener to rate, using five-point scales, the distortion of speech alone (SIG), the background noise alone (BACK) and the overall quality (OVRL). The sentences were presented to 24 listeners amongst which 8 English natives.

The experimental results are those of tables I and II. The perceptual weighting WPF1 performs generally better than the other filters. An analysis of variance (ANOVA) on the data of tables I and II with level of significance α = 0.05 shows the following. WPF1 achieves a statistically significant smaller noise distortion (higher BACK scores) than the other algorithms in every experimental condition (p-values< 0.05 for BACK). This is natural since WPF1 basically accentuates the denoising to reduce the background noise distortion (residual noise) and avoid the MAN phenomenon. As far as the signal quality is concerned, the results are not statistically The fact that WPF1 does not perform statistically better than the other filters with respect to signal quality can be explained as follows: the Wiener filtering alters speech quality and the perceptual filtering does not correct this distortion even though it reduces the residual noise intrusiveness. Now, regarding the overall quality, WPF1 performs significantly better in every situation except for babble noise at 5 dB and for non-native listeners.

According to these experimental results, the MAN phenomenon has a significant impact on the speech enhancement performance since such a simple method as WPF1 makes it possible to avoid it and to generally yield better performance than standard perceptual filters. In the next section, another perceptual weighting is proposed to overcome the limitation of WPF1 pointed out by the subjective ratings above.

IV. ATTENUATION AT FREQUENCY CANDIDATES TO THE MAN PHENOMENON

A. Principle

As noticed in the previous section, the main drawback of WPF1 is the following one: the Wiener filtering is applied to each frequency and this entails some signal distortion that cannot be compensated by the perceptual filtering. Instead of trying to avoid the MAN phenomenon, the new perceptual weighting WPF2 presented now intends to correct it. More specifically, the perceptual filtering of the noisy signal is now followed by a modified Wiener filter. This filter is a post-processing of the perceptual denoised signal only at frequencies ν where the value of the noise spectrum γ k (ν) lies between the absolute threshold of hearing ATH k (ν)and the masking threshold T k (ν). Because noise at these frequencies is candidate to the MAN phenomenon. The expression of the resulting filter where G k is given by Eq. ( 1) and

G WPF2 k (ν) is then G WPF2 k (ν) = G k (ν)W k (ν) (4)
W k (ν) = W k (ν) , if ATH k (ν) < γ k (ν) ≤ T k (ν) 1 , otherwise. (5) When G k (ν) is inactive (G k (ν) = 1), the weighting factor W k (ν) = W k (ν) attenuates noise within the MAN area (ATH k (ν) < γ k (ν) ≤ T k (ν)).

B. Objective and subjective tests

This section aims to compare WPF1 and WPF2 on the basis of objective and subjective tests. We assessed WPF2 with respect to WPF1 only because the latter performs better than the other methods tested above. Also, this allowed us to reduce the test burden for listeners. Only non-native English listeners were involved in these tests. Indeed, according to the results of table I and II, there was no significant difference between scores of English natives and scores of non-native English listeners. This is due to the fact that TIdigits database involves English digits only.

The objective criteria were the standard Segmental Signal to Noise Ratio (SSNR), the Modified Bark Spectral Distortion (MBSD) and the Perceptual Evaluation of Speech Quality (PESQ). The experiments were carried out by using 250 sentences, randomly chosen from the TIDIGITS database. These sentences were corrupted by additive babble and car noise with SNR ranging from -5dB to 20dB. According to table III, WPF2 achieves significant performance improvement in comparison with WPF1, whatever the criterion. On the other hand, the subjective tests carried out according to the same protocol as in section III-B show that generally the listeners better rate WPF1 than WPF2 (see table IV). However, it follows from a t-test analysis that the performance measurements in table IV are not statistically significant (p-values > 0.05). These results confirm the impact of the MAN phenomenon and the relevance of correcting it.

V. CONCLUSION

The goal of this letter was to describe and emphasise the importance of the MAN phenomenon in perceptual speech enhancement. In fact, the weighted perceptual filters proposed in this paper and which take into account this phenomenon perform better than standard perceptual ones that process audible noise only. The approach proposed in this letter can be generalised as follows: first, by extending it to other perceptual filters; second, by looking for more sophisticated weighting factors.

TABLE I THE

 I MEAN SCORES ON THE SIG, BACK AND OVRL SCALES FOR 4 METHODS EVALUATED IN CAR AND BABBLE NOISE AT SNR LEVELS OF 5DB AND 10DB (ENGLISH NATIVE LISTENERS)

	Car noise	Noisy	W k	G k	H k	WPF1
		SIG	4.53	4.5	4.53	4.39	4.61
	5 dB	BACK	2.44	3.96	3.90	3.96	4.54
		OVRL	3.09	4.28	3.88	3.85	4.39
		SIG	4.66	4.64	4.61	4.64	4.64
	10 dB	BACK	2.99	4.59	4.33	4.26	4.65
		OVRL	3.46	4.38	4.28	4.18	4.40
	Babble noise	Noisy	W k	G k	H k	WPF1
		SIG	3.96	3.08	2.31	2.31	3.01
	5 dB	BACK	1.41	3.33	2.88	2.79	3.88
		OVRL	2.09	2.70	1.99	1.81	2.71
		SIG	4.19	3.85	2.94	2.89	3.78
	10 dB	BACK	1.73	3.68	3.20	3.38	4.26
		OVRL	2.35	3.26	2.46	2.44	3.45

TABLE II THE

 II MEAN SCORES ON THE SIG, BACK AND OVRL SCALES FOR 4 METHODS EVALUATED IN CAR AND BABBLE NOISE AT SNR LEVELS OF 5DB AND 10DB (NON-NATIVE ENGLISH LISTENERS) different in car noise and the different methods are as good as each other (p-values > 0.05 for SIG); in contrast, in babble noise, the listeners statistically prefer the original noisy signal (p-values < 0.05 for SIG).

	Car noise	Noisy	W k	G k	H k	WPF1
		SIG	4.43	4.42	4.35	4.23	4.46
	5 dB	BACK	2.52	4.49	3.83	3.84	4.63
		OVRL	3.32	4.29	3.88	3.82	4.41
		SIG	4.40	4.47	4.30	4.42	4.49
	10 dB	BACK	2.78	4.50	4.18	4.20	4.54
		OVRL	3.41	4.30	4.08	4.22	4.37
	Babble noise	Noisy	W k	G k	H k	WPF1
		SIG	4.13	2.98	2.26	2.26	2.94
	5 dB	BACK	1.87	3.49	3.29	3.23	3.78
		OVRL	2.65	2.86	2.11	2.16	2.79
		SIG	4.39	3.54	2.91	2.86	3.51
	10 dB	BACK	2.08	3.83	3.42	3.52	3.98
		OVRL	2.92	3.24	2.63	2.61	3.34

TABLE III COMPARATIVE

 III PERFORMANCE MEASUREMENTS BETWEEN WPF1 AND WPF2, IN TERMS OF MEAN MBSD, PESQ AND SEGMENTAL SNR OVER 250 SENTENCES CORRUPTED BY BABBLE AND CAR NOISE AT DIFFERENT SNR LEVELS.

	Babble noise	-5dB	0dB	5dB	10dB	15dB
	MBSD	WPF1	0.054	0.045	0.037	0.031	0.027
		WPF2	0.038	0.033	0.030	0.027	0.024
	PESQ	WPF1	1.75	2.18	2.63	3.00	3.30
		WPF2	2.10	2.46	2.80	3.12	3.41
	SSNR	WPF1	1.15	2.44	4.29	6.45	8.80
		WPF2	1.93	3.54	5.50	7.60	9.89
	Car noise	-5dB	0dB	5dB	10dB	15dB
	MBSD	WPF1	0.025	0.022	0.020	0.016	0.013
		WPF2	0.023	0.020	0.017	0.013	0
	PESQ	WPF1	3.42	3.65	3.86	4.05	4.21
		WPF2	3.45	3.69	3.91	4.10	4.26
	SSNR	WPF1	6.50	8.47	10.41	12.25	14.13
		WPF2	7.96	10.13	12.29	14.5	16.61

TABLE IV MEAN

 IV SCORES FOR THE SIG, BACK, AND OVRL SCALES FOR WPF1 AND WPF2 EVALUATED IN CAR AND BABBLE NOISE AT SNR LEVELS OF 5DB AND 10DB WITH NON-NATIVE ENGLISH LISTENERS

	Car noise	WPF1	WPF2
		SIG	4.35	4.24
	5dB	BACK	4.08	3.86
		OVRL	4.01	3.89
		SIG	4.37	4.37
	10dB	BACK	4.21	4.03
		OVRL	4.12	4.02
	Babble noise	WPF1	WPF2
		SIG	2.81	2.654
	5dB	BACK	3.22	2.554
		OVRL	2.56	2.27
		SIG	3.46	3.49
	10dB	BACK	3.51	3.03
		OVRL	3.23	2.99


Julien Graffion

Dounia Dems

Mesut Demirelli

Thibaud Coradin

Nicolas Delsuc

Carole Aimé
email: carole.aime@upmc.frwww.labos.upmc.fr

An All-in-One Molecule for the One-Step Synthesis of Functional Hybrid Silica Particles with Tunable Sizes

Keywords: Nanoparticles, Self-assembly, Supramolecular chemistry, Alkoxysilanes, Hybrid materials

Spherical particles with well-defined diameters were obtained by self-assembly of trityl-based molecules. Thanks to the robustness of the organic scaffold, a variety of modifications could be covalently introduced into the network so as to stabilize the supramolecular structure by a sol-gel route. Using supramolecular chemistry, we showed that the synthesis of hybrid small molecules allowed engineering nanomaterials with tun-

Introduction

Antisolvent precipitation is an attractive technique for the elaboration of nanoparticles (NPs) and finds applications in various fields such as electro-optics, catalysis, and nanomedicine, notably with the elaboration of imaging agents and for improving drug formulation. Alternatively, the design and synthesis of nanoscale metal-organic frameworks have provided tunable nanoobjects to serve as functional platforms. [1][2][3][4][5][6] However, despite the size tunability of the resulting particles, their high interfacial area renders them prone to aggregation and Oswald ripening, which affect their stability. [7] Research efforts have been devoted to the selection of stabilizers to be added to the precipitating nanoparticles. [8] To this aim, amphiphilic diblock copolymers have been used, typically a hydrophilic poly(ethylene glycol) (PEG) block chemically linked to a hydrophobic block. [9] Thanks to its amphiphilic properties, the polymer adsorbs at the surface of the NPs, which limits further growth or dissolution and, thus, prevents aggregation. Polyelectrolytes such as polylysine, poly(ethylene imine), and chitosan have also been used to stabilize the NPs electrostatically and sterically and to ensure a narrow size distribution. [10] Alternatively, hybrid (organic-inorganic) materials have been intensively investigated for two decades owing to their unique properties that able size and functionality. The use of a combination of different characterization techniques, including dynamic light scattering, cryoTEM, and solid-state NMR spectroscopy, provided careful understanding of the relationship between the molecular and supramolecular structures for further chemical engineering of supramolecular hybrid materials. combine the hardness and stability of the inorganic phase with the flexibility and mild processing of polymer materials. [11,12] As recently reviewed, [13] silsesquioxane NPs have been synthesized from organotrialkoxysilane precursors by different routes, including hydrolytic [14][15][16][17] and nonhydrolytic [18] sol-gel routes and emulsion polymerization. [19][20][21][22] With such processes, spherical, [14,15,[23][24][25] hollow, [26][27][28][29][30] and Janus [31] NPs have been obtained. Their surface can be easily functionalized through silanization reactions of superficial silanol groups. [32,33] This offers a large number of possibilities to bind targeting moieties for imaging or drug-delivery systems [34] or to enhance cell compatibility and internalization. [35,36] Among polysilsesquioxanes, phenylsilsesquioxanes have drawn much attention, and the synthesis of PhSiO 3/2 particles has been reported. [16,37] In another study, Dirè et al. synthesized NPs from a series of alkyltriethoxysilane (methyl-, ethyl-, vinyl-, phenyl-, amyl-, and octyltriethoxysilane) with systematic assessment of the influence of the different organic groups on the synthesis and properties of the silsesquioxane NPs. [25] Such molecular design allowed the doping of particles, either by physical entrapment or covalent linking, with representative examples using fluorophores such as Nile red, [38] rhodamine B, [17,24,39] and rhodamine 6G. [25] Notably, for the abovecited alkyltriethoxysilanes, the formation of nanoparticles results from nucleation and growth induced by sol-gel polymerization, which is typically followed by a two-step hydrolysiscondensation process, as in the works of Dirè et al. On the other hand, the self-assembly of silsesquioxanes has been reported for the design of helical and lamellar structures with organogelator-like behavior. [40][41][42] In this work, we exploit the supramolecular interaction properties of the organic moiety of organotrialkoxysilanes as the only molecular key towards nanoparticles. We use different molecular designs based on the triphenylmethane (trityl) group (Figure 1) to take advantage of this moiety to drive the forma- tion of well-defined objects. The addition of a trialkoxysilane group (see compound 2, Figure 1) aims at allowing their stabilization through sol-gel chemistry by preventing aggregation and Oswald ripening. Finally, to investigate the scope of the tolerated modifications that preserve the structure of the objects, we introduced a chiral multifunctional backbone, that is, the amino acid L-cysteine. Its thiol group was functionalized with the trityl group, as in compounds 1 and 2, and the trialkoxysilane was introduced by dicyclohexylcarbodiimide (DCC)/ N-hydroxysuccinimide (NHS) coupling with (3-aminopropyl)triethoxysilane (APTES). A fluorescent anthracene moiety was conjugated to the amine function of the cysteine. This provided cues for functional implementation of the resulting supramolecular structures in an "all-in-one molecular structure" (see compound 3, Figure 1). The three molecules were synthesized as described in Schemes S1 and S2 (Supporting Information) and were obtained in high purity (> 95 %, no detectable impurities by 1 H NMR spectroscopy).

Results and Discussion

Antisolvent precipitation in an ethanol/water mixture was used to obtain supramolecular objects from 1. Precipitation was monitored by dynamic light scattering (DLS) in the range of 60 to 95 % water and led to well-defined objects with sizes from 210 to about 100 nm (Figure 2a). No supramolecular structure could be detected below this range owing to the solubility of 1 in ethanol, whereas macroscopic objects sedimented above 95 % water. Very interestingly, similar trends could be observed by precipitating 2 and 3, and a constant decrease in size from 60 to 80 % water gave rise to well-defined objects with a mean hydrodynamic diameter above 200 nm at 60 % (210 nm for 1, 260 nm for 3) and leveling around 90 nm from 80 to 95 % water. Notably, as the water content increased, the size distribution became narrower (Table S1). Playing with the tritylated molecule concentration at a given water content could also be used to tune the size of the resulting supramolecular structures. At 90 % water, upon increasing the alkoxysilane concentration from 100 to 400 μM, the hydrodynamic diameter of 2 steadily varied from 95 to 125 nm and that of 3 varied from 80 to 100 nm (Figure 2b). CryoTEM observations showed the formation of well-defined spherical objects in solution (Figure 2b, inset). Overall, the precipitation of 1 and the similar trends observed for 2 and 3 showed that the strong hydrophobic character and the nonplanar structure (one sp3 carbon atom) of the trityl group favored π-stacking interactions over a short range, and therefore, the trityl group is particularly interesting for antisolvent precipitation to obtain well-defined spherical objects. Given that trityl assembly is the driving force for the formation of spherical particles, its scaffolding role is strong enough to ensure particle formation upon the addition of even large organic molecules. Interestingly, neither the alkoxysilane nor the chiral center and anthracene moiety seemed to play an important role in obtaining well-defined nanostructures, as 1, 2, and 3 behaved similarly. In the presence of a chiral center, the preservation of a chiral environment within the hybrid particle could be ascertained by circular dichroism spectropolarimetry (Figure S3). Because this supramolecular chirality is expressed without modifying any of the previously described characteristics of the nanomaterials, this process is very tolerant to a wide range of modifications.

The introduction of the alkoxysilane group may offer a way to stabilize the supramolecular particles of 2 and 3. Indeed, the size of objects from 1 increased progressively up to (633 ± 40) nm after 1 month (Figure 3a). In the case of precursors 2 and 3, after nanoparticle formation, a sol-gel route was used by placing the samples under ammonia vapors for 24 h to catalyze their hydrolysis-condensation under basic conditions. For example, in the case of 3, after this treatment, the size of the spheres leveled out at 127 ± 3 nm after 1 day, whereas in the absence of ammonia, stable spherical structures exhibiting a size of 170 ± 7 nm could be obtained after 3 days. This showed, one, the role played by the alkoxysilane moiety in stabilizing the final structures and, two, that ammonia was not mandatory to stabilize 3 but that it sped up the kinetics of the sol-gel process. TEM observations after stabilization of 2 and 3 revealed that the spherical structures were preserved (Figure 3b). In terms of chemical structure, the hydrolysis-condensation of the alkoxysilane and the formation of the siloxane shell at the surface were ascertained by ζ-potential measurements and solid-state 29 Si NMR spectroscopy. Stabilized particles from 2 and 3 showed negative ζ potentials, which indicated the presence of silanolate groups, resulting from hydrolysis of the ethoxysilane groups, on their surface (Table S2). Probing the silicon chemical environment of 3 after stabilization by 29 Si NMR spectroscopy revealed broad signals around -45.4, -57.1, and -66.8 ppm that could be attributed to R-Si(OEt) 3 , R-Si(OSi) 2 (OH) (T2), and R-Si(OSi) 3 (T3), respectively; these signals confirmed partial alkoxysilane condensation (Figure S1). Further investigations by FTIR spectroscopy showed that the Si-O deformation region (900-1145 cm -1 ) did not change upon comparing 3 as the molecular precursor and as stabilized particles (Figure S2). In addition, no condensed siloxane bond [ν as (SiOSi)] could be identified owing to their low concentration, which was restricted to the surface. These results provided evidence that the precursor remained intact within the nanoparticles owing to the hydrophobic character of the core, which prevented internal diffusion of water and ammonia. Further investigations by using powder X-ray diffraction (PXRD) showed the presence of amorphous silica, which could be expected given the huge number of molecules involved in the formation of the 100 nm diameter particles (Figure S4).

As a means of imparting the particle with functionality, a fluorescent probe, namely, an anthracene moiety, was introduced to the molecular design (see compound 3). Covalent conjugation of the anthracene probe prevented leaching of the fluorophore. Figure 4a shows the emission spectra of 3 upon increasing the water content up to 90 %. From 0 to 50 % water, emission of the anthracene moiety (well-defined components, maximum intensity at 407 nm) indicates that the anthracene moiety exists as a monomer in solution. From 60 to 90 % water, the emission spectra drastically broaden and become redshifted (maximum intensity at ≈ 470 nm). These changes can be attributed to the transition of soluble 3 into its precipitated state, which is in good agreement with the DLS measurements, by which no supramolecular structure could be detected below 60 % water. As expected, the excitation spectrum was also slightly redshifted (7 nm) upon increasing the water content (Figure S5a).

Given that 3 was shown to be an efficient multifunctional scaffold to form hybrid nanoparticles with easily tunable diameters, we investigated the possibility to increase the functionality of the hybrid nanoparticles by physical encapsulation. As proof of concept, the encapsulation ability of the hybrid particles was assessed by the addition of rhodamine B (rhod-B) during the precipitation of 3. Rhod-B was chosen to obtain insight into the impact of adding sizeable molecules on the stability and size of the hybrid particles. In addition, rhod-B was selected because it provides a straightforward read out of the encapsulation by using Förster resonance energy transfer (FRET) with the anthra- cene moiety. Different concentrations of rhod-B were added during precipitation of 3 (200 μM, 90 % water). After stabilization and removal of nonencapsulated dye upon centrifugation, a pink precipitate consisting of hybrid particles incorporating rhod-B was recovered and further redispersed in water (Figure 4b, inset). DLS measurements showed an increase in size from 90 nm to ca. 160 nm upon rhod-B encapsulation, whereas no variation in the ζ potential could be detected; this was indicative of dye encapsulation rather than surface adsorption (Tables S3 andS4). A photoluminescence study was performed to provide additional evidence for the encapsulation of rhod-B in stabilized hybrid particles of 3. The emission spectra recorded at an excitation wavelength in the anthracene region (378 nm) allowed detection of the emission of both the anthracene (centered at 460 nm) and rhod-B (at 580 nm) moieties through FRET (Figure 4b). The excitation spectra recorded at the maximum intensity of rhod-B emission (λ em = 587 nm) also showed the contribution of the anthracene moiety in 3; this confirmed the FRET phenomenon between rhod-B and 3 and further evidenced the encapsulation of rhod-B (Figure S5b). Noteworthy, the emission intensity of rhod-B increased with its initial concentration, whereas the emission intensity of 3 remained constant (Figure 4b). The ratio of the areas of the emission bands of rhod-B and 3 rapidly increased with the rhod-B concentration before it started to level off around 2 mM, which may have indicated that saturation was approaching (Figure S6). The saturation regime could not be investigated further owing to destruction of the hybrid particles at higher rhod-B concentrations. Importantly, up to 2 mM rhod-B, the integrity of the particles could be preserved (Figure S7). These results show the robustness of the molecular design, for which stable spherical nanoparticles could still be obtained even after the addition of rhod-B.

Conclusion

Spherical nanoparticles with well-defined tunable diameters were obtained in a single step by antisolvent precipitation. Our "all-in-one molecule" was based on an organic moiety, a trityl group, which alone was responsible for nanoparticle formation.

The addition of an alkoxysilane group allowed stabilization of the hybrid particles by hydrolysis and condensation so that the final system could be frozen, and further shape and size evolution could be prevented. We showed that the propensity of the trityl group to generate supramolecular spheres provided this molecular design with a unique robustness that allowed the covalent introduction of large organic groups without changing the structural properties of the particles. This opens the door towards a broad range of functional implementations, as illustrated by the covalent conjugation of fluorophores. In addition, such nanoparticles allowed the encapsulation of large molecules without modifying their structure, which indeed opens up additional opportunities. Finally, we have to keep in mind that the functionality of those nanoparticles could be even broader given the versatility of silane chemistry for further functionalization of the nanoparticles' surface, in particular upon envisioning applications in biotechnologies and nanomedicine.

Experimental Section

Preparation of Hybrid Nanoparticles 1, 2, and 3: A solution of 1, 2, or 3 in ethanol was quickly added to freshly deionized water under vigorous stirring (always the same speed) at room temperature for 30 s. For each precursor, the water/ethanol ratio (from 0 to 95 % water) and concentration (100, 200, and 400 μM) were varied, and the final volume was always kept at 1.5 mL for reproducibility reasons. For nonstabilized systems, the suspension was left without further stirring. For stabilized systems, the suspension was placed in a sealed environment for 24 h in the presence of ammonia vapors, which was followed by aging for 48 h in the newly basified suspension through ammonia transfer. To recover the stabilized hybrid nanoparticles, the suspension was centrifuged and washed with a mixture of ethanol/deionized water (1:1) and then deionized water (4×).

Preparation of Hybrid Nanoparticles 3 Containing Rhodamine B:

A solution of 3 in ethanol was quickly added to a solution of freshly deionized water (90 % water content in all cases) containing 0.2 mg (42 μM), 0.5 mg (104 μM), 1 mg (208 μM), 2 mg (417 μM), 5 mg (834 μM), or 10 mg (2080 μM) of rhodamine B. The conditions of stirring, temperature, volume, stabilization, and washing were the same as those for the preparation of hybrid nanoparticles 3 (see above). Although in this case, the final concentration of 3 was kept at 200 μM only.

Dynamic Light Scattering (DLS) and ζ-Potential Measurements:

ζ-potential and size-distribution measurements were performed by using a Malvern Zetasizer Nano ZS90 instrument. The autocorrelation functions were recorded at a scattering angle θ = 90°and were analyzed by the non-negatively constrained least-squares technique (NNLS-Multiple Pass) for the determination of particle diameter. All the suspensions were sonicated for 5 s prior to measurement at 20 °C. For DLS measurements, all measurements were run in triplicate to obtain the size standard deviation.

Transmission Electron Microscopy:

A drop of sample in aqueous solution (7 μL) was deposited on carbon-coated copper grids (400 mesh, Electron Microscopy Sciences). After 3 min, the excess liquid was blotted with filter paper (Whatman #4). TEM was performed at room temperature by using a FEI Technai Spirit G2 operating at 120 kV. Images were recorded with a Gatan Orius CCD camera.

Fluorescence Spectroscopy: Fluorescence measurements were recorded with a Fluoromax-3 (Horiba Jobin Yvon Perkin-Elmer LS55) spectrometer, with a 150 W ozone-free xenon arc-lamp as excitation source. Emission spectra (350-650 nm for 3, and 390-725 nm for rhodamine B @ 3) were recorded with excitation wavelengths of 375 and 378 nm in the absence and presence of rhodamine B, respectively. The bandpass of the excitation and emission monochromators was set at 5 nm, and the scan speed was set at 10 nm s -1 . Samples were investigated at 25 °C in 1 cm path-length cells.

Figure 1 .

 1 Figure 1. Molecular structures synthesized from an alkyl-S-trityl moiety (see compound 1), in combination with a triethoxysilane group (see compound 2), and belonging to a cysteine functionalized with an anthracene (see compound 3).

Figure 2 .

 2 Figure 2. DLS of the supramolecular structures of 1, 2, and 3 as a function of (a) water content (100 μM), and (b) concentration (90 % water). The inset shows the cryoTEM image of 3 (200 μM, 90 % water).

Figure 3 .

 3 Figure 3. (a) DLS of the supramolecular structures obtained from 1 and 3 with (stabilized) or without (nonstabilized) ammonia vapors and (b) TEM images of the stabilized spheres obtained from 2 (b1) and 3 (b2) (100 μM, 90 % water).

Figure 4 .

 4 Figure 4. Emission spectra of 3 (a) as a function of the water content (100 μM, λ ex = 375 nm) and (b) after the encapsulation of rhod-B at different concentrations ([3] = 200 μM, 90 % water, λ ex = 378 nm).

Acknowledgments

The authors thank Clément Sanchez, François Ribot, and Cédric Boissière for fruitful discussions, Gervaise Mosser for cryoTEM observations and careful reading of the manuscript, Patrick le Griel for his help in electron microscopy, Cristina Coelho-Diogo for solid-state NMR spectroscopy measurements, Mohamed Selmane for PXRD experiments, and the Centre National de la Recherche Scientifique (CNRS) and the Idex PSL (ANR-10-IDEX-0001-02 PSL) for funding.