Modeling and analysis of adipocytes dynamic with a differentiation process

. We propose in this article a model describing the dynamic of a system of adipocytes, structured by their sizes. This model takes into account the diﬀerentiation of a population of mesenchymal cells into preadipocytes and of preadipocytes into adipocytes; the diﬀerentiation rates depend on the mean adipocyte radius. The considered equations are therefore ordinary diﬀerential equations, coupled with an advection equation, the growth rate of which depends on food availability and on the total surface of adipocytes. Since this velocity is discontinuous, we need to introduce a convenient notion of solutions coming from Filippov theory. We are consequently able to determine the stationary solutions of the system, to prove the existence and uniqueness of solutions and to describe the asymptotic behavior of solutions in some simple cases. Finally, the parameters of the model are ﬁtted thanks to some experimental data and numerical simulations are displayed; a spatial extension of the model is studied numerically

Abstract. We propose in this article a model describing the dynamic of a system of adipocytes, structured by their sizes. This model takes into account the differentiation of a population of mesenchymal cells into preadipocytes and of preadipocytes into adipocytes; the differentiation rates depend on the mean adipocyte radius. The considered equations are therefore ordinary differential equations, coupled with an advection equation, the growth rate of which depends on food availability and on the total surface of adipocytes. Since this velocity is discontinuous, we need to introduce a convenient notion of solutions coming from Filippov theory. We are consequently able to determine the stationary solutions of the system, to prove the existence and uniqueness of solutions and to describe the asymptotic behavior of solutions in some simple cases. Finally, the parameters of the model are fitted thanks to some experimental data and numerical simulations are displayed; a spatial extension of the model is studied numerically.

Introduction
Obesity is a worldwide major public health issue that doubled since 1980 and affects nowadays almost two billions of adults considered as overweight and 600 millions considered as obese [12]. Strikingly, obesity is the most prevalent cause for the development of cardio-metabolic diseases (cardiovascular diseases, type 2 diabetes, and liver diseases) as well as cancer, increasing mortality and morbidity and justifying the need for intensive research and intervention policy.
Obesity is characterized by an increase in adipose tissue (AT) mass. This expansion of AT is a complex process which requires succeeding steps of proliferation, differentiation and maturation of the cells from the adipocytes lineage [2]. Indeed, within AT, vascular-resident adipose progenitor cells (APCs) proliferate and, under specific signals, differentiate into pre-adipocytes. Pre-adipocytes also expand through proliferation before to differentiate into small mature adipocytes. Mature adipocytes have an impressive capacity of expanding their volume by more than 30-fold through triglyceride accumulation in lipid droplets [7,16]. Ultimately, one single large lipid droplet occupies most of the cytoplasm and stiffens locally the AT [20]. The mechanical forces generated by hypertrophic stiff adipocytes may both limit their size and stimulate the differentiation of APCs and pre-adipocytes to recruit new adipocytes [15]. Mature adipocyte size is not only critical for adipogenesis initiation but also for adipocyte functions. Indeed, hypertrophic adipocytes have less ability to properly store lipids, resulting in spillover of lipids and excessive fat deposition in other tissues, both favorizing cardio-metabolic diseases. Although there is growing evidence that impaired AT expandability plays a pivotal role in obesityrelated cardio-metabolic diseases, the molecular and cellular basis of this phenomenon is complex and far from being understood. Indeed, the expansion of AT depends on a large number of parameters including the rate of APCs and pre-adipocytes proliferation/differentiation/death, the mechanical feedback loop of adipogenesis stimulation, the size of the adipocytes and the kinetic of their death.
The pathological implication of these phenomena is a motivation for developing new approaches for a better understanding of the AT formation. Mathematical modeling can provide useful insights, in particular for identifying leading parameters. Noticeably, the mechanical feedback loop of adipogenesis stimulation is certainly a pivotal parameter that could control the ability of adipose tissue to expand through the recruitment of new adipocytes. We refer the reader to [17,18] for attempts in this direction. Here, we shall present a different modeling of the AT development, including some differentiation processes and a velocity growth depending on the total surface, and investigate, theoretically and numerically, the main features of the adopted model.

Description of the homogeneous in space model
We think of the process as a compartment model with three populations of cells: mesenchymal, pre-adipocytes and adipocytes, the latter population being structured by the size of the cells. Let t → m(t) and t → p(t) stand for the number of mesenchymal, and pre-adipocytes, respectively; in the first five sections, both quantities depend only on the time variable t ≥ 0. Adipocytes are described by their radius distribution function (t, r) → a(t, r): for given 0 ≤ r 1 ≤ r 2 , the integral r2 r1 a(t, s) ds gives the number of adipocytes with a radius r between r 1 and r 2 , at time t ≥ 0.
Mesenchymal and pre-adipocytes undergo proliferation and mortality. Furthermore, mesenchymal cells differentiate into pre-adipocytes, while pre-adipocytes differentiate into adipocytes. The mutation of pre-adipocytes gives rise to adipocytes with radius r ≥ 0, and r will be the minimal radius within the adipocytes population. Adipocytes can undergo mortality and dynamic change of their radius. Radius changes are modeled with a growth rate function (t, r) → V (t, r). The expression of the growth rate function V will be determined later on, through volume considerations. As it will be detailed below, we assume the existence of a critical size r c > 0, which stands for the maximal value of adipocyte radius.
The unknowns depend on the time variable t ∈ [0, +∞) and, for the adipocyte distribution, on the radius variable r ∈ [r , +∞). The evolution of the population of mesenchymal cells, pre-adipocytes and adipocytes is governed by the following set of equations, defined on the domain (t, r) ∈ [0, +∞) × [r , +∞) where all the parameters of the model are nonnegative and can be collected as follows: α, α proliferation rates for m and p, β, β differentiation rates for m and p, which depend on the mean radius of adipocytesr(t) at time t, γ, γ , γ mortality rates for m, p and a, r m radius of mesenchymal cells, r emergence radius of adipocytes (minimal radius) and radius of pre-adipocytes, r c maximal radius of adipocytes, V (t, r) growth rate of an adipocyte of radius r at time t. Denoting by δ = α − γ and δ = α − γ , we obtain the following system: Initial data and boundary conditions. The system is complemented by initial data Since V (t, r ) is positive, we also need to prescribe the boundary condition for a when r = r ; this is where we take into account the differentiation of the pre-adipocytes into adipocytes: For the largest adipocytes, the growth rate vanishes (see below), and we simply assume that V (t, r) a(t, r) = 0 for r > r c . It means that adipocytes beyond a certain size do not exist. Consistently, we also assume that the support of a 0 is included in [r , r c ]. The initial and boundary value problem under consideration is (1a), (1b), (1c), (4), (5), (6), (2a), (2b), (3).
Dependency of the mutation rates on the mean radius. To keep feedback mechanisms description tractable, we assume that the differentiation rates β and β are functions of the mean radiusr(t) at time t. As said above, rc r a(t, s) ds represents the total number of adipocytes, at time t. The mean radius of adipocytes at time t is therefore equal to:r (t) = rc r s a(t, s) ds rc r a(t, s) ds (4) and we also introduce for further purposes the total surface of the adipocytes at time t.
We define β and β as functions ofr, typically with a sigmoid shape that reproduces threshold effects: with all parameters positive and such that β M ≥ β m , β M ≥ β m . The image of β is [β m , β M ]; the sigmoid β is "centered" on r β and its slope at mid height is β M −βm 4R β . The behavior of β relatively to its parameters is similar, see Figure 1 for a typical shape. In particular, from a mathematical point of view, β and β are bounded  Description of the growth rate. The dynamics of the radius of an adipocyte t → R(t) with respect to time can be approximated with the ODE d dt where V is the growth rate. We assume that adipocytes capture all the excess of food, i.e. not used by the metabolism. Adipocytes gather the food through their membranes, thus the flux of food they receive is proportional to their surface. Given a radius r < R < r c the flux of food is proportional to the ratio of the surface of the considered adipocyte over the total surface of all the adipocytes, that is where k > 0 is the (excess of) available food. In the present paper, k is constant in time, but it could be relevant to consider it as evolving in time, for example with time periodic food input to mimic circadian cycle. Thus, in the time interval [t, t + dt], the volume variation of such an adipocyte with radius R(t) is governed by Letting dt go to 0, we deduce that the volume obeys the ODE Consequently, we have This relation holds as far as the radius is not too large. As already mentioned above, the size of the adipocytes is limited. There are several possibilities to model such a threshold, based either on phenomenological arguments, or on energetic considerations. In what follows, we make the equation as simple as possible, with the formula which assumes that the growth rate vanishes outside the domain r ∈ [r , r c ]. This simple definition will allow us to derive easily interesting formula for the stationary solutions of the model. However, for the numerical simulations, we adopt the following regularized growth rate : where > 0 is a small parameter. Typical values for the parameters are collected in Table 2.

Definition of the solutions to system (1)-(3)
The definition (6) of the growth rate V raises some slight technical difficulties. Indeed, the equation for the adipocytes concentration (t, r) → a(t, r) is a transport equation, which is classically understood by means of the characteristic curves associated with the growth rate V (t, r). However, V , as given by (6), does not fullfil the regularity required to apply the Cauchy-Lipschitz theorem. Nevertheless, assuming that S(t) ≥ > 0 for all t ∈ [0, T ), V satisfies the following one sided Lipschitz condition (OSLC): Hence, we can appeal to the generalized theory introduced by A. Filippov [10], see also [14] for an application to transport equations with discontinuous coefficients. It allows us to consider the family of continuous Filippov maps satisfying X s (s, r ) = r for all s ∈ [0, t], X 0 (0, r) = r for all r ∈ [r , +∞) and, for fixed (s, r) ∈ (0, ∞) × (r , ∞), the function t → X t (s, r) is absolutely continuous on [s, ∞) and satisfies the differential equation d dt X t = V (t, X t ) for almost every t with Cauchy datum X s (s, r) = r. Owing to the OSLC, X t is the unique flow of the differential equation for t ≥ s.
In our particular case, the characteristic maps associated with a given S can be explicitly computed and are defined by: Note that any trajectory issued from r ∈ [r , r c ] do not cross the line r = r c : sizes larger than r c cannot be reached, which corresponds to the intuition since the velocity vanishes of (r c , ∞). Moreover, there is no need of boundary condition for the transport equation for the adipocytes at r = r c since there is nothing entering from this end into the domain. Let a 0 be a positive finite measure supported on [r , r c ]. We denote by M − w * the set of finite measures on [r , +∞) endowed with the weak star topology.
We generalize Eq.(3) by replacing its right-hand side with a given positive continuous function f ∈ C 0 (R + ). In that context, we will call a solution to the problem (1c)-(2b)-(3) with initial condition, i.e.
any measure-valued function a ∈ C 0 (R + ; M − w * ) such that the following Duhamel formula holds: for all φ ∈ C 0 b ([r , +∞)) and for all t ∈ [0, T ), where T = inf{t : S(t) = 0}. Note that here, a t denotes the measure a at time t and the condition a ∈ C 0 (R + ; M − w * ) means that lim h→0 φ(r) da t+h (r) = φ(r) da t (r) holds for any continuous and bounded trial function φ. In the following, we shall use the definition (9) with f (t) = β (r(t))p(t). Moreover, we bear in mind that the problem is non linear, since the growth rate V depends on the total surface S(t) of the adipocytes, see (5). It turns out that the OSLC on which the construction is based relies on the positivity of S, that we are going to discuss now, showing that the lifespan of solution is infinite (for positive times). Bearing in mind the physical meaning of the unknown, the data f and a 0 are non negative, with, furthermore, supp(a 0 ) ⊂ [r , r c ]. Formula (9) then tells us that a t is a non negative measure too. We also remark that the characteristics issued from [r , r c ] cannot exceed r c (this is a consequence of the well-posedness of this Filippovtype Cauchy problem) and Accordingly, the support of a t remains in [r , r c ] for all t ∈ [0, T ). Finally, let us discuss formally that the model does not produce a shrinking of the surface, that would be an obstacle to the global existence of solutions. Since a t is compactly supported, we can use (9) with φ(r) = 4πr 2 to compute the surface as: This formula enables to bound from below the value of S(t) as The continuity of S, which comes from the continuity of t → a t , enables to conclude that T = +∞.

Outline of the article
In this article, a complete study of the model is proposed: in section 1, we begin with a computation of the stationary solutions of the model, followed by a proof of the existence and uniqueness of solutions in Section 2. The asymptotic in time behavior of the system is described in Section 3 and some numerical simulations are displayed in Section 4. Finally, this model is extended by considering space heterogeneities through a coupling with a fluid environment and numerical simulations are presented in Section 5.

Stationary solutions
We will use the previous framework to exhibit stationary solutions of the system (1)-(3). Finding a stationary solution (m, p, a 0 ) is equivalent to finding a solution for the following system: for all φ ∈ C 0 b ([r , +∞)) and for all t ∈ [0, +∞), with unknowns a 0 a positive finite measure supported on [r , r c ] and m, p ∈ [0, +∞). The last condition expresses the fact that a t = a 0 for all t ∈ [0, +∞) where a is a solution in the sense given in the previous section.
Let (m, p, a 0 ) be a solution to (11) and V = k S . The characteristics X S t are computed using the constant function S = 4π rc r r 2 da 0 (r) and are well defined since a 0 is positive, that is to say: We can hence rewrite the last equation of system (11) Using the equality a(r ) = β (r)p V , we deduce the form of a 0 as: where µ can be determined from (12) through the relation: and finally we find Therefore the stationary solutions satisfy with unknowns p, m ∈ R + , a(r ) ∈ R * + . Note that the expression for a 0 enables to computer and S as functions of V , i.e.: Proposition 1. System (13) has a non-trivial (i.e. with a non-zero a 0 ) solution if and only if one of the following two conditions hold: Proof. We begin with noticing that p cannot be zero, since it would imply that a(r ) = 0 and therefore a 0 = 0. Now, remark that the function Ψ such thatr = Ψ(V ) and defined by Ψ : Assume that there exists a non-trivial solution. If m = 0, we should impose that δ ∈ β ((r , r c )). If m = 0, then m > 0 (we are interested in physical solutions) and it implies that δ ∈ β((r , r c )), namely δ = β(r) with r such that δ − β (r) = −β(r)m/p < 0. thus if there exists a non-trivial solution, one of the two conditions is satisfied.
Conversely, if δ ∈ β ((r , r c )) , we find a stationary solution defined bȳ From Eq. (14), we can deduce the value of a(r ) knowing S and V , and therefore a 0 and p = V a(r ) β (r) .
Similarly, if δ ∈ β ((r , r c )) and δ − β β −1 (δ) < 0, then a solution is given bȳ Remark 1. The stationnary solution may not be unique. Precisely, if exactly one condition is fulfilled, the stationary solution is unique. However, if both conditions are fulfilled, we may have two different values ofr and thus two different stationary solutions: one with m = 0 and another one with m > 0. If one sees stationary solutions as asymptotic in time unsteady solutions, the non-uniqueness of stationary solutions is understood as the existence of different limit profiles with different initial conditions. We do not address here the problem of the stability of the stationary solutions.

Existence and uniqueness of unsteady solutions
Equations of system (1)-(3) are coupled altogether as follows: the ODE part of system (1), that is to say equations (1a) -(1b) is a simple linear ODE system, but which coefficients are defined through the mean radius functionr of adipocytes, defined at Eq.(4). We therefore need to know function a to computer. Conversely, the transport part of system (1), that is to say equation (1c), depends on the total surface S of a through the velocity V defined at Eq. (6) and depends on function p solution to the ODE part, through the boundary condition (3).
In this section, we will prove the existence and uniqueness of solutions to system (1)-(3). To do so, we will proceed in three steps: first we prove, in Subsection 2.1, the existence and uniqueness of solutions to the ODE part (1a) -(1b)-(2a). Then, in Subsection 2.2, we prove the existence and uniqueness of solutions to the transport equation (1c)-(3)-(2b) for a given non negative flux by solving a fixed point equation. We also prove some stability property with respect to the flux that will be used later. Finally, in Subsection 2.3, we prove the existence and uniqueness of a solution to the general system thanks to a fixed-point theorem, coupling the results obtained at the foregoing subsections.
2.1. Some preliminary results on the solution of the ODE part (1a) -(1b)-(2a) First we will give some existence and stability results on the solution of the ODE part (1a) -(1b) of system (1) with initial conditions (2a), assuming that the mean radiusr ∈ C 0 ([0, T ], [r , r c ]) is given.
Since β and β are bounded, the following function: is Lipschitz-continuous. We call L its Lipschitz modulus. Therefore, the system satisfies the hypotheses of Cauchy-Lipschitz Theorem and we have the following bound: Moreover, we can prove the following stability property of the solution with respect tor, where p(r) denotes the solution p computed with the mean radius functionr: for allr 1 Indeed, let us denote by M the following Lipschitz-continuous function of Lipschitz modulus k: We obtain the following inequality, using (15): and we conclude thanks to Grönwall's Lemma.

Solution to the transport equation (1c)-(3)-(2b) with a given non negative flux
In this section, the first theorem gives an existence and uniqueness result for the transport part (1c) of system (1), complemented with boundary condition (3) and initial condition (2b).
Remark that the velocity V , involved in the transport equation (8) and defined at Eq. (6), depends on the total surface S, which is computed thanks to the solution a of Eq. (8). We are therefore led to consider the following fixed-point problem, coming from the expression of the surface S as a function of the characteristic curves X S t , see Eq. (10): where It is easy to prove that solutions to Eq. (17) and solutions to Eq. (8) are the same, defining a by the following formula: for all t ∈ [0, T ] and for all φ ∈ C 0 b ([r , r c ]). Now, let us prove the existence and uniqueness of solutions to system (8).
, has a unique solution.
Proof. As explained before, it is enough to prove the existence and uniqueness of solutions to the fixed-point equation (17).
To do so, we define the following operator Γ on and we prove that Γ is a contraction.
We can prove easily that the range of Γ is and, for any µ > 0, we define the following norm on E: Then, the following inequality holds, for all S 1 , S 2 ∈ E and for all t ∈ [0, T ]: which comes from the fact that Indeed, for any λ, µ, ν ∈ R, one has Thus, if we take µ large enough, the operator Γ is a contraction of the complete metric space (E, . E ) and the Banach fixed point Theorem proves the existence and uniqueness of the solution to equation (17).
We will denote by S(f ) the unique solution to equation (17), defined at Theorem 1, with the positive flux condition f . We prove now a stability result, that is to say an estimate of the quantity Theorem 2. Let M ∈ R + . We have the following estimate: there exists a continuous function Proof. We denote by Γ 1 (respectively Γ 2 ) the operator defined at Eq. (18) with the flux f 1 (respectively f 2 ). We decompose the difference S(f 1 ) − S(f 2 ) as follows: Using inequality (20), we can bound from above the first term of the right-hand side as follows: We now control the second term of the right-hand side Γ 1 (S(f 2 )) − Γ 2 (S(f 2 )) E , i.e.: and therefore, using the value of µ and the definition (19) of the norm, For any t ∈ [0, T ], applying the previous inequality to the restriction of the solutions to [0, t], we obtain (1)-(3). Now let us prove the following theorem, that states the existence and uniqueness for the full system (1)-(3), coupling the results of the previous two subsections.

Existence and uniqueness of solutions to system
Once again, this is equivalent to a fixed-point problem. Indeed, the resolution of the ODE part (1a) -(1b) requires the knowledge of the mean radiusr, which is computed from the solution a of the transport equation (1c). In turns, the boundary condition (3), which is necessary to compute a, involves the function p, solution of the ODE part.
More precisely, we consider a solution (m, p, a) to system (1)-(3). Knowing a, we can compute the function S thanks to formula (5) and therefore the characteristics X S t associated to system (1c)-(2b)-(3) thanks to formula (7).
We can consequently writer(t) in the form of a fixed point equation, using formula (4) and Eq.(9) with f (t) = β(r(t))p(t), that is to saȳ Conversely, if we find a functionr, solution to Eq.(21), we can then solve the ODE part (1a) -(1b)-(2a) of system (1) in order to obtain functions m and p, see Subsection 2.1. We can then deduce surface S from the fixed point equation (17) with f (t) = β (r(t)) p(t), that is to say, following the notation of Sec. 2.2, S = S(β (r) p). Solution a is finally given by Eq. (9).
We are therefore reduced to find a solution to the fixed point problem (21) with p solution to Eq. (1b) and S = S(β (r) p).
Now, let us prove the existence and uniqueness of solutions to system (1)-(3). Proof. As explained before, we are reduced to find a solution to the fixed point problem (21) with p solution to Eq. (1b) and S = S(β (r) p). To do so, let 0 < T < ∞ and define the following operator: Λ : where S = S(β (r) p).
To prove the contraction of operator Λ, we will use the same idea and the same norm (19) than previously. We will also need the stability results demonstrated at subsections 2.1 and 2.2.
Let F = C 0 ([0, T ], [r , r c ]). We will prove in the following that there exists A > 0 (which may depend on T ) such that for allr 1 ,r 2 ∈ F and for all t ∈ [0, T ], Therefore, and finally, if we take µ large enough, Λ is a contraction of the complete metric space (F, . F ). Thus the Banach fixed-point theorem gives the existence and uniqueness of the solution to the fixed point equation (21). Now let us prove inequality (23). Considering the definition (22) of operator Λ, it is sufficient to prove that the following three functions are Lipschitz-continuous and bounded with respect tor in L 1 ([0, t]): . First, using the fact that β is a bounded Lipschitz-continuous function, the bound (15) and the stability property (16) regarding the solution p of the ODE part (1b), it is straightforward to prove that for allr 1 where p(r) denotes the solution p of Eq. (1b) computed with the given mean radius functionr.
Since β (r 1 )p(r 1 ) and β (r 2 )p(r 2 ) are bounded, see the properties of β and Eq. (15), we can use Theorem 2 and inequality (24) to obtain that which implies that functionr → e −γ t rc r X S t (0, r) da 0 (r) is Lipschitz-continuous. Function (3): A similar proof enables to prove that functionr → t 0 e −γ (t−s) X S t (s, r )β(r(s))p(s) ds is bounded from above and Lipschitz-continuous.

Asymptotic behavior of the solutions in some simple degenerate cases
In this section, the behavior of the solutions for large times is studied. We first prove that, under some conditions involving the growth rate and the differentiation rate, the number of mesenchymal cells (resp. of preadipocytes) goes either to zero or to infinity. We distinguish two situations depending on whether or not the death rate γ of the adipocytes vanishes. In the case γ = 0, we prove that the mean radius converges towards the minimal radius as time goes to infinity. In the case γ = 0, we prove that the number of adipocytes converges towards a Dirac mass at r c .
In some simple cases, we can therefore study the asymptotic behavior of the solution, solving explicitly Eq. (1a) and (1b) as: For the proof of the asymptotic properties, we also introduce the zeroth and the first momentum of a t , that is to say: Note thatr(t) = M 1 (t) M 0 (t) . We also have the following formula linking the zeroth and the first momenta of a: We have therefore the following proposition regarding the asymptotic behavior of m and p: Proposition 2. Let m 0 ≥ 0, p 0 ≥ 0, a 0 ∈ M − w * with a 0 positive. We denote by m, p, a t the solution given by Theorem 3. We have the following limits when t → +∞: Moreover, If δ − β(r ) < 0 and δ − β (r ) < 0, then there exist C > 0 and > 0 such that Proof. First of all, as r ≤r(s) ≤ r c for any s ∈ [0, +∞), we have the following inequality on m: which gives a limit for m in the two cases δ − β(r c ) > 0 and δ − β(r ) < 0. In the same manner, when m → 0, we can use the explicit expression to bound p as follows: Now, we denote by ρ = −(δ − β(r )) and ρ = −(δ − β (r )) and we assume that ρ > 0, ρ > 0. Thanks to (28), we obtain that m(t) ≤ m 0 e −ρt , ∀t ∈ [0, +∞), and we use this bound in the formula for p, leading to: Table 1. Summary of the asymptotic behaviors of functions t → (m(t), p(t)) and therefore The previous proposition can be summarized in Table 1, where the couple in each cell represents lim t→+∞ m, lim t→+∞ p . .

Case γ = 0
In the case when γ = 0, we can make precise the result on the mean radiusr. More precisely, if the growth of p is exponential, which holds in particular if δ − β (r c ) > 0, the total surface increases and the velocity decreases; meanwhile, the flux of mass in r increases and therefore a large proportion of the mass of adipocytes stays around r . Proposition 3. We assume that γ = 0. We consider m, p, a t the solution given by Theorem 3 and M 0 the zeroth momentum of a defined at Eq.(25). Then, we have the following limits when t → +∞: Moreover, if we assume that there exist > 0 and C > 0 such that for all t ∈ [0, +∞), p(t) ≥ Ce t , then lim t→+∞r (t) = r .
In particular, for t large enough, there exists C > 0 such that S(t) ≥ Ce t/2 and therefore, since r < r c , we have for t large enough that: Now, from the expression of the mean radius (21) and of the characteristics (7), we get the following equivalent when t → +∞:r

Case γ = 0
Now, in this section, we consider the case γ = 0. At first sight this assumption on the death rate might appear biologically questionable. However, it is very relevant. Indeed, adipocytes have usually a very low mortality rate (for humans, the renewal of half of the population takes about 8 years), see [19]. This could be different for certain specific experiments performed on mice, due to experimental conditions that can enhance the adipocytes mortality.
The results on p and m of the previous section still hold but the results on M 0 are different. We prove that in the case when δ − β(r ) < 0 and δ − β (r ) < 0, the number of adipocytes converges towards a Dirac mass at r c . Proposition 4. We assume that γ = 0, δ − β(r ) < 0 and δ − β (r ) < 0. We consider m, p, a t the solution given by Theorem 3 and M 0 the zeroth momentum of a defined at Eq.(25). Then, we have the following limits when t → +∞: Proof. The integrability of β (r(s))p(s) comes from the exponential decay of p, see Eq. (27), and the limit for M 0 (t) when t → +∞ comes directly from the expression of M 0 at Eq.(26) with γ = 0.
Let us denote this limit by M 0,∞ = M 0 (0) + +∞ 0 β (r(s))p(s) ds. Since the total surface S(t) ≤ r 2 c M 0,∞ , we can deduce that S is bounded from above and therefore that +∞ 0 k S(u) du = +∞. Therefore, for t large enough, we can define s c (t) as the unique real such that Since S is bounded from above, we can deduce that lim t→+∞ s c (t) = +∞.
Let φ ∈ C 0 ([r , r c ]). For t large enough , Eq. (9) can be written as: β (r(s))p(s) ds, Figure 2. Discretization of the interval. The variable a is computed in the middle of each interval, as well as the velocity.
we can compute the limit and we get: which ends the demonstration.

Numerical scheme
To find a numerical solution to the system of equations, we apply an explicit Euler scheme for the time approximation and an upwind scheme for the discretization of the transport-like equation of the adipocytes population. We consider a grid with a uniform radius-step ∆r > 0 (see Figure 2) and J intervals. The time step is denoted by ∆t. The time step is not constant and it is updated at each iteration according to the Courant-Friedrichs-Lewy stability condition, that will be described later on. Nevertheless, for the sake of simplicity, we denote it by ∆t instead of ∆t n . Let us denote by m n , p n the approximations of the solutions m(t n ) and p(t n ) at time t n . Moreover, a n j stands for the approximation of a(t n , r j ) at time t n and point r j and V n j+1/2,j∈{0,...,J} is the approximation of the velocity at time t n at the cells boundaries. Given the solutions (m n , p n , a n j∈{1,...,J} ) and the velocity V n j+1/2,j∈{0,...,J} at time t n , the solution at time t n+1 is updated by m n+1 = m n + ∆t(α(r n ) − γ)m n − ∆tβ(r n )m n , p n+1 = p n + ∆t(α (r n ) − γ )p n + ∆tβ(r n )m n − ∆tβ (r n )p n , a n+1 j = a n j − ∆t ∆r (V n j+1/2 ) + a n j − (V n j+1/2 ) − a n j+1 − (V n j−1/2 ) + a n j−1 + (V n j−1/2 ) − a n j − ∆tγ a n j , j ∈ {1, ..., J}, V n 1/2 a n 0 = β (r n )p n .
Here above, for any real number V , (V ) + and (V ) − stand for the positive and negative parts of V : (V ) + = max(V, 0) and (V ) − = max(−V, 0). In this scheme, both integrals involved in the computation ofr n and the V n j+1/2,j∈{0,...,J} are approximated via a centered rectangle formula (thus considering that the numerical approximation of a is piecewise constant in radius). We here decide to smooth out the velocity field by replacing the discontinuous one with the following Lipschitz-continuous one where > 0 is a small parameter, which will be equal to = 0.005r c in the simulations. Note that, in the case of a discontinuous velocity, such as (6), we would have used instead an upwind typed scheme, that would write (see [5,6,9]), a n+1 j = a n j − ∆t ∆r (V n j ) + a n j − (V n j+1 ) − a n j+1 − (V n j−1 ) + a n j−1 + (V n j ) − a n j − ∆tγ a n j , j ∈ {1, ..., J}.
We refer the reader to [9] for further details on the pros and cons of several schemes for transport with discontinuous velocities. Stability analysis: That the numerical unknowns remain non negative, as required by the modeling, imposes constraints on the time step.
Let us consider first the ODE part of the system. We are interested in situations where the parameters satisfy Therefore, we expect the mesenchymal cells to decrease in time as well as the preadipocytes. For the mesenchymal cells, the discrete relation casts as Thus m n+1 remains non negative (if m n is) as far as 1 + ∆t(α − γ − β(r n )) ≥ 0. A sufficient condition to satisfy this is that ∆t|α − γ − β(r)| ≤ 1 for any r, thus the condition we keep is In the same way we derive the condition for the preadipocytes equation. In the worst situation, m vanishes and the population of preadipocytes decreases with respect to time. Thus, p n+1 remains non negative provided Now let us consider the PDE in size space for the adipocytes. The upwind scheme is stable (in the linear case where the velocity field is given, see [5]) if the following Courant-Friedrichs-Lewy condition is satisfied: Therefore, at each iteration the time step is updated in order to fulfill all three conditions identified so far.

Simulations
The parameters used in the simulations are given in Table 2, according to unpublished experimental data from C3M. A fit was performed with these data in order to calibrate the parameters; a larger number of available data will be necessary to improve this calibration. In this initial configuration, the adipocytes have a small radius. We study the behavior of the solution in a large time scale, namely at the scale of a year (t ∈ (0, 350) days).   Table 2 and initial datum in Eq. (30), namely a Gaussian distribution of adipocytes with small radii. The mesenchymal cells differentiate into preadipocytes, which in turn differentiate into mature adipocytes. For large times, the number of mesenchymal cells and of preadipocytes go to zero, while the total number of adipocytes stabilizes to a stationary value and the mean radius of adipocytes stabilizes towards the critical value r c . Note that the mean radius grows quickly at the beginning, then decreases when adipocytes with a small radius are created from preadipocytes, and finally increases again. Figure 3 shows the solution to the system of equations. Since we imposed that the proliferation rate is smaller than the sum of the differentiation and of the death rates, either the mesenchymal cells or the preadipocytes tend to zero for large time. At the beginning (first 40 days) the number of preadipocytes increases. This is due to the differentiation of the mesenchymal cells into preadipocytes. The number of adipocytes increases in time.
Since we have set a very small death rate for the adipocytes (γ ∼ 10 −6 ), the population reaches a plateau and the death of the adipocytes is not sensitive on this time scale of observation. We observe that the mean radius of the adipocytes grows rapidly at the beginning (first 10 days). Then, there is a drop due to the differentiation of preadipocytes into adipocytes: indeed, this process produces a large number of new adipocytes which have a radius equal to r . Next, these small adipocytes grow until their radius reaches the critical value r c . The adipocytes growth can equally be observed in Figure 4. Initially, the adipocytes are small and they are subjected to a strong growth rate, since the growth rate is proportional to the inverse of the total surface of adipocytes. When the cells grow, the rate decreases because of an increasing of the surface of adipocytes. Figure  5 shows the distribution of adipocytes at the last time step. We observe the accumulation of the adipocytes population to the critical radius and the formation of the Dirac mass, according to the theoretical results.  . Adipocytes distribution (left) and growth rate (right) with respect to the radius at three different time steps (t = 0 day, t = 1.263 day and t = 3.103 days). Parameters are the same as in Fig.3. We can notice that there is indeed no adipocytes with a radius larger than the critical value r c . Initially, the growth rate is very high as the total surface is small; it then decreases when the total surface increases.   Fig.3. Adipocytes accumulate to the critical radius r c and form a Dirac mass at this value. The growth rate is quite small compared to its initial value, see Fig. 4, and null for radii larger than the critical radius r c .

Differentiation rates
We now observe the behavior of the proliferations rates β(r(t)) and β (r(t)). Figure 6 shows the dependency of the two functions with respect to the radius and to time (with parameters as given in Table 2). Since the mean radius increases rapidly, also the differentiation rates reach rapidly the maximum value. (r(t)) (r(t)) Figure 6. Differentiation rates of the mesenchymal cells (blue continuous line) and the preadipocytes (orange dotted line) with respect to the radius (left) and to time (right). Parameters are the same as in Fig.3. The differentiation rates are stiff and reach their maximum values for small radii. Since the mean radius increases quickly initially, the differentiation rates reach rapidly their maximum values.
The qualitative behavior does not change significantly when working with constant differentiation rates (see Figure 7). It likely means that the quantities of interest remain in the low or high regions of the sigmoid in this example.   Fig.3. We can notice that there are no differences between the two cases.
However, changing the parameters of the the sigmoid functions β and β can lead to different behaviors of the solutions. We consider three cases, detailed in Table 3: we fix the slope of the sigmoid functions and we make the inflection points vary. Figure 8 shows the comparison of the populations dynamics in the different cases. The larger the inflection point of β , the more preadipocytes are generated, and therefore differentiate. In fact, since the differentiation rate is slower, they have time to duplicate. However, the dynamics of the adipocytes and of the mean radius in Figure 8 suggest that the parameters considered in Table 2 are in good agreement with the experimental data. r ip,β r ip,β case 1 0.5r c 0.5r c case 2 0.5r c r c case 3 0.5r c 0.7r c Table 3. Inflection points relative to β and β in the different cases.  Tables 2 and 3 and initial datum in Eq. (30). The blue line refers to the parameters set in Table 2 (which correspond to a fit with the experimental data), the red dotted line to the case where the inflexion points for the differentiation rates are equal to r ip,β = r ip,β = 0.5r c , the dashed green line to r ip,β = 0.5r c , r ip,β = r c and the orange line to r ip,β = 0.5r c , r ip,β = 0.7r c . The black circles correspond to the experimental data. Behaviors of the solutions are significantly different: the larger the inflection point of β , the more preadipocytes are created, and therefore the more adipocytes are generated, leading to a smaller value of the mean radius in large times. In the critical case when r ip,β = r c , the preadipocyte number and the adipocyte number tend to infinity asymptotically, while the mean radius decreases.

Comparison with experimental data
We now compare the models with the experimental data. We test the model structured in size (1) (referred to as 0D), and the models structured in both size and space, which will be detailed in the next sections: we refer to these models as 1D and 2D, depending on the space dimension used in the simulations. We observe in Figure 9 that the solutions to the 0D, the 1D and the 2D models have similar behavior. Since the 0D model has been calibrated with the data, it is in good agreement with the experimental data.  Table 2 and initial datum in Eq. (30), whereas initial datum for the 1D model is given in Eq. (38), see also Fig.11. Initial datum for the 2D model is given by an equivalent formula. We observe that the three models show comparable tendency, but with different values. The 0D model has been calibrated with the experimental data.

5.
Extension to a spatial model

Reinterpretation of the unknowns
We are now going to extend the previous model by incorporating a spatial dependence intended to describe inhomogeneities and displacement of the adipose tissue, following some ideas of mixture theory, see [1,3,8]. The functions (t, x) → m(t, x), (t, x) → p(t, x) now depend also on the space variable x ∈ Ω ⊂ R n (n = 1 or 2, and Ω is a bounded domain), and they are defined on [0, ∞) × Ω with values in [0, ∞). For the adipocytes, (t, x, r) → a(t, x, r) is a function defined on [0, ∞) × Ω × [r , +∞) with values in [0, ∞). We also introduce a new function (t, x) → s(t, x) which is the volume fraction of the surrounding material (tissue) that is carrying the cells. The motion of all these species is driven by a velocity field (t, x) → u(t, x).
In order to derive the model, let us go back to the mass balance relations. Given r 2 > r 1 ≥ r , the integral a(t, x, r)r 3 dr dx, gives the volume occupied in O at time t by the adipocytes with a radius r ∈ (r 1 , r 2 ). Accordingly, 4π 3 ∞ r a(t, x, r)r 3 dr defines the volume fraction of the adipocytes, and the total volume of functional adipocytes at time t is given by Schematically, we can split the volume of an adipocyte with radius r ≥ r into two parts, see Fig. 10: the center, with radius r , is mainly made of water (volumetric mass density ρ w ), like the pre-adipocytes and mesenchymal cells, while the outer domain is made of lipids (volumetric mass density ρ l ). Therefore, the mass of such an adipocyte reads 4π 3 ρ l (r 3 − r 3 ) + ρ w r 3 = 4π 3 ρ a (r)r 3 r water ρ w r > r lipids ρ l Figure 10. Geometrical description of a growing adipocyte.
Accordingly, the mass density of the adipocytes is given by x, r)r 3 dr, and the total mass of adipocytes reads Ω M a (t, x) dx. The evolution of the mass density obeys γ ρ a (r)a(t, x, r)r 3 dr + 4πr 3 3 ρ w β (r(t, x))p(t, x) + 4π ∞ r ρ l V a(t, x, r)r 2 dr.
The last two terms describe two mechanisms of gain of mass: the transformation of pre-adipocytes into adipocytes with radius r and volumetric mass density ρ w and the input of lipids from the surface of the existing adipocytes. Similarly, the volume fractions of mesenchymal cells and pre-adipocytes are given by • 4π 3 r 3 m m(t, x) for the mesenchymal cells, • 4π 3 r 3 p(t, x) for the preadipocytes, respectively, with r m and r the typical radius of these cells. According to the biological data, we will use from now on that r m = r .
Given O ⊂ Ω, the integrals O The mass of each species is transported by a velocity field (t, x) ∈ [0, ∞) × Ω → u(t, x) ∈ R n , which amounts to say For the adipocytes, we get ∂ t (ρ a r 3 a) + ∇ x · (ρ a r 3 ua) + ρ a r 3 ∂ r (V a) = −γ ρ a r 3 a which can also be written in the conservative form In fact, we can get rid of the densities in all these equations, and (32) can be simplified as: Beyond the transport by the velocity u, the modeling also uses a pressure field, hereafter denoted (t, x) → q(t, x); we shall assume that the parameters β and β are now functions of bothr and this quantity q. It incorporates another source of space inhomogeneities. We assume that β, β are non decreasing with respect to q. The pressure acts as a mechanical constraint that limits the expansion of the adipocytes: the weaker the pressure, the easier the transformation of mesenchymal cells and pre-adipocytes to pre-adipocytes and adipocytes, respectively.
Let us now discuss the equations for the pair (u, q). The following constraint on the volume fractions holds, (using that the radius of pre-adipocytes and of mesenchymal cells coincides with the radius r of the smallest adipocytes) for a. e. x ∈ Ω, t ≥ 0. Therefore, adding the equations in (31) and (33), we find Using the same boundary condition as (3), we obtain Again, we observe that the last term can be integrated by parts to make two contributions appear, the former from the passage of pre-adipocytes to adipocytes with size r , the latter due to surface fluxes of lipids. This constraint is related to the pressure q, which may be seen as the corresponding Lagrange multiplier. Next, we use Darcy's equation: which have different frames of application [4,11]. In (36), x → K(x) takes values in the set of symmetric positive matrices. With (36) we are directly led to an equation for q, by using equation (35) We complement the whole system with initial conditions: satisfying constraint (34). Finally, we need to impose boundary conditions. For r = r , we make use of the same birth condition as in the homogeneous case V (t, r ) a(t, x, r ) = β (r(t, x), q(t, x))) p(t, x). On ∂Ω, we shall find boundary conditions. For example, we can use • for cells, the no-incoming-flux condition: if u · n < 0, then a = p = m = 0 and s = 1 where n stands for the unit outward normal vector on ∂Ω, we also set q = 0 on ∂Ω.
• or the wall conditions. These boundary conditions have to be compatible with condition (35).

Description of the numerical discretization
In order to discretize the system of equations, we apply the finite volume method. First, we find the spatial velocity u solving the Darcy's law. We formulate the Darcy's equation as a Laplace problem: From (37) we compute the pressure q. Then, we find the velocity u thanks to in Ω.
Applying the finite volume method, we obtain the following discretization on each element K i of the domain.
Finally, we apply an explicit upwind scheme in order to solve the system of hyperbolic equations. Regarding the adipocytes, we use a "directional" splitting and solve the transport equation either in space or with respect to the radius dimension applying the finite volume upwind scheme.

Numerical simulations for the spatial case
In one dimension, we consider an interval (0, L) where L = 1 cm, therefore x ∈ Ω = (0, 0.01). We discretize the interval in J = 300 subintervals with a homogeneous spatial step ∆x. Moreover, we have to consider the "radius" dimension of the adipocytes. To approximate their growth, we consider the radius of the adipocytes r to be in the interval I r = (r , R), and we set R = 90µm. We divide I r into 100 homogeneous intervals of size ∆r. Then, we solve the problem that has one dimension with respect to the mesenchymal cells m and the preadipocytes p and two dimensions with respect to the adipocytes a.

Parameters setting and initial conditions
The parameters used for the simulations are the ones introduced in Table 2. Regarding the initial condition, we consider the three populations to be distributed as gaussians centered in different points of Ω with various standard deviations (with respect to space and radius), namely For the numerical simulations, we take the following values for the standard variations : θ = 0.1, θ r = 0.1, θ a = 0.8. Figure 11. Initial conditions for m (left), p (center ) and a (right) relative to the 1D spatial model described at Section 5.1. The figure on the right represents the distribution of a in each point of the space x and with a radius r. Figure 11 shows the initial conditions of the three unknowns. At the beginning we consider the highest number of adipocytes to have the minimal radius r . However, since the growth rate is proportional to the inverse of the surface of adipocytes in each spatial point, we consider minimum number of cellsā 0 in each spatial point.

Simulations
We study the behavior of the populations and of the mean radius first looking at a small time scale. In particular, we focus on the growth of the initial adipocytes. Figure 12 shows the adipocyte dynamic in the first day. We observe a slow growth of cells at the center of the initial gaussian, where the total surface of the cells is the highest. If we move from this point, the total surface of the cells decreases and consequently the radius of the cells increases faster. In Figure 13 we observe the global growth of adipocytes at different time steps. Like with the space homogeneous equations, we observe a high growth rate at the beginning, leading to a fast growth of the adipocytes at the initial configuration.
Since the spatial velocity is very small, we do not observe a displacement of the different cells in the interval. The spatial dependency indeed introduces two different timescales relative to the growth of adipocytes and to the spatial displacement.
We observe on numerical simulations, see Fig.13, that dealing with the space dependent framework favors the apparition of bimodal size-distributions, that have been recorded in experiments, see [17,18]. We were unable to reproduce such bimodal shape with the space homogeneous model. It is likely that such distributions correspond only to transient states of the model; nevertheless they can be relevant on the time scale of observation. Going further in this direction requires a better knowledge of the parameters of the model and deserves a thorough investigation.

Conclusion
This work aims at setting up the basis of a multi-species model for adipose tissue growth, with the aim of improving our understanding of obesity. Our model is based on detailed behaviors and interactions of adipose tissue cells. It accounts for the adipocyte maturation, starting from mesenchymal cells differentiating into pre-adipocytes, which in turn differentiate into adipocytes whose sizes grow with food supply. The core of our model is the coupled regulation of cell differentiation and proliferation using phenomenological laws with sigmoid shapes that depend on the mean size of adipocytes.
The model in its non-spatial version is able to reproduce with reasonable precision experimental data after a fitting process, indicating that the main biological phenomena are probably accounted for in our model hypotheses. The time dynamics of the adipocyte radii exhibits an interesting early overshoot that cannot be seen in the data because of its sampling. The overshoot is easily explained by the model dynamics and its existence could be investigated in future experimental work. A shift in early adipocyte numbers is also observed, probably because of the difficult choice of an initial condition in the absence of a more detailed set of data.
The spatial model also keeps a reasonable agreement with the data, at least in order of magnitude and allows to exhibit a well known behavior for adipocyte distribution: a bimodal distribution of adipocytes in space [17,18]. This last result can be seen as a first step for validation of our model and makes it very promising. The fact that we do not observe spatial displacements of the cells needs more investigations and we can consider the mechanical properties of adipose tissues described in [13] to improve the spatial description in our model. The next step is now to gather more complete data sets from experimental studies and to reach a proper validation. Once validated, such a model could prove invaluable to identify and quantify the biological mechanisms involved in adipose tissue growth and to understand potential dysfunctions linked to obesity.   Table 2. Growth of cells is heterogeneous in space: it is slow at the center of the initial gaussian, where the total surface of the cells is the highest and faster away from the center. Spatial displacement of the cells is small, since the spatial velocity is very small.  . Adipocytes distribution with respect to the radius at different time steps (t = 0.1 day, t = 0.2 day, t = 0.5 day, and t = 1 day). Parameters are the same as in Fig. 12. Adipocytes grow with respect to time, quickly at early times. We observe in that case the apparition of bimodal size-distributions, with a peak at the crtitical value r c .