Short time existence of a quasi-stationary fluid-structure interaction problem for plaque growth

We address a quasi-stationary fluid-structure interaction problem coupled with cell reactions and growth, which comes from the plaque formation during the stage of the atherosclerotic lesion in human arteries. The blood is modeled by the incompressible Navier--Stokes equation, while the motion of vessels is captured by a quasi-stationary equation of nonlinear elasticity. The growth happens when both cells in fluid and solid react, diffuse and transport across the interface, resulting in the accumulation of foam cells, which are exactly seen as the plaques. Via a fixed-point argument, we derive the local well-posedness of the nonlinear system, which is sustained by the analysis of decoupled linear systems.


Model description.
In this article, we consider a quasi-stationary fluid-structure interaction problem for plaque growth, which describes the formation of plaque during the reaction-diffusion and transport of different cells in human blood and vessels. The problem is set up in a smooth domain Ω t ⊂ R 3 , with three disjoint parts Ω t = Ω t f ∪ Ω t s ∪ Γ t , where Γ t = ∂Ω t f , Ω t f ⊂ Ω t and Ω t f , Ω t s denote the domains for the fluid and solid, respectively. Γ t s = ∂Ω t stands for the outer boundary of Ω t , which is also a free boundary. Given T > 0, let As in [2,3,43], the flood is described by the incompressible Navier-Stokes equation where T f := −π f I + ν s (∇v f + ∇ ⊤ v f ) denotes the Cauchy stress tensor. v f : R 3 × R + → R 3 , π f : R 3 × R + → R are the unknown velocity and pressure of the fluid. ρ f > 0 stands for the fluid density and ν s represents the viscosity of the fluid.
Compared to the problems in [2,3,43], where the evolutionary neo-Hookean material was employed, the vessel in this manuscript is assumed to be quasi-stationary since it moves far slower than the blood from a macro point of view. Thus, we model the blood vessel by the equilibrium of a nonlinear elastic equation. To mathematically describe the elasticity conveniently, the Lagrangian coordinate was commonly used, see e.g. [18,19]. Thus, we set reference configuration as the initial domain which is defined by Ω := Ω t | t=0 , as well as Ω f = Ω 0 f , Ω s = Ω 0 s and Γ = Γ 0 , then Ω = Ω f ∪ Ω s ∪ Γ. Let X be the spatial variable in the reference configuration. Now we introduce the Lagrangian flow map ϕ : Ω × (0, T ) → Q T , with x(X, t) = ϕ(X, t) = X + u(x(X, t), t) (1.3) for all X ∈ Ω and x(X, 0) = X, where if X ∈ Ω s , denotes the displacement for fluid or solid. In the sequel, without special statement, the quantities with a hat will indicate those in Lagrangian reference configuration, e.g., u(X, t) = u(x(X, t), t), while the operators with a hat means those act on the quantities in Lagrangian coordinate. Then the tensor field F(x(X, t), t) =F(X, t) := ∂ ∂X ϕ(X, t) =∇ϕ(X, t) = I +∇û(X, t), ∀ X ∈ Ω, (1.4) withF f (X, 0) = I andF s (X, 0) = I +∇û 0 s is referred to be the deformation gradient and J =Ĵ := det(F) denotes its determinant. For the blood vessels, since the growth is taken into account, we impose the so-called multiplicative decomposition for the solid deformation gradientF s asF s =F s,eFs,g withĴ s =Ĵ s,eĴs,g , whereF s,e is the pure elastic deformation tensor andF s,g denotes the growth tensor, which will be specified later. For more details about the decomposition, see e.g. [18,22,35,43]. Inspired by Goriely [18, and Jones-Chapman [22,Section 3.2], a general incompressible hyperelastic material is considered for solid as − div T s = 0, in Q T s , (1.5) where T s := −π s I + J −1 s,e DW (F s,e )F ⊤ s,e stands for the Cauchy stress tensor. π s : R 3 × R + → R is the unknown pressure of the solid. The scalar function W : R 3×3 → R + is called the strain energy density function (also known as the stored energy density), which needs some general assumptions for the sake of analysis. for some constant C 1 > 0, where symF := 1 2 (F + F ⊤ ). Then for any a, b ∈ R 3 , one can derive the so-called Legendre-Hadamard condition by the Taylor expansion and the polar decomposition, which ensures that the operator − div D 2 W (I)∇ is strongly normally elliptic, see e.g. [33,Page 271].
Besides (1.5), the balance of mass should hold as well together with the growth, namely, (∂ t + v s · ∇) ρ s + ρ s div v s = f g s , in Q T s , (1.6) where ρ s > 0 is the solid density. v s := ∂ t u s denotes the velocity of the solid. The function f g s := γf r s with γ > 0 on the right-hand side of (1.6) stands for the growth rate, which comes from the cells reaction assigned later.
There are two essential assumptions when either the density or the volume does not change with growth. Generally, a constant-density type growth is assumed for the incompressible tissue, see e.g. [22,35], in this paper we assume this kind of growth as in [2,3,43]. Then by [18, (13.10)], one obtainŝ ρ s tr(F −1 s,g ∂ tFs,g ) = f g s , in Ω s × (0, T ). (1.7) As in [2,3,43], the plaque grows when the macrophages in the vessels accumulate a lot from the monocytes in the blood flow and turn to be the foam cells. Thus, denoting by c f , c s , c * s the monocytes, the macrophages and the foam cells respectively, we introduce where D f /s > 0 are the diffusion coefficients in the blood and the vessel, respectively. f r s := βc s with β > 0 stands for the reaction function, modeling the rate of transformation from macrophages into foam cells. Noticing that the foam cells is supposed to only move with vessels motion, one has (1.10) above. To close the system, one still needs to impose suitable boundary and initial conditions. For this purpose, we follow a similar setting as in authors' previous work [2]. On the free interface Γ t , one has the continuity of velocity and normal stress tensor for the fluid-structure part, while for the cells, the concentration flux is continuous and the jump of the concentration is determined by the permeability and flux across the vessel wall.
where n Γ t represents the outer unit normal vertor on Γ t pointing from Ω t f to Ω t s . For a quantity f , f denotes the jump defined on Ω t f and Ω t s across Γ t , namely, ζ denotes the permeability of the sharp interface Γ t with respect to the monocytes, which basically should depend on the hemodynamic stress T f n Γ t . It is supposed to be a constant for simplicity. Moreover, Γ t s is assumed to be a free boundary as well due to the physical compatibility, see e.g. [2,Remark 1.3]. Then we give the boundary conditions on Γ t s , T s n Γ t s = 0, on S T s , (1.13) D s ∇c s · n Γ t s = 0, on S T s . (1.14) Finally, the initial values are prescribed as In the case of 3d-3d fluid-structure interaction problem with a free interface, the strong solution results can be tracked back to Coutand-Shkoller [12], who addressed the interaction problem between the incompressible Navier-Stokes equation and a linear Kirchhoff elastic material. The results were extended to fluid-elastodynamics case by them in [13], where they regularized the hyperbolic elastic equation by a particular parabolic artificial viscosity and then obtained the existence of strong solutions by delicate a priori estimates. Thereafter, systems consisting of an incompressible Navier-Stokes equation and (damped) wave/Lamé equations were continuously investigated in e.g. [20,21,23,34] with further developments. Besides the models above, one refers to [8] for a compressible barotropic fluid coupled with elastic bodies and to [37] for the magnetohydrodynamics (MHD)-structure interaction system, where the fluid is described by the incompressible viscous non-resistive MHD equation and the structure is modeled by the wave equation of a superconductor material.
In the context of 3d-2d/2d-1d models, various kinds of models and results were established during the last twenty years. The widely focused case is the fluid-beam/plate systems where the beam/plate equations were imposed with different mechanical mechanism (rigidity, stretching, friction, rotation, etc.), readers are refer to e.g. [9,16,31,42] for weak solution results and to e.g. [5,14,17,25,28,30] for strong solutions respectively. Besides, the fluid-shells interaction problems were studied as well in e.g. [6,24] for weak solutions and in e.g. [10,11,27] for strong solutions. It is worth mentioning that in recent works [14,28], a maximal regularity framework, which requires lower initial regularity and less compatibility conditions compared to the energy method, was employed.
Recently inspired by [43], the authors considered the local well-posedness of a fluidstructure interaction problem for plaque growth in a smooth domain [2] and a cylindrical domain [3] respectively. By the maximal regularity theory, one obtains the wellposedness of linearized systems and then derives the existence of a unique strong solution to the nonlinear system via a fixed-point argument. However, concerning the origin of the model, which consists of blood flow and blood vessels, it is reasonable that the movement of vessels is supposed to be far slower than the blood, namely, the time scales of motions of the fluid and the solid are distinguished and the kinetic energy of the vessel is neglected. Thus, in this paper, we take an incompressible quasi-stationary hyperelastic equation into account, i.e., for every time t the elastic stresses are in equilibrium. To the best of our knowledge, this is the first result concerning on the strong solutions to a quasi-stationary fluid-structure interaction problem with growth.
1.3. Technical discussions. Under the above setting, the fluid-structure part is of parabolic-elliptic type, while the cells part is similar to the ones in [2,3]. To solve the nonlinear problem, our basic strategy is still a fixed-point argument in the framework of maximal regularity theory, while more issues come up when we consider the linearized systems. More precisely, for the linearization of the fluid-structure part, it is hard to solve it directly by the maximal regularity theory due to the parabolic-elliptic type coupling, which results in the unmatched regularity betweenv f andû s on the sharp interface Γ. To overcome the problem, one tries to decouple the system to a nonstationary Stokes equation with respect to fluid velocityv s and a quasi-stationary Stokes-type equation with regard to the solid displacementû s . Note that one key point is to separate the kinetic and dynamic condition on the interface correctly. Specifically, we impose the Neumann boundary condition forv s and a Dirichlet boundary condition forû s , see Section 4 below. Otherwise ifv f | Γ =v s , one may face the problem that there is no regularity information about the velocity of solidv s = ∂ tûs , since the solid equation is quasi-stationary without any damping.
Another issue is the choice of function spaces for the elastic equation, i.e., how to assemble suitable function spaces for the data in the linearized solid equation so that the regularity of∇û s matches with∇v s of fluid on the interface in the nonlinear system. Our choice is . This space is motivated by the observation that the nonstationary Stokes equation (4.1) is uniquely solvable if the Neumann boundary data holds, which implies that as well. Here H s q and W s q are Bessel potential space and Sobolev-Slobodeckij space respectively, which will be given later in Section 2.2. In fact, the anisotropic Bessel potential space we assigned is a sharp regularity if one goes back to the anisotropic trace operator, see Lemma 2.7 below, and it is natural to equip f with the regularity above. Because of this sharp regularity setting, one can not expect the Lipschitz estimates of nonlinear terms only with small time, which gives birth to an additional smallness assumption (3.7) on the initial solid displacement and pressure. Detailed discussion can be found in Remark 3.2 and Proposition 5.7 later.
To solve the quasi-stationary (linearized) elastic equation, one treated it as a Stokestype problem with respect to the displacementû s and the pressureπ s , due to the incompressibility. However, as we assigned the certain regularity space for it as above and the Stokes operator is not a standard one (namely, div(D 2 W (I)∇·)), one needs to consider a generalized stationary Stokes equation with f in L q and W −1 q,Γ 1 respectively, for which the maximal regularity of analytic C 0 semigroups is applied, as well as a complex interpolation method with a very weak solution in L q of a mixed-boundary Stokes-type equation, which can be solved by a duality argument.
For the cells part there is a problem of the positivity for the concentrations, compared to [2]. The idea in [2] to prove it is to apply the maximum principle to the original equation and deduce a contradiction with the help of Hopf's Lemma. However, due to the lack of regularity of v s = ∂ t u s , one can not expect it to be Hölder continuous in space-time, even continuous. To deal with this trouble, we make use of the idea of mollification, i.e., approximating v s by sufficient smooth functions v ǫ s such that t 0 v ǫ s → u s in certain space. Then arguing by a similar procedure in [2], we obtain an approximate nonnegative solution c ǫ . Finally one can show that it converges to a nonnegative function c, which exactly satisfies the original equations of cell concentrations.
1.4. Structure of the paper. The paper is organized as follows. In Section 2 we briefly introduce some notations and function spaces together with some corresponding properties. Moreover, a reformulation of the system is done in Section 3.1 and later we give the main result for the reformulated system. Section 4 is devoted to three linearized systems in Section 4.1, 4.2 and 4.3 respectively. The main results of this section are the L q -solvability for these linear problems, for which a careful analysis is carried out. In Section 5, we first introduce some preliminary lemmas in Section 5.1, which will be frequently used in proving the Lipschitz estimates later in Section 5.2. Then by the Banach fixed-point Theorem, we derive the short time existence of strong solutions to the nonlinear system in Section 5.3. Moreover, the cell concentrations are shown to be nonnegative, provided that the initial concentration is nonnegative. In addition, we establish the solvability of a Stokes-type resolvent problem with mixed boundary condition in Appendix A.

2.1.
Notations. Let us introduce some notations. For a vector u ∈ R d , d ∈ N + , we define the gradient as (∇u) ij = ∂ i u j . For a matrix A = (A ij ) d i,j=1 ∈ R d×d , we introduce the divergence of a matrix row by row as (div Then we denote by A : B := tr(B ⊤ A) = d j,k=1 B jk A jk the Frobenius product and |A| := √ A : A the induced modulus. For a differentiable and invertible mapping A : see e.g. [18,19]. Moreover, we define R d×d sym+ as the space of all positive-definite symmetric matrices and as the set of all proper orthogonal tensors.
Generally, B X (x, r) denotes the open ball with radius r > 0 around x in a metric space X. For normed spaces X, Y over K = R or C, the set of bounded, linear operators T : X → Y is denoted by L(X, Y ) and in particular, L(X) = L(X, X).
As usual, the letter C in the present paper denotes a generic positive constant which may change its value from line to line, even in the same line, unless we give a special declaration.

Function spaces and properties.
For an open set M ⊂ R d with Lipschitz boundary ∂M, d ∈ N + , we recall the standard Lebesgue space L q (M) and Sobolev space W m q (M), m ∈ N, 1 ≤ q ≤ ∞, as well as where q ′ is the conjugate exponent of q satisfying 1/q + 1/q ′ = 1. Furthermore, we set the Sobolev spaces associated with Γ ⊂ ∂M as , where X is a Banach space. In particular, L q (M) = W 0 q (M), L q (M; X) = W 0 q (M; X). Moreover, for 1 < q < ∞ and 1 ≤ p ≤ ∞ the standard Besov space B s q,p and Bessel potential space H s q coincide with the real and complex interpolation of Sobolev spaces respectively (see e.g. Lunardi [26], Runst-Sickel [36], Triebel [41]) In particular, W l q = H l q for l ∈ N and setting p = q above yields the Sobolev-Slobodeckij space Then we define the seminorm of W s q (M) as for f ∈ W s q (M) with 0 < s < 1, 1 < q < ∞. For an interval I ⊂ R, 0 < s < 1, 1 < q < ∞, we recall the Banach space-valued spaces K s q (I; X), K ∈ {W, H} with norm For convenience, with T > 0 we define the corresponding space with vanishing initial trace at t = 0 as Now we recall several embedding results and properties for Banach space-valued spaces that will be frequently used later.
Proof. The case K = W was shown in Simon [39,Corollary 17]. Then taking t satisfying r < t < s, one has , where the last inequality holds in virtue of H r q (I; X) ֒→ W t q (I; X) for r > t > 0. The second assertion can be easily derived by means of the observation and the definition of Sobolev-Slobodeckij space, so we omit it here. Lemma 2.2. Let 0 < s < 1, 1 < q < ∞ satisfying sq > 1, X be a Banach space and I = (0, T ) ⊂ R be a bounded interval for 0 < T < ∞. Then Proof. By Meyries-Schnaubelt [29, Proposition 2.10] with µ = 2 there, one has the first assertion and for K = W , 1/q < r < s, where C is independent of I. Then it follows from Lemma 2.1 that Adapting from Meyries-Schnaubelt [29, Proposition 3.2] and Prüss-Simonett [33, Section 4.5.5], we use the following time-space embedding lemma. Lemma 2.3. Let 1 < q < ∞, 0 < α, s < 2 and 0 < r < s , we have the embeddings In particular, . All these assertions remain true if one replaces W -and H-spaces by 0 W -and 0 Hspaces respectively, and the embedding constants in this case does not depend on T > 0. Now we include results on multiplication and composition.
Proof. The first part follows from Runst-Sickel [36, Theorem 5.5.1/1]. We note that in [36], the function spaces act on the full space R d . Here we just need to employ suitable extensions for Ω so that we can recover the case of full space. For the second part, let u, v be arbitrary two functions in K s , which completes the proof with (2.2) and the multiplication property Lemma 2.4 with s > d/p. Remark 2.6. We comment that for the case s = 1, the lemma above holds true as well due to [36].
In the following, an anisotropic trace lemma is introduced for a fractional power space. Lemma 2.7 (Anisotropic trace on the boundary). Let 1 < q < ∞ and Ω ⊂ R d , d ∈ N + , be a bounded domain with Γ := ∂Ω of class C 1 , T > 0, and where C > 0 is independent of T . Moreover, it is surjective and has a continuous right-inverse.
Proof. By means of a coordinate transformation and a partition of unity of Ω, one can easily reduce it to case of a half-space R d−1 × R + . Then thanks to [29,Theorem 4.5] with s = 1/2, m = 1, µ = 1 (see also [33,Proposition 6.2.4] with m = 1, µ = 1 and p = q there), one completes the proof.

Reformulation and main result
3.1. System in Lagrangian coordinates. In this section, we transform (1.1)-(1. 16) in deformed domain Ω t to the reference domain Ω, whose definition are given in Section 1.1. Let φ/φ be any scalar function in Ω t /Ω and w/ŵ be any vector-valued function in Ω t /Ω. Then one can easily derive the relations between derivatives in different configurations as whereF is defined as in (1.4).
Before the reformulation, we assume an isotropic growth, which is the simplest nontrivial form for the growth tensor. It is taken as a multiple of the identity, namely, F s,g =ĝI, in Ω s , whereĝ =ĝ(X, t) is the metric of growth, a scalar function depending on the concentration of macrophages. Note that there are other possibilities for growth, see e.g. [18,22], the isotropic one is taken for the sake of analysis. Then we haveĴ s,g =ĝ 3 indicating the isotropic change of a volume element and (1.7) becomes withĝ(X, 0) =ĝ 0 . Since it follows from (2.1) that we haveĴ By the decomposition ofF s and the incompressibility of the solid, we know thatĴ s,e = 1 andĴ Similar to [2], the reformulated system now reads as: s}, denotes the first Piola-Kirchhoff stress tensor associated with the Cauchy stress tensor T i defined in Section 1.1.

3.2.
Compatibility condition and well-posedness. Before stating our main theorem, one still needs to impose suitable function spaces and compatibility conditions. Following the general setting of maximal regularity, e.g. [2,3,33], where the basic space Analogous to [2], the compatibility conditions forv 0 f andĉ 0 read as Generally speaking, one does not need to assign any initial pressure for the Stokes equation. However, in this manuscript the coupling on the interface does lead to a condition on the initial fluid pressure since the solid equation is quasi-stationary and holds at t = 0. More specifically, we assume that there existsπ 0 for sufficiently small κ > 0, such that − div(DW (I +∇û 0 s )) +∇π 0 Here, the regularity forπ 0 f on the interface Γ is initiated from the matched regularity of∇v 0 f ,π 0 s and DW (I+∇û 0 s ) on Γ. Moreover, it coincides with the regularity ofπ f by the trace method of interpolation (see e.g. [33, Example 3.4.9]), i.e., Remark 3.2. In this paper, we need the smallness assumption of initial displacement to guarantee the estimates with respect to the deformation gradient, e.g. (5.7), which is a key element to derive the final contraction property of the certain operator. This is because we consider the general case ofû s | t=0 and linearize the elastic equation around the identity I, not the initial deformation gradient I +∇û 0 s . Specifically, one can not control (F s − I) by a small constant only with a short time. In particular, for the casê u 0 s = 0, one knowsF s | t=0 = I and hence the estimates later is uniform with respect to time T > 0. Moreover, for initial pressure it does also need the smallness due to the sharp regularity of pressure, see e.g. (5.15). Theorem 3.3. Let 5 < q < ∞ and κ > 0 be a sufficiently small constant. Ω ⊂ R 3 is the domain defined above with Γ, Γ s hypersurfaces of class C 3 . Assume that (v 0 f ,ĉ 0 ) ∈ D q satisfying the compatibility condition (3.6),π 0 .7) and (3.8). Then there is a positive [2,33], we prove Theorem 3.3 via the Banach fixed point theorem. To be more precise, we are going to linearize (3.5) in the first step, show the well-posedness of the linear system, estimate the nonlinear terms in suitable function spaces with small time and then constract a contraction mapping.
Remark 3.4. In fact, Theorem 3.3 still holds true in even more general dimensional case n ≥ 2 as long as q has a adapted restriction with respect to n. This is also an advantage of making use of maximal regularity theory.

3.3.
Linearization. Now following the linearization procedure in [2], we linearize (3.5) first, equate all the lower-order terms to the right-hand side and then arrive at the equivalent system:ρ  (2) The linearization is similar to the one in [2], but with several modifications, one of which is deduced above. It is possible to have other kinds of linearizations but we remark that in the present paper a divergence structure of divK s plays an essential role when we prove the linear theory and estimate it in a particular function space, see Corollary 4.4 and Proposition 5.7 later. Moreover, for the solid mass balance equation (1.6) (equivalently div v s = f g s /ρ s since the density is constant), we integrate it over (0, t) as (3.9e) to keep the Stokes-type structure for the elastic equation with respect to the displacementû s .
(3) Noticing that the continuity conditions (3.5f) on the interface are separated to (3.9c) and (3.9f) formally after the linearization, we remark here that this is for the sake of analysis due to the mismatch of the regularity on Γ. For instant, if one replaces (3.9c) with the boundary conditionv f = ∂ tûs , it has no chance to solve the fluid part since we have no first-order temporal derivative information for the solid displacementû s .

Analysis of the linear systems
In this section, we are devoted to solve the linear systems associated with (3.9).

Nonstationary Stokes equation.
Let Ω be a bounded domain with a boundary ∂Ω of class C 3− , T > 0. We consider the nonstationary Stokes equation where S µ (u, π) = −πI + µ(∇u + ∇ ⊤ u). ρ, µ > 0 are the constant density and viscosity. n denotes the unit outer normal vector on ∂Ω. Then we have the following solvability and regularity result, which can be adapted directly from e.g. Abels [ where P n := I−n⊗n denotes the tangential projection onto ∂Ω. For given data (f, g, h) with Moreover, there is a constant C > 0 independent of f, g, h, u 0 , T 0 , such that for 0 < T ≤ T 0

Quasi-stationary Stokes equation with mixed boundary conditions.
Let Ω be a bounded domain with a boundary ∂Ω of class C 3− , ∂Ω = Γ 1 ∪ Γ 2 consisting of two closed, disjoint, nonempty components. Consider the generalized stationary Stokes-type equation where n denotes the unit outer normal vector on ∂Ω. W : R 3×3 → R + is the elastic energy density such that H holds. Before going to the quasi-stationary case, we first investigate the weak solution and strong solution in L q -class of the stationary problem (4.2). (Ω), h 1 ∈ W 2+s−1/q q (Γ 1 ) 3 and h 2 ∈ W 1+s−1/q q (Γ 2 ) 3 . Then problem (4.2) admits a unique solution (u, π) ∈ W 2+s q (Ω) 3 × W 1+s q (Ω). Moreover, there is a constant C > 0 such that if s = 0 and when s = −1, Proof. First let s = 0, we reduce the system (4.2) to the case (g, h 1 , h 2 ) = 0. To this end, take a cutoff funtion ψ ∈ C ∞ 0 ((0, T )) such that Then ψ(t)g ∈ L p (0, T ; W 1 q (Ω)) ∩ W 1 p (0, T ; W −1 q,Γ 2 (Ω)), In view of Remark 1.1 and the maximal L q -regularity result for the generalized Stokes problems (e.g., [ in Ω, with 3 < p < ∞, 1 < q < ∞ to get a pair of solution (ũ,π) fulfilling Subtracting the solution to (4.2) with (ū,π), we are in the position to solve (4.2) with (g, h 1 , h 2 ) = 0, which can be referred to Theorem A.1 with λ = 0. Note that the case λ = 0 is applicable due to Remark A.2. Now we consider s = −1, namely the weak solution. In this case we only reduce (g, h 1 ) to zero since the Neumann boundary trace need to make sense on Γ 2 correctly.
Then one can subtract the solution of (4.2) with (ū,π) and solve (4.2) with reduced data (g, h 1 ) = 0 and modified (f, h 2 ) (not to be relabeled). The idea of the proof is to introduce a L q -class of very weak solution (see e.g. [15,38]), so that one can derive a solution with certain regularity in W 1 q (Ω) by complex interpolation, see e.g. Schumacher [38] for the stationary Stokes equation in fractional Bessel potential spaces.
For 1 < q, q ′ < ∞ satisfying 1/q +1/q ′ = 1, we define a generalized Stokes-type operator with respect to (4.2) as A q (u) := P q − div(D 2 W (I)∇u) for all u ∈ D(A q ), with D(A q ) = u ∈ W 2 q (Ω) 3 ∩ L q σ (Ω) : u| Γ 1 = 0, P n ((D 2 W (I)∇u)n) Γ 2 = 0 , where P q denotes the Helmholtz-Weyl projection onto L q σ (Ω), see e.g. [1, Appendix A] for the existence of the projection with mixed boundary conditions. P n := I − n ⊗ n is the tangential projection onto ∂Ω. By the result of s = 0 we see that A q : D(A q ) → L q σ (Ω) is well-defined and bijective. Then one knows that its dual operator ′ is bijective as well, which gives rise to the very weak solution. Note that A q and A * q ′ are consistent, namely, ij kl , i, j, k, l = 1, 2, 3. Then by the complex interpolation of operators, e.g. [38, Theorem 2.6], we record that Consequently, , for all w ∈ W 1,q ′ σ,Γ 1 (Ω). Moreover, by means of the open mapping theorem, one immediately deduce the estimate Up to now, one still needs to recover the pressure in the very weak sense, i.e., solving it is easy to verify that functional F defined above is well-defined in D(∆ q ′ ,DN ) ′ . For every u ∈ L q ′ (Ω), it follows from [33,Corollary 7.4.5] that there exists a unique solution ϕ(u) ∈ D(∆ q ′ ,DN ) satisfying ∆ϕ = u. Now we define π ∈ L q (Ω) by duality as a linear functional on L q ′ (Ω) acting for every u as π, u = F, ϕ .
(4.4) Indeed π is the very weak solution we are looking for, since for all ϕ ∈ D(∆ q ′ ,DN ), we have π, ∆ϕ = π, u = F, ϕ . The uniqueness can be showed by letting F = 0 in (4.4) so that for all u ∈ L q ′ (Ω), π, u = 0, which implies π = 0 a.e. in Ω. Then we have the estimate This completes the proof.
In Theorem 4.3, we consider the general case of data. In fact, for applications the right-hand side terms sometimes have special structure, which is of much help to derive the estimate in a concise form.
makes sense, which contributes to the nonlinear estimate for the fluid part. where D > 0 is a constant. u : Ω × (0, T ) → R stands for the system unknown, for example, the temperature, the concentration, etc. By e.g. [2,Proposition A.2], the existence and uniqueness result reads as Theorem 4.9. Let 3 < q < ∞ and T 0 > 0. Assume that u 0 ∈ W 2−2/q q (Ω) with the compatibility condition D∇u 0 | ∂Ω = g| t=0 holds. Given known functions (f, g) with regularity f ∈ F f (T ) := L q (0, T ; L q (Ω)), Then the parabolic equation (4.7) admits a unique strong solution u ∈ E(T ) where

Nonlinear well-posdeness
We denote by δ a universal positive function The most common example in the present paper is δ(t) = t θ for different θ > 0 from line to line.

Auxiliary lemmas.
In this section, we give some lemmas which we shall use later on.
Proof. First let us recall the equivalent definition of Bessel potential space f H s q (R;L q (Ω)) = f L q (R;L q (Ω)) + (−∆) s 2 f L q (R;L q (Ω)) , where the fractional Laplace operator is represented by the singular integral , see e.g. [40, Chapter V, Section 6.10]. For f, g ∈ H 1/2 q (R; L q (Ω)) ∩ L ∞ (R; L ∞ (Ω)) ∩ W α 2q (R; L 2q (Ω)) with q > 1 and 1/4 < α < 1/2, we see that the integrand has an algebraic decay rate greater than one with respect to |h|, which means the singular integral is actually integrable and one can omit the "lim" in the following. Then we have .
Dividing the region R into a neighborhood of the origin and its complement, we have where for every ε > 0 and 1/q + 1/q ′ = 1. By the Hölder's inequality, .
Lemma 5.6. Let T > 0, R > 0 and q ∈ (1, ∞). Ω s is the domain defined in Section Proof. By calculus, Then where we choose T R > 0 sufficiently small such that T Then we have the following Lipschitz estimates for the lower-order terms defined in (3.9).
Proposition 5.7. Let q > 5 and R > 0. There exist constants C, κ > 0 and a finite time T R > 0 both depending on R such that for 0 < T < T R and ∇û 0 Proof. For the estimates related to the fluid (with a subscript f ) except H 1 f and F j , j = 1, ..., 5, we refer to [ .

(5.15)
Estimate of G s . By means of Lemma 5.2 and 5.4, estimate in L q (0, T ; W 1 q (Ω s )) is clear. By the definition of G s , :∇û 2 s . Then with Lemma 5.2, 5.4 and the regularity (5.12), we have for q > 5 that It follows from the trace theorem tr Γ : whose Lipschitz estimate in Z 6 T can be controlled by δ(T ) w 1 − w 2 Y T thanks to Lemma 2.1 and 2.7, namely, . Now let us recall H 1 f = −K fnΓ +K snΓ . The first part can be addressed in the same way as in [2,Proposition 4.2]. By Lemma 2.7 the anisotropic trace theorem and C 3 interface that ensures an Γ of class C 2 , the second term can be estimated by Then together with Lemma 5.2-5. 6 and Assumption H, one gets Estimate of F j . We can estimate F 1 f and F j , j = 4, 5 analogously as in [2]. For others, since Y 5 , one can apply Lemma 5.2 and 5.4 to derive the corresponding estimates with q > 5, combining Lemma 2.7, the regularity ofF ⊤ s ,ĉ * s andĝ. This completes the proof.

5.3.
Nonlinear well-posedness. For a w defined in (5.13), define , where the elements are given by (3.9). Then the following proposition holds for M (w) : Y T → Z T , where Y T , Z T are given in Section 3.1, 5.2 respectively. Proposition 5.8. Let q > 5 and R > 0. Let w ∈ Y T be the function as in (5.13) with the associated initial data as w 0 as in (5.14). Then there exist constants C, κ > 0, a finite time T R > 0 both depending on R and δ(T ) as in (5.1) such that for 0 < T < T R , M : Y T → Z T is well-defined and bounded together with the estimates for w Y T ≤ R and ∇û 0 there exist a constant C > 0, a finite time T R > 0 depending on R and a function δ(T ) as in (5.1) such that for 0 < T < T R , Proof. The second part follows directly from Proposition 5.7. Then by setting w 2 = (0, 0, 0, 0, 0, 0, 1) in (5.17), one derives (5.16) immediately in view of the fact that M (0, 0, 0, 0, 0, 0, 1) = 0.
Now recalling the definition of solution and initial spaces in Section 3.1, we rewrite (3.9) in the abstract form where L (w) denotes the left-hand side of (3.9) and N (w, w 0 ) is the right-hand side. It follows from the linear theory in Section 4 that L : Y T → Z T ×D q is an isomorphism.
Proof of Theorem 3.3. For (v 0 f ,ĉ 0 ) ∈ D q satisfying the compatibility conditions, we may solve L (w) = N (0, w 0 ) by somew ∈ Y T . Then one can reduce the system to the case of trivial initial data by eliminatingw and we are able to set a well-defined constant , which can be verified to stay boundeded as T → 0 by the linear theories in Section 4 and the estimate (5.16), as in [2]. Choose R > 0 large such that R ≥ 2C L (v 0 f ,ĉ 0 ) Dq . Then For w i Y T ≤ R, i = 1, 2, we take T R > 0 and κ > 0 small enough such that C L C(R)(δ(T R ) + κ) ≤ 1/2, where C(R) is the constant in (5.17). Then for 0 < T < T R , we infer from Theorem 5.8 that  The uniqueness in Y T , 0 < T < T 0 , follows easily by repeating the continuity argument in [2, Proof of Theorem 2.1], so we omit it here. In summary, (3.9) admits a unique solution in Y T , equivalently, (3.5) admits a unique solution in Y T . Now we are in the position to prove the positivity of cells concentrations. Since the regularity ofv s = ∂ tûs is much lower than that in [2], we can not proceed as in [2] for c. To overcome this problem, we take a smooth mollificationv ǫ s ofv s for ǫ > 0 such that t 0v ǫ s (·, τ )dτ →û s , in Y 2 T , as ǫ → 0. Consider the problem in Q T f , ∂ t c ǫ s + div (c ǫ s v ǫ s ) − D s ∆c ǫ s = −f r s , in Q T s , with boundary and initial values as in Section 1.1. Then with same argument in [2, Proof of Theorem 2.1], one obtains 0 ≤ c ǫ (x, t) ∈ W 1 q (0, T ; L q (Ω t ) 3 ) ∩ L q (0, T ; W 2 q (Ω t \Γ t ) 3 ), which means there is a subsequent still denoted by c ǫ and a function c such that c ǫ ⇀ c weakly in W 1 q (0, T ; L q (Ω t ) 3 ) ∩ L q (0, T ; W 2 q (Ω t \Γ t ) 3 ), and c ǫ → c in D ′ (Q T \S T ), where D ′ (U) denotes the space of distributions on U, Q T , S T are defined in Section 1.1.
It is standard to verify c solves the same equation with v ǫ s replaced by v s . We only give the sketch of proof with respect to div(c ǫ s v ǫ s ) as an example.
for all φ ∈ D(Q T \S T ), φ ≥ 0, one concludes that c ≥ 0, a.e. in Q T \S T . The positivity ofĉ * s andĝ then follows automatically, as showed in [2], which completes the proof.
Appendix A. Stokes resolvent problem In this section, we give a short proof the solvability of the following Stokes resolvent problem with mixed boundary conditions. Let Ω ⊂ R 3 be a bounded domain of class C 3− with boundary ∂Ω = Γ 1 ∪ Γ 2 consisting of two closed, disjoint, nonempty components. Consider the resolvent problem λu − div(D 2 W (I)∇u) + ∇π = f, in Ω, div u = 0, in Ω, u = 0, on Γ 1 , (D 2 W (I)∇u − πI)n = 0, on Γ 2 , where n is the outer unit normal on the boundary, W : R 3×3 → R + is a scalar function with Assumption H holding.
Remark A.2. As D(A q ) embeds compactly into L q σ (Ω), the Stokes-type operator A q has compact resolvent. Therefore, its spectrum consists only of eigenvalues of finite algebraic multiplicity by spectral theory for compact operators (see e.g. Alt [4]), and is independent of q by Sobolev embeddings. So it is enough to investigate these eigenvalues for the case q = 2. Let ω be the eigenvalue of −A 2 . Employing the energy method, Lemma 1.1 and the Korn's inequality, we have which shows that ω is real and nonpositive. Since Ω is bounded and Γ 1 is assumed to be nonempty, the Poincaré's inequality is valid and one obtains ω Ω |u| 2 dx ≤ −C Ω |u| 2 dx, which implies w ≤ −C, where C > 0 does not depend on λ. Then one conclude that λ 0 defined in Theorem A.1 is negative and hence if λ = 0, the Theorem still holds true.