Laurent Boudin
email: laurent.boudin@upmc.fr

AND Céline Grandmont
email: celine.grandmont@inria.fr

Ayman Moussa
email: ayman.moussa@upmc.fr

Global existence of solutions to the incompressible Navier-Stokes-Vlasov equations in a time-dependent domain

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

Introduction

The collective motion of a dispersed phase of small particles inside a fluid is often described using the so-called spray or fluid-kinetic models, first introduced in the combustion theory framework [START_REF] Williams | Combustion theory[END_REF] (see also [START_REF] O'rourke | Collective drop effects on vaporizing liquid sprays[END_REF][START_REF] Dufour | Modélisation multi-fluide eulérienne pour les écoulements diphasiques à inclusions dispersées[END_REF][START_REF] Desvillettes | Some aspects of the modeling at different scales of multiphase flows[END_REF]). In such models, one couples a kinetic equation with fluid mechanics equations. The fluid unknowns are the standard macroscopic quantities (mass density, velocity, for instance), and the dispersed phase is represented thanks to a distribution function. When one considers thin sprays, the fluid and kinetic equations are coupled through a drag term. This term depends on the fluid unknowns, the distribution function and their variables, and allows momentum and energy exchanges between both continuous and dispersed phases. Note that, with more dense (thick) sprays, the fluid volume fraction also appears in the equations, see [START_REF] Boudin | A modeling of compressible droplets in a fluid[END_REF] for instance.

The model investigated in this article, the incompressible Navier-Stokes-Vlasov system in a time-dependent domain, arises to study the transport and deposition of a therapeutic aerosol in a Newtonian, viscous, incompressible airflow, inside the human upper airways. The blueprints of the model were presented, in the aerosol therapy context, in [START_REF] Baranger | A modeling of biospray for the upper airways[END_REF][START_REF] Gemci | A numerical and experimental study of spray dynamics in a simple throat model[END_REF], and the extensive model was eventually written in [START_REF] Boudin | Modelling and numerics for respiratory aerosols[END_REF].

The fluid-kinetic models have been studied from a mathematical point of view for two decades or so. They depend on the physical phenomena we take into account, such as the fluid compressibility or viscosity, the particle transport, the interactions of the particles with the fluid or the wall, the time dependence of the domain, etc. We focus on systems where only mass and momentum are exchanged. The case with energy exchanges is, for instance, discussed in [START_REF] Boudin | Fluid-particles flows: a thin spray model with energy exchanges[END_REF].

In the compressible case, up to our knowledge, [START_REF] Baranger | Coupling Euler and Vlasov equations in the context of sprays: the local-in-time, classical solutions[END_REF] is the first contribution, for the coupled compressible Euler-Vlasov system, where Baranger and Desvillettes obtain the local-in-time existence of classical solutions. Mathiaud [START_REF] Mathiaud | Local smooth solutions of a thin spray model with collisions[END_REF] obtains the same kind of result for the Euler-Vlasov-Boltzmann system. For the compressible Navier-Stokes-Vlasov-Fokker-Planck system, in [START_REF] Mellet | Global weak solutions for a Vlasov-Fokker-Planck/Navier-Stokes system of equations[END_REF], Mellet and Vasseur prove global existence of weak solutions, and provide an asymptotic analysis for the system in [START_REF] Mellet | Asymptotic analysis for a Vlasov-Fokker-Planck/compressible Navier-Stokes system of equations[END_REF]. More recently, Chae, Kang and Lee [START_REF] Chae | Global classical solutions for a compressible fluid-particle interaction model[END_REF] and Li, Mu and Wang [START_REF] Li | Global well-posedness and large time behavior of strong solution to a kinetic-fluid model[END_REF] investigate, among other topics, the existence of global strong solutions close to the equilibrium.

There has been more mathematical contributions in the incompressible case, starting with [START_REF] Anoshchenko | The existence of the global generalized solution of the system of equations describing suspension motion[END_REF][START_REF] Hamdache | Global existence and large time behaviour of solutions for the Vlasov-Stokes equations[END_REF]. In [START_REF] Boudin | Global existence of solutions for the coupled Vlasov and Navier-Stokes equations[END_REF], the authors of the present article prove, with Desvillettes, the global existence of weak solutions to the incompressible Navier-Stokes-Vlasov equations in a periodic framework. The result is extended to the bounded domain case in [START_REF] Yu | Global weak solutions to the incompressible Navier-Stokes-Vlasov equations[END_REF]. The present work then appears as a natural continuation of those previous articles, with the additional difficulty of the time-dependent domain, which was not tackled yet, up to our knowledge. It also involves absorption boundary conditions for the distribution function, which are expected when dealing with depositing aerosols. Those boundary conditions induce extra difficulties to deal with, when compared to the periodic or whole-space cases. Note that Goudon, Jabin and Vasseur [START_REF] Goudon | Hydrodynamic limit for the Vlasov-Navier-Stokes equations. I. Light particles regime[END_REF][START_REF] Goudon | Hydrodynamic limit for the Vlasov-Navier-Stokes equations[END_REF] investigate the hydrodynamic limit of the Vlasov-Navier-Stokes system, for some particular regimes of the dispersed phase. We also mention a particular model [START_REF] Moussa | On a Vlasov-Euler system for 2D sprays with gyroscopic effects[END_REF] of inviscid fluid-particles, where the drag force is replaced by a lift force.

Various other problems in the incompressible case have also been investigated. Goudon, He, Moussa and Zhang [START_REF] Goudon | The Navier-Stokes-Vlasov-Fokker-Planck system near equilibrium[END_REF] establish the global existence of strong solutions near the equilibrium for the incompressible Navier-Stokes-Vlasov-Fokker-Planck system, whereas Carrillo, Duan and Moussa [START_REF] Carrillo | Global classical solutions close to equilibrium to the Vlasov-Fokker-Planck-Euler system[END_REF] studied the corresponding inviscid case. The global existence of weak and classical solutions for the Navier-Stokes-Vlasov-Fokker-Planck equations in a torus is investigated by Chae, Kang and Lee in [START_REF] Chae | Global existence of weak and classical solutions for the Navier-Stokes-Vlasov-Fokker-Planck equations[END_REF]. If the fluid is inhomogeneous, Wang and Yu [START_REF] Wang | Global weak solution to the inhomogeneous Navier-Stokes-Vlasov equations[END_REF] prove the existence of global weak solutions to the Navier-Stokes-Vlasov equations. Eventually, Benjelloun, Desvillettes and Moussa [START_REF] Benjelloun | Existence theory for the kinetic-fluid coupling when small droplets are treated as part of the fluid[END_REF] consider the formal asymptotic limit of the incompressible Vlasov-Navier-Stokes system with a fragmentation kernel, and prove the existence of global weak solutions for the resulting system, which shares some ties with the one considered in [START_REF] Wang | Global weak solution to the inhomogeneous Navier-Stokes-Vlasov equations[END_REF].

The article is organized as follows. In the next section, we present our problem in the time-dependent domain, state the main result, and discuss the proof strategy such as how we build an approximated problem to first deal with. Then, in Section 3, we recall some standard results on transport equation of Vlasov type in bounded domains and extend a result from [START_REF] Boyer | Mathematical tools for the study of the incompressible Navier-Stokes equations and related models[END_REF] to the case of a moving and unbounded domain. Section 4 is dedicated to the proof of existence of solutions to the approximated problem. Finally, in the last section, we go back to the whole problem by passing to the limit in the approximated one as the regularization parameters vanish, using, in particular, an Aubin-Lions lemma-like result from [START_REF] Moussa | Some variants of the classical Aubin-Lions Lemma[END_REF].

Presentation of the problem

We investigate the interaction of a size-monodispersed aerosol with a Newtonian, viscous and incompressible fluid. The evolution of this system is studied in a time-dependent domain.

Consider a finite constant T > 0, bounded an open set Ω ⊂ R 3 , and an open ball B such that Ω ⊂ B. We assume that Ω has a Lipschitz boundary.

Define

a mapping A ∈ C 2 (R + × R 3 ; R 3 ), (t, x) → A(t, x) = A t (x) such that, for each t ≥ 0, A t is a C 1 -diffeomorphism, A t = Id R 3 on B c , and A 0 = Id R 3 . The time-dependent bounded open domains Ω t , 0 ≤ t ≤ T , are obtained as Ω t = A t (Ω), and clearly satisfy Ω t ⊂ B for all t ∈ [0, T ]. We denote Ω = 0<t<T {t} × Ω t ⊂ B = (0, T ) × B, Γ = 0<t<T {t} × ∂Ω t .
Moreover, for all t ∈ [0, T ], we denote by n t the outgoing unit normal vector field of ∂Ω t . Finally, we introduce w, the Eulerian velocity associated to the flow t → A t , characterized by

∂ t A(t, x) = w(t, A t (x)), (t, x) ∈ [0, T ] × R 3 .
Note that the assumptions on A t imply that w ≡ 0 on (0, T ) × B c . Without loss of generality, we assume that div x w = 0, which is equivalent to assume that the Jacobian of the transformation A t does not depend on t. Besides, we need to introduce the following phase-space boundaries for the aerosol

Σ = Γ × R 3 , Σ ± = {(t, x, ξ) ∈ Σ | ± (ξ -w(t, x)) • n t (x) > 0}, Σ 0 = {(t, x, ξ) ∈ Σ | (ξ -w(t, x)) • n t (x) = 0},
and, for any t,

Σ ± t = {(x, ξ) ∈ ∂Ω t × R 3 | ± (ξ -w(t, x)) • n t (x) > 0}.
The fluid is described via macroscopic quantities, pressure p(t, x) and velocity u(t, x), thanks to the classical incompressible Navier-Stokes equation, with constant density and viscosity both chosen equal to 1,

(2)

∂ t u + u • ∇ x u + ∇ x p -∆ x u = F in Ω, div x u = 0 in Ω,
where F is the action of the aerosol on the fluid that will be defined later on. We prescribe the following Dirichlet boundary conditions

(3) u = w in Γ.
The aerosol is described through a density function f (t, x, ξ), which solves a Vlasov-like equation ( 4)

∂ t f + ξ • ∇ x f + ∇ ξ • [Af ] = 0 in Ω × R 3 ,
where A is the drag acceleration exerted by the fluid on the particles. We prescribe the following absorption boundary conditions for the aerosol:

(5)

f = 0 in Σ -.
The coupling terms A and F, which respectively depend on (t, x, ξ) and (t, x), are defined by ( 6)

A = u -ξ, F = - R 3 f A.
The system is eventually supplemented with initial conditions for u and f which are

(7) (8)
u(0, •) = u in , in Ω, f (0, •, •) = f in , in Ω × R 3 .
As usual in the framework of fluid-kinetic coupling, we may compute, at least formally, an energy equality reflecting the dissipation of energy and the exchange between the two phases. Recall Reynolds' formula for any real-valued function k : Ω → R:

d dt Ωt k = Ωt ∂ t k + ∂Ωt k w • n t .
We multiply equation ( 1) by uw and integrate over Ω t (t is fixed). After integration by parts, using (2), (3), [START_REF] Boudin | Fluid-particles flows: a thin spray model with energy exchanges[END_REF] and the Reynolds formula, we get

1 2 d dt Ωt |u| 2 + Ωt |∇ x u| 2 = d dt Ωt u • w - Ωt u • ∂ t w - Ωt (u • ∇ x w) • u + Ωt ∇ x u : ∇ x w - Ωt×R 3 (u -w) • (u -ξ)f.
We then multiply (4) by |ξ| 2 /2 and integrate on Ω t × R 3 to get after integration by parts, using (6),

1 2 d dt Ωt×R 3 f |ξ| 2 = Ωt×R 3 ξ • (u -ξ)f + 1 2 ∂Ωt×R 3 f |ξ| 2 (w -ξ) • n t .
Summing both previous equalities, we obtain the following energy balance

dE u,f dt + Ωt |∇ x u| 2 + Ωt×R 3 f |u -ξ| 2 = d dt Ωt u • w - Ωt u • ∂ t w - Ωt (u • ∇ x w) • u + Ωt ∇ x u : ∇ x w + Ωt×R 3 w • (u -ξ)f + 1 2 ∂Ωt×R 3 f |ξ| 2 (w -ξ) • n t ,
where E u,f denotes the total kinetic energy of the system and writes

E u,f (t) = 1 2 Ωt |u| 2 + 1 2 Ωt×R 3 f |ξ| 2 .
Assuming f to be nonnegative, we note that

Ωt×R 3 w • (u -ξ)f ≤ 1 2 Ωt×R 3 f |u -ξ| 2 + 1 2 Ωt×R 3 f |w| 2 .
Thanks to the absorption boundary condition (5), we may eventually write

dE u,f dt + Ωt |∇ x u| 2 + 1 2 Ωt×R 3 f |u -ξ| 2 ≤ d dt Ωt u • w - Ωt u • ∂ t w - Ωt (u • ∇ x w) • u + Ωt ∇ x u : ∇ x w + 1 2 Ωt×R 3 f |w| 2 .
We deduce, thanks to Young's inequality, using the regularity of w, and integrating over [0, t],

(9) E u,f (t) + t 0 Ωs |∇ x u| 2 + t 0 Ωs×R 3 f |u -ξ| 2 ≤ C E u,f (0) + C w + C w t 0 Ωs |u| 2 + C w t 0 Ωs×R 3 f,
where C > 0 is a universal constant and C w > 0 is a constant depending on the domain velocity w. Moreover, integrating (4) on Ω t × R 3 , we get, using once again [START_REF] Benjelloun | Existence theory for the kinetic-fluid coupling when small droplets are treated as part of the fluid[END_REF],

d dt Ωt×R 3 f = ∂Ωt×R 3 f (w -ξ) • n ≤ 0,
and consequently,

Ωt×R 3 f ≤ Ω×R 3 f in .
All in all, going back to [START_REF] Boudin | Modelling and numerics for respiratory aerosols[END_REF] we have, for some constant C w,f in > 0 depending on w and f in ,

E u,f (t) + t 0 Ωs |∇ x u| 2 + t 0 Ωs×R 3 f |u -ξ| 2 ≤ C E u,f (0) + C w,f in T + t 0 Ωs |u| 2 ,
so that, thanks to Gronwall's lemma, we infer

(10) E u,f (t) + t 0 Ωs |∇ x u| 2 + t 0 Ωs×R 3 f |u -ξ| 2 ≤ (C E u,f (0) + C w,f in T )e T 2 C w,f in .
The previous energy inequality, together with equation (2) motivates the introduction of several function spaces. Set, for the fluid part,

L 2 (0, T ; H 1 (Ω t )) = {ψ ∈ L 2 ( Ω) | ∇ x ψ ∈ L 2 ( Ω)}, V = {ψ ∈ L 2 (0, T ; H 1 (Ω t )) | div x ψ = 0}, V 0 = {ψ ∈ V | ψ = 0 on Γ}, V = {ψ ∈ C 1 ( Ω ) | div x ψ = 0 on Ω, ψ = 0 on Γ}.
Note that V 0 is the closure of V in L 2 (0, T ; H 1 (Ω t )). For any function (scalar or vectorvalued) g defined on Ω, we denote by g the extension (in the space variable) by 0 of g on (0, T ) × R 3 . For the kinetic part, we now set, for any p ∈ [1, ∞],

L ∞ (0, T ; L p (Ω t × R 3 )) = {f measurable | f ∈ L ∞ (0, T ; L p (R 3 × R 3 ))}, W = {φ ∈ C 1 c ( Ω × R 3 ) | φ = 0 on Σ + ∪ Σ 0 }. For any function h : Ω × R 3 → R + and α ∈ R + ,
we further introduce the following notations on the moments of h which will be useful in the study of the Vlasov equation. For a.e. (t, x), we set

m α h(t, x) = R 3 h |ξ| α , M α h(t) = Ωt×R 3 h |ξ| α = Ωt m α h.
In particular, we recall a standard interpolation estimate, see [START_REF] Hamdache | Global existence and large time behaviour of solutions for the Vlasov-Stokes equations[END_REF] for instance:

Lemma 2.1. Let β > 0, and h be a nonnegative function in L ∞ ( Ω × R 3 ) such that m β h(t,
m α h(t, x) ≤ 4 3 π h L ∞ ( Ω×R 3 ) + 1 m β h(t, x) α+3 β+3 .
Let us now give the assumptions on the initial data and the domain motion:

Assumption 2.1. u in ∈ L 2 (Ω), div x u in = 0, and u in • n 0 = w(0, •) • n 0 on ∂Ω. Assumption 2.2. f in is nonnegative, f in ∈ L ∞ (Ω × R 3 ), f in (1 + |ξ| 2 ) ∈ L 1 (Ω × R 3 ).
Weak solutions of the coupled problem are defined as follows.

Definition 2.1. We say that a couple (u, f ) is a weak solution of system (1)-( 5) with initial datum (u in , f in ) if the following conditions are satisfied:

• u ∈ L ∞ (0, T ; L 2 (R 3 )), • u -w ∈ V 0 , • f ∈ L ∞ (0, T ; L ∞ (R 3 × R 3 )) ∩ C 0 ([0, T ]; L p (R 3 × R 3 )), for all p ∈ [1, ∞), • f (1 + |ξ| 2 ) ∈ L ∞ (0, T ; L 1 (R 3 × R 3 )), • for all ψ ∈ V , for a.e. t ∈ [0, T ], Ωt u(t) • ψ(t) - t 0 Ωs u • (∂ t ψ + u • ∇ x ψ) - t 0 Ωs ∇ x u : ∇ x ψ - t 0 Ωs u • R 3 (ξ -u)f = Ω u in • ψ(0),
• for all φ ∈ W , for all t ∈ [0, T ],

Ωt×R 3 f (t)φ(t) - t 0 Ωs×R 3 f (∂ t φ + ξ • ∇ x φ + (u -ξ) • ∇ ξ φ) = Ω×R 3 f in φ(0).
The main result of this paper is the following Theorem 2.1. Under Assumptions 2.1-2.2, there exists at least one weak solution (u, f ) to (1)-( 8) in the sense of Definition 2.1.

The proof follows a standard scheme: introduction of an approximated system, itself solved through a fixed-point procedure.

Let us explain how we obtain the approximated system. First, we regularize the convection term in the Navier-Stokes equations. Then, in order to work in the cylindrical domain B containing the non-cylindrical domain Ω, we add a penalty term in (1), as in [START_REF] Fujita | On existence of weak solutions of the Navier-Stokes equations in regions with moving boundaries[END_REF]. Since the relative velocity term uξ governs the coupling of the system through [START_REF] Boudin | Fluid-particles flows: a thin spray model with energy exchanges[END_REF], to be able to apply standard existence results for the Navier-Stokes equations, we are led to truncate the right hand side of [START_REF] Anoshchenko | The existence of the global generalized solution of the system of equations describing suspension motion[END_REF]. Moreover, to preserve the energy estimate of the coupled approximated system, the same truncation is performed in [START_REF] Bardos | Problèmes aux limites pour les équations aux dérivées partielles du premier ordre à coefficients réels; théorèmes d'approximation; application à l'équation de transport[END_REF]. We thus introduce an odd, increasing and bounded function χ ∈ C ∞ (R), satisfying 0 ≤ χ(z) ≤ z for any z ≥ 0, and we use the abuse of notation χ(z) = (χ(z 1 ), χ(z 2 ), χ(z 3 )), when z ∈ R 3 . Hence, the resulting problem system writes [START_REF] Boyer | Mathematical tools for the study of the incompressible Navier-Stokes equations and related models[END_REF] (12) (13)

∂ t u + (u ⋆ ϕ) • ∇ x u + ∇ x p -∆ x u + λ1 B\ Ω (u -w) = R 3 f χ(ξ -u) in B, div x u = 0 in B, ∂ t f + ξ • ∇ x f + ∇ ξ • [χ(u -ξ)f ] = 0 in Ω × R 3 .
In [START_REF] Boyer | Mathematical tools for the study of the incompressible Navier-Stokes equations and related models[END_REF], λ > 0 is the penalty parameter and ϕ denotes an element of D(R 3 ), the convolution product of the convection velocity being made in x only. The previous system is completed with the following boundary conditions: homogeneous Dirichlet boundary conditions for u on (0, T ) × ∂B and ( 5) for f . The initial data are also extended and/or regularized, to ensure that, at least, they satisfy the following assumptions.

Assumption 2.3. u in ∈ H 1 0 (B) and div x u in = 0, Assumption 2.4. f in is nonnegative, has a compact support in ξ, f in ∈ L ∞ (Ω × R 3 ) ∩ L 1 (Ω × R 3 ).
The existence of weak solutions to ( 11)-( 13) is obtained thanks to Schaefer's fixed point theorem (for the proof, see, for instance, [21, Theorem 11.6, p.286]) which we here recall, for the sake of completeness.

Theorem 2.2. Let E be a Banach space and

Θ : E × [0, 1] → E a continuous mapping sending bounded subsets of E × [0, 1] on relatively compact subsets of E. Denoting Θ σ = Θ(•, σ), if Θ 0 =
0 and the set of all the fixed points of the family

(Θ σ ) σ∈[0,1] is bounded in E, then Θ 1 has at least one fixed point in E.
We choose E = L 2 (0, T ; H 1 0 (B)), and define the mapping Θ : E × [0, 1] → E in the following way. Starting from u ∈ E, we define f u as the unique weak solution of (13), with absorption boundary conditions and initial datum f in satisfying Assumption 2.4, and we set [START_REF] Chae | Global classical solutions for a compressible fluid-particle interaction model[END_REF] S χ (u) =

R 3 f u χ(ξ -u).
Then, for each σ ∈ [0, 1], consider u σ the unique solution of (

(18)

∂ t u σ + ( u σ ⋆ ϕ) • ∇ x u σ + ∇ x p -∆ x u σ = σS χ (u) -λσ1 B\ Ω (u -w) in B, ∇ x • u σ = 0 in B, u σ = 0 on (0, T ) × ∂B, u σ (0) = σu in in B,
where u in satisfies Assumption 2.3. Eventually, Theorem 2.2 is used for the mapping Θ :

E × [0, 1] → E, (u, σ) → u σ .
Now that we have introduced the approximated problems, let us describe more precisely the main steps and the difficulties in the existence proof. We first deal with the Vlasov equation ( 13) with absorption boundary conditions, and a non smooth drag acceleration in a time-dependent domain. Consequently, the solution to [START_REF] Chae | Global existence of weak and classical solutions for the Navier-Stokes-Vlasov-Fokker-Planck equations[END_REF], together with its trace on the boundary Σ -, can be defined in a weak sense, following the arguments developped in [START_REF] Boyer | Trace theorems and spatial continuity properties for the solutions of the transport equation[END_REF][START_REF] Boyer | Mathematical tools for the study of the incompressible Navier-Stokes equations and related models[END_REF] and extended to our time-dependent domain setting. In particular, the use of renormalization theory is crucial to ensure uniqueness of solutions, and strongly relies on the regularity of χ(u-ξ). The study of the incompressible Navier-Stokes system in a timedependent domain is more standard, see for instance [START_REF] Fujita | On existence of weak solutions of the Navier-Stokes equations in regions with moving boundaries[END_REF][START_REF] Salvi | On the existence of weak solutions of a nonlinear mixed problem for the Navier-Stokes equations in a time dependent domain[END_REF]. Nevertheless, because of the time dependence of the domain, a special care to apply or adapt the Aubin-Lions lemma is required to prove the needed compactness on the fluid velocity field and to take the limit in the sequence of approximated problems as the regularization parameters vanish.

Here we choose to apply a result obtained in [START_REF] Moussa | Some variants of the classical Aubin-Lions Lemma[END_REF]. Note finally that, as in [START_REF] Anoshchenko | The existence of the global generalized solution of the system of equations describing suspension motion[END_REF][START_REF] Chae | Global existence of weak and classical solutions for the Navier-Stokes-Vlasov-Fokker-Planck equations[END_REF], the approximation strategy preserves the energy estimate satisfied by the full coupled system. In particular, it ensures that the set of fixed points of the mappings Θ σ is bounded in E.

Vlasov's equation in a moving domain

We here consider the following Vlasov equation

(19) (20) (21) ∂ t f + ξ • ∇ x f + ∇ ξ • (f G) = 0 in Ω × R 3 , f = 0 in Σ -, f (0, •) = f 0 in Ω × R 3 ,
where the vector field G(t, x, ξ) and the initial datum f 0 (x, ξ) are given. In this section, we aim to obtain an appropriate functional setting in order to define weak solutions to ( 19)-( 21) and prove their existence and uniqueness. For uniqueness and to give a weak sense to the trace on Σ, we closely follow the approach developped by Boyer in [START_REF] Boyer | Trace theorems and spatial continuity properties for the solutions of the transport equation[END_REF] and adapt it to our time-dependent domain. As for the existence, we follow the standard characteristic method together with a regularization procedure of the field G.

The first subsection is devoted to the definition and existence of weak solutions to ( 19)- [START_REF] Gilbarg | Elliptic partial differential equations of second order, volume 224 of Grundlehren der Mathematischen Wissenschaften[END_REF]. Next, we prove that the trace on Σ of such solutions can properly be defined together with a renormalized weak formulation satisfied by any bounded weak solution. Uniqueness is then obtained, strongly relying on the renormalized weak formulation satisfied by the trace.

Existence of solution. We only consider the case of bounded solutions and assume that

G ∈ L 1 (0, T ; W 1,1 loc (B × R 3 )), div ξ G ∈ L ∞ (0, T ; L ∞ (B × R 3 )) and f 0 ∈ L 1 (Ω × R 3 ) ∩ L ∞ (Ω × R 3
). We introduce the following definition:

Definition 3.1. A function f ∈ L ∞ ( Ω × R 3
) is a weak solution of (19)-( 21) if, for any test function φ ∈ W such that φ(T, •) = 0, one has [START_REF] Goudon | The Navier-Stokes-Vlasov-Fokker-Planck system near equilibrium[END_REF] -

T 0 Ωt×R 3 f (∂ t φ + ξ • ∇ x φ + G • ∇ ξ φ) = Ω×R 3 f 0 φ(0).
Remark 3.1. Note that the previous kind of solution is weaker than the one introduced for f in Definition 2.1. Indeed, the latter included some continuity in time and a larger set of test functions. Nevertheless, the arguments of Subsection 3.2 enable to prove that both definitions are in fact equivalent.

Let us first consider a vector field G defined in the whole space and regular enough, namely G ∈ C 0 (R; W 1,∞ (R 3 × R 3 )). With this regularity, the Cauchy-Lipschitz theorem allows us to define global characteristic curves associated to the field (ξ, G). More precisely, for any s ∈ R and any (x, ξ) ∈ R 3 × R 3 , there exists a unique global solution of the differential system

             dx s dt (t) = ξ s (t), dξ s dt (t) = G(t, x s (t), ξ s (t)), (x s (s), ξ s (s)) = (x, ξ).
We denote the associated curve by t → T s,t (x, ξ) = (x s (t), ξ s (t)). For (t, x, ξ) ∈ Ω × R 3 , let us then consider the retrograde outgoing time

τ -(t, x, ξ) = inf{s ≤ t | ∀σ ∈ [s, t], T t,σ (x, ξ) ∈ Ω σ × R 3 }.
In that case, we can define, for (t, x, ξ)

∈ Ω × R 3 , ( 23
) f (t, x, ξ) = e h(t,x,ξ) f 0 (T t,0 (x, ξ))1 τ -<0 (t, x, ξ),
where h is the solution in (0, T )×R 3 ×R 3 of the following transport equation with a source term

∂ t h + ξ • ∇ x h + G • ∇ ξ h = div ξ G
with initial condition h(0, •) = 0. One can check, thanks to quite straightforward though tedious computations, that the function defined in ( 23) is a weak solution of ( 19)- [START_REF] Gilbarg | Elliptic partial differential equations of second order, volume 224 of Grundlehren der Mathematischen Wissenschaften[END_REF] in the sense of Definition 3.1. The existence of weak solutions to the transport equation with absorption boundary conditions based on formula ( 23) is investigated in [START_REF] Bardos | Problèmes aux limites pour les équations aux dérivées partielles du premier ordre à coefficients réels; théorèmes d'approximation; application à l'équation de transport[END_REF] when the domain is bounded and the field G is regular using the semigroup theory, whereas the Vlasov equation case is studied in [START_REF] Mischler | On the trace problem for solutions of the Vlasov equation[END_REF] when G is not regular, but only depends on t and x. Note that Boyer and Fabrie [START_REF] Boyer | Mathematical tools for the study of the incompressible Navier-Stokes equations and related models[END_REF] follow another strategy to obtain the existence of weak solutions, which is based on a parabolic regularization of the transport equation.

The function defined in ( 23) furthermore satisfies

f L ∞ (0,T ;L 1 (R 3 ×R 3 )∩L ∞ (R 3 ×R 3 )) ≤ exp T div ξ G L ∞ (0,T ;L ∞ (B×R 3 )) f 0 L 1 (Ω×R 3 )∩L ∞ (Ω×R 3 ) . If G ∈ L 1 (0, T ; W 1,1 loc (B×R 3 )), with div ξ G ∈ L ∞ (0, T ; L ∞ (B×R 3 )), it can be approximated in L 1 loc ( B × R 3 ) by a sequence (G k ) lying in D([0, T ] × R 3 × R 3 ). Furthermore, (G k ) can be chosen such that div ξ G k L ∞ (0,T ;L ∞ (R 3 ×R 3 )) ≤ div ξ G L ∞ (0,T ;L ∞ (B×R 3 )) .
Letting k go to +∞, we obtain the following

Proposition 3.1. If G ∈ L 1 (0, T ; W 1,1 loc (B × R 3 )) is such that div ξ G ∈ L ∞ (0, T ; L ∞ (B × R 3 )), then, for any f 0 ∈ L 1 (Ω × R 3 ) ∩ L ∞ (Ω × R 3 )
, there exists at least one weak solution to the system (19)-( 21) in the sense of Definition 3.1. Furthermore, this solution satisfies

(24) f L ∞ (0,T ;L 1 (R 3 ×R 3 )∩L ∞ (R 3 ×R 3 )) ≤ exp T div ξ G L ∞ (0,T ;L ∞ (B×R 3 )) f 0 L 1 (Ω×R 3 )∩L ∞ (Ω×R 3 ) .
Remark 3.2. It is clear from the previous proof that the assumption G ∈ L 1 (0, T ; W 1,1 loc (B× R 3 )) can be replaced by G ∈ L 1 (0, T ; L 1 loc (B×R 3 )). However, the stronger assumption on G is crucial in the next subsection, where we define the trace of the solution on the boundary.

Remark 3.3. If G ∈ L 1 (0, T ; L ∞ (B × R 3 )
) and if f 0 has a compact support in ξ, then the weak solution f to (19)-( 21) also has a compact support in ξ. It is clear in the regularized case, when f is given by [START_REF] Goudon | Hydrodynamic limit for the Vlasov-Navier-Stokes equations. I. Light particles regime[END_REF]. Moreover, the support of f in ξ can be chosen independently from the regularization parameter thanks to the assumption G ∈ L 1 (0, T ; L ∞ (B × R 3 )). Consequently, the property still holds at the limit. It will be useful to obtain a priori estimates for the regularized coupled problem involving moments of f .

Trace theorem.

Let us first state a trace theorem which holds for bounded solutions of transport equations in bounded domains. A first version of this theorem has been established by Boyer [10, Theorem 3.1], we here use the more general version presented in [START_REF] Boyer | Mathematical tools for the study of the incompressible Navier-Stokes equations and related models[END_REF]. Note nevertheless that we do not state that theorem in its whole generality. We choose to adapt it to our framework.

Theorem 3.1. Consider a Lipschitz bounded open set

O ⊂ R d . Assume that a ∈ L 1 (0, T ; W 1,1 (O)) such that div a ∈ L 1 (0, T ; L ∞ (O)). Let ρ ∈ L ∞ ((0, T ) × O) be a

distributional solution of the following transport equation

∂ t ρ + div(aρ) = 0, in (0, T ) × O.
Then ρ ∈ C 0 ([0, T ]; L p (O)) for all p ∈ [1, ∞), and admits a well-defined trace γρ on (0, T ) × ∂O: it is the only element of L ∞ ((0, T ) × ∂O, |dµ a |) with dµ a = a • n dσ dt, which satisfies, for any test function ϕ ∈ D(R × R d ), any real-valued function β ∈ C 1 (R), and all t 0 , t 1 such that 0 ≤ t 0 ≤ t 1 ≤ T , (25)

t 1 t 0 O β(ρ)(∂ t ϕ + a • ∇ϕ)) - t 1 t 0 O (ρβ ′ (ρ) -β(ρ))ϕ div a = O β(ρ(t 1 ))ϕ(t 1 ) - O β(ρ(t 0 ))ϕ(t 0 ) + t 1 t 0 ∂O β(γρ) ϕ a • n.
Here we need to generalize this result to our non-cylindrical and unbounded setting. In this case, the theorem we obtain writes

Theorem 3.2. Assume that f ∈ L ∞ ( Ω × R 3 ) is a solution of (26) ∂ t f + ξ • ∇ x f + div ξ (Gf ) = 0 in D ′ ( Ω × R 3 ), where G ∈ L 1 (0, T ; W 1,1 loc (B × R 3 )). Then, for all p ∈ [1, ∞), f ∈ C 0 ([0, T ]; L p loc (R 3 × R 3 )). Moreover the trace of f on Γ is well-defined: it is the unique element γf ∈ L ∞ ( Γ) satisfying for any test function ψ ∈ D(R × R 3 × R 3 ), all real-valued function β ∈ C 1 (R), and all 0 ≤ t 0 ≤ t 1 ≤ T : (27) t 1 t 0 Ωt×R 3 β(f ) (∂ t ψ + ξ • ∇ x ψ + G • ∇ ξ ψ) - t 1 t 0 Ωt×R 3 (f β ′ (f ) -β(f ))ψ div ξ G = Ωt 1 ×R 3 β(f (t 1 ))ψ(t 1 ) - Ωt 0 ×R 3 β(f (t 0 ))ψ(t 0 ) + t 1 t 0 ∂Ωt×R 3 β( γf ) ψ (ξ -w) • n t .
Remark 3.4. Following [11, Corollary VI.1.5], β can be chosen continuous and piecewise C 1 . Indeed, from [START_REF] Li | Global well-posedness and large time behavior of strong solution to a kinetic-fluid model[END_REF], we can prove that, for any α and for any test function

ψ ∈ D(R × R 3 × R 3 ), Ωt×R 3 ψ 1 f =α f div ξ G = 0.
Remark 3.5. Note that, in the boundary term [START_REF] Lions | Quelques méthodes de résolution des problèmes aux limites non linéaires[END_REF], the velocity w appears because of our time-dependent domain setting. However, the case w = 0, which corresponds to a cylindrical setting, does not fall directly into the scope of Boyer's theorem because Ω × R 3 is not bounded. Consequently, one of the main step of the proof will be to generalize Theorem 3.

1 for O = Ω × R 3 .
Proof. Since we deal with a time-dependent domain, we perform a change of variables in order to work in a cylindrical domain. We associate to f a function g defined in the reference domain as follows:

(28) g(t, y, ξ) = f (t, A t (y), ξ), a.e. (t, y, ξ)

∈ (0, T ) × Ω × R 3 . The fact that f ∈ L ∞ ( Ω × R 3 ) implies that g ∈ L ∞ ((0, T ) × Ω × R 3
). In the same way as in (28), we set

K(t, y, ξ) = G(t, A t (y), ξ), a.e. (t, y, ξ) ∈ (0, T ) × Ω × R 3 .
Since f solves (26), we have

∂ t g + (ξ -v) • (Cof(∇ y A t )∇ y ) g + div ξ (Kg) = 0 in D ′ ((0, T ) × Ω × R 3 ),
where v is the domain Lagrangian velocity, i.e.

v(t, y) = w(t,

A t (y)), for any (t, y) ∈ (0, T ) × Ω, and Cof(∇ y A t ) denotes the cofactor matrix of ∇ y A t . Thanks to the additional assumption div x w = 0 on the domain Eulerian velocity w, we have div y (Cof(∇ y A t ) ⊺ v) = 0. Moreover, Piola's identity ensures that div y (Cof(∇ y A t )) = 0 (see [15]). Consequently, introducing the new variable z = (y, ξ) ∈ Ω × R 3 and considering the vector field ( 29)

a(t, z) = ((Cof(∇ y A t )) ⊺ (ξ -v), K) ⊺ ,
the previous equation may be written, in

D ′ ((0, T ) × Ω × R 3 ), ( 30
) ∂ t g + div z (ag) = 0.
We cannot directly apply Theorem 3.

1 because Ω × R 3 is not bounded. Let us use Boyer's theorem for O = O R , where R > 0, O R = Ω × B R and B R
is the open ball of radius R > 0 centred at 0. It is possible because g also solves [START_REF] Mellet | Asymptotic analysis for a Vlasov-Fokker-Planck/compressible Navier-Stokes system of equations[END_REF] ∞) and for all R > 0, so that g ∈ C 0 ([0, T ]; L p loc (Ω × R 3 )). Consequently, going back to the non-cylindrical domain, we get f ∈ C 0 ([0, T ]; L p loc (R 3 × R 3 )). Next, we would like to define the trace of g on (0, T ) × ∂Ω × R 3 . Using, once again, Theorem 3.1 in (0, T )×O R , we get existence and uniqueness of the trace γ R g in L ∞ ((0,

in D ′ ((0, T ) × O R ) and g ∈ L ∞ ((0, T ) × O R ). We thus infer that g ∈ C 0 ([0, T ]; L p (O R )) for p ∈ [1,
T )× ∂O R , |dµ a,R |) such that (25) holds with ρ = g and O = O R . Note that ∂O R = (∂Ω × B R ) ∪ (Ω × ∂B R ).
We restrict ourselves to the space part of the boundary ∂Ω × B R , on which dµ a,R does not depend on R and equals [START_REF] Mischler | On the trace problem for solutions of the Vlasov equation[END_REF] dµ a = (n 0 , 0

) ⊺ • a dσ dξ dt = (ξ -v) • Cof(∇ y A t )n 0 dσ dξ dt.
The uniqueness of the trace implies that, for any R ′ > R,

(γ R g) |(0,T )×∂Ω×B R = (γ R ′ g) |(0,T )×∂Ω×B R in L ∞ ((0, T ) × ∂Ω × B R , |dµ a |).
Since dσ dξ dt-almost everywhere, we have (ξ

-v) • Cof(∇ y A t )n 0 = 0 on (0, T ) × ∂Ω × B R , we have L ∞ ((0, T ) × ∂Ω × B R , |dµ a |) ⊂ L ∞ ((0, T ) × ∂Ω × B R
) and the previous equality also holds in the latter space. It is hence possible to define, without ambiguity, an element

γg ∈ L ∞ loc ([0, T ] × ∂Ω × R 3 ) by γg = γ R g on [0, T ] × ∂Ω × B R , ∀R > 0.
From the trace of g on the boundary of the reference domain, we define the element γf of L ∞ loc ( Γ) by γg(t, y, ξ) = γf (t, A t (y), ξ).

The next step consists in proving the weak renormalized formulation [START_REF] Lions | Quelques méthodes de résolution des problèmes aux limites non linéaires[END_REF], which ensures that γf is, in fact, the trace of f on Γ. We first check that g satisfies the same type of equality as ( 25) on (0, T ) × Ω × R 3 . Note that, for any ϕ ∈ D(R × R 3 × R 3 ), the support of ϕ lies in fact in R 3 × R 3 × B R , and thus we get [START_REF] Moussa | Some variants of the classical Aubin-Lions Lemma[END_REF]

t 1 t 0 Ω×R 3 β(g)(∂ t ϕ + a • ∇ z ϕ) - t 1 t 0 Ω×R 3 (gβ ′ (g) -β(g))ϕ div z a = Ω×R 3 β(g(t 1 ))ϕ(t 1 ) - Ω×R 3 β(g(t 0 ))ϕ(t 0 ) + t 1 t 0 ∂Ω×R 3 β(γg)ϕ (ξ -v) • a • (n 0 , 0) ⊺ .
From the definition (29) of a, we obtain div z a = div ξ K. Thus, taking (31) into account, (32) writes ( 33)

t 1 t 0 Ω×R 3 β(g)(∂ t ϕ + ξ • (Cof(∇ y A t )∇ y )ϕ + K • ∇ ξ ϕ) - t 1 t 0 Ω×R 3 (gβ ′ (g) -β(g))ϕ div ξ K = Ω×R 3 β(g(t 1 ))ϕ(t 1 ) - Ω×R 3 β(g(t 0 ))ϕ(t 0 ) + t 1 t 0 ∂Ω×R 3 β(γg)ϕ (ξ -v) • Cof(∇ y A t )n 0 .
Thanks to a density argument, ( 32) is still satisfied if ϕ is only

C 1 c (R × R 3 × R 3
). The final step consists in going back to the time-dependent domain in order to define the trace of f and altogether obtain the weak formulation [START_REF] Lions | Quelques méthodes de résolution des problèmes aux limites non linéaires[END_REF]. For any fixed t ≥ 0 and any given function ζ defined on ∂Ω t , thanks to [15], and with obvious notations, we have (whenever the integrals make sense)

∂Ω (ζ • A t ) Cof(∇ y A t )n 0 dσ = ∂Ωt ζ n t dσ t .
Consequently, since g is the image of f by the mapping A t , from [START_REF] Moussa | On a Vlasov-Euler system for 2D sprays with gyroscopic effects[END_REF], we easily obtain [START_REF] Lions | Quelques méthodes de résolution des problèmes aux limites non linéaires[END_REF] by applying the change of variables x = A t (y). At this stage, this weak formulation ( 27) is satisfied for any function ψ such that ϕ(t, y, ξ) = ψ(t, A t (y, ξ)),

with ϕ ∈ C 1 c (R×R 3 ×R 3 ). Since A t is a C 1 -diffeomorphism of R 3 , (27) actually holds for any test function ψ ∈ C 1 c (R × R 3 × R 3 ).
To end the proof of the theorem, we just have to check that the trace γf belongs to L ∞ ( Γ). In the same way as in [START_REF] Boyer | Trace theorems and spatial continuity properties for the solutions of the transport equation[END_REF], we choose an appropriate renormalization function

β. Since f ∈ L ∞ ( Ω×R 3 ), there exist two real numbers c 1 , c 2 such that c 1 ≤ f ≤ c 2 a.e. on Ω × R 3 . We choose β ∈ C 1 (R) such that β(x) = 0 if and only if x ∈ [c 1 , c 2 ]. For such a β, the weak formulation (27) leads to Γ β( γf ) ψ (ξ -w) • n t = 0, ∀ψ ∈ C 1 c (R × R 3 × R 3 ).
This implies that β( γf )(ξw) • n t = 0 a.e. on Γ. But (ξw) • n t = 0 a.e. (this is the counterpart of the already used property (ξv) • Cof(∇ y A t )n 0 = 0 a.e.), so that we get β( γf ) = 0 a.e., and thus c 1 ≤ γf ≤ c 2 a.e. Finally, the uniqueness of γf in L ∞ ( Ω × R 3 ) is straightforward.

Remark 3.6. If f ≥ 0 almost everywhere, following the above argument and taking c 1 = 0, we can prove that γf ≥ 0 almost everywhere.

Remark 3.7. Note that, thanks to the just proved renormalized weak formulation [START_REF] Lions | Quelques méthodes de résolution des problèmes aux limites non linéaires[END_REF], any weak solution f of ( 19)-( 21) in the sense of Definition 3.1 satisfies γf = 0 on Σ -.

Uniqueness.

In this subsection, we prove uniqueness of the weak solution of ( 19)- [START_REF] Gilbarg | Elliptic partial differential equations of second order, volume 224 of Grundlehren der Mathematischen Wissenschaften[END_REF]. The following result holds.

Proposition 3.2. Consider a weak solution

f ∈ L ∞ (0, T ; L 1 (R 3 ×R 3 ))∩L ∞ (0, T ; L ∞ (R 3 × R 3
)) to the system (19)-( 21) in the sense of Definition 3.1 with

f 0 ≡ 0. Assume that G ∈ L 1 (0, T ; W 1,1 loc (B × R 3 )), div ξ G ∈ L ∞ (0, T ; L ∞ (B × R 3 )), and
(34) G 1 + |ξ| ∈ L 1 (0, T ; L 1 (B × R 3 )) + L 1 (0, T ; L ∞ (B × R 3 )).
Then f ≡ 0.

Proof. The proof follows the strategy developped in [START_REF] Diperna | Ordinary differential equations, transport theory and Sobolev spaces[END_REF] and later used for instance in [START_REF] Mischler | On the trace problem for solutions of the Vlasov equation[END_REF][START_REF] Boyer | Trace theorems and spatial continuity properties for the solutions of the transport equation[END_REF]. The main idea is to use the renormalized weak formulation obtained in Theorem 3.2.

Let f be a bounded weak solution of ( 19)- [START_REF] Gilbarg | Elliptic partial differential equations of second order, volume 224 of Grundlehren der Mathematischen Wissenschaften[END_REF] in the sense of Definition 3.1 with f 0 ≡ 0. Following Remark 3.7, γf = 0 on Σ -. Next, we choose appropriate test and renormalized functions in formulation [START_REF] Lions | Quelques méthodes de résolution des problèmes aux limites non linéaires[END_REF]. Thanks to Remark 3.4, we can choose β(z) = |z|. Let us also consider a nonnegative test function ϕ ∈ D(R 3 ), only depending on the velocity unknown ξ, such that ϕ ≡ 1 if |ξ| ≤ 1 and ϕ ≡ 0 if |ξ| ≥ 2. Then we set ϕ R (ξ) = ϕ(ξ/R) for any R ≥ 1, and use the weak formulation [START_REF] Lions | Quelques méthodes de résolution des problèmes aux limites non linéaires[END_REF] to write, for any t 1 ∈ [0, T ],

Ωt 1 ×R 3 |f (t 1 )|ϕ R + t 1 0 Σ + t |γf |ϕ R (ξ -w) • n t = t 1 0 Ωt×R 3 |f | G • ∇ ξ ϕ R .
The boundary term is clearly nonnegative. We can then successively write

Ωt 1 ×R 3 |f (t 1 )|ϕ R ≤ t 1 0 Ωt×R 3 |f | G • ∇ ξ ϕ R ≤ C R t 1 0 Ωt×R 3 |f | |G| 1 R≤|ξ|≤2R ≤ C t 1 0 Ωt×R 3 |f | |G| 1 + |ξ| 1 R≤|ξ| .
The last integral vanishes when R goes to +∞, thanks to the assumptions on both f and G. That eventually implies that f (t 1 , •) = 0 for any t 1 , and consequently f ≡ 0.

In the next section, we apply the previous results to Eqn. [START_REF] Chae | Global existence of weak and classical solutions for the Navier-Stokes-Vlasov-Fokker-Planck equations[END_REF] with G(t, x, ξ) = χ(u(t, x)ξ).

Existence of solutions to the approximated problem

We want to prove that the approximated problem ( 11)-( 13), supplemented with homogeneous Dirichlet boundary conditions for u on (0, T ) × ∂B and (5) for f , and initial data satisfying Assumptions 2.3-2.4, has at least one weak solution.

We must check that the mapping Θ σ is well defined for any σ ∈ [0, 1]. Choose u ∈ E. We first prove that there exists a unique f u solution of [START_REF] Chae | Global existence of weak and classical solutions for the Navier-Stokes-Vlasov-Fokker-Planck equations[END_REF], with absorption boundary conditions and initial datum satisfying Assumption 2.4, such that f u ∈ L ∞ (0, T ; L 1 (R 3 × R 3 ))∩L ∞ (0, T ; L ∞ (R 3 ×R 3 )). To do so, we just have to verify that G(t, x, ξ) = χ(u(t, x)ξ) has the required regularity properties. Since u ∈ E and thanks to the properties of χ, it is clear that G ∈ L 1 (0, T ; W 1,1 loc (B × R 3 )). Because of the asymptotical properties of χ, χ ′ is bounded on R, which obviously implies that div ξ G ∈ L ∞ (0, T ; L ∞ (B × R 3 )). Eventually, since E ֒→ L 2 (0, T ; L 6 (B)), we can write

|χ(u -ξ)| 1 + |ξ| ≤ |u| 1 + |ξ| + 1 ≤ |u| 6 (1 + |ξ|) 6 + C,
for some constant C. The previous estimate clearly allows to obtain [START_REF] O'rourke | Collective drop effects on vaporizing liquid sprays[END_REF]. Having existence and uniqueness of f u allows to focus on the approximated Navier-Stokes equations ( 15)- [START_REF] Dufour | Modélisation multi-fluide eulérienne pour les écoulements diphasiques à inclusions dispersées[END_REF].

In (15), we need to prove that the right-hand side term σS χ (u) -λσ1 B\ Ω (uw) lies in L 2 ((0, T ) × B)). The penalty term clearly belongs to L 2 ((0, T ) × B)), because of the regularity of u and w, and satisfies

(35) λσ1 B\ Ω (u -w) L 2 ((0,T )×B) ≤ |λ|( u L 2 (0,T ;L 2 (B)) + w L 2 (0,T ;L 2 (B)) ).
The regularity of S χ (u) is also clear. Indeed, f u has a compact support in ξ, see Remark 3.3, and lies in L ∞ (0, T ; L

1 (R 3 × R 3 )) ∩ L ∞ (0, T ; L ∞ (R 3 × R 3 )), with (36) f u L ∞ (0,T ;L 1 (R 3 ×R 3 ))∩L ∞ (0,T ;L ∞ (R 3 ×R 3 )) ≤ exp(3 χ ′ L ∞ (R) T ) f in L 1 (Ω×R 3 )∩L ∞ (Ω×R 3
) . Using [START_REF] Goudon | Hydrodynamic limit for the Vlasov-Navier-Stokes equations[END_REF] with G = χ(uξ), and since χ is a bounded function, that allows to state that S χ (u) ∈ L 2 ((0, T ) × B)), and

(37) S χ (u) L 2 ((0,T )×B) ≤ C T χ L ∞ (R) exp(3 χ ′ L ∞ (R) T ) f in L 1 (Ω×R 3 )∩L ∞ (Ω×R 3
) . Note that the proof of the previous estimate strongly uses the compact support assumption for f in . Consequently, since the convection velocity in ( 15) is regularized, we can apply standard existence and regularity results which hold for the Stokes system [START_REF] Lions | Quelques méthodes de résolution des problèmes aux limites non linéaires[END_REF][START_REF] Temam | Theory and numerical analysis[END_REF], namely, there exists a unique solution u σ ∈ H 1 (0, T ; L 2 (B)) ∩ L 2 (0, T ; H 2 (B)) ∩ E to ( 15)-( 18), which continuously depends on the data of the problem. That allows to define Θ σ for any σ ∈ [0, 1], and is useful to check the required properties to apply Schaefer's fixed point theorem to the mapping Θ.

Let us first prove that Θ sends bounded subsets on relatively compact subsets. Indeed, using estimates ( 35)- [START_REF] Wang | Global weak solution to the inhomogeneous Navier-Stokes-Vlasov equations[END_REF] on the right-hand side of the regularized Navier-Stokes equation (15), we have [START_REF] Williams | Combustion theory[END_REF]

u σ L 2 (0,T ;H 2 (B))∩H 1 (0,T ;L 2 (B)) ≤ C u L 2 (0,T ;L 2 (B)) , u in H 1 (B) , f in L 1 (Ω×R 3 )∩L ∞ (Ω×R 3 ) ,
where C also depends on χ, ϕ, λ, w and T . Consequently, the previous estimate and the compact embedding of H 1 (0, T ; L 2 (B)) ∩ L 2 (0, T ; H 2 (B)) ∩ E in E allow to obtain the required relative compactness property.

Let us now focus on the continuity of Θ. Consider a sequence (u n , σ n ) n∈N converging towards (u, σ) in E × [0, 1] when n goes to +∞. Set u n = Θ σn (u n ) for any n. We want to prove that ( u n ) strongly converges in E towards u = Θ σ (u). Thanks to [START_REF] Temam | Theory and numerical analysis[END_REF]

, (f un ) is uniformly bounded in L ∞ (0, T ; L 1 (R 3 × R 3 )) ∩ L ∞ (0, T ; L ∞ (R 3 × R 3 )).
Up to a subsequence, (f un ) weakly- * converges in L ∞ (0, T ; L ∞ (R 3 × R 3 )), and (f un ) weakly converges in L p loc ( Ω × R 3 ), for any p ∈ (1, +∞). That allows us to pass to the limit in the weak formulation [START_REF] Goudon | The Navier-Stokes-Vlasov-Fokker-Planck system near equilibrium[END_REF] of Definition 3.1, which, as stated in Remark 3.1, is equivalent to the one appearing in Definition 2.1. By uniqueness of the weak solution to the Vlasov equation, the weak limit is f u , and the whole sequence (f un ) converges. Thanks to estimate [START_REF] Williams | Combustion theory[END_REF], ( u n ) is uniformly bounded in H 1 (0, T ; L 2 (B)) ∩ L 2 (0, T ; H 2 (B)) ∩ E, and thus compact in E. Hence, up to a subsequence, ( u n ) strongly converges in E. We can then pass to the limit in the regularized Navier-Stokes system, also using the convergences of (u n ), (σ n ) and (f un ). We obtain the continuity of Θ, again thanks to a uniqueness argument.

Eventually, let us prove that the set of all the fixed points of the family Θ σ is bounded in E. Consider σ ∈ [0, 1], and u ∈ E such that u = Θ σ (u). To obtain the required boundedness property, we derive a priori bounds for both f u and u σ = u. Since f u is currently compactly supported in ξ, we can rigorously apply [START_REF] Lions | Quelques méthodes de résolution des problèmes aux limites non linéaires[END_REF] with G = χ(uξ), β(z) = z and ψ = |ξ| 2 . We then obtain, also recalling the homogeneous boundary condition on each Σ - s ,

2 t 0 Ωs×R 3 f u ξ • χ(u -ξ) = Ωt×R 3 f u |ξ| 2 - Ωt×R 3 f in |ξ| 2 + t 0 Σ + s γf u |ξ| 2 (ξ -w) • n s ,
where the last term is clearly nonnegative. Consequently, we can write

(39) 1 2 M 2 f u (t) ≤ 1 2 M 2 f in + t 0 Ωs×R 3 f u ξ • χ(u -ξ).
Besides, multiplying (15) by u σ = u, we get

(40) 1 2 u(t) 2 L 2 (B) + t 0 ∇ x u 2 L 2 (B) = σ 2 u in 2 L 2 (B) + σ t 0 B×R 3 χ(ξ -u) • uf u + σλ t 0 B 1 B\ Ω (w -u) • u.
Then, multiplying (39) by σ, and adding (40), we obtain, since σ ∈ [0, 1], and thanks to the properties of χ and Young's inequality,

(41) σ 2 M 2 f u (t) + 1 2 u(t) 2 L 2 (B) + t 0 ∇ x u 2 L 2 (B) + σ t 0 B×R 3 χ(u -ξ) • (u -ξ)f u ≤ σ 2 M 2 f in + 1 2 u in 2 L 2 (B) + σλ t 0 B 1 B\ Ω (w -u) • u ≤ 1 2 M 2 f in + 1 2 u in 2 L 2 (B) + λ t 0 B |w • u| ≤ 1 2 M 2 f in + 1 2 u in 2 L 2 (B) + λ 2 w 2 L 2 (0,T ;L 2 (B)) + t 0 u 2 L 2 (B) .
Gronwall's lemma then ensures that u is bounded in E by a quantity which only depends on the data of the problem and the regularization parameters.

Consequently, thanks to Schaefer's theorem, we can state that Θ 1 has a at least one fixed point, which leads to the required existence result.

Cauchy-Schwarz and Young's inequalities imply, together with the nonnegativity of f n , that

-n t 0 R 3 1 B\ Ω (u n -w) • w ≤ 1 4 u n (t) 2 L 2 (B) + 1 2 u n in 2 L 2 (B) + C w t 0 u n 2 L 2 (B) + 1 2 t 0 ∇ x u n 2 L 2 (B) + C w + C w t 0 M 0 f n + 1 2 t 0 B×R 3 f n |χ n (u n -ξ)| 2 .
If we add (45) to the previous estimate, we obtain, thanks to the inequality χ n (z) 2 ≤ χ n (z)z for any z ∈ R, and the nonnegativity of f n ,

1 2 M 2 f n (t) + 1 4 u n (t) 2 L 2 (B) + 1 2 t 0 ∇ x u n 2 L 2 (B) + 1 2 t 0 B×R 3 f n |χ n (u n -ξ)| 2 + n t 0 B 1 B\ Ω (u n -w) 2 ≤ 1 2 M 2 f n in + u n in 2 L 2 (B) + C w t 0 u n 2 L 2 (B) + C w t 0 M 0 f n + C w .
Since, for a given n, f n is compactly supported, we straightforwardly have, from the weak formulation of the Vlasov equation in Definition 2.1,

t 0 M 0 f n ≤ T M 0 f n in .
Taking into account the boundedness and convergence properties of the sequences of initial data, we deduce

(46) 1 2 M 2 f n (t) + 1 4 u n (t) 2 L 2 (B) + 1 2 t 0 ∇ x u n 2 L 2 (B) + n t 0 B 1 B\ Ω (u n -w) 2 ≤ 1 2 M 2 f in + u in 2 L 2 (B) + C w t 0 u n 2 L 2 (B) + T C w M 0 f in + C w .
Hence, thanks to Gronwall's lemma, we get the following bound for (u n ) and (f n ): The previous bounds also imply that there exists (u, f ) such that u ∈ L ∞ (0, T ; L 2 (B)) ∩ L 2 (0, T ; H 1 (B)), f ∈ L ∞ (0, T ; L ∞ (R 3 × R 3 ) ∩ L 1 (R 3 × R 3 )), and that the following weak convergences hold, up to a subsequence: Besides, to perform the asymptotics n → +∞ in the weak formulations of (42)-(44), in particular to handle the convection and coupling terms, we need some strong convergence on the fluid velocity. A standard way to obtain compactness is the Aubin-Lions lemma, which requires uniform bounds on ∂ t u n . Nevertheless, because of the time dependence of the domain and our approximation penalty strategy, this lemma cannot directly be applied to (u n ). Various compactness proofs of the fluid velocity can be found, in the very same context, for instance, in [START_REF] Fujita | On existence of weak solutions of the Navier-Stokes equations in regions with moving boundaries[END_REF][START_REF] Salvi | On the existence of weak solutions of a nonlinear mixed problem for the Navier-Stokes equations in a time dependent domain[END_REF]. Here we choose to apply [32, Theorem 4 and Remark 4.5].

u n ⇀ u, ∇ x u n ⇀ ∇ x u
More precisely, we know, thanks to the energy estimate, that (u n ) is bounded in L ∞ (0, T ; L 2 (B)) ∩ L 2 (0, T ; H 1 0 (B)). In order to apply [START_REF] Moussa | Some variants of the classical Aubin-Lions Lemma[END_REF]Theorem 4], we have to prove two properties.

The first one is the compactness of the normal trace of (u n ) in L 2 (0, T ; H -1/2 (∂Ω t )). This control is obtained since (u n ) strongly converges towards w in L 2 ( B \ Ω), and the normal trace operator on Γ is continuous from {ψ ∈ L 2 ( B \ Ω) | div x ψ ∈ L 2 ( B \ Ω)} to L 2 (0, T ; H -1/2 (B \ Ω t )).

The second property is the following: for any test function ψ ∈ D( Ω) such that div x ψ = 0 on Ω, (48)

T 0 B ∂ t u n ψ ≤ C |α|≤2 ∂ α x ψ L 2 ( Ω) ,
where C is a constant only depending on the data of the problem. This bound is obtained using the weak form of the Navier-Stokes equations, and relies on the already derived energy estimate, and on a uniform bound of (S χn (u n )). Since ψ is compactly-supported in Ω, the integral involving the penalty term cancels. As for (S χn (u n )), we proceed as follows. Recalling that |χ n (z)| ≤ |z|, we have, almost everywhere,

|S χn (u n )| ≤ m 1 f n + |u n |m 0 f n .
We know that (M 2 f n ) is bounded in L ∞ (0, T ) thanks to the energy estimate (47). Hence, Lemma 2.1 ensures that (m 1 f n ) and (m 0 f n ) are respectively bounded in L ∞ (0, T ; L 5/4 (B)) and L ∞ (0, T ; L 5/3 (B)). Since (u n ) is bounded in L 2 (0, T ; H 1 (B)), which is continuously embedded in L 2 (0, T ; L 6 (B)), we eventually get the boundedness of (S χn (u n )) in the space L 2 (0, T ; L 1 (B)). Estimate (48) follows from where the last inequality is a consequence of the embedding of H 2 in L ∞ in three dimensions.

This allows to let n go to +∞ in the weak formulation of the regularized coupled system to obtain a weak solution to (1)-( 5) as given in Definition 2.1. In particular, the strong convergence of (u n ) enables to pass to the limit in the nonlinear terms (convection and coupling terms). Moreover, we prove that M 2 f lies in L ∞ (0, T ) by performing the asymptotics in the energy estimate (47), following the very same method as in [7, p. 1268].

( 47 )

 47 u n L ∞ (0,T ;L 2 (B))∩L 2 (0,T ;H 1 (B)) + M 2 f n L ∞ (0,T ) ≤ C w (T, u in L 2 (B) , M 0 f in , M 2 f in ),where the constant depends on w, T , u in L 2 (B) , M 0 f in and M 2 f in .Estimates (46)-(47) ensures that (u n ) → w strongly in L 2 ( B \ Ω).

  weakly in L 2 (0, T ; L 2 (B)),u n ⇀ u weakly * in L ∞ (0, T ; L 2 (B)), and f n ⇀ f weakly * in L ∞ (0, T ; L ∞ (R 3 × R 3 )).

S

  χn (u n ) • ψ ≤ S χn (u n ) L 2 (0,T ;L 1 (B)) ψ L 2 (0,T ;L ∞ (B)) ≤ C |α|≤2 ∂ α x ψ L 2 ( Ω) ,

This work was partially funded by the French ANR-13-BS01-0004 project Kibord headed by L. Desvillettes.

Back to the whole system

The results obtained in Section 4 allows us to state the existence, for any n ∈ N * , of (u n , f n ) weakly solving the following system (42) (43) (44)

where (ϕ n ) is a mollifying sequence. The truncation functions χ n satisfy the same assumptions as χ in the previous sections and are furthermore chosen such that (χ n ) converges to the identity mapping on R when n goes to +∞. However, note that χ n is not uniformly bounded with respect to n. Moreover, the system is supplemented with regularized initial conditions (see Assumptions 2.3-2.4)

The sequences (u n in ) and (f n in ) are chosen such that (u n in|Ω ) strongly converges to u in in L 2 (Ω), that (f n in ) strongly converges to f in in L p (Ω × R 3 ) for all 1 ≤ p < +∞, weakly * in L ∞ (Ω × R 3 ), and that (M 2 f n in ) is uniformly bounded with respect to n and strongly converges to M 2 f in in L ∞ (Ω).

Remark 5.1. Note that, for a given n, u n is in fact a strong solution to the Navier-Stokes equations, satisfied in L 2 (0, T ; L 2 (B)), whereas f n is only a weak solution of the Vlasov equation in the sense of Definition 2.1.

We need to obtain, for the previous approximated coupled problem, an energy estimate similar to (41), but which does not depend on the regularization parameter n. It relies on a proper treatment of the penalty term. Instead of having

as in the first inequality of (41), we bring out

which is nonnegative. To recover the previous term, we rewrite (41) as

Then we multiply (42) by w and integrate on B. After integration by parts on both variables t and x, remembering that w = 0 on ∂B and div x (u n ⋆ ϕ n ) = 0, we get


Adrien Tessier

Julie Bertrand

Marylore Chenel

Emmanuelle Comets

Combined analysis of phase I and phase II data to enhance the power of pharmacogenetic tests

Keywords: Pharmacogenetics, Pharmacokinetics, Nonlinear mixed effects models, Design, Clinical drug development

We show through a simulation study how the joint analysis of data from phase I and phase II studies enhances the power of pharmacogenetic tests in pharmacokinetic (PK) studies. PK profiles were simulated under different designs along with 176 genetic markers. The null scenarios assumed no genetic effect, while under the alternative scenarios, drug clearance was associated to 6 genetic markers randomly sampled in each simulated dataset. We compared penalised regression Lasso and stepwise procedures to detect the associations between empirical Bayes estimates of clearance, estimated by nonlinear mixed effects models, and genetic variants. Combining data from phase I and phase II studies, even sparse, increases the power to identify the associations between genetics and PK due to the larger sample size. Design optimisation brings a further improvement, and we highlight a direct relationship between η-shrinkage and loss of genetic signal.

INTRODUCTION

Studying the sources of the variability observed in drug response facilitates individualisation of prescription. One of the sources of variability in drugs' pharmacokinetics (PK) [START_REF] Aarons | Population pharmacokinetics: theory and practice[END_REF] is the variation in activity of enzymes and transporters involved in the drug absorption, distribution, metabolism or elimination. Pharmacogenetics 2 studies the genetic component of interindividual variability (IIV) observed in PK to identify populations at risk of treatment inefficacy or adverse effects [START_REF] Guo | Pharmacogenomics and drug development[END_REF] . Single nucleotide polymorphisms (SNPs) are the genetic variants most frequently studied in pharmacogenetics and screened more and more often in clinical studies.

Genetic data offer some unique challenges, in particular because they may lead to very unbalanced number of subjects, which impacts the power of tests in pharmacogenetic analyses [START_REF] Bertrand | Pharmacogenetics and population pharmacokinetics: impact of the design on three tests using the SAEM algorithm[END_REF][START_REF] Bertrand | Some alternatives to asymptotic tests for the analysis of pharmacogenetic data using nonlinear mixed effects models[END_REF] . In a previous simulation work, we showed that typical phase I studies have low power to detect genetic effects because of the limited sample size [START_REF] Tessier | Comparison of nonlinear mixed effects models and noncompartmental approaches in detecting pharmacogenetic covariates[END_REF] . On the other hand, phase I studies generally provide good quality PK information allowing to characterise the PK profile of the drug. We showed that from the different approaches used at this stage to estimate PK parameters, nonlinear mixed effects models (NLMEM) [START_REF] Sheiner | Modelling of individual pharmacokinetics for computer-aided drug dosage[END_REF] could be considerably more powerful than noncompartmental analyses (NCA) [START_REF] Rowland | Clinical pharmacokinetics and pharmacodynamics: concepts and applications[END_REF] for complex PK models [START_REF] Tessier | Comparison of nonlinear mixed effects models and noncompartmental approaches in detecting pharmacogenetic covariates[END_REF] . Our simulations also showed that increasing the sample size, as in phase II studies, would improve the power to detect genetic variants. However sparse designs typically used in phase II may result in biased estimations for empirical Bayes estimates (EBE) [START_REF] Savic | Importance of shrinkage in empirical bayes estimates for diagnostics: problems and solutions[END_REF] used in GAM covariate analysis procedure (generalised additive models) [START_REF] Mandema | Building population pharmacokinetic-pharmacodynamic models. I. Models for covariate effects[END_REF] .

To increase the detection of genetic covariates, a way could be to combine for the analysis data from a study collected with a rich design, as expected in phase I, with sparser, but still informative, data from a phase II study.

In the present work we propose practical designs involving phase I and phase II data, and we quantify through simulations their ability to detect genetic associations with PK. A motivating example was provided by IRIS, a pharmaceutical industry, to generate realistic genetic and PK data. We compared two association methods, a penalised regression method and a stepwise procedure [START_REF] Tessier | Comparison of nonlinear mixed effects models and noncompartmental approaches in detecting pharmacogenetic covariates[END_REF] .

MATERIALS AND METHODS

Simulation study

Figure 1 presents the framework of the simulation study, which was designed based on PK data from drug S (IRIS) collected in 78 subjects from 3 phase I clinical studies [START_REF] Tessier | High-throughput genetic screening and pharmacokinetic population modeling in drug development[END_REF] . All subjects were genotyped at baseline using a DNA microarray developed by IRIS of 176 SNPs known for being involved in the PK of drugs. These 176 polymorphisms were matched to a reference Hapmap panel (Hapmap 3 release 2) for a Caucasian population [START_REF]International HapMap Consortium The International HapMap Project[END_REF] and we used the Hapgen2 software [START_REF] Su | HAPGEN2: simulation of multiple disease SNPs[END_REF] to simulate genetic variants retaining their frequencies and the correlations between polymorphisms found in the human genome (see details in Supplementary Material, Supplementary Figures S1-S3).

PK profiles were simulated with a two-compartment model with dose-dependent double absorption (Supplementary Figure S4), with the parameters in

𝛽 𝑘 = √ 𝑅 𝐺𝐶 𝑘 × 𝜔 2 𝐶𝐿 2𝑝 𝑘 (1 -𝑝 𝑘 ) -𝑅 𝐺𝐶 𝑘 × 2𝑝 𝑘 (1 -𝑝 𝑘 ) (2)
where 𝜔 2 𝐶𝐿 is the variance of interindividual random effects on CL due to non-genetic sources. 𝑅 𝐺𝐶 𝑘 was respectively equal to 1, 2, 3, 5, 7 and 12% for the 6 causal variants [START_REF] Bertrand | Multiple single nucleotide polymorphism analysis using penalized regression in nonlinear mixed-effect pharmacokinetic models[END_REF] to mimick a multifactorial genetic effect. Then under H 0, 𝜔 2 𝐶𝐿 = 0.06 (as in Table I), while under H 1 30% of the variance is explained by the genetics so that 𝜔 2 𝐶𝐿 = 0.04 (example on the magnitude of simulated effect sizes is available in Supplementary Table S1).

The simulated datasets were then fitted with the base model without genetic covariates.

Individual clearance estimates (𝐸𝐵𝐸 𝐶𝐿 𝑖 ) were estimated and all associations with the 176 simulated polymorphisms were tested assuming a linear relation without re-estimating model parameters, as in a GAM analysis [START_REF] Mandema | Building population pharmacokinetic-pharmacodynamic models. I. Models for covariate effects[END_REF] .

We compared two association methods to detect gene effects. Lasso 15 is a multivariate penalised regression which simultaneously estimates effect size coefficients and selects variants by setting a large number of coefficients to 0. The penalty is set by a tuning parameter (𝜉) which depends on 𝛼, the type I error per test, and the number of subjects [START_REF] Hoggart | Simultaneous analysis of all SNPs in genome-wide and re-sequencing association studies[END_REF][START_REF] Bertrand | Integrating dynamic mixed-effect modelling and penalized regression to explore genetic association with pharmacokinetics[END_REF] (Figure 1). Alternatively in practice the penalty can be determined through permutation or cross-validation methods, which are more time consuming. Stepwise procedure includes relationships one by one depending on the significance of a Wald test compared to a threshold 𝛼. The correlation between two significant SNPs, due to linkage disequilibrium, is computed through the Pearson's correlation coefficient 𝑟 and if two significant SNPs are strongly correlated (|𝑟| > 0.89), only the most significant is kept. Finally the most significant variant among selected SNPs is kept in the final model and steps are repeated until no more association is significant [START_REF] Tessier | Comparison of nonlinear mixed effects models and noncompartmental approaches in detecting pharmacogenetic covariates[END_REF] (Figure 1).

In both approaches, we control the Family Wise Error Rate (FWER), representing the percentage of datasets where at least one variant is selected under H 0 , by correcting the nominal 𝛼 by the number of tests performed (Sidak correction) corresponding to the number of polymorphic SNPs considered 𝑁 𝑡 (Figure 1). The FWER was set to 20% (with a prediction interval for 200 datasets equal to [14.5-25.5]) for an exploratory analysis. The prediction interval determined when to adjust 𝛼 to control the FWER under H 0 .

Simulated designs and analysis scenarios

We simulated a phase I study corresponding to the motivating example, including 78 subjects (N1) receiving 8 different single doses (5, 10, 20, 50, 100, 200, 400 or 800 units, for respectively 6, 6, 24, 12, 12, 6, 6 and 6 subjects per dose) and sampled at 16 times. Three designs of phase II study were simulated. They included 306 subjects (N2), receiving 3 doses (20, 50 or 100 units, 102 subjects per dose), sampled at steady state. Two phase II studies included three samples per subject, optimised using the PFIM software [START_REF] Bazzoli | Design evaluation and optimisation in multiple response nonlinear mixed effect models: PFIM 3.0[END_REF] to ensure a reasonable precision of CL estimates. The last sampling time was limited to 24h in one, while a late sample was allowed after the last dose administration in the other. The third study included only one trough concentration (24h). We considered 4 analysis scenarios (Figure 1), three combining the phase I and one of the phase II study (respectively SPI/II 3s.24h , SPI/II 3s.96h

and SPI/II 1s.24h ), and one, for comparison, with only the phase I subjects (SPI).

We also investigated the impact of a higher variability on phenotype on the results. For this, we simulated the same four scenarios increasing the IIV on CL to 60% (instead of 25% in previous settings).

Evaluation

For each analysis scenario, 200 datasets were simulated under H 0 and H 1 .

The ability of the designs to estimate the population and individual parameters under H 0 was first evaluated through estimation bias and η-shrinkage (see details in Supplementary Material, Supplementary Figures S5-S6).

Under H 1 we evaluated the performance of each scenario in terms of true and false positives counts (TP and FP) and rates (TPR, the proportion of TP detected among the causal variants;

and FPR, the proportion of FP detected among all potential false associations) for parameter CL, as well as the probability to detect genetic variants. Assuming that SNPs located on genes coded for metabolism enzymes and transporters affect mostly the drug distribution and elimination, we also applied association tests on Q, the intercompartmental clearance and V2, the peripheral volume, separately. Any variants associated to Q and V2 were counted as false positives. The central volume V1 was not considered because it had no random effects.

We also evaluated the loss of genetic signal between simulated and estimated individual clearances, comparing slopes 𝑏 of univariate linear regressions on log(𝐶𝐿 𝑠𝑖𝑚 𝑖 ) or log(𝐸𝐵𝐸 𝐶𝐿 𝑖 ) for each causal variant. A relative deviation of the genetic signal 𝑅𝐷 𝑠𝑖𝑔𝑛𝑎𝑙 was computed as follows:

𝑅𝐷 𝑠𝑖𝑔𝑛𝑎𝑙 (%) = 𝑏 𝐸𝐵𝐸 𝐶𝐿 𝑖 -𝑏 𝐶𝐿 𝑠𝑖𝑚 𝑖 𝑏 𝐶𝐿 𝑠𝑖𝑚 𝑖 × 100 (3)
𝑅𝐷 𝑠𝑖𝑔𝑛𝑎𝑙 quantifies the departure of the estimated genetic signal (𝑏 𝐸𝐵𝐸 𝐶𝐿 𝑖 ) from the one simulated (𝑏 𝐶𝐿 𝑠𝑖𝑚 𝑖 , see details in Supplementary Material).

RESULTS

Control of FWER under H 0

Lasso and stepwise procedure both tended to be too conservative, as the FWER was lower than expected in some scenarios (Table 2). After an empirical correction by increasing the type I error per SNP α, FWER was properly controlled around 20%. This correction was applied in the corresponding simulation under H1. Previous simulations suggested that this decrease in FWER is influenced by correlations between polymorphisms [START_REF] Tessier | Comparison of nonlinear mixed effects models and noncompartmental approaches in detecting pharmacogenetic covariates[END_REF] .

Detection of genetic effects

Under H 1 the TPR (Figure 2, top left) was higher in scenarios including phase II data (from 22 to 32%) compared to scenario with only phase I data (SPI, 4%) and was the highest in scenario SPI/II 3s.96h . The FPR was lowest (0.2%) in scenario SPI where a limited number of SNPs was selected, and only slightly higher in scenarios including phase II data, ranging from 0.6 to 0.8% for both methods. Very few TP were effectively detected in scenario SPI (around 44 for both methods) where the number of subjects was limited (N1 = 78) (Supplementary

Tables S2-S3

). By adding more subjects (N2 = 306) to the analysis, the number of TP increased sharply. Scenario SPI/II 3s.96h allowed detecting the largest number of TP (380 or more), while in SPI/II 3s.24h around 326 TP were detected. In SPI/II 1s.24h the number of TP was lower (around 270 TP), but remained much higher than scenario SPI with only phase I data.

In the same way, the number of FP increased when including phase II data to the analysis, but to a much lesser extent.

With only phase I data, the probability to detect at least one genetic variant on CL was low (Figure 2, bottom left), around 20% (SPI). This probability decreased quickly when trying to detect more polymorphisms and reached 0 for 3 variants or more. Adding phase II data to the analysis increased the probability to detect at least one variant about 85% in scenario SPI/II 1s.24h , and up to 95% in scenario SPI/II 3s.96h . Scenarios including phase II data showed good detection of 1 to 3 SNPs and SPI/II 3s.96h had always the higher detection. This shows that the major determinant of power is the number of subjects, and that optimising the design for more informativeness can bring a smaller further improvement. The low probability to detect 4 SNPs or more (≤ 4%) in scenarios combining phase I and phase II data can be explained by those variants having a very weak impact; polymorphisms only explaining 1, 2 or 3% of the variability of CL.

In Supplementary Table S4, the TPR was computed separately for each causal SNP. The variants associated to the lowest R GC had low TPR, close to the FPR. Thus the signal associated to these variants was close to the noise created by the non-causal variants

Shrinkage

Two η-shrinkage estimates were computed using metric proposed by Bertrand et al. [START_REF] Bertrand | Pharmacogenetics and population pharmacokinetics: impact of the design on three tests using the SAEM algorithm[END_REF] based on estimated variances, with respect to the estimate of 𝜔 ̂2 in the dataset; one over the 𝜂î from phase I subjects and one over the 𝜂̂𝑖 from phase II subjects (Figure 3). The η-shrinkage for phase I subjects was low (median = 23%) thanks to the large number of observations per subject. A large range of η-shrinkage estimates for phase II data was observed across analysis scenarios, but was below 50% in scenario SPI/II 3s.96h .

Loss of signal

𝑅𝐷 𝑠𝑖𝑔𝑛𝑎𝑙 was always negative for the 6 SNPs, indicating that part of the signal was lost during the estimation step (Figure 4, top). This loss was smaller in the scenario with phase I data alone (SPI) than in scenarios combining phase I and phase II data. In each scenario the signal loss was of the same magnitude for the 6 SNPs, regardless of the value of associated R GC . For phase I data (Figure 4, bottom), the loss was of a constant magnitude across scenarios (median = -30%). For phase II data, in the most informative scenario (SPI/II 3s.96h )

the loss was of a similar magnitude (median = -41%) than the loss in phase I data, where subjects were extensively sampled. The loss was higher in scenario SPI/II 3s.24h (median = -56%), and even more when only one time was sampled (SPI/II 1s.24h , median = -70%).

The signal loss and η-shrinkage values changed accordingly across phase II scenarios, while the probability of detection changed in the opposite direction.

Influence of the phenotype variance

Increasing IIV for the CL parameter to 60% led to a sharp increase in the number of TP (Supplementary Tables S5-S6), resulting in higher TPR and higher probabilities to detect the causal variants (Figure 2, right), compared to when individual CL were simulated with a moderate IIV. This higher number of TP is explained first and foremost by the increase in simulated effect sizes which depended on the variance of interindividual random effects on CL due to non-genetic sources (equation 2). A second consequence of the larger IIV was that the estimated η-shrinkages became much smaller. Lower η-shrinkages resulted in lower signal losses in all scenarios for phase I and phase II data (Supplementary Figure S7-S8), which again favoured a higher probability to detect the genetic effects.

DISCUSSION

In this work, we show and evaluate practical designs to combine data from studies occurring in phase I and II of a drug development. We assess through a simulation study, inspired by a real example, the probability to detect genetic variants and the influence of the phase II study design. We considered phenotypes estimated by NLMEM, which can handle the analysis of heterogeneous data involving sparsely sampled subjects.

Genetic variants are unbalanced and so the amount of information they provide is directly related to the variant allele frequency and the study sample size. On the other hand, PK information depends also on the number and times of sampling which drives the precision of the PK model parameter estimates. A limited number of samples, as in phase II studies, may lead to missing a true association when EBE are used as phenotypes [START_REF] Savic | Importance of shrinkage in empirical bayes estimates for diagnostics: problems and solutions[END_REF] . Savic and Karlsson suggested a more extensive use of the likelihood ratio test (LRT) for covariate selection when η-shrinkage is large, but Combes et al. showed that the power to detect a covariate effect is the same with a LRT or a simple correlation test on EBE [START_REF] Combes | Powers of the Likelihood Ratio Test and the Correlation Test Using Empirical Bayes Estimates for Various Shrinkages in Population Pharmacokinetics[END_REF] .

The effect of sample size can be distinctly observed in our simulations. In the context of phase I studies, where the number of subjects is limited, the probability to detect the genetic effects was low, in line with our previous results [START_REF] Tessier | Comparison of nonlinear mixed effects models and noncompartmental approaches in detecting pharmacogenetic covariates[END_REF] . The combined analysis of phase I and phase II data allowed a marked improvement in this detection probability, irrespective of the phase II study design. By modifying the design of the phase II data, we highlighted a direct link between η-shrinkage, loss of genetic signal and probability to detect genetic variants. Our results showed that poor PK information due to the phase II study design results in higher η-shrinkage, which increases the loss of genetic signal at the estimation step and translates to a lower probability to detect genetic variant. The dilution of the individual information by adding subjects with sparse designs to subjects with rich designs increases as expected the loss of genetic signal. But this is accompanied with a sharp increase in detection power thanks to a larger sample size. η-shrinkage may also modify the EBE-EBE relationship, falsely inducing or masking correlations between model parameters [START_REF] Savic | Importance of shrinkage in empirical bayes estimates for diagnostics: problems and solutions[END_REF] . This could result in an increased number of false positives associated to other parameters than CL, although in our simulations the number of FP on CL, V2 and Q remained of a similar magnitude across scenarios (Supplementary Table S3), showing no systematic effect.

We assume homogeneity of the PK between subjects simulated for the phase I and the phase II study. In practice healthy volunteers are often included in phase I while phase II studies focus on patients. A difference in typical values, for example of CL, between the two populations should not impact the detection power by combination of data, as the association tests use the phenotype variance, provided that the genetic effect is the same and that the model accounts for the systematic difference between clearances. It is more difficult to predict what would happen if the variability of clearance is different in the two populations, as the magnitude of the shrinkage in each subpopulation could affect the signal detection. When the assumption that the two populations are similar breaks down, we would suggest instead to combine rich and sparse data within the phase II study.

Pharmacogenetic studies including a large number of subjects combining sparse and rich designs have already been published [START_REF] Chou | Population pharmacokinetic-pharmacogenetic study of nevirapine in HIVinfected Cambodian patients[END_REF][START_REF] Bertrand | Multiple genetic variants predict steady-state nevirapine clearance in HIV-infected Cambodians[END_REF] , showing that the combination of different sampling designs is feasible within the same study to assure more homogeneity.

Situations where pharmacogenetic analyses in PK studies are recommended are described by health authorities [START_REF]EMA Guideline on the use of pharmacogenetic methodologies in the pharmacokinetic evaluation of medicinal products[END_REF] . In our work, we simulated a blinded pharmacogenetic analysis, exploring a large number of genetic markers. In real applications, other considerations than the statistical significance of genetic variants such as their physiological and clinical relevance could be factored in the analysis and its interpretation. Lehr et al. [START_REF] Lehr | Integration of high-throughput genotyping data into pharmacometric analyses using nonlinear mixed effects modeling[END_REF] proposed in their stepwise procedure to select only significant polymorphisms having a physiologic relevance in the final model, and the same constraint could be integrated in penalised regression approaches. The probability to detect genetic variants could also be increased through the targeted inclusion of subjects for a few polymorphisms of interest, but this approach requires hypotheses on which polymorphisms to test, with a risk to miss important associations. We focused in this work on PK variability, which is a part of the variability in drug response. But the conclusions from the simulation study could be extended to pharmacodynamics. A previous survey indicated that most pharmacogenetic analyses in clinical PK studies used a phenotype estimated by NCA and furthermore included a limited number of subjects (lower than 50 subjects in two thirds) [START_REF] Tessier | Comparison of nonlinear mixed effects models and noncompartmental approaches in detecting pharmacogenetic covariates[END_REF] . Authorities in fact recommend to study pharmacogenetics in phase I [START_REF]EMA Guideline on the use of pharmacogenetic methodologies in the pharmacokinetic evaluation of medicinal products[END_REF] , where the number of subjects is limited. Our work

shows that such analyses do not have the power to detect polymorphisms efficiently but can generate hypotheses to assess in later studies. A recent simulation work [START_REF] Kloprogge | Statistical Power Calculations for Mixed Pharmacokinetic Study Designs Using a Population Approach[END_REF] studied the sample size required to detect a binary covariate. They conclude that around 60 subjects combining rich or sparse designs was sufficient to detect the covariate with at least 80% power. Again our simulations showed that genetic covariates require higher sample sizes because they are highly unbalanced.

In the first series of simulations a moderate IIV on CL (25%) was used, resulting in a low impact of the genetics on PK, since overall 30% of the moderate CL variability was explained by genetic variants. This setting represented a realistic case to challenge the detection of genetic variants through modelling. We also evaluated the same scenarios with a higher IIV for CL, set to 60%. The η-shrinkage was much lower, as a higher IIV downweighs the population priors in the combined criterion used to compute EBE. This decrease in CL ηshrinkages resulted in lower signal loss because of the direct relationship between the two.

Associated with larger simulated effect sizes, the number of TP and the probability to detect genetic variants increased in these scenarios. The effect of η-shrinkage on the probability to detect genetic effects was in these simulations higher than the one we observed with the main settings, because the decrease of η-shrinkage was associated with a sharp increase in the number of TP. This shows that our conclusions don't depend on the level of IIV.

This simulation study also confirms the results of our previous work concerning the relative performance of the different association methods [START_REF] Tessier | Comparison of nonlinear mixed effects models and noncompartmental approaches in detecting pharmacogenetic covariates[END_REF] . The penalised regression method Lasso and the stepwise procedure showed a similar probability to detect genetic variants in all scenarios. However, the Lasso is a slightly more complex method which requires computing the penalty in a first step before testing the associations. In this work we assessed methods to detect genetic effects on EBE, after an initial fit. An algorithm proposed by Lehr et al. [START_REF] Lehr | Integration of high-throughput genotyping data into pharmacometric analyses using nonlinear mixed effects modeling[END_REF] uses univariate regressions to select variants to test in the PK model through LRT. This approach is easy to implement but run-times depend on the number of iterations leading to the full covariates model. An alternative is to use an integrated approach where effect sizes are estimated and significant variants selected using a penalised regression in the same step [START_REF] Bertrand | Integrating dynamic mixed-effect modelling and penalized regression to explore genetic association with pharmacokinetics[END_REF] ; it showed similar performance than the stepwise procedure proposed by Lehr et al., but with longer computing times [START_REF] Bertrand | Integrating dynamic mixed-effect modelling and penalized regression to explore genetic association with pharmacokinetics[END_REF] . The results for the two other penalised regression methods tested in the previous work, ridge regression and HyperLasso, were similar (Supplementary Tables S7-S9, Figures S9-S10). None of the methods detected the 6 SNPs simultaneously, as 3 of the polymorphisms only explained 1 to 3% of the clearance variability, making them difficult to detect. Because association methods relate the polymorphisms to the phenotype variance, we fixed the variance explained by the causal variants (through the parameter R GC ) and computed the effect sizes as a function of their allelic frequencies. For a given R GC an infrequent polymorphism was therefore associated with higher effect sizes. This reflects that a clinically relevant polymorphism (with a high impact on PK), present in few subjects because of its low frequency, will explain a limited proportion of the phenotype variance. Detecting such polymorphisms is crucial to identify subpopulation at risk but required much larger sample sizes, as in genome-wide studies [START_REF] Takeuchi | A genome-wide association study confirms VKORC1, CYP2C9, and CYP4F2 as principal genetic determinants of warfarin dose[END_REF] . As an example, the rs3918290 polymorphism from gene DPYD has a frequency lower than 1%, but results in a deficient dihydropyrimidine dehydrogenase activity associated with a 40% decrease of the maximum conversion capacity of the chemotherapeutic drug 5fluorouracil [START_REF] Kuilenburg | Evaluation of 5-fluorouracil pharmacokinetics in cancer patients with a c.1905+1G>A mutation in DPYD by means of a Bayesian limited sampling strategy[END_REF] , resulting in severe toxicities.

The power to detect polymorphisms is also closely related to the type I error chosen for the analysis. In a context of exploratory analyses, we fixed the global type I error to 20%. But using the Sidak correction the significance thresholds were finally lower than 0.1% for each test, so that only strong effects of causal variants will be detected, and our simulations show that polymorphisms explaining a limited part of the phenotype variance are not detected.

Approaches based on FWER and corrections as Bonferroni or Sidak are easy to implement but are conservative and may reduce the power of analyses, but limit the number of polymorphisms to test in later confirmatory trials. In practice other corrections for type I error could be considered, as permutation methods which are more time consuming but less conservative.

Although this correction was conservative and was calibrated under H0 to control the FWER, the proportion of FP under H1 amongst selected variants was higher than the expected 20%.

This could reflect the correlations between polymorphisms we simulated.

To make more specific recommendations for study designs is difficult because it is closely related to the developed drug. In our simulations a late sample allowed larger information on the elimination phase to estimate CL. This result can be generalised to pharmacogenetic studies involving clearance and drugs with a long half-life. Taking a late sample requires to suspend treatment long enough to observe a decrease in concentrations, which may not be possible in patients from phase II trials.

In any case, it is essential that the sampling protocol, although limited, is as informative as possible to minimise the estimation error and shrinkage in individual parameters estimation.

The detection of genetic polymorphisms could highly benefit from the use of larger sample sizes through combined analysis and optimised design [START_REF] Bazzoli | Design evaluation and optimisation in multiple response nonlinear mixed effect models: PFIM 3.0[END_REF][START_REF] Combes | Prediction of shrinkage of individual parameters using the bayesian information matrix in non-linear mixed effect models with evaluation in pharmacokinetics[END_REF] .

In conclusion, this work confirmed the very limited likelihood that weak genetic effects can be detected in a typical phase I study, due to the small sample size. Such studies have to be considered only as hypothesis generating [START_REF] Bromley | Designing pharmacogenetic projects in industry: practical design perspectives from the Industry Pharmacogenomics Working Group[END_REF] . On the basis of our results in term of detection probability when analysing together data from phase I and phase II study, we claim that phase II is the best moment to identify the impact of genetic variants on drug response. It would be less efficient to start the study of the pharmacogenetics of a new drug in phase III trials or in post-marketing, because these take place too late in drug development [START_REF] O'donnell | Pharmacogenomics in early-phase oncology clinical trials: is there a sweet spot in phase II?[END_REF] and the new treatment could be administered in non-responders or expose subjects to high toxicities. Furthermore genetic subpopulations can be better targeted and potentially some subjects excluded from the study to increase the efficacy and reduce the risk of toxicity of the drug in these phase III studies.

STUDY HIGHLIGHTS

What is the current knowledge on the topic?

Most pharmacogenetic analyses in pharmacokinetic studies recently published included a limited number of subjects (less than 50). Previous simulations showed that such sample sizes result in low probability to detect polymorphisms. But with large number of subjects, extensive pharmacokinetic information is difficult to obtain in drug development.

What question did this study address?

This simulation study explored realistic ways to increase the amount of information by combining rich phase I data and sparse phase II data, and optimising such sparse designs.

What this study adds to our knowledge?

This study shows that even sparse data from phase II allow a marked improvement in the probability to detect genetic variants when combined with rich data from phase I, even more when sparse designs are optimised.

How this might change clinical pharmacology and therapeutics?

The pharmacogenetic analyses should be planned later in the drug development to take advantage of larger sample sizes by combining data which would increase the power to detect genetic effects.

TABLES

  of empirical family wise error rates (FWER) obtained without correction. b. Set of empirical FWER obtained after correction of type I error per tests. The 95% prediction interval around 20 for 200 simulated datasets is [14.5-25.5].

FIGURESFigure 1 .

 1 FIGURES

Figure 2 .

 2 Figure 2. True Positive Rate (TPR) versus False Positive Rate (FPR) under H 1 (top) and probability estimates (points) and 95% confidence interval (bars) to detect at least 𝑥 causal variants explaining the interindividual variability of CL (𝑥 = 1, … ,6) under H 1 (bottom) for main scenarios simulated with IIV CL = 25% (left) or modified scenarios simulated with IIV CL = 60% (right). Different symbols are used for each scenario, and colours denote the Lasso (grey) and the stepwise procedure (light blue).

Figure 3 .

 3 Figure 3. Distribution of the η-shrinkages on clearance for subjects in the phase I dataset (blue) and for subjects in the phase II dataset (brown), for each main scenario simulated under H 0 with IIV CL = 25% .

Figure 4 .

 4 Figure 4. Boxplot showing the loss of the signal for genetic effect in the overall population (top), as well as separately for the phase I data (blue borders) and for the phase II data (brown borders) (bottom). A boxplot is shown separately for each main scenario simulated under H 1 with IIV CL = 25% as a function of increasing 𝑅 𝐺𝐶 (boxplots colour).

Table 1 ,

 1 under two conditions: (i) no gene effect (H 0 ); (ii) gene effect on clearance CL (H 1 ). Under H 1 , 6 SNPs were drawn randomly without physiological assumptions or prior knowledge, and assumed to explain in total 30% of the IIV on CL through the following additive genetic model on the

	log-transformed CL:
	6
	𝑙𝑜𝑔(𝐶𝐿 𝑠𝑖𝑚 𝑖 ) = 𝑙𝑜𝑔(𝜇 𝐶𝐿 ) + ∑ 𝛽 𝑘 × 𝑆𝑁𝑃 𝑖𝑘	+ 𝜂 𝑖 𝐶𝐿	(1)
	𝑘=1

where 𝐶𝐿 𝑠𝑖𝑚 𝑖 is the simulated individual clearance, 𝜇 𝐶𝐿 the typical clearance, 𝛽 𝑘 the effect size associated to the variant allele of 𝑆𝑁𝑃 𝑖𝑘 and 𝜂 𝑖 𝐶𝐿 the interindividual random effect for clearance of subject 𝑖. Causal SNPs were different from one dataset to another. Assuming an additive genetic model, genotypes take values 0, 1 or 2, reflecting the number of mutated alleles. We chose this model to simplify the simulations but dominant or recessive genetic models could be easily simulated by changing genotype values. 𝛽 𝑘 was computed as a function of the coefficient of genetic component (𝑅 𝐺𝐶 𝑘 , the percentage of the interindividual variability in CL explained by the SNP) and the minor allele fraction (𝑝 𝑘 ), as follows:

Table 1 .

 1 Population values (µ) and interindividual variability (ω) for the model parameters of drug S used in the simulation study.

		Parameters	µ	ω (%)
		F a	Imax F D50 F	0.8 41.7	32.9
		FRAC b	Emax FRAC 0.45 D50 FRAC 18.6	-
		Tlag 1		0.401 35.1
		Tk0		1.59	31.6
		Tlag 2		22.7	-
		Ka		0.203	-
		V1		1520	-
		Q		147	89.9
		V2		2130	44.2
		CL		94.9	25.1
		σ slope (%)		20	-
					𝐼𝑚𝑎𝑥 𝐹 (𝑑𝑜𝑠𝑒 -20) 𝐷50 𝐹 + 𝑑𝑜𝑠𝑒 -20 , where dose is the
	amount administered.
	b. 𝐹𝑅𝐴𝐶 =	𝐸𝑚𝑎𝑥 𝐹𝑅𝐴𝐶 𝑑𝑜𝑠𝑒 𝐷50 𝐹𝑅𝐴𝐶 + 𝑑𝑜𝑠𝑒
	F: bioavailability ; FRAC: fraction of dose ; Tk0: zero order absorption duration ; Tlag 1 : lag time of zero
	order absorption ; Ka: first order absorption constant rate ; Tlag 2 : lag time of first order absorption ;
	V1: central compartment volume ; V2: peripheral compartment volume ; Q: intercompartmental
	clearance ; CL: linear elimination clearance

a. For doses < 20 units 𝐹 = 1, for doses ≥ 20 units 𝐹 = 1 -

Table 2 .

 2 Empirical estimates of Family-Wise Error Rate under H 0 for both association tests.

FWER (

%) Method SPI SPI/II 3s.96h SPI/II 3s.24h SPI/II 1s.24h

	Lasso	Without correction a	14	17.5	21.5	13.5
	Stepwise procedure Without correction a	20	18.5	22.5	15.5
	Lasso	After empirical correction b	20	19.5	21.5	19.5
	Stepwise procedure After empirical correction b

Acknowledgments

Adrien Tessier received funding from Institut de Recherches Internationales Servier. The authors thank Laurent Ripoll and Bernard Walther from Institut de Recherches Internationales Servier for their advices in pharmacogenetics. The authors would also like to thank Hervé Le Nagard for the use of the computer cluster services hosted on the "Centre de Biomodélisation UMR1137".