Catherine Doss
email: doss@ann.jussieu.fr

Michèle Thieullen
email: michele.thieullen@upmc.fr

Oscillations and Random Perturbations of a FitzHugh-Nagumo System

Keywords: FitzHugh-Nagumo system, fast-slow system, excitability, equilibrium points, bifurcation parameter, limit cycle, bistable system, random perturbation, large deviations, metastability, stochastic resonance

We consider a stochastic perturbation of a FitzHugh-Nagumo system. We show that it is possible to generate oscillations for values of parameters which do not allow oscillations for the deterministic system. We also study the appearance of a new equilibrium point and new bifurcation parameters due to the noisy component.

1 Introduction.

Let us consider the following family of deterministic systems indexed by the parameters a ∈ R and δ > 0.

δ ẋt = -y t + f (x t ), X 0 = x (1.1) ẏt = x t -a, Y 0 = y (1.2)
and their stochastic perturbation by a one dimensional Wiener process (W t ) as follows

δdX t = (-Y t + f (X t ))dt + √ ǫdW t , X 0 = x (1.3) dY t = (X t -a)dt, Y 0 = y (1.4)
The function f is a cubic polynomial: f (x) = -x(x -α)(x -β) with α < 0 < β. The parameter δ is small. The deterministic system (1.1)-(1.2) is an example of a slow-fast system: the two variables x, y have different time scales, x t evolves rapidly while y t evolves slowly. This system is one version of the so called FitzHugh-Nagumo system and plays an important role in neuronal modelling. In this context x t denotes the voltage or action potential of the membrane of a single neuron. It was first proposed by FitzHugh and Nagumo (cf. [START_REF] Fitzhugh | Impulses and physiological states in theoretical models of nerve membrane[END_REF], [START_REF] Nagumo | An active pulse transmission line simulating nerve axon[END_REF]). One interest of this model is that it reproduces periodic oscillations observed experimentally. Indeed FitzHugh-Nagumo system finds its origin in the nonlinear oscillator model proposed by van der Pol. It is also a simplification of the Hodgkin-Huxley model which describes the coupled evolution of the membrane potential and the different ionic currents: existence of different time scales enable to pass from a four dimensional model to a two dimensional one. Oscillations can take place because the deterministic system (1.1)-(1.2) exhibits bifurcations; more details will be given in section 3. Let us mention that oscillations in this system (1.1)-(1.2) can only occur when a ∈]a 0 , a 1 [ where a 0 < a 1 are two particular values of parameter a namely the bifurcation parameters.

Our main interest in the present paper is to generate oscillations even for a < a 0 (symmetrically a > a 1 ) )by adding a stochastic perturbation to the deterministic system. What may be interpreted as some resonance type effect (cf. [START_REF] Herrmann | The exit problem for diffusion with time periodic drift and stochastic resonance[END_REF], [START_REF] Gammaitoni | Stochastic Resonance[END_REF]). We will therefore investigate possible oscillations for system (1.3)- (1.4). The presence of parameter ǫ introduces a third scale in the system and the relative strength of δ and ǫ measured by the ratio ǫ| log δ| δ will determine its evolution. Our study was inspired by reference [START_REF] Freidlin | On stable oscillations and Equilibrium Induced by Small Noise[END_REF] where M. Freidlin considers a random perturbation of the second order equation δ d 2 yt dt 2 = g( dyt dt , y t ) and performs the study of its solution using the theory of large deviations (cf. [START_REF] Freidlin | Random perturbation of dynamical systems[END_REF]). See also [START_REF] Freidlin | On stochastic perturbation of dynamical systems with fast and slow components[END_REF] for the study of a more general situation. In our case g( ẏ, y) = y -f ( ẏ + a). Although our argument is close to M. Freidlin's, the presence of parameter a leads to a richer behaviour. We prove the existence of equilibrium point and limit cycles different from the deterministic ones; as in [START_REF] Freidlin | On stable oscillations and Equilibrium Induced by Small Noise[END_REF], as well as a new bifurcation point which did not exist for the deterministic system. Our study relies on transitions between basins of attraction of stable equilibrium points due to noise. Relying on some estimation of a family of exit times(propositions 3.3 and 3.4) we study conditions on the parameters under which a convenient stochastic dynamic system approach its main state(proposition 3.5 ),(which corresponds to the equilibrium point exhibited in main theorem 2.2), or approach a metastable state (which corresponds to the limit cycle and to the new bifurcation parameters exhibited in main theorem 2.1). A general study of slow-fast systems perturbed by noise can be found in [START_REF] Berglund | Noise-Induced Phenomena in Slow-Fast Dynamical Systems: A Sample Paths Approach[END_REF]. Bursting oscillations in which a system alternates periodically between phases of quiescence and phases of repetitive spiking has been studied for stochastically perturbed systems in [START_REF] Hitczenko | Bursting Oscillations induced by small noise[END_REF] and may be studied later in our stochastic setting. We recall that in the deterministic one a bursting-type behaviour has been generated in [START_REF] Doss | it Bursting oscillations in two coupled Fitzhugh-Nagumo system[END_REF].

The paper is organized as follows. In section 2 we recall basic facts about (1.1)-(1.2) and we state the two main theorems (2.1) and (2.2). Section 3 is devoted to the application of large deviation theory to (1.3)-(1.4) and section 4 to the proof of the main theorems.

2 Some Basic Results.

Deterministic FitzHugh-Nagumo System.

In (1.1)-(1.2) let us consider α < 0 < β and f (x) = -x(x -α)(x -β). In order to investigate the asymptotic behaviour of (x t , y t ) one first looks for equilibrium points of the system and their stability. The equilibrium points are defined as the points (x, y) where the right-hand sides of both equations of the system vanish. For any value of a there is therefore a unique equilibrium point for (1.1)-(1.2) which is (a, f (a)). Moreover let a 0 , a 1 , a 0 < a 1 be the two points where f ′ vanishes. The stability of the equilibrium point (a, f (a)) changes when a passes through the value a 0 (resp. a 1 ); a 0 and a 1 are called the bifurcation parameters of the system. Let us focus on a 0 ; an analogous argument holds for a 1 . By linearizing system (1.1) at (a 0 + η, f (a 0 + η)) for η small, we obtain system Ż = AZ with

A = f ′ (a 0 +η) δ -1 δ 1 0
A admits the two eigenvalues

λ ± = 1 2δ (f ′ (a 0 + η) ± i 4δ -f ′ 2 (a 0 + η))
. The sign of f ′ (a 0 + η) is the same as that of η since f ′ is increasing in the neighbourhood of a 0 . The point (a 0 + η, f (a 0 + η)) is therefore an attracting (resp. repulsive) focus when η < 0 (resp. η > 0). In particular, (a, f (a)) is stable when a < a 0 , unstable when a > a 0 . For a ∈]a 0 , a 1 [ the system admits a limit cycle. The bifurcation is of Hopf type [START_REF] Françoise | Oscillations en biologie. Analyse qualitative et modèles[END_REF] It can be verified numerically that if δ < 0.01 the limit cycle is very close to the loop made up with the two attracting branches of the curve y = f (x) where x → f (x) is decreasing and y ∈ [f (a 0 ), f (a 1 )], and the portions of the two horizontal segments y = f (a 0 ), y = f (a 1 ) connecting them. When δ → 0 the period of this limit cycle is 0(1); by example for f (x) = x(4 -x 2 ) it is equal to 2 (cf.

[14]).

Main Theorems

Consider (X t , Y t ) the solution of (1.3)-(1.4) and S > 0 given in definition 3.3. Let us assume that ǫ > 0 and δ > 0 go to zero in such a way that for some constant c > 0.

ǫ|

log δ| δ → c (2.1)
Let us denote by lim * any limit on ǫ and δ going to 0 under condition (2.1).

Theorem 2.1 Let c ∈]0, S[; then: 1.If a ∈]x -(c), x + (c

2, s.t for all

A, h > 0, y ∈]f (a 0 ), f (a 1 )[, lim * P (x,y) ( A 0 |X t -Φ a c (t)| 2 dt > h) = 0 (2.2) lim * P (x,y) (sup [0,A] |Y t -Ψ a c (t)| > h) = 0 (2.3) 2. If a < x -(c) or a > x + (c) then for all y ∈]f (a 0 ), f (a 1 )
[, for all h > 0 there exists t(y, h) such that for all A > t(y, h),

lim * P (x,y) ( sup [ t(y),A] |X t -a| + |Y t -f (a)| > h) = 0 (2.4) Figure 1: Limit cycle when f (x) = x(4 -x 2 ) and ǫ| log δ| δ → c

Remark

In the first case the solution stabilizes when δ → 0 and ǫ| log δ| δ → c around a limit cycle defined by c and different from the one obtained in the deterministic case when δ → 0 and ǫ = 0 (see figure 1). Moreover x -(c) and by symmetry x + (c) play the role of bifurcation parameters for the stochastic FitzHugh-Nagumo system (1.3)- (1.4). Indeed for a in the neighborhood of x -(c) but smaller the limit of (X t , Y t ) is a unique equilibrium point, whereas for a in the neighborhood of x -(c) but greater it is the graph of a periodic function. These bifurcation parameters are different from those of the deterministic system (1.1). This theorem is a new result w.r.t. [START_REF] Freidlin | On stable oscillations and Equilibrium Induced by Small Noise[END_REF]. It is made possible by the freedom on parameter a. The new limit cycle (Φ a c (t), Ψ a c (t)) is defined in the same way as in [START_REF] Freidlin | On stable oscillations and Equilibrium Induced by Small Noise[END_REF] Theorem 1 Part 3., provided we take into account the presence of a in our system. Other regimes are considered in the work of Berglund and Gentz (cf. [START_REF] Berglund | Noise-Induced Phenomena in Slow-Fast Dynamical Systems: A Sample Paths Approach[END_REF] and references therein).

Theorem 2.2 Let c > S. Consider x * -(y * ) and x * + (y * ) defined in proposition 3.2 and definition 3.3; then for all y ∈]f (a 0 ), f (a 1 )[ and for all h > 0 there exists t(y, h) such that for all A > t(y, h),

1. If a ∈]x * -(y * ), x * + (y * )[ lim * P (x,y) ( sup [ t(y),A] |Y t -y * | > h) = 0 (2.5) 2. If a < x * -(y * ) or a > x * + (y * ), lim * P (x,y) ( sup [ t(y),A] (|X t -a| + |Y t -f (a)|) > h) = 0 (2.

6)

Remark: Case 1 may be considered as a degenerate version of the limit cycle of case 1 theorem 2.1. In fact y * is a fixed point but X t oscillates between x * -(y * ) and x * + (y * ).

3 Exit Time, Main State and Metastable State

dx y t = (-y + f (x y t ))dt, x y 0 = x (3.1)
So we are led to consider the real valued deterministic system

dx t = b(x t )dt, x 0 = x (3.2)
and its perturbation by a brownian motion

dx t = b(x t )dt + √ ǫdW t , x0 = x (3.3)
Let us briefly recall some results from [START_REF] Freidlin | Random perturbation of dynamical systems[END_REF]. The process (x t ) describes the movement of a particle on the real line submitted to the force field b(x) and to a stationnary Gaussian noise of amplitude √ ǫ. When ǫ → 0, (x t ) converges to the solution (x t ) of (3.2):

∀η > 0 ∀T > 0 lim ǫ→0 P(sup [0,T ] |x t -x t | > η) = 0 (3.4)
However because of diffusion due to the presence of noise, some trajectories of the process (x t ) may present large deviations from those of the deterministic system (x t ). Such deviations are measured by means of the action functional S T 2 T 1 (ϕ) independant of ǫ and defined by

S T 2 T 1 (ϕ) = 1 2 T 2 T 1 | φu -b(ϕ u )| 2 du (3.5)
when ϕ is absolutely continuous, by

S T 2 T 1 (ϕ) = +∞ otherwise Theorem 3.1 ( [5], Lemma 2.1, Chap. 4) Let η > 0. Then P x (sup [0,T ] |x t -x t | ≥ η) ≤ exp(- 1 ǫ [inf ∆ S T 0 (ϕ) + o(1)]) (3.6)
when ǫ → 0 and where ∆ := {ϕ;

ϕ 0 = x, sup [0,T ] |ϕ t -x t | ≥ η}.
Large deviations theory also provides estimates on the first exit time of (x t ) from a domain (cf. [START_REF] Freidlin | Random perturbation of dynamical systems[END_REF], Theorem 4.2, Chap. 4). Domains of interest are basins of attraction of stable equilibrium points of (3.2). The key tools are quasipotentials.

Definition 3.1 The quasi potential of the deterministic system (3.2) w.r.t. a point x (also called transition rate) is defined as the function

u → V (u) := inf{S T 2 T 1 (ϕ); 0 ≤ T 1 < T 2 , ϕ(T 1 ) = x, ϕ(T 2 ) = u} (3.
=]α 1 , α 2 [ with V (α 1 ) < V (α 2 )
. Then for all x ∈ D and h > 0,

lim ǫ→0 P x (x τ = α 1 ) = 1 (3.9) lim ǫ→0 P x (e V (α 1 )-h ǫ < τ < e V (α 1 )+h ǫ ) = 1 (3.10)

Remark:

With the notations of Theorem 3.2, with great probability when ǫ → 0, the behaviour of the process (x t ) on an interval [0, T (ǫ)], in particular whether the process has jumped out of the basin of attraction D or not, depends on the value of ǫ log T (ǫ) compared to V (α 1 ).

Two Families of Quasipotentials.

Let us now apply these results to the following family of one dimensional dynamical systems indexed by parameter y introduced in the preceeding subsection. V ± (y) = -2

dx y t = (-y + f (x y t ))dt, x y 0 = x (3.11) Proposition 3.2 Let y ∈]f (a 0 ), f (a 1 )
x * 0 (y)

x * ± (y)

(-y + f (u))du.

(3.12)

From Proposition 3.1 we see that V ± (y) is the value at x * 0 (y) of the quasipotential of (3.2) w.r.t. x * ± (y). Both functions V ± (y) are strictly monotone: V -is strictly increasing, V + is strictly decreasing. Therefore their graphs restricted to ]f (a 0 ), f (a 1 )[ intersect at a unique point. For any function U satisfying y -f (u) ≡ -∂ x U /2 the following identities hold:

V -(y) = U(x * 0 (y)) -U(x * -(y)), V + (y) = U(x * 0 (y)) -U(x * + (y)) (3.13)
In our case such a function U is a polynomial of degree 4 which admits x * 0 (y) as relative maximum and x * ± (y) as relative minima. The graph of the function U has two wells with respective bottoms at x * ± (y) and one top at x * 0 (y). Identities (3.13) express that V ± (y) are the respective depths of these wells. Therefore, on E 1 the well with bottom x * -(y) is the deepest one while it is the contrary on E 2 . Second, the portion of the curve z = f (x) connecting (x * -(y), y) to (x * 0 (y), y) is situated below the horizontal line L y := {(x, z); z = y}; thus the positive quantity 1 2 V -(y) measures the surface of the area limited by the curve z = f (x) and the segment of L y on which x ∈ [x * -(y), x * 0 (y)]. In the same way 1 2 V + (y) measures the surface of the area limited by the curve z = f (x) and the segment of L y on which x ∈ [x * 0 (y), x * + (y)] but in this case the portion of the curve is above the line segment. Moreover as we will see in the following section V ± is connected to exit times of diffusions

d Zy t = (-y + f ( Zy t ))dt + √ ǫdW t , Zy 0 = x (3.14)
from the basins of attraction of x * ± (y) (cf.Theorem 3.2).

Exit Times, Main state and Metastable States

We refer the reader to [START_REF] Freidlin | On stochastic perturbation of dynamical systems with fast and slow components[END_REF] for the present section. Let us recall the fundamental difference between the two parameters ǫ and δ. Parameter δ is already present in the deterministic system (1.1)-(1.2) where it measures the difference between the time scale of the slow variable y t and the time scale of the fast variable x t . In particular after the time change s := δt, the trajectory (x s , ỹs ) := (x δt , y δt ) satisfies

ẋs = -ỹ s + f (x s ) (3.15) ẏs = δ(x s -a) (3.16)
Since δ is small the component ỹs may be considered as constant equal to ỹ0 . Let us define Z y t to be the family of solutions of

δdZ y t = (-y + f (Z y t ))dt + √ ǫdW t , Z y 0 = x. (3.17)
We will also consider the time change ( Zy t ) of (Z y t ) under s := δt

d Zy s = (-y + f ( Zy s ))ds + √ ǫd Ws , Zy 0 = x (3.18)
with ǫ = ǫ δ . ( Zy t ) is the stochastic perturbation of (3.11). In the two following propositions we give crucial estimates on exit time of the respective solutions of (3.17) and (3.18).

Proposition 3.3 Let τ y 1 (resp. τ y 2 ) denote the exit time of Zy from D y 1 (resp. D y 2 ) which is the basin of attraction of x * -(y) (resp. x * + (y)). Let us recall that D y 1 =] -∞, x * 0 (y)[ (resp. D y 2 =]x * 0 (y), +∞[). From Theorem 3.2 identity (3.10), for x ∈ D y 1 and h > 0 we obtain ∀x ∈] -∞, x * 0 (y)[, ∀h > 0 lim ǫ→0 P x (exp( V -(y) -h ǫ ) < τ y 1 < exp( V -(y)) + h ǫ )) = 1 (3.19)
An analogous result holds for x ∈ D y 2 by replacing τ y 1 by τ y 2 and V -(y) by V + (y): ∀x ∈]x * 0 (y), +∞[, ∀h > 0:

lim ǫ→0 P x (exp( V + (y) -h ǫ ) < τ y 2 < exp( V + (y)) + h ǫ )) = 1 (3.20) Proposition 3.4 Let τ y 1 (resp. τ y 2 ) denote the first exit time of Z y from D y 1 (resp. D y 2 )
. The law of τ y i is the same as the law of δτ y i for i = 1, 2. Let us assume that ǫ and δ go to 0 in such a way that

ǫ δ | log δ| → c ∈]0, +∞[. (3.21)
In this case we obtain:

∀x ∈] -∞, x * 0 (y)[, ∀h > 0 lim P x (δ c -1 (c-V -(y)+h) < τ y 1 < δ c -1 (c-V -(y)-h) ) = 1 (3.22)
An analogous result holds for x ∈ D y 2 by replacing τ y 1 by τ y 2 and V -(y) by

V + (y): ∀x ∈]x * 0 (y), +∞[, ∀h > 0 lim P x (δ c -1 (c-V + (y)+h) < τ y 2 < δ c -1 (c-V + (y)-h) ) = 1 (3.23)
These two properties enable us to introduce some remarks about the main theorems stated in the previous section. These remarks are linked to the notions of mainstate and metastable state introduced in [START_REF] Freidlin | Random perturbation of dynamical systems[END_REF]. In a general framework ( cf. [START_REF] Freidlin | On stochastic perturbation of dynamical systems with fast and slow components[END_REF]) the main state is the point towards which the cost of moving, or the transition rate, is minimum. It is not always unique: for instance in our bistable case there are two main states when y = y * . Before reaching the main state, the process may reach metastable ones accessible for shorter time lenghts. Main states may be considered as stable states, whereas metastable states are only stable in some time scale.

To study transitions of (3.18) between the two basins of attraction during [0, T (ǫ)], the relevant quantity to consider is ǫ log T (ǫ) which we must compare to the transition rates V ± (y) given by (3.12). For the time scale T (ǫ) = e c ǫ ; ǫ log T (ǫ) = c. So we must compare c to the transition rates V ± (y). Actually this amounts to compare first c to S defined in definition 3.3. We refer again to figure 1.

More precisely we can state using definition 3.

Remark:

1. When c > S the time interval [0, T (ǫ)] is long enough so that the process Z y reaches with great probability a small neighborhood of its main states. And as we can find an open interval I containing y * such that:

∀y ∈ I c > max(V -(y), V + (y)).
(3.24)

following proposition both exit time τ y 1 and τ y 2 tends to 0 in probability so we can only expect results on the slow component and the result depends on whether the boundary between the two main states is attracting or notas is shown in theorem(2.2). 2. On the countrary when c < S the time interval [0, T (ǫ)] is too short and the process only reaches with great probability a neighborhood of a metastable state. In this case we have y and it depends on whether the boundary {y * } between E 1 := {y > y * } and E 2 := {y < y * } is attractive for this system or not(cf. [START_REF] Freidlin | On stochastic perturbation of dynamical systems with fast and slow components[END_REF]). If a ∈ ]x * -(y * ), x * + (y * )[ the boundary is attractive. This is not the case when a < x * + (y * nor when a > x * + (y * ). 1. Let x * -(y) < a < x * + (y * ). The point y * is an attracting boundary since x * -(y * ) -a < 0 < x * + (y * ) -a. Assume for example that Y 0 = y < y * that is y ∈ E 2 . The evolution of (Y t ) is close to the solution of Ẏ t = x * + (Y t ) -a. Therefore since x * + (y) > a, process (Y t ) starts increasing until it reaches y * . After it has reached y * , the evolution of (Y t ) is close to the solution of Y 0 = y * , Ẏ t = B(Y t ) where B is a vector field tangent to boundary (cf. [START_REF] Freidlin | On stochastic perturbation of dynamical systems with fast and slow components[END_REF] Theorem 1) which is the 0-dimensional manifold {y * }; therefore B is zero and Y t ≡ y * . The time t(y, h) is the time necessary to reach a small ball B(y * , h). If y ∈ E 1 the evolution of (Y t ) is close to the solution of Ẏt = x * -(Y t ) -a. In this case process (Y t ) starts decreasing since x * -(y) < a until it reaches y * . The argument and the conclusion are then identical to the preceeding case.

x * -(z) < x -(c) < a < x + (c) < x * + (z) ( 4

  )[ where x -(c) and x + (c) are introduced in definition 3.4; then there exist two periodic functions Φ a c and Ψ a c given in definition 4.

3. 1

 1 Basic results on Large Deviations TheoryBecause of the slow-fast property of FitzHugh-Nagumo systems, the slow variable Y t of system (1.3)-(1.4) may be in a first approximation frozen at the value y which leads us to the study of the family of one dimensional dynamical systems indexed by parameter y,which plays a basic role in the study of FitzHugh-Nagumo systems (1.1)-(1.2) and (1.3)-(1.4):

7 ) 3 . 1

 731 Proposition The quasi potential of(3.2) w.r.t. x coincides with the function u → V (u) = -2 The above statement holds since (3.2) is one dimensional. It also holds in the multidimensional case when the drift b of (3.2) is a gradient. Theorem 3.2 Let x * be a stable equilibrium point of (3.2) such that b(r) < 0 for all r > x * , and b(r) > 0 for all r < x * . Let D be the basin of attraction of x * and τ denote the first exit time of x from D. Let us assume that D

[ with a 0 and a 1 2 Definition 3 . 2

 1232 the two points where f ′ vanishes. (i) The set {x ∈ R; f (x) = y} consists of three points x * -(y) < x * 0 (y) < x * + (y) (the equilibrium points of (3.11)) each being a continuous function of y with bounded first and second derivatives. (ii) The two points x * ± (y) are stable. x * 0 (y) is unstable.Proof of Proposition 3.2: Left to the reader; we refer to figure1-section Let us define the two functions V ± on ]f (a 0 ), f (a 1 )[ as follows:

Definition 3 . 3 Definition 3 . 4

 3334 We denote by (y * , S) the intersection point of the graphs ofV -and V + . Let E 1 := {y > y * } and E 2 := {y < y * }. For c ∈]0, S[ we denote by y ± (c) the points of ]f (a 0 ), f (a 1 )[ which satisfy y -(c) < y * < y + (c) and V -(y -(c)) = c = V + (y + (c)). Let us also define x -(c) := x * -(y -(c)) and x + (c) := x * + (y + (c)) (cf. figure1).Remark:By definition V -(y * ) = V + (y * ) = S. For f (x) = 4x -x 3 , y * = 0 and S = 4.

3 and 3. 4 :

 4 Proposition 3.5 ; 1. When c > S the main state of Zy is equal to x * + (y) (resp. x * -(y)) if y < y * (resp. y > y * ). When y = y * the two points x * ± (y * ) are both main state. 2. When c < S for y ∈]y -(c), y + (c)[ and x ∈ D y 1 (resp. x ∈ D y 2 ) the metastable state of Zy is equal to x * -(y) (resp. x * + (y)). Proof of Proposition 3.5: Direct consequence of the estimation of the time of exit given in proposition 3.4.

  -(c) < y + (c) and if we consider the interval ]y -(c), y + (c)[ when c < V -(y) < V + (y) and x ∈ D y 1 (resp. c < V -(y) < V + -(y) and x ∈ D y 2 ), for all y ∈]y -(c), y + (c)[ exit time τ y 1 (resp. τ y 2 )tends to infinity. So the process Z y remain in a neighborhood of one metastable state and switch to the other one as soon as Y t gets out of ]y -(c), y + (c)[; this give rise to a limit cycle as is shown in theorem (2.1). 4 Proof of the Main Theorems Definition 4.1 For a ∈]x -(c), x + (c)[ and y ∈]y -(c), y + (c)[ as we can check on figure 1 we have:

. 1 )

 1 Then for c ∈]0, S[, a ∈]x -(c), x + (c)[ and y ∈]y -(c), y + (c)[,the following definitions of time duration makes sense * -(v t ) -a provided v t > y *

2 .

 2 Let a > x * + (y). The point y * is not an attracting boundary since x * -(y * )a =< 0 and x * + (y * ) -a < 0. Let us assume for example that Y 0 = y ∈ E 2 . The evolution of (Y t ) is close to the solution of Ẏ t = x * + (Y t ) -a. The value f (a) is a stable equilibrium for (Y t ). If x * + (y) > a, process (Y t ) starts increasing until it reaches f (a) to which it is attracted. If now x * + (y) < a the argument and the conclusion are the same except that (Y t ) starts decreasing. If y ∈ E 1 the evolution of (Y t ) is close to the solution of Ẏt = x * -(Y t ) -a. Since x * -(y * ) -a < 0 process(Y t ) starts decreasing until it crosses y * (which is not attractive) towards E 2 . Then we are back to the previous case.

3 .

 3 The case a < x * -(y) is treated analogously.

Acknowledgment. We thank Prof. Mark Freidlin for sending us a copy of his paper [7] which was not available to us.

where y ± (c) and x ± (c) have been introduced in Proposition 3.2 and Definition 3.4.

On the interval [0, T a 1 + T a 2 ]we can introduce the solution of the following ordinary differential equations used to define the limit cycle in Theorem

Existence and regularity of solutions to odes (4.3)-(4.4) follow from the C 2 regularity of the functions x ± .. Under the assumptions we recall that we have (see figure 1)

so Ψ a c (t) increases from y -(c) to y + (c) in a duration of time equal to T a 1 (c) and decreases from y + (c) to y -(c) in a duration of time equal to T a 2 (c). The periodic function Ψ a c is obtained by continuously sticking together at y ± (c) the solutions y a ± of the following odes ẏa

Remark:

The function Ψ a c is continuous on R differentiable with continuous derivatives except at points in the set {kT a (c),

remember that V -is strictly increasing). Therefore the 13 process (X t ) leaves D y 1 instantaneously and is attracted to a neighborhood of x * + (y) > x + (c). The point (X t , Y t ) remains in a neighbourhood of the branch v = f (u) containing (x * + (y), y) and Y t increases as long as X t > x + (c) since dY t = (X t -a)dt. However identity (3.23) implies that lim * P x (τ z 2 = +∞) = 1 for z < y + (c) (resp. lim * P x (τ z 2 = 0) = 1 for z > y + (c) ). Therefore the point (X t , Y t ) is instantaneously attracted to a neighborhood of (x * -(y + (c)), y + (c)) after Y t has crossed y + (c) since the speed of Y t is strictly positive in a neighbourhood of (x + (c), y + (c)). The argument is then the same as before. We detail it for completeness. Since x * -(y + (c)) < x -(c) < a, the second coordinate Y t decreases as long as X t -a < 0. However identity (3.22) implies that lim * P x (τ z

Hence (X, Y ) converges in probability to a limit cycle of period T a (c) = T a 1 (c) + T a 2 (c). b) Let y > y -(c) and x ∈ D y 1 . By identity 3.22 lim * P x (τ y 1 = +∞) = 1, the fast process (X t ) is attracted to x * -(y). After this first phase, one can apply the same argument as in a).

2.Case a < x -(c). As already mentioned in section 3, when c < S, the relevant states are the metastable ones. The assumption a < x -(c) implies f (a) > y -(c). Let V be a small neighborhood V of (a, f (a)). Let us first assume that V is so small that all (u, v) ∈ V satisfies v > y -(c) and accordingly u < x -(c). Let (X t , Y t ) start from (u, v). In particular u ∈ D v 1 and τ v 1 = +∞ since V -(v) > c. Then (Y t ) will evolve like the solution of vt = x * -(v t ) -a, v 0 = v for which f (a) is a stable equilibrium point. Let us now assume that there are points (u, v) ∈ V such that v < y -(c). Then (X t ) instantaneously jumps to the right stable branch; (Y t ) becomes close to the solution of vt = x * + (v t ) -a and therefore increases until it reaches y + (c). At that time it jumps to the left stable branch and we are back to the previous argument since then (Y t ) becomes close to the solution of vt = x * -(v t ) -a. Let us notice that for the slow variable Y t we get a result using the uniform norm by large deviation estimates but for the quick variable X t the result is formulated in L 1 norm thanks to equation (1.4).

Proof of Theorem 2.2. When c > S, as pointed out in subsection (3) Z y has time to reach its main state so the evolution of Y t is close to the solution of vt = x * + (v t ) -a provided v t < y *

[14] Piquet, C., Personnal communication.


Sylvain Caubet

Yves Le Guer
email: yves.leguer@univ-pau.fr

Bruno Grassl

E L Kamal

Omari

Eric Normandin

A LOW-ENERGY EMULSIFICATION BATCH MIXER FOR CONCENTRATED OIL-IN-WATER EMULSIONS

Keywords: Batch mixer, Low energy emulsification, Concentrated emulsions, Droplet size distribution control

The presented work shows the formation of a monodisperse High Internal Phase Ratio (HIPR) Oil-in-Water (O/W) emulsion using a new type of a two-rod batch mixer. The mixture components have sharply different viscosities (1/3400 for W/O) and similar densities (1/0.974 for W/O). The oil ratio in the mixture remains constant at 91% [w/w] for each experiment. The simple design of the new mixer leads to a low-energy process, characterized by very low rotational speed, laminar flow and no need for fluid heating. The droplet size distribution during the emulsification was investigated according to different physical and formulation parameters such as stirring time, rotational speed, surfactant type, concentration, and salt addition. We show that all studied parameters except salt addition allow a precise control of the concentrated emulsion droplet size distribution, which governs the rheological behavior of the emulsion and its stability. The formed HIPR emulsion shows two shear-thinning domains.

Introduction

To disperse a fluid into another immiscible fluid, some mechanical energy, e.g., shear, is applied, together with a surfactant, to produce a dispersion with stable properties, including final droplet size, polydispersity and rheology. The nature of the two fluids, the surfactant and the process conditions (design of the mixer, mixing rate and time) all have a critical effect on the properties of the final emulsion. During the process of emulsification, external shear energy is used to break large drops into smaller ones. For an internal volume ratio above 0.74, emulsions are called High Internal Phase Ratio (HIPR) emulsions, bi-liquid foam systems or gel emulsions [START_REF] Kizling | On the formation and stability of high internal phase O/W emulsions[END_REF]. HIPR emulsions are encountered in a wide variety of industrial applications: petrochemical for bitumen or heavy crude oils [START_REF] Fournanty | Laminar Flow Emulsification Process to Control the Viscosity Reduction of Heavy Crude Oils[END_REF]; agrochemicals; pharmaceuticals; cosmetics; etc. The main works of Princen [START_REF] Princen | Highly concentrated emulsions, II. Real systems. The effect of film thickness and contact angle on the volume fraction in creamed emulsions[END_REF][START_REF] Princen | Rheology of foams and highly concentrated emulsions. II: Experimental study of the yield stress and wall effects for concentrated oil-in-water emulsions[END_REF] have shown the rheological and structural properties of gel emulsions. Most of the proprieties of gel emulsions, particularly the rheological behavior, are affected by the droplet size distribution, the volume fraction and the interfacial tension between oil and water [START_REF] Princen | Rheology of foams and highly concentrated emulsions. IV: An experimental study of the shear viscosity and yield stress of concentrated emulsions[END_REF][START_REF] Welch | Rheology of high internal phase emulsions[END_REF]. Thus, during the emulsification process, the droplet size distribution is an important parameter to control to obtain a HIPR emulsion with the desired properties. This depends on several parameters, such as the process, the surfactant, the water/oil ratio, the type of oil (vegetable, synthetic, bitumen, etc.) and other formulation and physical parameters used for the mixing. Typically, the droplet diameter can vary from one to hundreds of micrometers. Rheological studies concerning emulsions have only begun recently and are mostly concerned with "dispersed" emulsions (with an internal volume ratio lower than 0.74), while most of HIPR emulsion studies are made on Water-in-Oil (W/O) emulsions [START_REF] Malkin | Effect of droplet size on the rheological properties of highly-concentrated w/o emulsions[END_REF][START_REF] Masalova | Evolution of rheological properties of highly concentrated emulsions with aging-Emulsion-to-suspension transition[END_REF]. HIPR inverse emulsions (O/W), which are considered in this work, generally have a shear thinning behavior with viscosity ranging from a few to several thousand Pa.s as a function of the shear rate, the droplet size distribution and the internal volume ratio [START_REF] Cuellar | More efficient preparation of parenteral emulsions or how to improve a pharmaceutical recipe by formulation engineering[END_REF][START_REF] Otsubo | Effect of drop size distribution on the flow behavior of oil-in-water emulsions[END_REF][START_REF] Yaghi | Rheology of oil-in-water emulsions containing fine particles[END_REF]. Several processes exist to produce HIPR Oil-in-Water emulsions, for example: Coaxial Mixers [START_REF] Gingras | Concentrated bitumen-in-water emulsification in coaxial mixers[END_REF], Colloid Mill [START_REF] Karbstein | Developments in the continuous mechanical production of oil-in-water macro-emulsions[END_REF], High-Pressure Homogenizers [START_REF] Stang | Emulsification in high-pressure homogenizers[END_REF], Static Mixers [START_REF] Gingras | Inline bitumen emulsification using static mixers[END_REF] or the use of a complex mechanism of the inversion phenomenon [START_REF] Salager | Using emulsion inversion in industrial processes[END_REF]. Of course, this list is not exhaustive; other methods and apparata exist for making HIPR O/W emulsions. Here, the issue of mixing effects on the size, distribution and stability of concentrated O/W emulsions with contrasting viscosities (1/3, 400 for W/O) and similar densities (1/0.974 for W/O) is addressed for a new type of laminar batch mixer. The HIPR emulsions are made in a batch Two-Rod Mixer (TRM) at a high volume ratio (0.91) of castor oil in water. We investigate how the droplet size distribution is affected by the variation of physical and formulation parameters such as stirring time, rotational speed, surfactant type and its concentration and addition of salt. Rheological behaviors are also presented.

Materials and Methods

Materials

The oil used is blown castor oil from Seatons (Ricin oil). It was chosen mainly for its high density (974 kg.m -3 ) and high viscosity (3.4 P a.s at 20 • C). Two well-known non-ionic surfactants were used and compared: Triton X-405 (HLB = 17.9, CMC in water of 0.159 [wt/wt], product of Sigma-Adrich) and Triton X-100 (HLB = 13.5, CMC in water of 0.018 [wt/wt], product of Sigma-Aldrich) [START_REF] Paillet | Rapid and quantitative determination of critical micelle concentration by automatic continuous mixing and static light scattering[END_REF][START_REF] Hait | Determination of critical micelle concentration (CMC) of nonionic surfactants by donor-acceptor interaction with lodine and correlation of CMC with hydrophile-lipophile balance and other parameters of the surfactants[END_REF]. Deionized water, with an electrical conductivity of 18.3 M Ωcm, was filtered through a 0.22 µm Millipore filter. Sodium chloride (NaCl) was purchased from Sigma-Aldrich.

Emulsion preparation

The experimental setup used to produce the HIPR emulsion was a batch Two-Rod Mixer (TRM, Figure 1 (a) and (b)). The larger rod was not coaxial with the tank, and its rotational speed was controlled by a computer. A second, smaller, rod, not rotating, was located next to the first one. The geometrical parameters were as follows: tank diameter d ′ = 150 mm, larger rod radius r = 50 mm, smaller rod diameter d = 15 mm, rod-tank gap e = 3 mm, rod-rod gap e ′ = 2 mm and the angular position of the smaller rod θ = 30 • (Figure 1 (b)). These geometrical parameters remained the same for all of the experiments considered in this study. A system with a strength sensor was adapted to the motor of the larger rod in order to follow the evolution of the resistive torque during the emulsification process. The calculation of the torque and thus of the power input to keep the rotational speed constant can be calculated using the following relations:

C = F • a (1)
P = C • Ω ( 2
)
where F is the measured strength in newtons; a, the distance between the sensor and the axis of the motor (a = 0.185 m); C, the torque (N.m); Ω, the rotational speed of the larger rod, and P , the power (W ).

The protocol used to produce the HIPR Oil-in-Water emulsion is as follows. First, 400 mL of castor oil was added to the tank. Then, the larger rod was rotated at the selected rotational speed (Ω = 60, 90 or 120 rpm). Finally, after one minute, 38.6 mL of a mixture of distilled water and surfactant was rapidly added into the TRM batch process. The emulsification began at time t = 0s when the water-surfactant mixture is added in the tank. The emulsification can be decomposed into three steps. First, the water-surfactant mixture is introduced into the tank. Rapidly thereafter (after a few seconds), the emulsification begins, droplets of oil are formed in the water phase, and all of the Castor oil is gradually dispersed in the water. A crude emulsion is then created, which leads to the formation of the first monodisperse emulsion. The duration t 0 , which corresponds to the creation of this primary monodisperse emulsion, ranges between one and several minutes and depends in particular on the rotational speed of the larger TRM rod (see section 3.2). This initial time t 0 is an important characteristic time because it represents the limit between two regimes in the evolution of the mean droplet diameter as a function of the mixing time [START_REF] Mabille | Monodisperse fragmentation in emulsions: Mechanisms and kinetics[END_REF]. In the last step, the droplets of the initial monodisperse emulsion are gradually broken into smaller ones. This mechanism induces a reduction of the mean diameter until a limit value is reached (see section 3.1). During the process, the temperature was measured in situ and its variation (in Celcius) did not exceed 10% of its initial value for all studied rotational speeds.

Determination of droplet size distribution

The droplet size distribution was determined via optical microscopy according to the method described by Tcholakova et al. [START_REF] Tcholakova | Role of surfactant type and concentration for the mean drop size during emulsification in turbulent flow[END_REF]. The samples, carefully taken from the TRM with a syringe in order to not disturb the mixing, were immediately transferred to 4 wt% Triton (X-405 or X-100) solutions to prevent further drop coalescence and to reduce the droplet concentration to about 3 vol% because the original emulsion was too concentrated for droplet size analysis (Figure 2). This dilution was transferred for optical examination into microcapillaries with rectangular cross sections (depth 0.2 mm, width 0.7 mm, length 1 mm) no more than 1 h after their preparation. Before loading the sample into the microcapillary, the vial containing the Triton-stabilized emulsion was gently rotated to homogenize the emulsion. The oil droplets were observed via transmitted light microscopy with a Leica DMLM microscope (magnification ×20) connected to a charge-coupled device camera (Leica DFC280) and video recorder software (Leica IM4.0). The diameters of individual recorded oil droplets were measured using the publicly available software ImageJ 1.40, released by the National Institute of Health (NIH). These data were numerically processed to obtain droplet size histograms. The diameters of at least 800 droplets were measured for each system. Illustrative droplet size histograms are shown in Figure 3.

The mean diameter d m , the mean volume-surface diameter (the Sauter diameter) d 32 and the uniformity factor U were calculated from the following relations: The surfactant used was 15.9 wt% Triton X-405 in water.

d m = n i d i n i (3)
d 32 = n i d 3 i n i d 2 i (4) U = 1 d * n i d 3 i | d * -d i | n i d 3 i (5)
where n i is the number of droplets of diameter d i and is the median diameter (the diameter for which the cumulative undersize volume fraction is equal to 0.5). Emulsions are considered to be monodisperse if U is smaller than 0.25 [START_REF] Mabille | Monodisperse fragmentation in emulsions: Mechanisms and kinetics[END_REF].

Rheological measurements

The steady shear rheological properties of the HIPR emulsion were measured using a Malvern Bohlin C-VOR 150 rheometer in controlled stress mode with a cone-plate geometry (diameter of 60 mm, cone angle of 2 • ). Stepped shear stress values were selected to study shear rates ranging from 0.1 to 100 s -1 . The system was allowed to reach steady state at each shear stress prior to registering the measured values. The measurements were performed by increasing and decreasing the shear stress. Temperature was controlled with a Peltier plate system and was set at 20.0 ± 0.1 • C.

Experimental results and discussion

Effect of the mixing time on the mean diameter

In this study, the emulsion was made using a 15.9 wt% solution of Triton X-405 in distilled water as described in section 2.2. For this experiment, the rotational speed of the rod remained constant and equal to 90 rpm throughout the process. During the emulsification, samples were taken at times t = 300, 420, 780, 1140, 2400 and 3600 s. The droplet size distribution was then determined, and its probability distribution was plotted for each sampling time (Figure 2). The emulsion can be considered to be monodisperse because its uniformity factor U, calculated using equation [START_REF] Princen | Rheology of foams and highly concentrated emulsions. IV: An experimental study of the shear viscosity and yield stress of concentrated emulsions[END_REF], is less than or close to 0.25 (Figure 4). The evolution of the mean diameter during the mixing over one experiment is shown in Figure 3. Two periods can be distinguished. In the first period, the mean diameter decreases rapidly from 23 µm (t = 180 s) to 15 µm (t = 780 s); i.e., in 600 seconds, d m is reduced by 8 µm. In the second period, a slower mechanism takes place, and the mean diameter is reduced from 15 to 11 µm over 1620 seconds.

A crude emulsion can be transformed into a monodisperse one, as shown by Mason and Bibette [START_REF] Mason | Emulsification in Viscoelastic Media[END_REF], who have experimentally discovered this phenomenon by applying a shear step to a crude emulsion. Moreover, Mabille et al. [START_REF] Mabille | Monodisperse fragmentation in emulsions: Mechanisms and kinetics[END_REF] have shown that the fragmentation kinetics of a calibrated emulsion as a function of the initial droplet size involves two distinct regimes. At short timepoints, the droplet diameter decreases abruptly. The droplets deform into long threads that undergo Rayleigh instability. During this regime, the obtained diameter is mainly determined by the applied stress and weakly depends on the viscosity ratio between the dispersed and continuous phases. After this first step, the resulting droplets can, once again, break up into daughter droplets. This second mechanism is much slower, with a characteristic duration of several hundred seconds. Depending on the initial size, the first step can vanish, and only the second slow step persists. In our process, if the crude emulsion was mainly composed of large drops, even after a few seconds of shear, we obtained a wellcalibrated emulsion (with a mean diameter close to 24 µm); all drops deformed into threads of different lengths, but the Rayleigh instability occurs when the same critical radius is reached by each thread. Subsequently, when the first monodisperse emulsion was obtained, the second slow mechanism took place, and since the smaller droplets (between d 0 and d s ) also broke up, the size distribution became even narrower in the second regime. However, it is clear that the most efficient process to get narrowsize distributions is the initial Rayleigh instability. To fit their experimental data, Mabille et al. [START_REF] Mabille | Monodisperse fragmentation in emulsions: Mechanisms and kinetics[END_REF] Figure 4: Evolution of the mean droplet diameter during mixing. Points are experimental data, and the line is fitted to experimental data by the use of Eq. ( 7). used the following equation:

d m (t) = d s 1 - d 0 -d s d 0 exp -ω 0 t 3 (6)
In this equation, d 0 is the diameter of the initial monodisperse emulsion in the slow regime, d s is the saturating diameter, and ω 0 is the rupture rate. To examine the kinetics of this fragmentation in our process, we must account for the start of the primary monodisperse emulsion. At the time t = 0 s, we do not have an emulsion; the primary monodisperse emulsion is created during the process at a time defined by t = t 0 and corresponding to a drop diameter d 0 . This primary monodisperse emulsion is, in fact, the equivalent of the mother emulsion of Mabille et al. [START_REF] Mabille | Monodisperse fragmentation in emulsions: Mechanisms and kinetics[END_REF] Thus, to fit our data, we have to adapt equation ( 6) by introducing t 0 to take into account the time shift and to determine the rupture rate ω 0 :

d m (t) = d s 1 - d 0 -d s d 0 exp -ω 0 (t -t 0 ) 3 (7)
t 0 , d 0 , ω 0 and ds are the free parameters used to fit this model (Eq. 7) with our data, using appropriate limit values to fit this equation into the monodisperse regime. Thus, we determine t 0 = 723 s, d 0 = 15.4 µm, d s = 10.9 µm, and the rupture rate, ω 0 , is equal to 0.0039 s -1 . This value of ω 0 is five times smaller than the one found by Mabille et al. for their process. The inverse of ω 0 corresponds to the characteristic fragmentation time of 256 s. Compared to the Rayleigh instability, this regime is extremely long; therefore, the two regimes can be decoupled. This high value of the characteristic fragmentation time is an advantage because it allows the production of different emulsions solely by varying the emulsification time while keeping the same process and the same components. After t = 2400 s of mixing, emulsification is achieved, and the mean drop diameter remains constant.

Effect of rotational speed on the final diameter

In this section, the surfactant used was also Triton X-405. The composition, the surfactant concentration, the batch TRM geometry and the protocol used were the same as those previously described.

The only variable parameter is the rotational speed of the rod. For each speed studied (ω = 60, 90 and 120 rpm), the emulsification was performed over 2400 seconds (40 minutes), and the power necessary to maintain constant rotational speed was calculated (Eqs. ( 1) and (2), Figure 5). The mean droplet diameter d m and the uniformity factor U (Eqs (3) and ( 5)) were determined at different mixing times for each obtained emulsion and for the different rotational speeds (Figure 6). At the beginning of the emulsification, the power rapidly increased until a certain limit, at which time, its value remained nearly constant. This type of behavior has been previously observed by Sanchez et al. [START_REF] Sanchez | Emulsification rheokinetics of nonionic surfactant-stabilized oil-in-water emulsions[END_REF] for two types of impellers (anchor and helical ribbon) at different rotational speeds. The main difference observed with their results concerns the early stages of mixing, due to the presence in our process of a macro-mixing step that does not exist in theirs because it represents a pre-emulsification step. This first step can be explained by a lubrication of the rod, leading to a decrease of the measured power. During this macro-mixing phase, the power remains constant until emulsification begins.

(a) (b)

Figure 5: (a) Evolution of power during the emulsification for different rotational speeds of the rod (circles: ω = 120 rpm, squares: ω = 90 rpm, triangles: ω = 60 rpm). Curve fits were obtained using equation ( 9). (b) Expansion of the early time of mixing for ω = 90 rpm (the same behavior was observed at all rotational speeds). The surfactant used in each case was 15.9 wt% Triton X-405 in water.

To fit the experimental data, we adapted the first-order kinetics equation given by Sanchez et al. [START_REF] Sanchez | Emulsification rheokinetics of nonionic surfactant-stabilized oil-in-water emulsions[END_REF]:

P (t) = P e -P e exp[-k(t -t m )] (8)
To take into account the fact that in our process, we have, contrary to Sanchez et al., a first step in which the power P has an initial value, we adapted equation ( 8) to obtain:

P (t) = P e + (P m -P e ) exp[-k(t -t m )] (9)
with P e the final power, k the kinetic constant, t m the macro-mixing time, and P m the power at the end of the macro-mixing phase. The values of k, P e , P m and P 0 (the initial value of the power before addition of the surfactant mixture) are given in Table 1.

The macro-mixing time t m is shorter for the HIPR emulsion obtained with the rotational speed of 120 rpm and is up to 6 times longer than at 90 and 60 rpm, for which it remains roughly equal. The kinetic constant k and the final power P e decreased with the rotational speed, as was also observed by Sanchez et al. It is important to note that the power input necessary for the formation of the emulsion was always lower than the initial power P 0 given when only oil is present in the mixer. This is notable 7) and ( 9) and of t 5% as functions of the rotational speed.

Ω 1 (
because for all rotational speeds, the shear rate in the TRM was between 1 s -1 (in the larger tank space) and 200 s -1 (in the rod-rod space). Furthermore, in this range of shear rates, the viscosity of the final HIPR emulsion was higher than that of the castor oil (as shown hereafter). Thus, the strength applied to the rod to keep the rotational speed constant should be higher than that required by the castor oil alone; however, this was not the case. This could be explained by a thin layer of water adhering to and lubricating the rod. For each rotational speed, we determined the mean drop diameter dm at different mixing times (Figure 6). The experimental data were fitted using equation ( 7), and the different values of ω 0 , d s , d 0 and t 0 obtained for each rotational speed are presented in Table 1.

Figure 6: of the mean droplet diameter during mixing for different rotational speeds. ω = 60 rpm (triangles), ω = 90 rpm (squares) and ω = 120 rpm (circles). Experimental points were fitted using equation [START_REF] Malkin | Effect of droplet size on the rheological properties of highly-concentrated w/o emulsions[END_REF].

The mean droplet radius decreases with increasing rotational speed of the rod. For all rotational speeds, the evolution of the mean drop diameter can be modeled using equation [START_REF] Malkin | Effect of droplet size on the rheological properties of highly-concentrated w/o emulsions[END_REF], and the rupture rate ω 0 remains small regardless of the rotational speed of the rod. The initial emulsification time t 0 also decreases with increasing rotational speed. In each case, the mean droplet diameter remains constant after about 2400 seconds of mixing, and the value of this final droplet diameter is dependent on the rotational speed. A rotational speed of 60 rpm can be used to obtain an emulsion with a mean diameter value ranging between 15 and 30 µm. In contrast, the rotational speed of 120 rpm only permits the production of a mean droplet diameter lower than 15 µm.

Here, we can introduce another characteristic time: t 5% , corresponding to the time for which the relative deviation between the mean droplet diameter d 5% and the final one is equal to 5% of ds (d 5% = 1.05 d s ). Table 1 gives all of the characteristic times for each rotational speed studied. These characteristic times correspond to the three periods of our emulsification process explained in section 2.2. First, during the macro-mixing time (between t = 0 s and t = t m ), the introduced watersurfactant mixture lubricates the rod but does not disperse in the oil. Next, from t m to t 0 , the droplets of oil are gradually generated in the water to create the primary monodisperse emulsion. Finally, the fragmentation of the emulsion takes place by reducing the mean droplet diameter until the limit value. The last two periods correspond to the ones observed by Mabille et al [START_REF] Mabille | Monodisperse fragmentation in emulsions: Mechanisms and kinetics[END_REF]. By integrating equation ( 9) between the macro-mixing time tm and a mixing time t, we could determine the evolution of the amount of energy stored in the TRM per unit volume of dispersed phase, E V d (energy density in J.m -3 , Eq. ( 11)) for each rotational speed of the rod:

V d .E V d (t) = t tm P e -(P e -P m ) exp[-k(t -t m )]dt (10)
E V d (t) = P e (t -t m ) - (P e -P m ) k (exp[-k(t -t m )] -1) V d (11)
In this equation, V d is the volume of the dispersed phase (here, the volume of oil). The others parameters have been previously defined. In Figure 7, we plotted the Sauter diameter as a function of the energy density input for the different rotational speeds studied and for emulsification times ranging between tm and t 5 %. We have selected these times to account only for the energy used to reach the final mean droplet diameter (about 2400 s). During this effective period, The Sauter diameter can be expressed as a function of the energy density using equation [START_REF] Gingras | Concentrated bitumen-in-water emulsification in coaxial mixers[END_REF] to fit the experimental data [START_REF] Stang | Emulsification in high-pressure homogenizers[END_REF]:

d 32 = C.(E V ) -b (12)
For each rotational speed, the Sauter diameter decreases with increasing energy density input into the mixer. However, the energy density input necessary to obtain the same value of the Sauter diameter depends on the rotational speed of the rod (Figure 7); likewise, for the same energy density input, the obtained Sauter diameter is dependent upon the rotational speed (Figure 7). Karbstein and Schubert [START_REF] Karbstein | Developments in the continuous mechanical production of oil-in-water macro-emulsions[END_REF] have obtained HIPR emulsions composed of 80% rapeseed oil (viscosity: 6.1 × 10 -3 P a.s at 20 • C) in water using egg yolk as an emulsifier using different mixers: toothed colloid mill, toothed disc dispersing machine and high-pressure homogenizer. We found that to obtain emulsions in the TRM having Sauter diameters equivalent to those obtained by Karbstein and Schubert [START_REF] Karbstein | Developments in the continuous mechanical production of oil-in-water macro-emulsions[END_REF], the necessary energy density is ten times smaller in our process. However, the relationship between the Sauter diameter and the energy density also depends upon the composition of the emulsion; thus, this comparison is only indicative. Effectively, to compare the energy efficiency of the TRM to those of other mixers, we must make emulsions with the same components and the same proportions. Furthermore, to calculate the energy density, Karbstein and Schubert used the relation E V = P V , where V corresponds to the flow rate of the emulsion. For high-pressure homogenizers, Stang et al. [START_REF] Stang | Emulsification in high-pressure homogenizers[END_REF] have obtained energy densities ten times larger than those obtained by our process for different oil/water viscosity ratios (however, their Sauter diameter is smaller, ranging between 1 and 15 µm).

Effect of surfactant type and concentration

In these experiments, the only variable parameter was the concentration of surfactant used for the emulsification. Two non-ionic surfactants are compared: Triton X-100 and Triton X-405. Emulsions were made with a rotational speed of 90 rpm. The emulsification time chosen for the comparison was 780 seconds (13 minutes), which corresponds to the end of the first period of evolution of the mean droplet diameter during the mixing (for 15.9 wt% of Triton X-405, see section 3.1) and to the time at which the power is stabilized (see section 3.2). For this timepoint, the mean droplet diameters were determined and compared (Figure 8).

For the two studied surfactants, the evolution of the mean droplet diameters presented two regions. In the first time period, the mean droplet diameter decreased rapidly, and in a second one, this diameter decreased slowly towards an asymptote. These two regions have also been observed by Tcholakova et al. [START_REF] Tcholakova | Role of surfactant type and concentration for the mean drop size during emulsification in turbulent flow[END_REF] for less concentrated emulsions made using turbulent flow with another non-ionic surfactant (Brij 58). For each concentration, emulsions made with Triton X-100 have a mean droplet size smaller than that for the ones made with Triton X-405. The uniformity factor U has also been determined for each emulsion via equation [START_REF] Princen | Rheology of foams and highly concentrated emulsions. IV: An experimental study of the shear viscosity and yield stress of concentrated emulsions[END_REF]. At low and high concentrations of surfactant, emulsions are monodisperse; however, the ones produced at low concentration (below 4 wt% of surfactant in water) are unstable. Effectively, coalescence takes place after one day, whereas for emulsions made with higher concentrations of surfactant, the breaking occurred only after several months. All of the HIPR emulsions seemed to be monodisperse, with U values lower than 0.25.

Effect of the mean droplet diameter on the rheological behavior

We know that the rheological behavior of the emulsions with the same composition is mainly dependent upon the droplet size distribution. Thus, to compare the droplet size distributions between two emulsions, we can compare their rheological behaviors [START_REF] Langenfeld | Rheological behavior of fluorinated highly concentrated reverse emulsions with temperature[END_REF]. Figure 9 confirms the influence of the mean droplets radius on the rheological behavior of the HIPR castor Oil-in-Water emulsion obtained in the TRM. At all values of the mean droplet diameter, the viscosity behaviors of these HIPR emulsions (Figure 9) are similar and do not present a Newtonian plateau at low shear rates. These emulsions have a shear thinning behavior in the entire range of the shear rate study. However, there appears to be a break (change in the slope) in this shear thinning behavior for a shear rate value between 0.1 and 1 s -1 . This viscosity behavior is the result of the rearrangement of the droplets, which evolve from a state of equilibrium to a state of flow [START_REF] Saiki | Droplet deformability and emulsion rheology: steady and dynamic behavior[END_REF]. Moreover, Figure 9 also illustrates the dependence of the mean droplet size on the viscosity value, and, as was already observed [START_REF] Malkin | Effect of droplet size on the rheological properties of highly-concentrated w/o emulsions[END_REF], the viscosity of the HIPR emulsion increases with decreasing mean droplet radius. The evolution of the storage modulus G ′ and of the loss modulus G ′′ as functions of the frequency observed in Figure 9 was also observed by Mason [START_REF] Mason | New fundamental concepts in emulsion rheology[END_REF] for HIPR emulsions. At an early timepoint, the storage modulus increases until it exhibits a plateau, whereas the loss modulus presents a minimum in the frequency band where the storage modulus plateau appears. This minimum is also a manifestation of the rearrangement of the droplets in the HIPR emulsion. In accordance with the results observed by Langfeld et al. [START_REF] Langenfeld | Rheological behavior of fluorinated highly concentrated reverse emulsions with temperature[END_REF], we note a decrease in the elasticity with increasing mean droplet diameter of the HIPR emulsion.

Effect of salt

Within the context of a targeted application, such as petroleum, biomaterial engineering, pharmaceuticals, cosmetics, etc., it was interesting to evaluate the effect of added salt on the emulsification performance of HIPR emulsion via rheological behavior. An HIPR emulsion made with 30 g/L of salt in the initial surfactant mixture was compared to one produced classically, both with a rotational speed equal to 90 rpm and with Triton X-405 at a concentration of 15.9 wt% in water. In section 3.4, we discussed the influence of the mean droplet diameter on HIPR rheological behavior. Here, we compare the effect of shear rate on the viscosity and the evolution of storage and the loss modulus as functions of frequency between the two HIPR emulsions made with and without salt (Figure 10). As was explained previously, these comparisons permit to indirectly compare the mean droplet diameter of these HIPR emulsions.

After 40 minutes of mixing, the evolutions of the viscosity and the stress as functions of the shear rate present similar behaviors for the produced emulsions with and without salt. Likewise, the storage modulus G ′ and the loss modulus G ′′ present the same variation as functions of the frequency for both of the emulsions. Thus, we can affirm that, for these mixing conditions, added salt does not affect the emulsification.

Conclusion

We have studied and characterized several parameters affecting the droplet size distributions of HIPR emulsions produced in a Two-Rod Mixer batch process operating in the laminar flow regime. By adjusting certain parameters (stirring time, rotational speed, surfactant type and concentration and salt addition), we can control the droplet size distribution and thus the rheological and stability characteristics of the created emulsion. By comparison with other mixing techniques, the droplet size decreases slowly during the mixing time; thus, it is possible to easily obtain emulsions with targeted droplet size distributions. Here, a first comparison has been made of the droplet size as a function of the energy density input, and the Two-Rod mixer batch process seems to consume less energy than do other processes used to obtain similar values of the droplet size. Furthermore, the different periods of HIPR emulsification in the TRM have been highlighted, as well as their characteristics times. First, we have a macro-mixing regime from t = 0 s (when the water and surfactant mixture is introduced into the tank) to t = t m . Next, the emulsification begins, and we obtain a first monodisperse emulsion at the time t = t 0 . Finally, in the last regime, droplets of this initial monodisperse emulsion are gradually broken, and the mean droplet diameter decreases progressively to reach a limit value. Another advantage of the TRM is that the relatively slow rotational speed of the rod prevents unsuitable heating of the fluids during the emulsification process. This eliminates the need for a final step of cooling of the emulsion encountered in more classical processes that use turbulent flows with high rotational velocities.

Figure 1 :

 1 Figure 1: (a) Photo of the Two-Rod Mixer (TRM) batch device; (b) Schematic of the TRM batch with the dimensions of the elements; (c) 3D view of the TRM equipped with a strength sensor.

Figure 2 :Figure 3 :

 23 Figure 2: Optical photomicrograph of concentrated emulsion obtained using the batch Two-Rod Mixer process after 780 seconds of mixing with a rotational speed of the rod equal to 90 rpm. The surfactant used was 15.9 wt% Triton X-405 in water. Inset: Optical photomicrograph after dilution. In both cases, the scale bar is 20µm.

Figure 7 :

 7 Figure 7: Sauter diameter as a function of energy density EVd for the different rotational speeds: ω = 60 rpm (triangles), ω = 90 rpm (squares) and ω = 120 rpm (circles). The values of the couple (C, b) used in equation (10) are, respectively: (18085, 0.5); (6916, 0.4); (553, 0.3).

Figure 8 :

 8 Figure 8: Evolution of (a) the mean droplet diameter and (b) the uniformity factor U as functions of the surfactant concentration for emulsions stabilized by Triton X-100 (squares) and Triton X-405 (circles)

Figure 9 :

 9 Figure 9: Comparison of the rheological behaviors between HIPR Castor oil in water emulsions made with Triton X-405 (15.9 wt% in water) for different mean droplet radii: d m = 8.6 µm (circles), d m = 11µm (squares) and d m = 15µm (triangles). (a) Viscosity as a function of shear rate. (b) Evolution of the storage modulus G ′ (closed symbols) and the loss modulus G ′′ (open symbols) as functions of frequency.

Figure 10 :

 10 Figure 10: (a) Viscosity (closed symbols) and stress (open symbols) as functions of shear rate for emulsions made with (triangles) and without (circles) salt. (b) Evolution of the storage modulus G ′ (closed symbols) and the loss modulus G ′′ (open symbols) as functions of frequency for emulsions made with (triangles) and without (circles) salt.

Table 1 :

 1 Values of the different parameters appearing in equations (

	rpm) 120 (W) 90 (W) 60 (W)
	P 0 (W)	0.6	0.4	0.17
	P m (W)	0.10	0.14	0.05
	P e (W)	0.39	0.28	0.13
	k(s -1 )	6	38	35
	t m (s)	0.03	0.01	0.005
	ω 0	0.0033	0.0039 0.0030
	d 0 (µm)	10.5	15.4	20
	d s (µm)	8.3	10.9	14.5
	t 0 (s)	546	723	883
	t 5% (s)	1539	2318	1977

Acknowledgments

We would like to gratefully acknowledge Marc Rivalleto and Daniel Champier from the UPPA technical center "Innov'Adour" for their help during the building and the instrumentation of the Two-Rod Mixer batch process. This work was supported by a Graduate Fellowship (S.C.) from the Communauté d'Agglomération Pau-Pyrénées (CDAPP).


Céline Marie

Elisabeth Leclerc

Pauline Cavagna

Cécile Chung

Victoria Gauthier

Filomena Marques-Tavares

Christine Fernandez

Ariel Cohen

Marie Antignac

Medication reconciliation: Assessment of 4 years' experience in a cardiology department

Keywords: Medication reconciliation, Cardiology department, Unintentional medication discrepancy, Clinical pharmacy ATC, Anatomical, Therapeutic and Chemical, BPMH, best possible medication history, ICU, intensive care unit, IMD, intentional medication discrepancy, UMD, unintentional medication discrepancy, WHO, World Health Organization

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

In 2017 the World Health Organization (WHO) set up a global initiative to reduce serious preventable effects of medication errors by 50% in all countries over the subsequent 5 years.

Medication reconciliation is key to ensuring patient safety and reducing the risk of prescription error [START_REF] Mekonnen | Pharmacy-led medication reconciliation programmes at hospital transitions: a systematic review and meta-analysis[END_REF]. The aims of medication reconciliation are prevention and detection of prescription discrepancies at transition stages and at every step of medication management process. The WHO defines medication reconciliation as "a formal process in which healthcare professionals partner with patients to ensure accurate and complete medication information transfer at interfaces of care" [START_REF]The High5s Project -Standard Operating Protocol for Medication Reconciliation[END_REF].

The first step in medication reconciliation is to obtain the best possible medication history (BPMH), which is a review of all medications that the patient was taking before hospitalization. To obtain this information, three sources are required: the patient and/or family; medication packages; and medical records and/or the patient's physicians or pharmacists. The second step consists of comparing the BPMH with the first hospital prescription. Then, medication discrepancies between the BPMH and the hospital prescription can be identified by pharmacists, and can lead to a pharmaceutical intervention to propose a modification of medication to prescribers. This intervention allows a collaborative discussion between pharmacists and physicians to classify the type of discrepancy as intentional or unintentional, and to correct the prescription, if necessary. The French National Authority for Health defines intentional medication discrepancies (IMDs) as voluntary deviation, documented or not documented in the patient record; whereas unintentional medication discrepancies (UMDs) are defined as involuntary, and are considered as medication errors. IMDs are physicians' choices, i.e. there is a medical reason to modify the patient's prescription. UMDs are medication errors, such as overdosing, underdosing or omission or substitution of the patient's medication, and must be corrected.

Medication reconciliation was implemented in the cardiology department of Saint-Antoine Hospital in November 2017, after analysis of the CACTUS study [START_REF] Chung | Medication reconciliation: Predictors of risk of unintentional medication discrepancies in the cardiology department[END_REF], which identified patients' prioritization criteria: eight or more prescribed medications; French not spoken; low educational level; admission to a hospitalization unit rather than an intensive care unit (ICU); and two or more co-morbidities.

To assess the benefit of medication reconciliation on the evolution of medication prescriptions, we performed quantitative and qualitative retrospective analyses of the medication reconciliation process in the cardiology department from November 2017 to June 2021. A trained pharmaceutical team, composed of pharmacy students and a pharmacy resident, performed medication reconciliation within 24 to 72 hours after admission. Collected covariates were sociodemographic characteristics (sex, age), hospitalization data (hospitalization unit or ICU) and reconciliation data (number of prescribed medications, type of discrepancies, number of pharmaceutical interventions and corrected number of UMDs). As mentioned previously, UMDs are considered to be medication errors that can lead to adverse events for the patient, and need to be corrected; therefore, this analysis focused on unintentional discrepancies.

Continuous data are expressed as means ± standard deviations, and categorical variables are expressed as numbers and percentages. A P value < 0.001 was considered significant. A  2 test of independence was used to compare adjusted corrected UMD rates, and a  2 test of trend was used to study rate evolution, with R software (R Foundation for Statistical Computing, Vienna, Austria).

From November 2017 to June 2021, 6820 patients were hospitalized in the cardiology department. Following the prioritization criteria from the CACTUS study, 1164 patients were included in the medication reconciliation process, and 1361 medication reconciliations were realized: n = 89 in 2017 (2 months); n = 669 in 2018; n = 367 in 2019; n = 156 in 2020; and n = 80 in 2021 (6 months).

Among them, 1291 (94.9%) were successful, defined as achievement of: review of the BPMH based at least on three sources; comparison with the first hospital prescription; and discussion between physicians and pharmacists to classify discrepancies and correct them, if necessary. There were 70 (5.1%) unsuccessful reconciliations because of: discharge before BPMH achievement; no known medication at admission; fewer than three sources available.

The mean age was 73 ± 15.0 years, 679/1164 (58.3%) were male and 73% of patients were admitted to a hospitalization unit rather than an ICU. The average length of hospitalization was 13.6 ± 7.5 days. Hospitalization was mainly related to circulatory system diseases (81.2%): ischaemic heart disease; cerebrovascular disease; cardiopulmonary circulation disease; or disease of the arteries or veins. Other reasons for hospitalization were: follow-up of medical device implantation or investigation after treatment for other diseases (9.4%); and abnormal results of clinical or biological investigations (5.4%).

A total of 10,283 prescribed medications were analysed in the reconciliation process, with a mean of 7.6 ± 4.2 medications/patient. Discrepancies were observed for 4803/10,283 (46.7%) prescribed medications. Of these, 4448 (92.6%) prescribed medications were flagged as IMDs, and 355 (7.4%) as UMDs. At least one UMD was identified for 191 (16.4%) patients, with 1-12 UMDs per patient. The median number of UMDs per patient was 1.0 (interquartile range 1-2), with a mean of 1.8 ± 1.4 UMDs per patient. Among the 1291 successful medication reconciliations, we observed no UMDs for 1099 reconciliations (85.1%), one UMD for 100 (7.7%), two UMDs for 54 (4.2%), three UMDs for 25 (1.9%) and more than three UMDs for 13 (0.2%). After pharmaceutical intervention, 209/355 (58.9%) of UMDs were corrected. The corrected UMD rate was significantly different by year (P < 0.001), and increased significantly from 2017 to 2021 (P < 0.001; Fig. 1).

A successful medication reconciliation process was achieved for 18.9% of hospitalization stays in the cardiology from 2017 to 2021. Almost half of the prescribed medications had at least one discrepancy, and 7.4% were flagged as UMDs.

The medication reconciliation process corrected almost 60% of UMDs as a result of pharmaceutical interventions. The corrected UMD rate increased with time, and reached 87% in 2020. increase probably reflects an improvement in the relationship between cardiologists and pharmacists after years of collaboration.

Our work had some limitations. The study was performed in only two cardiology units (the hospital ward and the ICU) in the same hospital. Medication reconciliation could not be performed and achieved successfully for all patients. We observed that medication reconciliation inclusions varied from year to year; we also noticed decreases caused by a temporary lack of pharmacists or related to the COVID-19 pandemic. Patient care is influenced by the evolution of resources allocated to the clinical pharmacy. Because of the retrospective analysis, we were unable to collect some data, such as classification of medications involved in discrepancies. Cardiovascular treatments are the most likely to be involved in UMDs [START_REF] Fernandes | Medication Reconciliation at hospital admission and discharge: Evaluation of fidelity and process outcomes in a real-world setting[END_REF]. It would be pertinent to use the Anatomical, Therapeutic and Chemical (ATC) classification system, which classifies drugs according to organs and/or systems targeted by their therapeutic pharmacological effects. Using the ATC classification system could help to identify the types of medications and therapeutic indications that are concerned by UMDs, and also to grade UMD severity [START_REF] Marinovic | Impact of an integrated medication reconciliation model led by a hospital clinical pharmacist on the reduction of post-discharge unintentional discrepancies[END_REF].

There were also strengths to this study. The medication reconciliation process is based on a standardized method described by French authorities, and set up after a feasibility study (CACTUS) [START_REF] Chung | Medication reconciliation: Predictors of risk of unintentional medication discrepancies in the cardiology department[END_REF], which allowed an important number of medication reconciliations.

In conclusion, the medication reconciliation process allowed the identification of a huge number of discrepancies in medication prescriptions, with 17% of patients concerned at cardiology department admission. Pharmacists have a major responsibility in preventing medication errors [START_REF] Mekonnen | Effectiveness of pharmacist-led medication reconciliation programmes on clinical outcomes at hospital transitions: a systematic review and meta-analysis[END_REF], and are an essential part of the cardiology care team. The multidisciplinary collaboration of health actors (city and hospital) is crucial for the proper conduct of medication reconciliation.

Figures :Figure 1 .

 1 Figures :

Disclosure of interest

The authors declare that they have no competing interest.