O Angulo

F Crauste
email: crauste@math.univ-lyon1.fr

J C López-Marcos

Numerical integration of an erythropoiesis model with explicit growth factor dynamics

Keywords: erythropoiesis model, nonlinear age-structured system, numerical scheme, long-time integration 2000 MSC: 65M25, 35B40, 35L60, 35Q92, 92D25

Erythropoiesis, the red blood cell production process, involves interactions between cell populations with different differentiation states, mainly immature progenitor cells and mature erythrocytes, and growth factors such as erythropoietin and glucocorticoids, known to respectively inhibit cell apoptosis, stimulate proliferation and differentiation, and stimulate self-renewal. The feedback regulation of this process allows a very fast and efficient recovery in the case of a severe anemia. We consider an age-structured model of red blood cell production accounting for these feedback regulations and the dynamics of growth factors. We theoretically show the existence of a unique positive steady state for the model and we propose a numerical method to obtain an approximation to its solution. Experiments are reported to show numerically, on one hand, the optimal convergence order of the numerical scheme and, on the other hand, a fine approximation to real experimental data, with a suitable selection of the parameters involved.

Introduction

Production of red blood cells occurs in the bone marrow, where undifferentiated immature hematopoietic stem cells differentiate, throughout a number of divisions, into more and more mature cells until fully differentiated mature red blood cells are produced. The process of production and regulation of red blood cells is called erythropoiesis. The first stage of differentiation of hematopoietic stem cells is the progenitor stage: erythroid progenitors are undifferentiated cells committed to the red blood cell -or erythroid -lineage. These cells share a lot of properties with stem cells, such as the ability to proliferate and differentiate but also to self-renew, that is to produce by division two daughter cells identical to the mother cell. Erythroid progenitors proliferate and differentiate in more mature cells, meanwhile losing their self-renewal ability, until they reach the stage of reticulocytes. These latter are immature red blood cells that will finish their differentiation process within few divisions, ejecting their nucleus, and finally differentiating in erythrocytes that leave the bone marrow and enter the bloodstream as red blood cells. The entire process is regulated by a variety of feedback controls, involving growth factors. Among them, erythropoietin has been identified in 1990 as the main regulator of erythropoiesis: released by the kidneys when the organism detects a lack of red blood cells, erythropoietin binds to the surface of erythroid progenitors and prevents their death, called apoptosis, hence increasing the amount of differentiated cells [START_REF] Koury | Erythropoietin retards DNA breakdown and prevents programmed death in erythroid progenitor cells[END_REF]. Glucocorticoids have also been shown to influence erythropoiesis in stress conditions by inducing erythroid progenitor self-renewal [START_REF] Bauer | The glucocorticoid receptor is required for stress erythropoiesis[END_REF][START_REF] Gandrillon | TGF-beta cooperates with TGF-alpha to induce the self-renewal of normal erythrocytic progenitors: evidence for an autocrine mechanism[END_REF].

We propose a mathematical model of stress erythropoiesis in mice based on age-structured nonlinear partial differential equations describing erythroid progenitor and erythrocyte dynamics in the next section. This model is based on previous models, proposed in [START_REF] Angulo | Investigating the Roles of the Experimental Protocol in Phenylhidrazine-Induced Anemia in Mice[END_REF][START_REF] Crauste | Adding self-renewal in committed erythroid progenitors improves the biological relevance of a mathematical model of erythropoiesis[END_REF], that were inspired by the initial works by Mackey [START_REF] Mackey | Unified hypothesis of the origin of aplastic anaemia and periodic hematopoïesis[END_REF] and co-authors, and Loeffler et al. [START_REF] Loeffler | A mathematical model of erythropoiesis in mice and rats. Part 1. Structure of the model[END_REF][START_REF] Roeder | A novel dynamic model of hematopoietic stem cell organization based on the concept of within-tissue plasticity[END_REF], to describe stress erythropoiesis in mice. Several mathematical models of erythropoiesis have been proposed over the last 30 years (see [START_REF] Roeder | Quantitative stem cell biology: computational studies in the hematopoietic system[END_REF] for a review), in order to address the mechanisms of regulation and their role in stress or pathological conditions. Mackey Mackey [START_REF] Mackey | Unified hypothesis of the origin of aplastic anaemia and periodic hematopoïesis[END_REF] published a pioneering paper on hematopoietic stem cell dynamics, which has then been used and modified by many authors, including Mackey and co-workers, to describe differentiation and maturation processes involved in hematopoiesis more precisely [START_REF] Adimy | On the stability of a maturity structured model of cellular proliferation[END_REF][START_REF] Adimy | A mathematical study of the hematopoiesis process with applications to chronic myelogenous leukemia[END_REF][START_REF] Colijn | A mathematical model of hematopoiesis -I. Periodic chronic myelogenous leukemia[END_REF][START_REF] Colijn | A mathematical model of hematopoiesis -II. Cyclical neutropenia[END_REF][START_REF] Pujo-Menjouet | Long period oscillations in a G 0 model of hematopoietic stem cells[END_REF]. Belair et al [START_REF] Bélair | Age-structured and two-delay models for erythropoiesis[END_REF] proposed a model of erythropoiesis considering the influence of growth factors on stem cell differentiation in erythroid progenitors, later improved by Mahaffy et al [START_REF] Mahaffy | Hematopoietic model with moving boundary condition and state dependent delay: applications in erythropoiesis[END_REF], and recently analyzed in more details [START_REF] Ackleh | A structured erythropoiesis model with nonlinear cell maturation velocity and hormone decay rate[END_REF][START_REF] Ackleh | Parameter Estimation in a Structured Erythropoiesis Model[END_REF]. Another erythropoiesis model, inspired by the same article, was introduced by Adimy et al. [START_REF] Adimy | Modelling hematopoiesis mediated by growth factors with applications to periodic hematological diseases[END_REF] in 2006, in which Epo is the only growth factor supposed to act during erythropoiesis. An important contribution to mathematical modeling of erythropoiesis is also due to Loeffler and his collaborators [START_REF] Loeffler | A mathematical model of erythropoiesis in mice and rats. Part 1. Structure of the model[END_REF][START_REF] Roeder | Quantitative stem cell biology: computational studies in the hematopoietic system[END_REF][START_REF] Wichman | Mathematical Modeling of Cell Proliferation[END_REF]. Their models consider feedback controls from progenitors on the stem cell level and from mature cells on progenitors, and were fitted to various experiments (including irradiations, bleeding, and phenylhydrazine treatments of mice). Most of these works were nevertheless performed before the role of Epo was definitely identified and long before erythroid progenitor self-renewal was hypothesized.

The equations we propose are coupled via age-dependent feedback functions associated with cell ability to die by apoptosis, to self-renew and to differentiate. Contrary to previous, simpler models on which it is based [START_REF] Angulo | Investigating the Roles of the Experimental Protocol in Phenylhidrazine-Induced Anemia in Mice[END_REF][START_REF] Crauste | Adding self-renewal in committed erythroid progenitors improves the biological relevance of a mathematical model of erythropoiesis[END_REF], that partially described regulation of erythropoiesis (for instance, none of them described agedependency of feedback functions, and only the model studied in [START_REF] Angulo | Investigating the Roles of the Experimental Protocol in Phenylhidrazine-Induced Anemia in Mice[END_REF] assumed a non-constant mortality rate of erythrocytes), our model additionally considers explicit growth factor dynamics, unlike all previously published models of erythropoiesis (except for Adimy et al. [START_REF] Adimy | Modelling hematopoiesis mediated by growth factors with applications to periodic hematological diseases[END_REF]). We first establish the theoretical existence of a unique positive stationary solution in section 3 then propose a new numerical scheme to simulate solutions of the system in section 4. In section 5 we illustrate the efficiency of our method by showing numerically that it is of second order of convergence, and that it allows to compare simulated solutions of the model with real biological data.

Model of Stress Erythropoiesis

Let consider a population of erythroid progenitors and a population of erythrocytes, defined by cell age a and the time of the observation t (see Figure 1). Cell age a is considered to be the time spent from the birth of a cell up to the time of observation, hence cells age with unitary velocity (see ( 4)). Among progenitors, a difference is made between self-renewing progenitors, whose population is denoted by s(t, a), and differentiating progenitors that are not self-renewing, denoted by p(t, a). The duration of the differentiating progenitor compartment is denoted by τ p , and the duration of one self-renewing cycle by τ c . One may note that τ c < τ p . All progenitors are supposed to be localized in the bone marrow, and all of them can die by apoptosis. The number of erythrocytes, circulating in blood, is denoted by e(t, a), and E(t) denotes the total number of erythrocytes at time t, defined by

E(t) = +∞ 0 e(t, a)da.
Two growth factors, erythropoietin (Epo) and glucocorticoids (GC), are assumed to be the main regulators of erythroid progenitor apoptosis and differentiation, and self-renewal respectively [START_REF] Crauste | Adding self-renewal in committed erythroid progenitors improves the biological relevance of a mathematical model of erythropoiesis[END_REF]. Denote by Epo(t) the concentration of Epo at time t in bloodstream, and by GC(t) the concentration of glucocorticoids at time t. These concentrations are regulated by the total number of erythrocytes E(t): the more erythrocytes, the less Epo and glucocorticoids [START_REF] Bauer | The glucocorticoid receptor is required for stress erythropoiesis[END_REF][START_REF] Koury | Erythropoietin retards DNA breakdown and prevents programmed death in erythroid progenitor cells[END_REF].

Denote by α(Epo) the erythroid progenitor apoptosis rate, and by σ(a, GC) the erythroid progenitor self-renewal rate. The function α is assumed to be decreasing with respect to Epo (negative feedback, [START_REF] Koury | Erythropoietin retards DNA breakdown and prevents programmed death in erythroid progenitor cells[END_REF]). The function σ is assumed to be increasing with respect to GC (positive feedback), and decreasing with respect to the age a of progenitor cells [START_REF] Bauer | The glucocorticoid receptor is required for stress erythropoiesis[END_REF] with σ(τ p , GC) = 0. This latter assumption illustrates a lack of sensitivity to GC of mature progenitors, contrary to early progenitors that are more sensitive to GC. These assumptions are summarized hereafter, with ∂ x denoting the partial derivative with respect to x:

α (Epo) ≤ 0, ∂ a σ(a, GC) ≤ 0, ∂ GC σ(a, GC) ≥ 0, σ(τ p , GC) = 0.
Denote by δ(a, Epo) the differentiation rate of progenitors into erythrocytes. This rate depends on the age of progenitors, older ones being more inclined to differentiate [START_REF] Hattangadi | Homeodomain-interacting protein kinase 2 plays an important role in normal terminal erythroid differentiation[END_REF], and on the concentration of Epo, the more Epo the more differentiation [START_REF] Hattangadi | Homeodomain-interacting protein kinase 2 plays an important role in normal terminal erythroid differentiation[END_REF]. The function δ is hence assumed to be increasing with respect to a and Epo, ∂ a δ(a, Epo) > 0, ∂ Epo δ(a, Epo) > 0.

In addition, since all progenitors which did not die in the differentiating compartment must differentiate in erythrocytes, we assume lim a→τp δ(a, Epo) = +∞ whatever the level of Epo.

The following functions α, σ and δ, will be used in section 5 to perform numerical experiments,

α(Epo) = C α θ nα α θ nα α + (log 10 (Epo)) nα , (1)
where C α > 0 is the maximal apoptosis rate, θ α > 0 is a threshold value, and n α > 1 is a sensitivity parameter;

σ(a, GC) = C σ (τ p -a) (log 10 (GC)) nσ θ nσ σ + (log 10 (GC)) nσ , (2)
where, similarly to the function α, θ σ > 0, n σ > 1, and C σ > 0 is a constant characterizing the maximum self-renewal rate given, for a fixed value of GC, by

C σ (τ p -a); δ(a, Epo) = C δ1 τ p τ p -a + C δ2 (log 10 (Epo)) n δ θ n δ δ + (log 10 (Epo)) n δ =: δ 1 (a) + δ 2 (Epo), (3)
where C δ1 is the minimum of the age-dependent differentiation rate δ 1 (a), reached when a = 0, C δ2 is the maximum of the Epo-dependent differentiation rate δ 2 (Epo), and C δ1 > 0, C δ2 > 0, θ δ > 0 is a threshold value for the Epo-dependent part of the differentiation rate, n δ > 1 a sensitivity parameter.

It must be noted that in the absence of GC, the self-renewal rate is supposed to vanish, whatever the age of cells, whereas differentiation occurs in the presence or the absence of Epo: Epo is assumed to positively act on erythroid progenitor differentiation, yet differentiation can be induced independently of Epo levels. This explains why the differentiation rate is composed of two additive terms, whereas the self-renewal rate is composed of only one term.

Denote by γ the mortality rate of erythrocytes. To our knowledge, no precise measurement of this rate is available in the literature, so it is usually supposed to be constant. It can however be modified under specific circumstances (see [START_REF] Crauste | Adding self-renewal in committed erythroid progenitors improves the biological relevance of a mathematical model of erythropoiesis[END_REF] and section 5 for instance), so without loss of generality we will consider that it depends on both erythrocyte age a and observation time t, and we will write γ(t, a) (see [START_REF] Angulo | Investigating the Roles of the Experimental Protocol in Phenylhidrazine-Induced Anemia in Mice[END_REF]).

Then, the quantities p, s and e satisfy the following system, for t > 0,

           ∂ t p(t, a) + ∂ a p(t, a)=-[α(Epo(t)) + δ(a, Epo(t))] p(t, a)
-σ(a, GC(t))p(t, a), a ∈ (0, τ p ),

∂ t s(t, a) + ∂ a s(t, a)=-α(Epo(t))s(t, a), a ∈ (0, τ c ),
∂ t e(t, a) + ∂ a e(t, a)=-γ(t, a)e(t, a), a > 0.

(

) 4
Boundary conditions associated with (4) describe cell flux between compartments,

             p(t, 0) = HSC + 2s(t, τ c ), s(t, 0) = τp 0 σ(a, GC(t))p(t, a)da, e(t, 0) = A τp 0 δ(a, Epo(t))p(t, a)da. (5)
New erythroid progenitors come both from the division of self-renewing progenitors in two newborn cells and the differentiation of hematopoietic stem cells: HSC denotes the flux of hematopoietic stem cells differentiating in erythroid progenitors, assumed to be constant over time. Self-renewing progenitors are produced at a rate σ from differentiating progenitors. Erythrocytes are produced from differentiating progenitors, and A denotes a constant amplification coefficient accounting for divisions of mature progenitors. For instance, A = 2 n where n is the number of differentiation stages during the reticulocyte stage [START_REF] Crauste | Adding self-renewal in committed erythroid progenitors improves the biological relevance of a mathematical model of erythropoiesis[END_REF]. Concentrations of growth factors Epo(t) and GC(t) satisfy the following ordinary differential equations [START_REF] Crauste | Mathematical study of feedback control roles and relevance in stress erythropoiesis[END_REF],

     d dt Epo(t) = f (E(t)) -k Epo Epo(t), d dt GC(t) = g(E(t)) -k GC GC(t), (6)
where k Epo > 0 and k GC > 0 are degradation rates, and the functions f and g describe feedback controls by the total number of erythrocytes on the production of growth factors. Both functions are assumed to be decreasing (see Figure 1),

f (E) < 0, g (E) < 0.
The model is completed with suitable initial conditions

                 p(0, a) = p 0 (a), a ∈ [0, τ p ], s(0, a) = s 0 (a), a ∈ [0, τ c ],
e(0, a) = e 0 (a), a > 0,

Epo(0) = Epo 0 , GC(0) = GC 0 . (7)
Existence and uniqueness of solutions of system ( 4)-( 7) follow from the classical theory of age-structured equations [START_REF] Webb | Theory of nonlinear age-dependent population dynamics, Monographs and textbooks in Pure[END_REF], under conditions of smoothness of the various nonlinear feedback functions. In section 3, we analyze the existence of stationary solutions for this model.

The effort made on increasing the realism in the model is however achieved at the expense of loss in mathematical tractability. We point out that, without other restrictive assumptions, this model cannot be solved analytically. Therefore, the use of efficient methods that provide a numerical approach is the most suitable mathematical tool for studying the problem and, indeed, it is often the only one available. Besides, numerical methods have been successfully applied to structured models to replicate available field and/or laboratory data for a variety of different systems (e.g. [START_REF] Angulo | Mass Structured Systems with Boundary Delay: Oscillations and the Effect of Selective Predation[END_REF][START_REF] Angulo | Numerical approximation of singular asymptotic states for a size-structured population model with a dynamical resource[END_REF][START_REF] Angulo | A semi-Lagrangian method for a cell population model in a dynamical environment[END_REF][START_REF] Angulo | Numerical analysis of a population model of marine invertebrates with different life stages[END_REF] and references therein). We completely describe, in section 4, the explicit second order numerical scheme, built "ad hoc", which has been used to obtain the solutions of ( 4)- [START_REF] Adimy | Numerical integration of a mathematical model of hematopoietic stem cell dynamics[END_REF]. And, finally, we perform a wide numerical experimentation which includes convergence, long-time integration with the simulation of the theoretical steady state and a comparison with real biological data in section 5.

Existence of a Positive Steady State

We focus in this section on the existence of steady states for system (4)-( 6) when γ does not depend on time t, that is in normal erythropoiesis. From a biological point of view, system (4)-( 6) is expected to possess in this case only one positive steady state, due to the constant positive flux HSC from the hematopoietic stem cell compartment.

Suppose (p * (a), s * (a), e * (a), Epo * , GC * ) is a steady state solution of system (4)-( 6). Then, it satisfies

             dp * da (a) = -[α(Epo * ) + σ(a, GC * ) + δ(a, Epo * )] p * (a), a ∈ (0, τ p ), ds * da (a) = -α(Epo * )s * (a), a ∈ (0, τ c ), de * da (a) = -γ(a)e * (a), a > 0, (8)
with

             p * (0) = HSC + 2s * (τ c ), s * (0) = τp 0 σ(a, GC * )p * (a)da, e * (0) = A τp 0 δ(a, Epo * )p * (a)da, (9)
and

k Epo Epo * = f (E * ), k GC GC * = g(E * ),
where

E * = +∞ 0 e * (a)da. ( 10
)
Let us introduce the following notations,

           ᾱ(E * ) := α(F (E * )), σ(a, E * ) := σ(a, G(E * )), a ∈ [0, τ p ], δ(a, E * ) := δ(a, F (E * )), a ∈ [0, τ p ], β(a, E * ) := ᾱ(E * ) + σ(a, E * ) + δ(a, E * ), a ∈ [0, τ p ],
where F (E * ) := f (E * )/k Epo and G(E * ) := g(E * )/k GC . Consider the following assumptions, (H1) Suppose for all θ ∈ [0, τ p ],

∂ Y β(θ, Y ) > 0, ∀Y ≥ 0.
(H2) Suppose, for all θ ∈ [0, τ p ], there exists Ȳ > 0 such that

∂ Y β(θ, Y ) < 0, ∀Y ∈ [0, Ȳ ), ∂ Y β(θ, Y ) > 0, ∀Y ∈ ( Ȳ , +∞).
From a biological point of view, Assumption (H1) describes a situation in which as the number of erythrocytes increases progenitor apoptosis increases. This is expected, otherwise the erythrocyte population could be unbounded. Assumption (H2) has a similar interpretation, but it also considers that for small erythrocytes counts, progenitor apoptosis rate becomes negligible compared to self-renewal and differentiation rates. This should theoretically allow a repopulation of the erythrocyte population. Both assumptions are in agreement with biological knowledge. The next proposition states the existence of a unique positive steady state.

Proposition 1. Assume inequality e -α(0)τc τp 0 σ(a, 0) exp - a 0 α(0) + σ(θ, 0) + δ(θ, 0)dθ da < 1 2 (11)
holds true. Then whether Assumption (H1) or (H2) is satisfied, system (4)-( 6) has a unique steady state.

Proof. This proof is in two parts. We first show that system (8)-( 10) has a solution if and only if condition (21) holds true. Then, as a second step, we show that under Assumption (H1) or (H2) condition ( 21) holds true as soon as [START_REF] Angulo | Numerical approximation of singular asymptotic states for a size-structured population model with a dynamical resource[END_REF] is satisfied.

Step 1. Solving system (8)- [START_REF] Angulo | Mass Structured Systems with Boundary Delay: Oscillations and the Effect of Selective Predation[END_REF] gives

p * (a) = exp - a 0 β(θ, E * )dθ [HSC + 2s * (τ c )] , a ∈ [0, τ p ], (12)
s * (a) = exp (-ᾱ(E * )a) τp 0 σ(θ, E * )p * (θ)dθ, a ∈ [0, τ c ], (13)
e * (a) = exp -

a 0 γ(θ)dθ A τp 0 δ(θ, E * )p * (θ)dθ, a > 0. ( 14
)
Using ( 13) in ( 12), we deduce

p * (a) = exp - a 0 β(θ, E * )dθ Ψ(E * , p * ), a ∈ [0, τ p ], (15)
where

Ψ(E * , p * ) := HSC + 2e -ᾱ(E * )τc τp 0 σ(θ, E * )p * (θ)dθ. ( 16
)
Using now [START_REF] Bauer | The glucocorticoid receptor is required for stress erythropoiesis[END_REF] in [START_REF] Angulo | A numerical method for nonlinear age-structured population models with finite maximum age[END_REF] one obtains

e * (a) = exp - a 0 γ(θ)dθ AΛ(E * )Ψ(E * , p * ), a > 0, with Λ(E * ) := τp 0 δ(θ, E * ) exp - θ 0 β(u, E * )du dθ. (17)
It follows that

E * := +∞ 0 e * (a)da = ξAΛ(E * )Ψ(E * , p * ), (18)
where

ξ := +∞ 0 exp - a 0 γ(θ)dθ da > 0.
From [START_REF] Bauer | The glucocorticoid receptor is required for stress erythropoiesis[END_REF], one gets for a ∈ [0,

τ p ] σ(a, E * )p * (a) = σ(a, E * ) exp - a 0 β(θ, E * )dθ Ψ(E * , p * ). (19)
Introducing the following notations,

     µ(E * ) := e -ᾱ(E * )τc , ν(E * ) := τp 0 σ(a, E * ) exp - a 0 β(θ, E * )dθ da, (20)
and using [START_REF] Bélair | Age-structured and two-delay models for erythropoiesis[END_REF], we finally get from the integration of ( 19)

[1 -2µ(E * )ν(E * )] τp 0 σ(θ, E * )p * (θ)dθ = ν(E * )HSC.
Hence, providing that 1 -2µ(E * )ν(E * ) = 0, we have from ( 16)

Ψ(E * , p * ) = HSC + 2µ(E * ) ν(E * )HSC 1 -2µ(E * )ν(E * ) = HSC 1 -2µ(E * )ν(E * ) ,
and from ( 18)

Ψ(E * , p * ) = E * AξΛ(E * )
. This yields that system (4)-( 6) has a steady state (p * (a), s * (a), e * (a), Epo * , GC * ), with

+∞ 0 e * (a)da := E * , if 1 -2µ(E * )ν(E * ) = 0 and E * = A.HSC.ξ Λ(E * ) 1 -2µ(E * )ν(E * ) . (21)
Conversely, suppose there exists a constant X > 0 such that 1 -2µ(X)ν(X) = 0 and X = A.HSC.ξ Λ(X) 1 -2µ(X)ν(X) .

Define

K := HSC ν(X) 1 -2µ(X)ν(X) ,
and set

             z(a) = exp - a 0 β(θ, X)dθ [HSC + 2µ(X)K] , a ∈ [0, τ p ], y(a) = exp (-ᾱ(X)a) K, a ∈ [0, τ c ],
x(a) = exp -

a 0 γ(θ)dθ A.Λ(X) [HSC + 2µ(X)K] , a > 0.
Then, it is straightforward, using [START_REF] Crauste | Mathematical study of feedback control roles and relevance in stress erythropoiesis[END_REF] and setting Epo * = F (X) and GC * = G(X), to check that (z, y, x) satisfies ( 8)-( 9) and [START_REF] Angulo | A numerical integrator for a model with a discontinuous sink term: the dynamics of the sexual phase of monogonont rotifera[END_REF]. In particular, using the definition of x and ξ one obtains

+∞ 0 x(a)da = AξΛ(X) [HSC + 2µ(X)K] ,
and using the definitions of K and X one finally gets

+∞ 0 x(a)da = A.HSC.ξΛ(X) 1 + 2µ(X) ν(X) 1 -2µ(X)ν(X) = X, so (10)
is satisfied. This defines a steady state solution of system ( 4)-( 6).

Step 2. We now prove that, under assumptions (H1) or (H2), there exists a unique E * > 0 satisfying [START_REF] Crauste | Adding self-renewal in committed erythroid progenitors improves the biological relevance of a mathematical model of erythropoiesis[END_REF]. We first set

χ(Y ) := 1 -2µ(Y )ν(Y ), Y ≥ 0,
and focus on the problem

Find Y ≥ 0 such that χ(Y ) > 0 and cY = Λ(Y ) χ(Y ) , (22)
where c := 1/(A.HSC.ξ) > 0. From ( 17) and [START_REF] Crauste | Mathematical study of feedback control roles and relevance in stress erythropoiesis[END_REF], and the previous assumptions on the functions α, σ and δ (see section 2), we obtain

d dY Λ(Y ) = τp 0 ∂ Y δ(θ, Y ) - θ 0 ∂ Y β(u, Y )du δ(θ, Y ) e -θ 0 β(u,Y )du dθ, and
d dY ν(Y ) = τp 0 ∂ Y σ(θ, Y ) - θ 0 ∂ Y β(u, Y )du σ(θ, Y ) e -θ 0 β(u,Y )du dθ, and δ(θ, Y ) > 0, σ(θ, Y ) > 0, ∂ Y δ(θ, Y ) < 0 and ∂ Y σ(θ, Y ) < 0.
First assume (H1) holds true. Then

d dY Λ(Y ) < 0 and d dY ν(Y ) < 0 for Y ≥ 0.
A first consequence, since by definition dµ(Y )/dY < 0, is that

d dY χ(Y ) = -2 ν(Y ) d dY µ(Y ) + µ(Y ) d dY ν(Y ) > 0.
Hence, the function χ is increasing on [0, +∞). Moreover, [START_REF] Angulo | Numerical approximation of singular asymptotic states for a size-structured population model with a dynamical resource[END_REF] implies that lim Y →+∞ χ(Y ) > 0. Hence, there exists Y 0 ≥ 0 such that χ(Y ) > 0 for Y ∈ (Y 0 , +∞) and χ(Y 0 ) ≥ 0. Moreover, from (H1) and as long as χ(Y ) > 0,

d dY Λ(Y ) χ(Y ) = χ(Y ) d dY Λ(Y ) -Λ(Y ) d dY χ(Y ) χ(Y ) 2 < 0.
Hence, the function Λ(Y )/χ(Y ) is decreasing on the domain (Y 0 , +∞), and Λ(Y 0 )/χ(Y 0 ) > 0. Consequently [START_REF] Gandrillon | TGF-beta cooperates with TGF-alpha to induce the self-renewal of normal erythrocytic progenitors: evidence for an autocrine mechanism[END_REF] has a unique solution Y > Y 0 . Assume (H2) holds true. Similarly to the above reasoning, Λ(Y ) and ν(Y ) are decreasing on ( Ȳ , +∞), and the function Λ(Y )/χ(Y ) is decreasing for Y > max{Y 0 , Ȳ }. We then conclude similarly to the previous case to the existence and uniqueness of Y > max{Y 0 , Ȳ } satisfying [START_REF] Gandrillon | TGF-beta cooperates with TGF-alpha to induce the self-renewal of normal erythrocytic progenitors: evidence for an autocrine mechanism[END_REF].

Under assumptions (H1) or (H2) we can then define a unique steady state to our problem. This concludes the proof.

Let us briefly comment on [START_REF] Angulo | Numerical approximation of singular asymptotic states for a size-structured population model with a dynamical resource[END_REF]. This condition implies that no more than half of the progenitor population can self-renew in order to obtain a stationary solution. Indeed, the term exp -a 0 α(0) + σ(θ, 0) + δ(θ, 0)dθ describes the survival rate of differentiating progenitors, σ(a, 0) refers to the maximal fraction of the population that will self-renew, and exp(-α(0)τ c ) is the maximal survival fraction of self-renewing progenitors. Inequality ( 11) is satisfied with our choices for α, σ and δ (see Equations ( 1), ( 2), (3)). If inequality [START_REF] Angulo | Numerical approximation of singular asymptotic states for a size-structured population model with a dynamical resource[END_REF] is not satisfied, then one may hypothesize the non-existence of a stationary solution and the unboundedness of the cell population.

Numerical Scheme

We introduce a new numerical method to simulate the system (4)- [START_REF] Adimy | Numerical integration of a mathematical model of hematopoietic stem cell dynamics[END_REF]. It is based on the integration along the characteristic curves with the use of a representation formula of the solution [START_REF] Abia | Age-structured population models and their numerical solution[END_REF].

The existence of a maximum age of differentiation, τ p , makes the differentiation rate δ unbounded. This fact has been circumvented in [START_REF] Adimy | Numerical integration of a mathematical model of hematopoietic stem cell dynamics[END_REF][START_REF] Angulo | A numerical method for nonlinear age-structured population models with finite maximum age[END_REF] with the use of an exact integration value at the end of the age-differentiation interval. However, this method is only possible when such an exact value is known, and this only happens with specific functions. In general, it is necessary to introduce a new different dependent function as in [START_REF] Iannelli | On the approximation of the Lotka-McKendrick equation with finite life-span[END_REF]. We assume, from now on, that age and Epo dependences of the function δ are independent, and we write, as in (3), δ(a, Epo) = δ 1 (a) + δ 2 (Epo).

Thus, let q(t, a) be defined as p(t, a) = Π(a) q(t, a), where

Π(a) = exp - a 0 δ 1 (x) dx ,
and, accordingly, p 0 (a) = Π(a) q 0 (a). The first equation of ( 4) is replaced with

∂ t q + ∂ a q = -(α(Epo(t)) + δ 2 (Epo(t)) + σ(a, GC(t))) q, a ∈ (0, τ p ). ( 23
)
This change also modifies the boundary condition [START_REF] Adimy | A mathematical study of the hematopoiesis process with applications to chronic myelogenous leukemia[END_REF], that writes

             q(t, 0) = HSC + 2s(t, τ c ), s(t, 0) = τp 0 σ(a, GC(t)) Π(a) q(t, a)da, e(t, 0) = A τp 0 δ(a, Epo(t)) Π(a) q(t, a)da. (24)
The unbounded age domain for erythrocytes is another circumstance that must be avoided in the numerical integration. Thus, we introduce a maximum age A max for erythrocytes as in [START_REF] Adimy | Numerical integration of a mathematical model of hematopoietic stem cell dynamics[END_REF], that would be large enough to simulate the dynamics of the model. This means that the numerical method approaches the evolution of a "truncated" version of problem ( 23)-( 24) and ( 6)- [START_REF] Adimy | Numerical integration of a mathematical model of hematopoietic stem cell dynamics[END_REF], which in turn is an approximation to the solution of the original problem (see [START_REF] Angulo | A numerical integrator for a model with a discontinuous sink term: the dynamics of the sexual phase of monogonont rotifera[END_REF]).

The formula we use in the numerical method is based on the theoretical integration along the characteristics, which gives us the next representation of the solution to the new system,

                                       q(t * + h, a * + h) = q(t * , a * ) exp - h 0 β * (a * + x, Epo(t * + x), GC(t * + x)) dx , a * + h ≤ τ p , s(t * + h, a * + h) = s(t * , a * ) exp - h 0 α(Epo(t * + x)) dx , a * + h ≤ τ c , e(t * + h, a * + h) = e(t * , a * ) exp - h 0 γ(t * + x) dx , a * + h ≤ A max , (25)
where β * (a, E, G) = α(E) + δ 2 (E) + σ(a, G). Note that, in this new layout, we have two different problems: solving [START_REF] Loeffler | A mathematical model of erythropoiesis in mice and rats. Part 1. Structure of the model[END_REF], which provides the solution to the population densities, and integrating [START_REF] Adimy | Modelling hematopoiesis mediated by growth factors with applications to periodic hematological diseases[END_REF], which gives the solution to the growth factor concentrations. We use discretization procedures in order to solve them.

We consider the numerical integration of model ( 23)-( 24) and ( 6)-( 7) along the time interval [0, T ]. Thus, given a positive integer J c , we define h = τ c /J c , J p = [τ p /h], J max = A max /h and N = [T /h] (it is usual that τ c and τ p would be proportional, and A max might be chosen to be proportional to τ c , these restrictions do not play an important role in the numerical integration). We introduce the discrete time levels t n = n h, 0 ≤ n ≤ N , and the grid nodes a j = j h, 0 ≤ j ≤ J max . Then, we denote the numerical approximations to the theoretical solutions by: Q n j to q n j = q(a j , t n ), 0 ≤ j ≤ J p ; S n j to s n j = s(a j , t n ), 0 ≤ j ≤ J c ; I n j to e n j = e(a j , t n ), 0 ≤ j ≤ J max ; O n to o n = Epo(t n ); and G n to g n = GC(t n ); 0 ≤ n ≤ N . We propose a new method to obtain them: (i) we start from an approximation of the initial data [START_REF] Adimy | Numerical integration of a mathematical model of hematopoietic stem cell dynamics[END_REF] of the problem, for instance the grid restriction of functions q 0 , s 0 , e 0 and the values Epo 0 and GC 0 , and (ii) the numerical solution at a new time level is described in terms of the previous one by means of a two steps scheme.

Such a general time level is obtained by means of the following second-order discretization of ( 24)-( 25) and [START_REF] Adimy | Modelling hematopoiesis mediated by growth factors with applications to periodic hematological diseases[END_REF]. For, 0 ≤ n ≤ N -1,

Q n+1 j+1 = Q n j exp (-h (β * (a j , O n , G n ) + β * (a j+1 , O n+1, * , G n+1, * ))/2), 0 ≤ j ≤ J p -1, S n+1 j+1 = S n j exp (-h (α(O n ) + α(O n+1, * ))/2), 0 ≤ j ≤ J c -1, I n+1 j+1 = I n j exp (-h (γ(a j , t n ) + γ(a j+1 , t n+1 ))/2), 0 ≤ j ≤ J max -1, with Q n+1 0 = HSC + 2 S n+1 Jc , S n+1 0 = Q h (σ n+1 • Π • Q n+1 ), I n+1 0 = A Q h (δ n+1 • Π • Q n+1 ),
and

O n+1 = O n + h (f (Q h (I n )) + f (Q h (I n+1, * )) -k Epo (O n + O n+1, * ))/2, G n+1 = G n + h (g(Q h (I n )) + g(Q h (I n+1, * )) -k GC (G n + G n+1, * ))/2.
Where the intermediate step is given by

Q n+1, * j+1 = Q n j exp (-h β * (a j , O n , G n )), 0 ≤ j ≤ J p -1, S n+1, * j+1 = S n j exp (-h α(O n )), 0 ≤ j ≤ J c -1, I n+1, * j+1 = I n j exp (-h γ(a j , t n )), 0 ≤ j ≤ J max -1, with Q n+1, * 0 = HSC + 2 S n+1, * Jc , I n+1, * 0 = A Q h (δ n+1, * • Π • Q n+1, * ),
and

O n+1, * = O n + h (f (Q h (I n )) -k Epo O n ), G n+1, * = G n + h (g(Q h (I n )) -k GC G n ).
One can note that initializing S n+1, * 0 is not necessary. In these formulae, Q h (V n ) represents a quadrature rule to approximate the integral on the interval [0, A max ] (approximation to the total number of erythrocytes E(t)) or [0, τ p ] (discretization of the boundary condition ( 24)) of the function with grid values V n , with different number of nodes each case, this fact is indicated by the vector so it does not cause confusion. In this case

Π j = Π(a j ), δ n+1 j = δ(a j , O n+1 ), δ n+1, * j = δ(a j , O n+1, * ), σ n+1 j = σ(a j , G n+1 ), 0 ≤ j ≤ J p , and δ n+1 • Π • Q n+1 , δ n+1, * • Π • Q n+1, * , σ n+1 • Π • Q n+1 ,
denote the componentwise product of the involved vectors. Here, a second order quadrature rule is appropriate, and we propose the composite trapezoidal quadrature rule. We point out that our method is completely explicit and also shows a good behaviour in the long-time integration to get the approximations to the theoretical steady state. We do not describe its convergence analysis because it is beyond the purposes of this work.

Numerical Experiments

We present in this section the analysis performed to experimentally check the behavior of the proposed numerical method with different test problems. The first test consists in confronting predictions of the theoretical analysis (steady state values) with the method's ability to generate a good approximation to the solution of the problem, even in the case of a long time integration allowing to approach steady state values. The second test aims at showing the capacity to obtain an accurate approximation when confronting the model to real biological data: we simulate an experimentally-induced anemia in mice and the response of the organism. The accuracy in the numerical integration is more important in this case than the long time integration.

In order to perform each test, we have used the functions α, σ and δ defined in ( 1)-( 3).

The average lifespan of an erythrocyte is known to be 40 days in mice [START_REF] Crauste | Adding self-renewal in committed erythroid progenitors improves the biological relevance of a mathematical model of erythropoiesis[END_REF], so we assume a constant mortality rate, γ(a, t) = γ = 1/40 d -1 .

Finally, the following Hill functions have been used in [START_REF] Adimy | Modelling hematopoiesis mediated by growth factors with applications to periodic hematological diseases[END_REF],

f (E) = C f θ n f f θ n f f + E n f , g(E) = C g θ ng g θ ng g + E ng , (26)
where (avoiding the subscripts) C > 0 denotes the maximum of the function, θ > 0 a threshold value, and n > 1 a sensitivity parameter.

Parameter values are listed in Tables 1 and3. Some values have been obtained from existing values in the literature, others have been estimated using a classical minimization of a weighted residual sum of squares in order to fit biological measurements (this is indicated by [NS] in Tables 1 and3).

Test Problem 1. The following numerical experiment shows the efficacy of the numerical method in simulating our model. To this end, we compare the behavior of the numerical solution with the solution provided by the theoretical analysis in section 3. In this section, existence of a unique positive stationary population density was shown, and we shall use it in the comparison of the numerical order of convergence of our scheme. We will show the optimal rate of convergence obtained with the numerical method.

We take the functions in (1)-( 3) and ( 26) and use the parameter values in Table 1.

In order to avoid discontinuities caused by an incompatible initial condition, we use

s 0 (a) = 1 7 τ p HSC σ(0, GC 0 ) (1 -a/τ c ), a ≤ τ c , e 0 (a) =    A HSC C δ1 5 + δ 2 (0, Epo 0 ) 6 (τ p -a), a ≤ τ p , 0, a > τ p ,
and for a ≤ τ p ,

p 0 (a) = HSC 1 - a τ p 5 , q 0 (a) = p(a, 0) Π(a) , (27)
= log 10 (2Epo * ), θ * σ = log 10 (3GC * ), C * f = kepoEpo * (θ n f f + (E * ) n f )/θ n f f , C * g = kgcGC * (θ ng g + (E * ) ng )/
C α 0.6 d -1 C σ 1 d -2 θ α θ * α mU.µL -1 θ σ θ * σ mU.µL -1 n α 2 N.U. n σ 2 N.U. Differentiation Rate [21], [NS] C δ1 1.15 d -1 θ δ 1 mU.µL -1 C δ2 0.1 d -1 n δ 2 N.U. Epo Production Rate [20], [NS] GC Production Rate [20], [NS] C f C * f mU.µL -1 .d -1 C g C * g mU.µL -1 .d -1 θ f 0.7E * cells.g -1 θ g 0.4E * cells.g -1 n f 7 N.U. n g 6 N.U.
which satisfy the first compatibility condition. Note that ( 27) is well defined. We complete the initial condition with Epo 0 = 5 mU.µL -1 and GC 0 = 44 mU.µL -1 . We do not know the analytical solution to the problem then, in order to compare with the numerical solution, we use the value of the steady state of the problem, given by the following formulae s * (a) = HSC σ(0, GC * ) i 1 1 -2 σ(0, GC * ) i 1 e -α(a,Epo * ) τc e -α(a,Epo * ) a , a ≤ τ c ,

e * (a) = A HSC

C δ1 i -1 + δ 2 (0, Epo * ) i 0 1 -2 σ(0, GC * ) i 1 e -α(a,Epo * ) τc e -γ a , a ≤ τ p , (29)
p * (a) = HSC e -a (α(a,Epo * )+δ2(a,Epo * )+(1-a/(2 τp)) σ(0,GC * ))

1 -2 σ(0, GC * ) i 1 e -α(a,Epo * ) τc 1 - a τ p τp C δ 1 , a ≤ τ p , (30)
q * (a) = HSC e -a (α(a,Epo * )+δ2(a,Epo * )+(1-a/(2 τp)) σ(0,GC * )) 1 -2 σ(0, GC * ) i 1 e -α(a,Epo * ) τc , a ≤ τ p , [START_REF] Pujo-Menjouet | Long period oscillations in a G 0 model of hematopoietic stem cells[END_REF] where, for r = -1, 0, 1,

i r = τp 0 1 - a τ p C δ 1 τp+r
e -[α(a,Epo * )+δ2(a,Epo * )+(1-a/(2 τp)σ(0,GC * )]a da. We have carried out an extensive numerical experimentation with different final-times , step-sizes h and maximum erythrocites age A max values. We observe that A max = 1600 is enough to obtain a good approximation to the nontruncated problem and T = 1600 produces a sufficiently long time integration in order to provide a numerical approximation to the steady state.

In Table 2, we present the results obtained with the method for different values of the step size. For each h, we compare at the final time T = 1600, the computed numerical solution U n h , with the steady state u n h , where U n h and u n h represent, respectively, P n h , S n h , I n h , O n and G n and p n h , s n h , e n h , Epo * and GC * , 0 ≤ n ≤ N . In each row of Table 2, the upper number shows the maximum error with different step sizes h (first column); that is

E u h = max 0≤n≤N U n h -u n ∞ ,
and the lower number the numerical order of convergence, which we compute with the formula [START_REF] Ackleh | A structured erythropoiesis model with nonlinear cell maturation velocity and hormone decay rate[END_REF] .

s u 2 h = log(E u 2 h /E u h ) log
Results in Table 2 clearly confirm the expected second-order of convergence.

Test Problem 2. The second test deals with the reproduction of biological measurements. Data come from an experiment of induced anemia in mice, introduced in [START_REF] Crauste | Adding self-renewal in committed erythroid progenitors improves the biological relevance of a mathematical model of erythropoiesis[END_REF]. We assume that the initial condition of the system is given by the steady state ( 28)-( 31), then we simulate the anemia (induced by a drug called phenylhydrazine that kills red blood cells upon injection on days 0 and 1, see Figure 2) and we observe the reaction of the organism that leads to full recovery within few days (see Figure 2). Thus, we are not interested on a long time integration but rather on an accurate description of transient dynamics.

Let first recall that "hematocrit" is a test that measures the volume of red blood cells in a blood sample. It gives a percentage of erythrocyte volume found in the whole blood volume. It can be considered [START_REF] Crauste | Adding self-renewal in committed erythroid progenitors improves the biological relevance of a mathematical model of erythropoiesis[END_REF] that a blood sample is mainly composed with erythrocytes and plasma, since platelet and white cell volumes can be easily neglected. The system (4)-( 6) is then numerically solved with the method proposed in section 4, and the corresponding simulated is given by the red curve. Mice are rendered anemic by two consecutive injections of phenylhydrazine on days 0 and 1. The hematocrit shows a strong fall following the anemia (the hematocrit reaches very low values), then it rapidly increases to reach a high value, and then returns to the equilibrium. The simulated hematocrit is able to correctly capture all these features. All parameter values are given by Tables 1 and3.

hematocrit HCT is computed using the formula in [START_REF] Crauste | Adding self-renewal in committed erythroid progenitors improves the biological relevance of a mathematical model of erythropoiesis[END_REF],

HCT (t) = E(t) E(t) + E * (1 -HCT * )/HCT * ,
where E(t) is the total number of erythrocytes and HCT * is the steady state value of the hematocrit, obtained through experimental data (in mice, HCT * is around 45% -50%, see Figure 2).

Again, we use the functions in ( 1)-( 3) and ( 26) to perform the simulations. We have carried out an extensive numerical experimentation with different finaltimes T , step-sizes h and values of the parameters involved in the problem. We observe that T = 200 and h = 7.8125e-3 provide sufficiently accurate approximations to the theoretical stationary state.

In [START_REF] Crauste | Adding self-renewal in committed erythroid progenitors improves the biological relevance of a mathematical model of erythropoiesis[END_REF] the authors noticed that in order to correctly reproduce experimental data (the rapid increase of the hematocrit following anemia and the damped oscillations observed after the peak of the response), a modification of the value of the mortality rate of erythrocytes should be accounted for. This has been identified as a consequence of using phenylhydrazine to induce anemia: one effect of this substance is to dramatically alter the lifespan of erythrocytes that were not killed by it (see [START_REF] Shimada | The maturation of reticulocytes. II. Life-span of red cells originating from stress reticulocytes[END_REF] for a study on chicken's erythrocytes lifespan, as well as [START_REF] Berlin | Life span of the red blood cell of the rat following acute hemorrhage[END_REF][START_REF] Nagai | Studies on the short-lived reticulocytes by use of the in vitro labeling method[END_REF][START_REF] Nagai | Studies on the erythroblast dynamics based on the production of fetal hemoglobin[END_REF][START_REF] Stohlman | Humoral regulation of erythropoiesis. VII. Shortened survival of erythrocytes by erythropoietin or severe anemia[END_REF]). We then introduce a modification of the erythrocyte with t i the time of the injection, t f the duration of the phenylhydrazine effect in the blood, R a nonnegative residual effect, K a positive constant which determines the maximum effect of the phenylhydrazine injection, and m an exponent associated with the clearance of the phenylhydrazine while in blood.

In order to compare the model simulations with the experimental data, we then reproduced the experimental protocol consisting in two injections of phenylhydrazine, with the initial injection occurring at time t = 0 and the second injection at time t 2 = 1 day. Since experimentally both injections are similar (same dose, same route of injection), we assumed that they both have the same effects and we used the same parameter values for the function phz, except for the value of t i . Additional parameter values are shown in Table 3.

Figure 2 shows the simulation performed numerically with γ = 1/ √ 40, parameter data in Tables 1, and 3, and associated experimental data. One can appreciate how the numerical solution reproduces the behavior of the phenylhydrazine action in inducing anemia and the recovery of the organism.

Conclusions

We considered an age-structured partial differential equation system, coupled to a system of ordinary differential equations, which generalizes a previous model published in [START_REF] Crauste | Adding self-renewal in committed erythroid progenitors improves the biological relevance of a mathematical model of erythropoiesis[END_REF] to describe erythropoiesis. In addition to cell population dynamics, this new model incorporates growth factor dynamics, and explicitly considers an age-dependency of feedback control functions describing immature cell self-renewal and differentiation. The complexity of the mathematical model we considered limited the theoretical analysis that could be performed, however a study about the existence of steady states has been performed.

In order to perform a more complex analysis of the model, we needed to introduce numerical approximation to the solution. Therefore, we presented an "ad hoc" new numerical method to generate a simulated solution of the system.

We have studied numerically both its accuracy and its potential to simulate real biological behaviors, based on experimental data. On one hand, we have shown numerically a second order of convergence. On the other hand, the method proved able to accurately replicate an experimental protocol, consisting in injecting mice with a substance that kills red blood cells and induces anemia, which has been modeled following [START_REF] Angulo | Investigating the Roles of the Experimental Protocol in Phenylhidrazine-Induced Anemia in Mice[END_REF]. Thus, the numerical scheme shows a second order convergence in a finite time integration, a good behavior in a longtime one and it is able to simulate and replicate biological data.

Figure 1 :

 1 Figure 1: Schematic Model of Erythropoiesis. Plain lines with arrows represent cell process (differentiation, death, proliferation, cell cycle), whereas dashed lines represent feedback controls.

Figure 2 :

 2 Figure 2: Experimental and Simulated Hematocrits. Experimental data are given by the blue dashed line, with error bars on every experimental time point. The simulatedis given by the red curve. Mice are rendered anemic by two consecutive injections of phenylhydrazine on days 0 and 1. The hematocrit shows a strong fall following the anemia (the hematocrit reaches very low values), then it rapidly increases to reach a high value, and then returns to the equilibrium. The simulated hematocrit is able to correctly capture all these features. All parameter values are given by Tables1 and 3.

Table 1 :

 1 Parameter values. Sources of parameter values are given between brackets, with [NS] = values obtained by performing numerical simulations. The values of θα, θσ, C f and Cg are respectively given by θ * α

Table 2 :

 2 Test problem 1. Error and numerical convergence order, for T = 1600.

	h	E p h	E s h	E e h	E Epo h	E GC h
	3.125e-2	5.6446e-3 3.4988e-3 5.1756e-1 2.3044e-6 1.6660e-5
	1.5625e-2	1.4187e-3 8.8050e-4 1.3073e-1 5.7095e-7 4.1277e-6 1.99 1.99 1.99 2.01 2.01
	7.8125e-3	3.5520e-4 2.2052e-4 3.2773e-2 1.4241e-7 1.0296e-6 2.00 2.00 2.00 2.00 2.00
	3.9063e-3	8.8846e-5 5.5164e-5 8.2001e-3 3.5589e-8 2.5729e-7 2.00 2.00 2.00 2.00 2.00
	1.9531e-3	2.2236e-5 1.3806e-5 2.0524e-3 8.9049e-9 6.4378e-8 2.00 2.00 2.00 2.00 2.00

Table 3 :

 3 Parameter values. Sources of parameter values are given between brackets, with [NS] = values obtained by performing numerical simulations. (N.U. = no unit) = γ (1 + phz(a, t, 0, t f ) + phz(a, t, t 2 , t f )) , where the function phz accounts for the effect of phenylhydrazine, as follows phz(a, t, t i , t f )=

	Par.		Value	Unit	Par.	Value Unit
	Steady States [21], [NS]
	E *	1.72 × 10 7	cells	HCT *	48%	N.U.
	Phenylhydrazine Injections [NS]
	R		2	d -1	K	9	d -1
	t f		1	d	m	4	N.U.
	t 2		1	d
	mortality rate γ(a, t), as follows
	γ(t, a) 	0,		t ≤ t
		   
		   

i or a ≤ t -t i , K (t i + t f -t) m t m f + R, t i ≤ t ≤ t i + t f and a ≥ t -t i , R, t i + t f ≤ t and a ≥ t -t i ,

Aknowledgements

FC work has been supported by ANR grant ProCell ANR-09-JCJC-0100-01. OA and JCLM have been supported in part by project MTM2014-56022-C2-2-P of the Spanish Ministerio de Economía y Competitividad and European FEDER Funds, and by Consejería de Educación, JCyL, project VA191U13.


Alain Barrat

Andrea Baronchelli

Luca Dall'asta

Vittorio Loreto

Agreement dynamics on interaction networks with diverse topologies

We review the behavior of a recently introduced model of agreement dynamics, called the Naming Game. This model describes the self-organized emergence of linguistic conventions and the establishment of simple communication systems in a population of agents with pairwise local interactions. The mechanisms of convergence towards agreement strongly depend on the network of possible interactions between the agents. In particular, the mean-field case in which all agents communicate with all the others is not efficient, since a large temporary memory is requested for the agents. On the other hand, regular lattice topologies lead to a fast local convergence but to a slow global dynamics similar to coarsening phenomena. The embedding of the agents in a small-world network represents an interesting trade-off: a local consensus is easily reached, while the long-range links allow to bypass coarsening-like convergence. We also consider alternative adaptive strategies which can lead to faster global convergence.

Numerous statistical physics models have been developed or adapted for the description and study of social behaviour and social phenomena, such as opinion dynamics and consensus formation. The Voter model is a well-known example: each agent can has one of two possible opinions; at each time step, an agent picks at random one of his neighbours and adopts his opinion. The possible convergence to a state in which all agents share the same opinion has been particularly studied for agents on the nodes of regular lattices (or in the case of all agents interacting with all others). This is also the case for most "sociophysics" models. Recently however, the study of complex networks has put in evidence the fact that social networks are typically very heterogeneous, and in all cases very different from regular lattices. The consequences of the topology of the interaction network on the dynamical behaviour of all these models has thus become of great relevance.

We consider a model for the decentralized formation of a communication system among interacting agents. The model has important differences with other usual opinion models, such as the existence of memory. We have studied the model's dynamics for agents interacting on various types of networks, and shown in particular that the formation of a common communication system is more easily obtained if the agents form a small-world network: in fact, and maybe counterintuitively, the possibility for the agents to all interact with each other allows as well to reach a consensus on the communication system, but in a less efficient way. We have also studied the convergence phenomenon for various networks, putting in evidence the main mechanisms and the effect of various parameters.

I. INTRODUCTION

The recent past has witnessed an important development of the activities of statistical physicists in the area of social sciences (for a recent collection of papers see [START_REF]Eighth Granada Lectures on Modeling cooperative behaviour in the social sciences[END_REF]). Indeed, the standard methods of statistical physics are very appropriate to study collective behaviors, neglecting details and retaining only few general ingredients observed in real social interactions. For this reason, physicists have put forward a large number of theoretical models of social dynamics, borrowing a suite of statistical methods from the theory of interacting particles systems [START_REF] Ligget | Interacting particle systems[END_REF][START_REF] De Oliveira | Evolution, Money, War and Computers -Non-Traditional Applications of Computational Statistical Physics[END_REF][START_REF] Durlauf | Statistical Mechanics Approaches to Socioeconomic Behavior[END_REF][START_REF] Blume | The statistical Mechanics of Social Interaction[END_REF]. In particular, many models have been proposed for the study of opinion formation, such as the Voter model [START_REF] Ligget | Interacting particle systems[END_REF][START_REF] Krapivsky | Kinetics of monomer-monomer surface catalytic reactions[END_REF][START_REF] Krapivsky | Dynamics of Majority Rule in Two-State Interacting Spin Systems[END_REF][START_REF] Ben -Naim | Coarsening and Persistence in the Voter Model[END_REF][START_REF] Sood | Voter Model on Heterogeneous Graphs[END_REF][START_REF] Frachebourg | Exact results for kinetics of catalytic reactions[END_REF], the Sznajd-Weron model [START_REF] Sznajd-Weron | Opinion evolution in closed community[END_REF], the Axelrod model for the dissemination of culture [START_REF] Axelrod | The Dissemination of Culture: A Model with Local Convergence and Global Polarization[END_REF], and their variants (e.g. the models proposed by Deffuant et al. [START_REF] Deffuant | Mixing beliefs among interacting agents[END_REF] and by Krause and Hegselmann [START_REF] Hegselmann | Opinion dynamics and bounded confidence models,. analysis, and simulation[END_REF]).

The behavior of these models has been much studied on regular topologies or in situations where each agent can interact with all the others. Recently however, network science [15][START_REF] Dorogovtsev | Evolution of Networks: from biological nets to the Internet and WWW[END_REF][START_REF] Pastor-Satorras | Evolution and Structure of the Internet: A statistical physics approach[END_REF] has led to a better knowledge of the topological properties of real social groups [START_REF] Granovetter | The strength of weak ties[END_REF], and in particular to show that the topology of the network on which agents interact is not regular. Models of social dynamics have thus been reconsidered, in order to integrate the new framework of complex networks, and to study the influence of various complex topologies on the corresponding dynamical behavior.

In this article, we review the behavior of a recently proposed model for the emergence of a communication system, called Naming Game, investigating its dynamics both on regular topologies and complex networks [START_REF] Baronchelli | Topology Induced Coarsening in Language Games[END_REF][START_REF] Baronchelli | Strategies for fast convergence in semiotic dynamics[END_REF][START_REF] Dall'asta | Agreement dynamics on small-world networks[END_REF][START_REF] Dall'asta | Non-equilibrium dynamics of language games on complex networks[END_REF].

Social interactions are indeed based on the existence of a communication system among the agents, who are able to understand each other by means of common linguistic patterns or, more generally, by means of a common vocabulary of symbols. Such a communication system is the result of a self-organized process in which individuals select specific symbols (words) and associate them to concepts and ideas (objects). The emergence of a shared lexicon inside social groups and communities of people is very likely to be driven by simple criteria, like popularity, imitation, negotiation, and agreement. When a new concept is introduced, people refer to it using several different names or words. These words start spreading among the population, competing one against the other, until the choice of one of them is taken (with a sudden transition or with a long process) and everyone uses the same word (or symbol, etc) [START_REF] Lass | Historical Linguistics and Language Change[END_REF][START_REF] Briscoe | Linguistic evolution through language acquisition: formal and computational models[END_REF][START_REF]Approaches to the Evolution of Human Language[END_REF]. This kind of dynamics has recently become of broad interest after the diffusion of a new generation of web-tools which enable human users to self-organize a system of tags in such a way to ensure a shared classification of information about different arguments (see, for instance, del.icio.us or www.flickr.com and Refs. [START_REF] Cattuto | Semiotic dynamics and collaborative tagging[END_REF][START_REF] Golder | The structure of collaborative tagging systems[END_REF]). Another application concerns global coordination problems in artificial intelligence, where a group of artificial embodied agents moving in an unknown environment have to exchange informations about the objects they gradually discover. The emergence of consensus about the objects names allows to establish a communication system. A practical example of this type of dynamics is provided by the well-known Talking Heads experiment [START_REF] Steels | The Origins of Ontologies and Communication Conventions in Multi-Agent Systems[END_REF][START_REF] Steels | The Talking Heads Experiment. Volume 1. Words and Meanings, Laboratorium[END_REF], in which embodied software agents managed to develop their vocabulary observing objects through digital cameras, assigning them randomly chosen names and negotiating these names with other agents.

The paper is organized as follows. In Sec. II we introduce and discuss the Naming Game. In Sec. III we describe the main phenomenology of the model in the meanfield case, i.e. when the population is completely unstructured and each agent can in principle interact with anybody else. In the following sections different interaction patterns are investigated: regular lattices (Sec. IV), small-world networks (Sec. V), and different kinds of complex networks (Sec. VI). In Sec. VII, the microscopic rules of the model are modified in order to improve the efficiency of the model. Finally, conclusions are drawn in Sec. VIII.

II. MODEL DEFINITION

A minimal model of Naming Game has been put forward by Baronchelli et al. in Ref. [START_REF] Baronchelli | [END_REF] to reproduce the main features of Semiotic Dynamics and the fundamental results of adaptive coordination observed in the Talking Heads experiment. The minimal Naming Game model consists of a population of N agents observing a single object, for which they invent names that they try to com-FIG. 1: Agents interaction rules. Each agent is described by its inventory, i.e. the repertoire of known words. The speaker picks up at random a name in its inventory and transmits it to the hearer. If the hearer does not know the selected word the interaction is a failure (top), and it adds the new name to its inventory. Otherwise (bottom), the interaction is a success and both agents delete all their words but the winning one. Note that if the speaker has an empty inventory (as it happens at the beginning of the game), it invents a new name and the interaction is a failure. municate to one another through pairwise interactions, in order to reach a global agreement. The agents are identical and dispose of an internal inventory, in which they can store an a priori unlimited number of names (or opinions). All agents start with empty inventories. At each time step, a pair of neighboring agents is chosen randomly, one playing as "speaker", the other as "hearer", and negotiate according to the following rules (see also Fig. 1):

• the speaker selects randomly one of its words and conveys it to the hearer;

• if the hearer's inventory contains such a word, the two agents update their inventories in order to keep only the word involved in the interaction (success);

• if the hearer does not possess the uttered word, the latter is added to those already stored in the hearer's inventory (failure), i.e. it learns the word.

Before entering in the detailed description of the dynamics, it is worth noting some visible differences of the Naming Game with other commonly studied models of social dynamics and, in particular, of opinion formation [START_REF] Ligget | Interacting particle systems[END_REF][START_REF] Sznajd-Weron | Opinion evolution in closed community[END_REF][START_REF] Axelrod | The Dissemination of Culture: A Model with Local Convergence and Global Polarization[END_REF][START_REF] Deffuant | Mixing beliefs among interacting agents[END_REF]. First of all, each agent can potentially be in an infinite number of possible discrete states (words, names, opinions), and the maximum number of states depends on the dynamical evolution itself (see [START_REF] Dall'asta | Microscopic activity patterns in the naming game[END_REF] for a detailed analysis of this point). The two-steps decision process is moreover rather realistic: an agent can accumulate in its memory different possible names for the object, waiting before reaching a decision. These points are strongly in contrast with traditional models (Voter, Potts, etc) in which the number of states is a fixed external parameter taking finite (and usually small) values [START_REF]Potts model with q states has a well defined q → ∞ limit[END_REF], and in which each agent has a unique well-defined opinion at each time step. Each dynamical step can be seen as a negotiation between speaker and hearer, with a certain degree of stochasticity, that is absent in deterministic models such as the Voter model. The stochastic component is however of a different nature compared to that of standard Glauber dynamics used in majority rule models [START_REF] Glauber | Time-dependent statistics of the Ising model[END_REF], since here it comes from an internal selection criterion, and involves only the speaker, without affecting the (deterministic) decision process of the hearer.

An important remark also concerns the random extraction of the word in the speaker's inventory. Most previously proposed models of semiotic dynamics attempted to give a more detailed representation of the negotiation interaction assigning weights to the words in the inventories. In such models, the word with largest weight is automatically chosen by the speaker and communicated to the hearer. Success and failures are translated into updates of the weights: the weight of a word involved in a successful interaction is increased to the detriment of those of the others (with no deletion of words); a failure leads to the decrease of the weight of the word not understood by the hearer. An example of a model including weights dynamics can be found in Ref. [START_REF] Lenaerts | The evolutionary language game: An orthogonal approach[END_REF] (and references therein). For the sake of simplicity the minimal Naming Game avoids the use of weights. Indeed, these are apparently more realistic, but their presence is not essential for the emergence of a global collective behavior of the system [START_REF] Baronchelli | Statistical mechanics approach to language games[END_REF].

Finally, we stress that, in the minimal Naming Game, all agents refer to the same single object, while in the original experiments the embodied agents could observe a set of different objects. This is actually possible only if we assume that homonymy is excluded, i.e. two distinct objects cannot have the same name. Consequently, in this model, all objects are independent and the general problem reduces to a set of independently evolving systems, each one described by the minimal model. In more realistic situations, however, homonymy should probably be taken into account.

III. MEAN-FIELD

Many studies of social dynamics have focused on populations of agents in which all pairwise interactions are allowed, i.e. the agents are placed on the vertices of a fullyconnected graph. In statistical mechanics, this topological structure is commonly referred as "mean-field" topology. The main quantities of interest which describe the system's evolution are [START_REF] Baronchelli | [END_REF] • the total number N w (t) of words in the system at the time t (i.e. the total size of the memory);

• the number of different words N d (t) in the system at the time t;

• the average success rate S(t), i.e. the probability, computed averaging over many simulation runs, that the chosen agent gets involved in a successful interaction at a given time t.

The consensus state is obtained when N d = 1 and N w = N (so that S = 1). The temporal evolution of the three main quantities is depicted in Fig. 2 (circles). At the beginning, many disjoint pairs of agents interact, with empty initial inventories: they invent a large number of different words (N/2, on average), that start spreading throughout the system, through failure events. Indeed, the number of words decreases only by means of successful interactions. In the early stages of the dynamics, the overlap between the inventories is very low, and successful interactions are limited to those pairs which have been chosen at least twice. Since the number of possible partners of an agent is of order N , it rarely interacts twice with the same partner, the probability of such an event growing as t/N 2 . The initial trend of S(t) (black circles) is indeed linear with a slope of order 1/N 2 . In this phase of uncorrelated proliferation of words, the number of different words N d invented by the agents grows, rapidly reaching a maximum that scales as O(N ). Then N d saturates, displaying a plateau, during which no new word is invented anymore (since every inventory contains at least one word). The total number N w of words stored in the system has a similar behavior, but it keeps growing after N d has saturated, since the words continue to propagate throughout the system even if no new one is introduced. The peak of N w has been shown to scale as O(N 1.5 ) [START_REF] Baronchelli | [END_REF], meaning that each agent stores O(N 0.5 ) words. This peak occurs after the system has evolved for a time t max ∼ O(N 1.5 ). In the subsequent dynamics, strong correlations between words and agents develop, driving the system to a final rather fast convergence to the absorbing state in a time t conv ∼ O(N 1.5 ) (note that throughout this review, we do not show data corresponding to the scaling of t conv with N , except in section VII, but refer the interested reader to our previously published papers [START_REF] Baronchelli | Topology Induced Coarsening in Language Games[END_REF][START_REF] Dall'asta | Agreement dynamics on small-world networks[END_REF][START_REF] Dall'asta | Non-equilibrium dynamics of language games on complex networks[END_REF][START_REF] Baronchelli | [END_REF]). The S-shaped curve of the success rate in Fig. 2 summarizes the dynamics: initially, agents hardly understand each others (S(t) is very low); then the inventories start to present significant overlaps, so that S(t) increases until it reaches 1, and the communication system is completely set in.

IV. COARSENING PHENOMENON ON LATTICES

A first study of the effects of topological embedding on the Naming Game dynamics is reported in Ref. [START_REF] Baronchelli | Topology Induced Coarsening in Language Games[END_REF]. When the interacting agents sit on the nodes of lowdimensional lattices, the long-time behavior is still characterized by the convergence to a homogeneous consensus state, but the evolution of the system changes considerably. In particular, the time required by the system to reach the global consensus displays a different scaling with the size N , and the effective size of the inventories is considerably diminished. Actually, the existence of different dynamical patterns are clearly visible in Fig. 2. Since each agent can interact only with a limited number of neighbors (2d in a d-dimensional lattice), at the local scale the dynamics is very fast: agents can rapidly interact two or more times with their neighbors, favoring the establishment of a local consensus with a high success rate, i.e. of small sets of neighboring agents sharing a common unique word. These "clusters" of neighboring agents with a common unique word are separated by individuals having a larger inventory with two or more words, playing the role of "interfaces". These interfaces then start a diffusion process, and the clusters of unique words grow in time with a law that is typical of coarsening phenomena [START_REF] Baronchelli | Topology Induced Coarsening in Language Games[END_REF]: the competition among the clusters is driven by the fluctuations of the interfaces. The easily reached local consensus thus leads to a slow dynamics, and the global consensus takes much longer to be reached than in mean-field: for example, O(N 3 ) in dimension 1 vs. O(N 1.5 ) in mean-field. However, another important aspect of the problem concerns the memory used by the agents. In mean-field indeed, each agent needs a memory capacity scaling as O(N 1/2 ), i.e. diverging with the system size. In contrast, the consequence of the embedding in a finite-dimensional lattice (with a finite number of neighbors), and of the subsequent coarsening like phenomena, with rapid local consensus, is that each agent uses only a finite capacity: the maximum total number of words in the system (maximal memory capacity) scales linearly with the system size N (as for the number of different words). In summary, low-dimensional lattice systems require more time to reach the consensus compared to mean-field, but a lower use of memory.

V. THE TRADE-OFF OF SMALL-WORLD NETWORKS

The precise knowledge of the dynamical behavior of the Naming Game model on low-dimensional lattices, and in particular on the one-dimensional ring, makes possible to understand, by means of simple arguments and numerical simulations, the effect of the small-world property, that is a relevant feature of real complex networks.

In the following, indeed, we investigate the effect of introducing long-range connections which link agents that are far from each other on the regular lattice. In other words, we study the Naming Game on the small-world model proposed by Watts and Strogatz [START_REF] Watts | Collective dynamics of small-world networks[END_REF]. Starting from a quasi-one-dimensional banded network in which each node has 2m neighbors, the edges are rewired with probability p, i.e. p represents the density of long-range connections introduced in the network. For p = 0 the network retains a one-dimensional topology, while the random network structure is approached as p goes to 1. At small but finite p (1/N p 1), a small-world structure with short distances between nodes, together with a large clustering, is obtained. When p = 0, the system is one-dimensional and the dynamics proceeds by slow coarsening. At small p, the typical distance between shortcuts is O(1/p), so that the early dynamics is not affected and proceeds as in one-dimensional systems. In particular, at very short times many new words are invented since the success rate is small. The maximum number of different words scales as O(N ), as in the other cases, while the average used memory per agent remains finite, since the number of neighbors of each site is bounded (the degree distribution decreases exponentially [START_REF] Barrat | On the properties of smallworld network models[END_REF]).

The typical cluster dynamics on a small-world network is graphically represented in Fig. 3. As long as the typical cluster size is smaller than 1/p, the clusters are essentially one-dimensional, and the system evolves by means of the usual coarsening dynamics. However, as the average cluster size reaches the typical distance between two shortcuts ∼ 1/p, a crossover phenomena toward an accelerated dynamics takes place. Since the cluster size grows as t/N, this corresponds to a crossover time t cross = O(N/p 2 ). For times much larger than this crossover, one expects that the dynamics is dominated by the existence of shortcuts, entering a mean-field like behavior. The convergence time is thus expected to scale as N 3/2 and not as N 3 . The condition in order for this picture to be possible is exactly the small-world condition; indeed, the crossover time N/p 2 has to be much larger than 1, and much smaller than the consensus time for the one-dimensional case N 3 , that together imply p 1/N . Figure 4 displays the evolution of the average number of words per agent as a function of time, for a smallworld network with average degree k = 8, and various values of the rewiring probability p. While N w (t) in all cases decays to N , after an initial peak whose height is proportional to N , the way in which this convergence is obtained depends on the parameters. At fixed N , for p = 0 a power-law behavior N w /N -1 ∝ 1/ √ t is observed due to the one-dimensional coarsening process. As soon as p 1/N however, we observe deviations from this power-law behavior, which get stronger as p is increased: the decrease of N w is first slowed down after the peak, but leads in the end to a very fast convergence.

As previously mentioned, a crossover phenomenon is expected when the one-dimensional clusters reach sizes of order 1/p, i.e. at a time of order N/p 2 . Since the agents with more than one word in memory are localized at the interfaces between clusters, their number is O(N p). The average excess memory per site (with respect to global consensus) is thus of order p, so that one expects N w /N -1 = pG(tp 2 /N ). Detailed numerical investigations have confirmed this picture [START_REF] Dall'asta | Agreement dynamics on small-world networks[END_REF], and allowed to show that the convergence towards consensus is reached on a timescale of order N βSW , with β SW ≈ 1.4±0.1 (not shown, see [START_REF] Dall'asta | Agreement dynamics on small-world networks[END_REF][START_REF] Dall'asta | Non-equilibrium dynamics of language games on complex networks[END_REF] for data). This behavior is close enough to the meanfield case N 3/2 to consider that they belong to the same universality class, as expected from the above arguments and other studies of dynamical phenomena on Watts-Strogatz networks [START_REF] Barrat | On the properties of smallworld network models[END_REF], and in strong contrast with the N 3 behavior of purely one-dimensional systems. Note that the time to converge scales as p -1.4±.1 (not shown), that is consistent with the fact that for p of order 1/N one should recover an essentially one-dimensional behavior with convergence times of order N 3 .

In summary, the small-world topology allows to combine advantages from both finite dimensional lattices and mean-field networks: on the one hand, only a finite memory per node is needed, in opposition to the O(N 1/2 ) in mean-field; on the other hand the convergence time is expected to be much shorter than in finite dimensions.

VI. THE NAMING GAME ON GENERAL COMPLEX NETWORKS

A. Pair selection strategies

Before describing the behavior of the Naming Game dynamics on general networks, it is worth noting that the definition of the model itself has in fact to be specified. Indeed, the two neighboring agents chosen to interact have different roles: one (the speaker) transmits a word and is thus more "active" than the other (the hearer). One should therefore specify whether, when choosing a pair, one chooses first a speaker and then a hearer among the speaker's neighbors, or the reverse order. If the agents sit on a fully connected graph or on a regular lattice, or even on a random graph with homogeneous degree distribution, they have an equivalent neighborhood so the order is not important. In the case of heterogeneous networks however, the degrees of the first and the second chosen nodes can have very different distributions (respectively P (k) and kP (k)/ k ). The asymmetry between speaker and hearer can couple to the asymmetry between a randomly chosen node and its randomly chosen neighbor, leading to different dynamical properties (this is the case for example in the Voter model, as studied by Castellano [START_REF] Castellano | Effect of network topology on the ordering dynamics of voter models[END_REF]). We therefore distinguish more possibilities for the definition of the Naming Game on generic networks.

• (i) A randomly chosen speaker selects (again randomly) a hearer among its neighbors. This is probably the most natural generalization of the original rule. We call this strategy direct Naming Game. In this case, larger degree nodes will preferentially act as hearers.

• (ii) The opposite strategy, here called reverse Naming Game, can also be carried out: we choose the hearer at random and one of its neighbors as speaker. In this case the hubs are preferentially selected as speakers.

• (iii) A neutral strategy to pick up pairs of nodes is that of considering the extremities of an edge taken uniformly at random. The role of speaker and hearer are then assigned randomly with equal probability among the two nodes.

As shown in [START_REF] Baronchelli | Bootstrapping communication in language games. Strategy, topology and all that[END_REF], a larger memory is used for the reverse rule, although the number of different words created is smaller, and a faster convergence is obtained. This corresponds to the fact that the hubs, playing principally as speakers, can spread their words to a larger fraction of the agents, and remain more stable than when playing as hearers, enhancing the possibility of convergence. Depending on the network under study, and similarly to the Voter model case [START_REF] Castellano | Effect of network topology on the ordering dynamics of voter models[END_REF], the scaling laws of the convergence time can even be modified. From the point of view of a realistic interaction among individuals or computer-based agents, the direct Naming Game in which the speaker chooses a hearer among its neighbors seems somehow more natural than the other ones. In the remainder of this section therefore, we will focus on the direct Naming Game.

B. Global quantities

Figure 5 reports, for homogeneous Erdös-Renyi (ER) networks (left) and heterogeneous Barabási-Albert (BA) networks (right), the temporal evolution of the three main global quantities: the total number N w (t) of words in the system, the number of different words N d (t), and the rate of success S(t). The curves for the average use of memory N w (t) show a rapid growth at short times, a peak and then a plateau whose length increases as the size of the system is increased. The time and the height of the peak, and the height of the plateau, are proportional to N . A systematic study of the scaling behavior shows that the convergence time t conv scales as N β with β 1.4 for both ER and BA (see [START_REF] Dall'asta | Non-equilibrium dynamics of language games on complex networks[END_REF]). The apparent plateau of N w , on the other hand, does not correspond to a steady state, as revealed by the continuous decrease of the number of different words N d in the system: in this re-organization phase, the system keeps evolving by elimination of words, although the total used memory does not change significantly. The observed scaling law for the convergence time is a general robust feature that is not affected by other topological details (average degree, clustering, etc), and more surprisingly it seems to be independent of the particular form of the degree distribution. We have indeed checked the value of the exponent β 1.4±0.1 for various k , clustering, and exponents γ of the degree distribution P (k) ∼ k -γ for scale-free networks constructed with the uncorrelated configuration model. These parameters have instead an effect on other quantities such as the time and the value of the maximum of memory, as shown in Fig. 6, which displays the effects of increasing the average degree on the behavior of the main global quantities. In both ER (left) and BA (right) models, increasing the average degree provokes an increase in the memory used, while the global convergence time is decreased: there is a trade-off between memory and rapidity. We have also investigated the effect of increasing the clustering at fixed average degree and degree distributions: the number of different words is not changed, but the average memory used is smaller and the convergence takes more time: it is more probable for a node to speak to 2 neighbors that share common words because they are themselves connected and have already interacted, so that it is less probable to learn new words. At fixed average degree, i.e. global number of links, less connections are available to transmit words from one part of the network to the other since many links are used in "local" triangles. The local cohesiveness is therefore in the long run an obstacle to the global convergence. This effect is similar to the observation of an increase in the percolation threshold in clustered networks, due to the fact that many links are "wasted" in redundant local connections [START_REF] Serrano | Percolation and epidemic thresholds in clustered networks[END_REF].

C. Cluster statistics

Further insights into the convergence mechanisms is obtained by the investigation of the behavior of clusters of words (recall that a cluster is any set of neighboring agents sharing a common unique word). In lowdimensional lattices indeed, the dynamics of the Naming Game proceeds by formation of such clusters, that grow through a coarsening phenomenon. As shown instead in Fig. 7 for ER networks (a similar behavior is observed in the BA case), the normalized average cluster size remains very close to zero (in fact, of order 1/N ) during the reorganization phase that follows the peak in the number of words, and converges to one with a sudden transition. The same behavior is shown also by the number of clusters N cl (t), that decreases to one very sharply. The emerging picture is not that of a coarsening or growth of clusters, but that of a slow process of correlations between inventories, followed by a multiplicative process of cluster growth triggered by a sort of symmetry breaking event in the success probability of the words (in favor of the word that will ultimately survive).

VII. STRATEGIES FOR FASTER CONVERGENCE

In all the investigated cases, the time to convergence grows quite fast as a function of the system size. A natural and important question is therefore whether it is possible to improve the performance of the system. More precisely, a major challenge would be to improve the population-scale performances of the process without loosing the simplicity of the microscopic rules, which is the precious ingredient that allows for in-depth investigations of global-scale dynamics. To investigate this issue, we come back for simplicity to the framework of meanfield dynamics, in which all agents can interact with all others.

Up to this point agents, when playing as speakers, extract randomly a word in their inventories. This feature, along with the drastic deletion rule that follows a successful game, is the distinctive trait of the model. Other previously proposed models use the possibility of associating a weight to each word in each inventory [START_REF] Lenaerts | The evolutionary language game: An orthogonal approach[END_REF]: weights are updated at each interaction, with rewards for winning words and penalties for the others. While such sophisticated structures could in principle lead to faster convergence, they also make the models more complicated, compromising the possibility of a clear global picture of the convergence process. We follow here a different route which maintains the simplicity of the dynamical rules. Among the words of a given agent, two words can be easily distinguished: the last recorded one and the last one that gave rise to a successful game, i.e. the first that was recorded in the new inventory generated after the successful interaction. Natural strategies to investigate consist therefore in choosing systematically one of these particular words. We shall refer to these strategies as "play-last" and "play-first" respectively.

A. "Play-last" strategy

When the "play-last" strategy is adopted, the peak time and height scale respectively as t max ∼ N α with α ≈ 1.3 and N max w ∼ N γ with γ ≈ 1.3, i.e. the used memory is reduced, while the convergence time scales as t conv ∼ N β with β ≈ 2.0 (not shown). At the beginning of the process, playing the last registered word creates a positive feedback that enhances the probability of a success. In particular a circulating word has more probabilities of being played than with the usual stochastic rule, thus creating a scenario in which less circulating words are known by more agents. On the other hand the "last in first out" approach is highly ineffective when agents start to win, i.e. after the peak. In fact, the scaling t conv ∼ N β can be explained through simple analytical arguments. Let us denote by N a the number of agents having the word "a" as last recorded one. This number can increase by one unit if one of these agents is chosen as speaker, and one of the other agents is chosen as hearer, i.e. with probability N a /N × (1 -N a /N ); the probability to decrease N a of one unit is equal to the probability that one of these agents is a hearer and one of the others is a speaker, i.e. (1 -N a /N )N a /N . These two probabilities are perfectly balanced so that the resulting process for the density ρ a = N a /N can be written as an unbiased random walk (with actually a diffusion coefficient ρ a (1ρ a )/N 2 ); it is then possible to show that the time necessary for one of the ρ a to reach 1 is of order N 2 . In summary, in this framework it is much more difficult to bring to convergence all the agents, since each residual competing word has a good probability of propagating to other individuals.

B. "Play-first" strategy

The "play-first" strategy, on the other hand, leads to a faster convergence. Due to a sort of arbitrariness in the strategy before the first success of the speaker, the peak related quantities keep scaling as in the usual model, so that t max ∼ N α and N max w ∼ N γ with α ≈ γ ≈ 1.5. This seems natural, since playing the first recorded word is essentially the same as extracting it randomly when most agents have only few words. In fact, in both cases no virtuous correlations or feedbacks are introduced between circulating and played words. However, playing the last word which gave rise to a successful interaction strongly improves the system-scale performances once the agents start to win. In particular it turns out that for the difference between the peak and convergence time we obtain (t conv -t max ) ∼ N δ with δ ≈ 1.15, so that the behavior of the convergence time is the result of the combination of two different power law regimes, i.e. t conv ∼ aN α + bN δ . On the other hand, the usual stochastic rule leads to (t convt max ) ∼ N 1.5 . This means that the "play-first" strategy is able to reduce the time that the system has to wait before reaching the convergence, after the peak region. This seems the natural consequence of the fact that successful words increase their chances to be played while suppressing the spreading of other competitors. FIG. 9: Success rate curves S(t) for the various strategies: stochastic, "play-last", "play-first" and "play-smart". At the beginning of the process the stochastic and "play-first" strategies yield similar success rates, but then the deterministic rule speeds up convergence. On the other hand also the "playsmart" and the "play last" evolve similarly at the beginning, but the latter reaches the final state much earlier through a steep jump. It is worth noting that for three strategies the S(t) curves present a characteristic S-shaped behavior, while in the "play-last" one the disorder-order transition is more continuous (see inset in the top figure). All curves, both for N = 10 3 and N = 10 4 , have been generated averaging over 3 × 10 3 simulation runs.

C. "Play-smart", an adaptive strategy

Compared to the usual random extraction of the played word, the "play-last" strategy is more performing at the beginning of the process, while the "play-first" one allows to fasten the convergence of the process, even if it ). The "play-smart" approach combines the advantages of "play-last" and "playfirst" strategies.

is effective only after the peak of the total number of words. We therefore define a third alternative strategy which results from the combination of the two. The new prescription, called "play-smart", is the following:

• If the speaker has never took part in a successful game, it plays the last word recorded;

• Else, if the speaker has won at least once, it plays the last word it had a communicative success with.

The first rule will thus be applied mostly at the beginning, and as the system evolves, the second rule will be progressively adopted by more and more agents. Since the change in strategy is not imposed at a given time, but takes place gradually, in a way depending of the evolution of the system, such a strategy has also the interest of being in some sense self-adapting to the system's actual state. In Fig. 8, the scaling behaviors relative to the "play-smart" strategy are reported. Both the height and time of the maximum follow the scaling of the "play-last" strategy: t max ∼ N α and N max w ∼ N γ with α ≈ γ ≈ 1.3. The convergence time, on the other hand, scales as a superposition of two power laws: t conv ∼ aN α + bN δ with α ≈ 1.3, δ ≈ 1.0. Thus, the global behavior determined by the "play-smart" modification is indeed less demanding in terms of both memory and time. In particular, while the lowering of the peak height yields in fact a slower convergence for the "play-last" strategy, the progressive self-driven change in strategy allows to fasten the convergence further than for the "play-first" strategy.

Finally, in order to have an immediate feeling of what different playing word selection strategies imply, we report in Figs. 9 and 10 the success rate S(t) and the total number of words, N w (t) relative to the four strategies described previously, for two different sizes. The TABLE I: Scaling with the system size N of the maximum number of words (memory) and time of convergence. Networks, thanks to the small-world property and the finite connectivity, ensure a trade-off between the fast convergence of mean-field topology and the small memory requirements of lattices.

"play-first" and "play-smart" curves exhibit the same "Sshaped" behavior for S(t) as in the case of the stochastic model, while the "play-last" rule affects qualitatively the way in which the final state is reached. Indeed, in this case the transition between the initial disordered state and the final ordered one is more continuous (see the inset in the top figure). Moreover, Fig. 10 illustrates that the choice of the strategy has substantial quantitative consequences for both necessary memory and time needed to reach convergence, even if the changes in scaling behavior could at first appear rather limited (from N 1.5 to N 1.3 ). In particular, the "play-smart" strategy, which adapts itself to the state of the system, leads to a drastic reduction of the memory and time costs and thus to a dramatic increase in efficiency.

VIII. CONCLUSIONS

The Naming Game describes how a shared communication system may emerge in a population of agents. Here we have reviewed some of the main aspects of the model, focusing in particular on the role of different interaction patterns on the global quantities of the model. We have also investigated alternative strategies that speed up the establishment of the final consensus state.

Different underlying topologies have the main effect of changing the scaling with the system size N of two crucial quantities of the model: the time required to reach a consensus, t conv , and the maximum memory demanded to the agents during the process, N max w . In particular, our analysis has pointed out that (i) finite connectivity implies finite memory requirements for the agents and (ii) the small-world property ensures a fast convergence (see Table I). Indeed, these properties affect the mechanism itself that leads to the final state. On the complete graph (mean-field case), the process starts with an initial spreading of words (linear with time) throughout the system followed by a longer period (O(N 3/2 )) in which words are exchanged among the agents. After the peak of the total number of words (whose height is O(N 3/2 )), the total number of words decreases till a final rapid convergence process leads the population to the adsorbing configuration. On low-dimensional lattices, on the other hand, the model converges very slowly (t conv /N ∼ N 2/d , where d is the dimension of the lattice), but the maxi-mum number of words scales as N , meaning that single agents require a finite memory that does not depends on N . The reason is related to the fast formation of many different local clusters of agents with the same unique word, that grow through a coarsening process. A tradeoff between lattices and fully connected graphs is offered by small-world networks, in which a finite memory per agent comes with a fast convergence N 1.4±0.1 ensured by the spreading of different words across the shortcuts that connect otherwise far-away regions of the graph. This picture has been confirmed also by extensive numerical studies on different kinds of complex networks, where in addition we have addressed also the role of other important parameters (e.g. average connectivity, clustering coefficient).

Finally, we have investigated whether the efficiency of the Naming Game can be improved acting on the microscopic rules of the model. Focusing on the mean-field case, we have checked the procedure describing which word has to be played by the speaker. We have substituted the usual random extraction with extremely simple deterministic rules, which exploit only the information regarding the time at which different words are inserted in the inventory of the agent. In particular, each agent can identify two special words, namely the oldest and the newest. We have shown that simple deterministic rules can capitalize on, and somehow tune, the correla-tion among the inventories, thus increasing the efficiency of the Naming Game in terms of both the individual use of memory and the convergence time.

In conclusion, we have shown that different topologies have major effects on the global dynamics of the Naming Game. We have also pointed out that simple modifications of the original rules can give rise to higher performances that allow the population to save memory and time to reach the final state. We believe that this findings could be extremely relevant as far as application are concerned, and in future more efforts should be made to improve the model further. As we have seen, however, each new rule should be tested on different topologies in order to understand its potential consequences in real world applications.

FIG. 2 :

 2 FIG.2: Evolution of the total number of words Nw (top), of the number of different words N d (center), and of the average success rate S(t) (bottom), for a mean-field system (black circles) and low-dimensional lattices (1D, red squares and 2D, blue triangles) with N = 1024 agents, averaged over 10 3 realizations. The inset in the top graph shows the very slow convergence in low-dimensional lattices.

FIG. 3 :FIG. 4 :

 34 FIG.3:A naive representation of clusters growth in the smallworld model of Watts and Strogatz. A cluster (in red) starts to expand locally by coarsening dynamics like in dimension one. When the size of the cluster is of the order of the average distance between shortcuts, long-range interactions take place. The effect of these long-range interactions is that of boosting up the dynamics.

3 FIG. 5 :

 35 FIG.5: ER random graph (left) and BA scale-free network (right) with k = 4 and sizes N = 10 3 , 10 4 , 5.10 4 . Top: evolution of the average memory per agent Nw/N versus rescaled time t/N . For increasing sizes a plateau develops in the reorganization phase preceding the convergence. The height of the peak and of the plateau collapse in this plot, showing that the total memory used scales with N . Bottom: evolution of the number of different words N d in the system. (N d -1)/N is plotted in order to emphasize the convergence to the consensus with N d = 1. A steady decrease is observed even if the memory Nw displays a plateau. The mean-field (MF) case is also shown (for N = 10 3 ) for comparison.

FIG. 6 :

 6 FIG.6: ER networks (left) and BA networks (right) with N = 10 4 agents and average degree k = 4, 8, 16. The increase of average degree leads to a larger memory used (Nw, top) but a faster convergence. The maximum in the number of different words is not affected by the change in the average degree (bottom).

FIG. 7 :

 7 FIG. 7: Number of clusters N cl and normalized average cluster size s /N vs. time for ER networks with N = 10 4 , k = 4 (circles), k = 8 (squares), k = 16 (crosses).

3 FIG. 8 :

 38 FIG.8: Play-smart strategy -scaling with the population size N . Top -For the time of the peak tmax ∼ N α , α ≈ 1.3, while for the convergence time we have tconv ∼ aN α + bN δ with δ ≈ 1.3, δ ≈ 1.0. Bottom -the maximum number of words scales as N max w ∼ N γ with γ ≈ 1.3. The "play-smart" rule gives rise to a more performing process, from the point of view of both convergence time and memory needed.

IX. ACKNOWLEDGMENTS

A. Baronchelli and V. L. are partially supported by the EU under contract IST-1940 "Embodied and Communicating Agents" (ECAgents). A. Barrat and L.D. are partially supported by the EU under contract 001907 "Dynamically Evolving, Large Scale Information Systems" (DELIS).