Thierry Colin

Marie-Christine Durrieu

Julie Joie

Yifeng Lei

Youcef Mammeri

Clair Poignard

Olivier Saut

Modelling of the migration of endothelial cells on bioactive micropatterned polymers

Keywords: Cell migration modelling, Chemotaxis, Non-linear P.D

In this paper a macroscopic model describing endothelial cells migration on bioactive micropatterned polymers is presented. It is based on a system of partial differential equations of Patlak-Keller-Segel type that describes the evolution of the cell densities. The model is studied mathematically and numerically. We prove existence and uniqueness results of the solution to the differential system and also that fondamental physical properties such as mass conservation, positivity and boundedness of the solution are satisfied. The numerical study allows us to show that the model behaves in good agreement with the experiments.

Modelling of the migration of endothelial cells on bioactive micropatterned polymers

Résumé : Dans cet article nous présentons un modèle macroscopique décrivant la migration de cellules endothéliales sur un micropattern de polymers bioactifs. Ce modèle est basé sur un système d'équations aux dérivées partielles du type Patlak-Keller-Segel. Les propriétés mathématiques et numériques du modèles sont présentées. Nous démontrons des résultats d'existence et d'uncité ainsi que les propriétés physiques, telles que la conservation de la masse, la positivité et le caractère borné de la solution. L'étude numérique nous permet de montrer que le modèle corrobore les résultats expérimentaux.

Mots-clés : Cell migration modelling, Chemotaxis, Non-linear P.D.E

Introduction

Tissue engineering is the use of combination of cells, engineering, materials, and suitable biochemical factors to improve or replace biological functions [START_REF] Nerem | Tissue engineering: the hope, the hype, and the future[END_REF]. Tissue engineering has being quickly developing since the 1990s [START_REF] Nerem | Tissue engineering: the hope, the hype, and the future[END_REF]. However, the major roadblock for engineering large tissues is the lack of functional microvasculature networks, which provide nutrients and oxygen for tissue survival and remove the waste product from metabolism [START_REF] Jain | Engineering vascularized tissue[END_REF]. The lack of vascularization has hampered the survival of engineered tissues after implantation [START_REF] Jain | Engineering vascularized tissue[END_REF]. Researchers rely on the increasing knowledge of angiogenic and vasculogenic processes to stimulate vascular network formation [START_REF] Santos | Vascularization in bone tissue engineering: physiology, current strategies, major hurdles and future challenges[END_REF][START_REF] Phelps | Engineering more than a cell: vascularization strategies in tissue engineering[END_REF]. This complex process of new blood vessel formation is orchestrated by the interaction between endothelial cells (ECs) and their neighboring mural cells via a complex network of intracellular signaling mechanisms [START_REF] Otrock | Understanding the biology of angiogenesis: review of the most important molecular mechanisms[END_REF][START_REF] Jain | Molecular regulation of vessel maturation[END_REF]. Ever since the introduction of the in vitro model of angiogenesis [START_REF] Folkman | Angiogenesis in vitro[END_REF], there has been great research interest to understand the intricate process of tube formation. Although many efforts have been made, the mechanism associated with angiogenesis and vascularization is poorly understood. A deeper comprehension of cellsbiomaterials interaction is required for basic understanding of angiogenesis and vascularization in tissue engineering [START_REF] Carmeliet | Common mechanisms of nerve and blood vessel wiring[END_REF].

One strategy in developing clinical implants is therefore to use bioactive materials which can either elicit a regenerative response at the site of damage in vivo or be used to grow tissue in vitro for subsequent implantation [START_REF] Anselme | The interaction of cells and bacteria with surfaces structured at the nanometre scale[END_REF][START_REF] Lei | Mimicking angiogenesis by SVVYGLR peptide micropatterning[END_REF]. Different bioactive ligands have been used to study their effects on cell functions for a better understanding of vascularization [START_REF] Phelps | Engineering more than a cell: vascularization strategies in tissue engineering[END_REF]. In the aim at promoting angiogenesis in the case of tissue engineering or inhibiting angiogenesis in the case of cancer, it is now important to understand the mechanisms that regulate lumen formation. Successful micropatterning of cells is becoming a key component of this field [START_REF] Ito | Surface micropatterning to regulate cell functions[END_REF]. Researchers are now interested in the behavior of cells on substrates that have been patterned by micro/ nanofabrication [START_REF] Folch | Microengineering of cellular interactions[END_REF][START_REF] Nicolau | Patterning neuronal and glia cells on light-assisted functionalized photoresists[END_REF]. It is known that cell positioning and physiology can be controlled by the substrate on which the cells adhere [START_REF] Chen | Geometric control of cell life and death[END_REF]. In our study, using cell adhesion peptides micropatterned onto material, we observed a tube-like formation in comparison to virgin material or homogeneously grafted materials [START_REF] Lei | Micropatterning of polyethylene terephthalate with RGD peptides for induction of endothelial cell alignment and morphogenesis[END_REF][START_REF] Lei | Mimicking angiogenesis by SVVYGLR peptide micropatterning[END_REF].

Experimental studies using micropatterned substrates revealed that the cell migration is governed by the geometry of the patterns (size of patterns and distance between patterns). Endothelial cells so cultured form extensive cell-cell interactions. Accumulation of endothelial cells junctions implies that some cells form tube-like structures. The goal of the present paper is to provide a model that describes such results.

Adhesive areas are composed of cell adhesion peptides or growth factors peptides that make the cells adhere. These areas are surrounded by non-adhesive areas [START_REF] Lei | Micropatterning of polyethylene terephthalate with RGD peptides for induction of endothelial cell alignment and morphogenesis[END_REF]. We assume (and this is actually confirmed by experiments) that active principles (cell adhesion peptides or growth factors) are not diffusing spatially. Therefore, endothelial cells outside adhesive areas have no mean to "feel" directly these areas. They find out these areas indirectly. We do not consider the influence of nutrients and assume that cells obtain enough nutrients from the material (due to grafted active principles onto material). Endothelial cells are seeded onto micropatterned bioactive materials during several hours, then they are washed out. Only the adhered endothelial cells remain on material. The initial density of cells is around 40 000 cells per cm 2 . At the beginning of the experiments, during the migration phase we observe that cells have a random motility and stop on adhesive areas. Moreover the attraction of endothelial cells on adhesive areas seems to be higher than the one of cells outside these areas. Experiments show that endothelial cells are grouping together along the micropatterns. According to pattern size, endothelial cells line their cytoskeleton to adjust it with the adhesive area. One can also notice that, with micropatterns of 10 µm thin strips, tubes containing a central lumen may appear [START_REF] Lei | Mimicking angiogenesis by SVVYGLR peptide micropatterning[END_REF][START_REF] Dike | Geometric control of switching between growth, apoptosis, and differentiation during angiogenesis using micropatterned substrates[END_REF]. In Inria other words, blood vessels are created from an initial random density of endothelial cells. Such phenomenon is not observed with larger strips [START_REF] Dike | Geometric control of switching between growth, apoptosis, and differentiation during angiogenesis using micropatterned substrates[END_REF][START_REF] Kamei | Endothelial tubes assemble from intracellular vacuoles in vivo[END_REF][START_REF] Lubarsky | Tube morphogenesis: making and shaping biological tubes[END_REF].

To illustrate these experiments, we present in Fig. 1 pictures of the micropatterned bioactive materials at the end of the migration step. Two different micropatterns are considered: on Fig. 1(a) thin adhesive areas (bioactive pattern size: 10 µm and distance between patterns: 100 µm) have beed used, whereas Fig. 1(b) shows the end of the migration on large strips (bioactive pattern size: 300 µm and distance between patterns: 100 µm). For a detailed description of such experiences, one may refer to [START_REF] Lei | Micropatterning of polyethylene terephthalate with RGD peptides for induction of endothelial cell alignment and morphogenesis[END_REF][START_REF] Lei | Mimicking angiogenesis by SVVYGLR peptide micropatterning[END_REF]. µm (B) stripes of SVVYGLR peptides) [START_REF] Lei | Micropatterning of polyethylene terephthalate with RGD peptides for induction of endothelial cell alignment and morphogenesis[END_REF]. The distance between bioactive patterns is 100 µm.

We observed that for the largest adhesive areas the adhered cells density is smaller than for thin strips. Therefore, the geometries of the micropatterns play a crucial role in the endothelial cells migration and then in the formation of new vessels.

In this paper, we are interested in understanding how these patterns (size and distance between microfeatures) do influence endothelial cells migration. The model we present here is a continuous Patlak-Keller-Segel type model [START_REF] Anderson | Continuous and discrete mathematical models of tumorinduced angiogenesis[END_REF][START_REF] Hillen | A user's guide to PDE models for chemotaxis[END_REF][START_REF] Keller | Traveling band of chemotactic bacteria: a theoretical analysis[END_REF][START_REF] Patlak | Random walk with persistence and external bias[END_REF]. The chemotaxis term takes the cellcell interactions into account instead of the cell-chemical attractant interactions. We show that this new model based on a system of coupled partial differential equations conserves the mass and existence and uniqueness results of weak solution hold. We also provide numerical results in accordance with the experiments, which ensures the validity of our model. Moreover, these simulations allow us to obtain informations on the influence of the geometry and of the initial concentration of cells on the migration of the cells.

The outline of the article is the following. In section 2, we describe the mathematical model and we state the main result of global existence and uniqueness of the weak solution to the P.D.E system. Section 3 is devoted to the proof of the main theorem. We then provide numerical results in section 4 in order to compare the simulations to the experiments.

Description of the model and main result

In this section, we describe the Patlak-Keller-Segel type continuous model we study throughout the paper. The model is composed of a diffusion term coupled with a reaction term, which describes the effect of the chemoattractants, which themselves statisfy a diffusion equation.

Various continuous models of Patlak and Keller-Segel type have been used to describe cell motility. Typically, the governing equations of these models are written in the following general form: on a domain

Ω ⊂ R n ∂ t u = ∇. (D 1 (u, v)∇u -χ(u, v)u∇v) + f (u, v) on Ω ∂ t v = ∇. (D 2 (u, v)∇v) + g(u, v) -h(u, v)v on Ω,
where u denotes the cells density, v is the chemical signal concentration. The diffusive terms take the random motility of the cells into account, whereas the advection describes the influence of the chemical signal on the motion of the cells. The two corresponding diffusion parameters are denoted by D 1 and D 2 while χ is the chemotaxis coefficient. The function f describes the growth and the death of cells, whereas g and h are the production and degradation of the chemotaxic signal. These equations have been theoretically studied [START_REF] Biler | Existence and nonexistence of solutions for a model of gravitational interaction of particles[END_REF][START_REF] Blanchet | Two dimensional Keller-Segel model : optimal critical mass and qualitative properties of the solution[END_REF][START_REF] Dolbeault | Optimal critical mass in the two dimensional Keller-Segel model in R 2[END_REF][START_REF] Gajewski | Global behavior of a reaction-diffusion system modelling chemotaxis[END_REF][START_REF] Horstmann | The nonsymmetric case of the Keller-Segel model in chemotaxis: some recent results[END_REF][START_REF] Senba | Chemotactic collapse in a parabolic-elliptic system of mathematical biology[END_REF] for several years. Based on this extensive literature we provide a slightly modified model to describe the cell migration on bioactive micropatterns.

Statement of the equations

According to the experiments, the behavior of the cells is drastically different on the adhesive area and outside this area. Actually, outside the adhesive strip the cells seem to attract each other (probably thanks to the chemoattractant they produce) and also diffuse in the domain, but as soon as they reach the adhesive strip they seem stuck on the strip and then they diffuse only on this area, ignoring the outer cells. Moreover it seems that the production of chemoattractant of the cells located on the adhesive strip is bigger than outside. Since there is no clear understanding of the way that endothelial cells communicate, we chose to consider the chemotaxis term as the attraction between endothelial cells (and we do not consider any gradient of concentration of the chemoattractant).

Based on these assumptions, we derive the following model. Consider a domain Ω composed by the adhesive areas, denoted Ω, and the non-adhesive areas Ω\ Ω, both domains being bounded domains with smooth boundary.

Two different types of endothelial cells are considered. We denote by u 1 (t, x, y) the density of endothelial cells, at any point (x, y) and time t, that can freely move (i.e. they have yet to move over adhesion proteins). Cells that are adhering on the substrate are tracked through their density u 2 . The function v represents the density of the chemoattractant. The equations governing the endothelial cells migration are given for t > 0 by

∂ t u 1 = d 1 ∆u 1 -λ1 Ω u 1 (1 -u 2 ) -∇. (χ(u 1 , v)u 1 ∇v) , in Ω, (1a)
∂ t u 2 = d 2 ∆u 2 + λ1 Ω u 1 (1 -u 2 ), in Ω, (1b)
∂ t v = ∆v -ηv + γ 1 u 1 + γ 2 u 2 , in Ω, (1c)
with the homogeneous boundary conditions on ∂Ω and ∂ Ω:

∂ n u 1 | ∂Ω = 0, ∂ n u 2 | ∂ Ω = 0, ∂ n v| ∂Ω = 0, (1d)
and with the initial conditions (u 0 1 , u 0 2 , 0)

u 1 | t=0 = u 0 1 , u 2 | t=0 = u 0 2 , v| t=0 = 0. ( 1e
)
We then denote by u the sum

u(t, x) = u 1 (t, x) + u 2 (t, x), t ≥ 0, x ∈ Ω,

Inria

where u 2 is extended by 0 in Ω \ Ω.

The parameters d 1 , d 2 , η, γ 1 , γ 2 and λ are strictly positive and they will be fitted by the experiments in a forthcoming work, but we consider here that they are given constants. The coefficients d 1 and d 2 denote the diffusion coefficients of the cells u 1 and u 2 respectively. The coefficient η > 0 is the self-degradation rate of the chemoattrant produced by the cells, while the coefficients γ 1 and γ 2 are the coefficients of the production of the chemoattractant respectively for the cell u 1 and u 2 . The parameter λ is the rate for the cell u 1 to become u 2 as soon as u 1 lies in the micropatterns. The two first equations describe the migration of the cells on Ω. Outside the strips, endothelial cells diffuse and attract the neighboring cells via the chemotaxis sensitivity function :

χ(u 1 , v) = χ 0 v 1 + |v| (1 -u 1 ), with χ 0 > 0.
Here above, χ 0 is a chemotaxis parameter, and the term (1 -u 1 ) is settled to prevent the overcrowding of the cells u 1 . Endothelial cells once they reach the adhesive area Ω are captured and then diffuse only in the strip. This is handled by the penalty term -λ1 Ω u 1 (1 -u 2 ). Cells on the strips still have a random motility and their concentration grows up as the term λ1 Ω u 1 (1-u 2 ), where 1 -u 2 prevents the blow-up of u 2 in equation (1b). The third equation describes the production of the chemoattractant by the cells. Since the cells on the strip seem to be more attractive, we suppose that the production coefficients checks 0 < γ 1 < γ 2 . We also add a degradation coefficient η > 0 describing the metabolization of the chemoattractant.

Main theoretical result

We have the following theorem which is a straightforward consequence of the results of the next section 3:

Theorem 2.1. Let d 1 , d 2 , η, γ 1 , γ 2 and λ be strictly positive constants. Suppose that the initial data

(u 0 1 , u 0 2 ) ∈ L ∞ (Ω) × L ∞ ( Ω) are such that ∀x ∈ Ω, 0 ≤ u 0 1 (x) ≤ 1, ∀x ∈ Ω, 0 ≤ u 0 2 (x) ≤ 1.
There exists a unique weak solution (u 1 , u 2 , v) to problem (1) such that

(u 1 , u 2 , v) ∈ L ∞ [0, +∞); L ∞ (Ω) × L ∞ [0, ∞); L ∞ ( Ω) × L ∞ [0, ∞); L ∞ (Ω) ,
and for almost any t > 0

0 ≤ u 1 (t, •) ≤ 1, 0 ≤ u 2 (t, •) ≤ 1, and 0 ≤ v(t, •) ≤ 1 η (γ 1 + γ 2 ) .
The next section is devoted to prove this theorem. The proof is based on Gaussian upper bounds for heat kernels [START_REF] Ouhabaz | Analysis of heat equations on domains[END_REF]- [START_REF] Wang | Gradient Estimate on Convex Domains and Application[END_REF].

Theoretical study of the model

In this section, we study the mathematical properties of the model. Throughout this section we suppose that Ω and Ω are smooth domains of R 2 . We remind that d 1 , d 2 and η are strictly positive coefficients.

Kernels of the operators

The aim of this paragraph is to define and provide estimates of the kernels of the operators ∂ t -∆ + η and ∂ t -d 1 ∆ on Ω and for the kernel of ∂ t -d 2 ∆ on Ω, with homogeneous Neumann conditions respectively on ∂Ω and ∂ Ω.

∀(x, y) ∈ Ω, lim t→0 + B(t, x, y) = δ y (x),
and for any (t, y) ∈ (0, ∞) × Ω,

∂ t B(t, x, y) = ∆B(t, x, y) -ηB(t, x, y), ∀x ∈ Ω, ∂ n B(t, x ∂Ω , y) = 0, ∀x ∂Ω ∈ ∂Ω, (2a)
for B, while G is given by

∀(x, y) ∈ Ω, lim t→0 + G(t, x, y) = δ y (x), (3a)
and for any (t, y) ∈ (0, ∞) × Ω,

∂ t G(t, x, y) = d 1 ∆G(t, x, y), ∀x ∈ Ω, ∂ n G(t, x ∂Ω , y) = 0, ∀x ∂Ω ∈ ∂Ω, (3b)
and G is the solution to

∀(x, y) ∈ Ω, lim t→0 + G(t, x, y) = δ y (x), (4a)
and for any (t, y) ∈ (0, ∞) × Ω,

∂ t G(t, x, y) = d 2 ∆ G(t, x, y), ∀x ∈ Ω, ∂ n G(t, x ∂ Ω , y) = 0, ∀x ∂ Ω ∈ ∂ Ω. (4b)
Note that the above kernels are symmetric in their second and third variables.

Proposition 1. For any y ∈ Ω (respectively for any y ∈ Ω), we have the following estimates for positive constants C Ω and C Ω , which depend on the domain Ω and Ω respectively:

G(t, •, y) L 1 (Ω) ≤ C Ω , (5a)
G(t, •, y) L 1 (Ω) ≤ C Ω , (5b)
B(t, •, y) L 1 (Ω) ≤ C Ω , (5c)
and gradient estimates hold too:

∇ x G(t, •, y) L 1 (Ω) ≤ C Ω max(1, t -3/4 ), (6a)
∇ x G(t, •, y) L 1 (Ω) ≤ C Ω max(1, t -3/4 ), (6b)
∇ x B(t, •, y) L 1 (Ω) ≤ C Ω max(1, t -3/4 ). ( 6c
)

Inria

In addition, due to the boundedness of Ω, we also have

∇ y G(t, •, y) L 1 (Ω) ≤ C Ω max(1, t -3/4 ), (7a)
∇ y G(t, •, y)

L 1 (Ω) ≤ C Ω max(1, t -3/4 ), (7b)
∇ y B(t, •, y) L 1 (Ω) ≤ C Ω max(1, t -3/4 ). (7c)
Proof. Obviously the coefficients of diffusion d 1 and d 2 , since they are strictly positive constants do not play a crucial role, and can be supposed to be 1, after an appropriate rescaling of the time variable t. Moreover it is sufficient to prove the above estimates for the heat kernel G, since B = e -ηt G.

For t ≥ 1, estimates (3.2)-(3.3) of [START_REF] Wang | Gradient Estimate on Convex Domains and Application[END_REF] straightforwardly provide the result. Suppose that 0 < t ≤ 1. Estimates (5a) easily come from Theorem 6.10 pp 171 of [START_REF] Ouhabaz | Analysis of heat equations on domains[END_REF], since for any x ∈ Ω,

0 ≤ 1 √ t Ω e -|x-y| 2 /t dy ≤ 2π.
Estimates (6a) are consequences of the section 6.6 entitled Weighted Gradient Estimates and in particular Theorem 6.19 p 185 [START_REF] Ouhabaz | Analysis of heat equations on domains[END_REF]. Actually by Cauchy-Schwarz inequality

∇ x G(t, •, y) 2 L 1 (Ω) ≤ Ω |∇ x G(t, •, y)| 2 e 2β|x-y| 2 /t dy Ω e -2β|x-y| 2 /t dy ≤ Ct -2 e ct Ω e -2β|x-y| 2 /t dy, ≤ 2πCt -3/2 e ct ≤ Ct -3/2 ,
hence the estimates (6a). Now let φ ∈ L ∞ (Ω), by estimates (6a) and since the measure |Ω| of Ω is bounded we infer

Ω |φ(y)| Ω |∇ y G(t, x, y)| dx dy = Ω×Ω |φ(y)| |∇ y G(t, x, y)| dy dx, ≤ φ L ∞ (Ω) Ω ∇ y G(t, x, •) L 1 (Ω) dx, ≤ |Ω|Ct -3/4 φ L ∞ (Ω) ,
hence estimates (7a), which ends the proof of the proposition.

Remark 1. The above estimates are probably not optimal, since for the half-plane the heat kernel write: 4t) , where y c = (y 1 , -y 2 ), and therefore the power t -3/4 has to be replaced by t -1/2 similarly to the heat kernel of the whole plane R 2 . However these results are sufficient to prove existence and uniqueness of the solution to problem (1).

G(t, x, y) = 1 4πt e |x-y| 2 /(4t) + e |x-y c | 2 /(
Corollary 1. In particular, for T > 0, and for any φ ∈ L ∞ ([0, T ]; L ∞ (Ω)) the solution to the following problem:

∂ t u = ∆u -ηu + φ(t, •), in Ω, ∂ n u| ∂Ω = 0, u| t=0 = 0, (8)
satisfies almost everywhere in (0, T ) × Ω:

|∇u(t, •)| ≤ C Ω t 1/4 sup s∈(0,T ) φ(s, •) L ∞ (Ω) . (9)
Proof. Since

u(t, •) = t 0 Ω B(t -s, •, y)φ(s, y) dy
and thus

|∇u(t, x)| ≤ t 0 Ω |∇ x B(t -s, x, y)φ(s, y)| dy ≤ sup s∈(0,T ) φ(s, •) L ∞ (Ω) t 0 ∇ x B(t -s, x, •) L 1 (Ω) ,
hence inequality (9).

Local existence

Using the above appropriate kernels, we deduce that a weak solution to problem (1) writes:

u 1 (t, x) = Ω G(t, x, y)u 0 1 (y) dy -λ t 0 Ω
G(t -s, x, y)u 1 (s, y)(1 -u 2 )(s, y) dy ds

+ t 0 Ω u 1 (s, y)χ(u 1 , v)(s, y)∇ y G(t -s, x, y) • ∇v(s, y) dy ds, (10a)
u 2 (t, x) = Ω G(t, x, y)u 0 2 (y) dy + λ t 0 Ω G(t -s, x, y)u 1 (s, y)(1 -u 2 )(s, y) dy ds, (10b)
v(t, x) = t 0 Ω B(t -s, x, y) (γ 1 u 1 (s, y) + γ 2 u 2 (s, y)) dy ds. (10c)
In this paragraph, we aim at proving a local-existence result.

Definition of the appropriate functional space X T M

Let M be a strictly positive constant, and let T > 0 that will be chosen later. We define the functional space X T M as

X T M = Λ ∈ L ∞ ([0, T ]; L ∞ (Ω)) : sup t∈[0,T ] Λ(t, •) L ∞ (Ω) ≤ M .
Let L be the linear operator defined on X T M × X T M by

L : (ν 1 , ν 2 ) → t 0 Ω B(t -s, •, y) (γ 1 ν 1 (s, y) + γ 2 ν 2 (s, y)) dy ds.

Inria

Using estimates (5a)-(6a) we infer that for any (ν 1 , ν 2 ) ∈ X T M × X T M :

L(ν 1 , ν 2 )(t, •) L ∞ (Ω) ≤ C Ω (γ 1 + γ 2 )M, ∇L(ν 1 , ν 2 )(t, •) L ∞ (Ω) ≤ C Ω (γ 1 + γ 2 )M t 1/4 . ( 11
)
Define now the operator T on X T M × X T M by

T : (ν 1 , ν 2 ) → T 1 -T 2 (ν 1 , ν 2 ) , T 2 (ν 1 , ν 2 ) ,
where T 1 is the operator defined on X M × X T M by

T 1 (ν 1 , ν 2 ) = t 0 Ω ν 1 χ (ν 1 , L(ν 1 , ν 2 )) ∇ y G(t -s, •, y) • ∇ y L(ν 1 , ν 2 ) dy ds,
and T 2 is defined by

T 2 (ν 1 , ν 2 ) = λ t 0 Ω G(t -s, x, y)ν 1 (s, y)(1 -ν 2 )(s, y) dy ds.
Remark 2. Proving that T is a contraction mapping from X T M ×X T M onto itself for small enough time T will then ensure the local existence of the weak solution given by [START_REF] Folch | Microengineering of cellular interactions[END_REF] to problem (1).

Contraction mappings

Proposition 2. The operators T is a contraction mapping from X T M × X T M onto itself for T small enough.

Proof. The proof is based on the properties of the kernel B, G and G given by Proposition 1. Thanks to estimate (5a) we deduce for any (ν 1 , ν 2 ) ∈ X T M × X T M :

T 2 (ν 1 , ν 2 ) L ∞ ( Ω) ≤ C Ω λM (1 + M )T,
hence for T small enough T 2 maps X T M × X T M onto X T M . Moreover using inequality:

|ν 1 (1 -ν 2 ) -µ 1 (1 -µ 2 )| ≤ (1 + |ν 2 |)|ν 1 -µ 1 | + |ν 2 ||µ 1 -µ 2 |,
we infer for T small enough the operator T 2 is a contraction mapping from

X T M × X T M onto X T M . Prove now that T 1 is a contraction mapping from X T M × X T M onto X T M . First observe that for any s ∈ R, |s|/(1 + |s|) ≤ 1 hence for any ν 1 ∈ X T M , for any s ∈ R, χ 1 (ν 1 , s)(t, •) L ∞ (Ω) ≤ χ 0 (1 + M ), for almost any t ∈ (0, T ), hence for any (ν 1 , ν 2 ) ∈ X T M × X T M χ 1 (ν 1 , L(ν 1 , ν 2 ))(t, •) L ∞ (Ω) ≤ χ 0 (1 + M ), for almost any t ∈ (0, T ),
and thanks to estimates (6a)-( 11)

|T 1 (ν 1 , ν 2 )(t, •)| ≤ C Ω (γ 1 + γ 2 ) χ 0 (1 + M )M 2 √ T .
This implies that for T small enough T 1 maps X T M × X T M onto X T M . In addition, observe that for two couples (ν 1 , ν 2 ) and (µ 1 , µ 2 ) belonging to X T M × X T M we have

T 1 (ν 1 , ν 2 ) -T 1 (µ 1 , µ 2 ) = t 0 Ω (ν 1 -µ 1 )χ ν1,ν2 ∇ y G(t -s, •, y) • ∇ y L(ν 1 , ν 2 ) dy ds, + t 0 Ω µ 1 χ ν1,ν2 ∇ y G(t -s, •, y) • ∇ y L(ν 1 -µ 1 , ν 2 -µ 2 ) dy ds, + t 0 Ω µ 1 (χ ν1,ν2 -χ µ1,µ2 ) ∇ y G(t -s, •, y) • ∇ y L(µ 1 , µ 2 ) dy ds,
where to simplify notations we have denoted by χ ν1,ν2 the function

χ ν1,ν2 = χ(ν 1 , L(ν 1 , ν 2 )),
and similarly for χ µ1,µ2 . According to estimate (6a), and thanks to the definition of L, we infer

∇ y L(ν 1 -µ 1 , ν 2 -µ 2 ) ≤ C Ω (γ 1 + γ 2 )t 1/4 ν 1 -µ 1 L ∞ (Ω) + ν 2 -µ 2 L ∞ (Ω) .
Moreover, observing that

χ ν1,ν2 -χ µ1,µ2 = χ 0 L(ν 1 , ν 2 ) 1 + |L(ν 1 , ν 2 )| (µ 1 -ν 1 ) + χ 0 (1 -µ 1 ) L(ν 1 , ν 2 ) 1 + |L(ν 1 , ν 2 )| - L(µ 1 , µ 2 ) 1 + |L(µ 1 , µ 2 )| ,
we deduce from estimates (5a)-(6a)-(7a) and ( 11) that there exists a constant C > 0 which depends on M , and on the parameters χ 0 , γ 1 , γ 2 , λ such that

T 1 (ν 1 , ν 2 ) -T 1 (µ 1 , µ 2 ) L ∞ (Ω) ≤ C √ T ,
which ensures the strict contractility of T 1 for T small enough, and therefore T is a strict contraction from X T M × X T M onto itself.

The Picard fixed point theorem straightforwardly implies the following theorem of existence and uniqueness for small time.

Theorem 3.2. Let (u 0 1 , u 0 2 ) ∈ L ∞ (Ω) × L ∞ ( Ω).
Then, for T small enough, there exists a unique weak solution (u 1 , u 2 , v) to (1) such that

(u 1 , u 2 , v) ∈ L ∞ ([0, T ]; L ∞ (Ω)) × L ∞ [0, T ]; L ∞ ( Ω) × L ∞ [0, T ]; L ∞ (Ω) .

Mass conservation and global existence

We first observe that the total mass of cells is conserved.

Proposition 3. Let (u 0 1 , u 0 2 ) ∈ L ∞ (Ω) × L ∞ ( Ω)
and let T small enough so that a weak solution (10) to (1) exists. Then, for any t ∈ [0, T ] we infer

Ω u(t, x) dx = Ω (u 1 + 1 Ω u 2 )(t, x) dx = Ω u 0 1 dx + Ω u 0 2 dx

Inria

Proof. Actually integrating (1a) and (1b) respectively and summing the integrands imply, since

∂ n u 1 | ∂Ω , ∂ n u 2 | ∂ Ω and ∂ n v| ∂Ω vanish ∂ t Ω u(x) dx = ∂ t Ω (u 1 (x) + u 2 (x)) dx = 0.
We now show that if u 0 1 and u 0 2 are positive and bounded by 1 then u 1 and u 2 stay positive and bounded by 1 on [0, T ].

Proposition 4. Let (u 0 1 , u 0 2 ) ∈ L ∞ (Ω) × L ∞ ( Ω)
and let T small enough so that a weak solution given by [START_REF] Folch | Microengineering of cellular interactions[END_REF] to problem (1) exists. If (u 0 1 , u 0 2 ) are such that

0 ≤ u 0 1 ≤ 1, 0 ≤ u 0 2 ≤ 1, then for almost any t ∈ [0, T ] 0 ≤ u 1 (t, x) ≤ 1, 0 ≤ u 2 (t, x) ≤ 1.
In addition,

0 ≤ v(t, x) ≤ 1 η (γ 1 + γ 2 ) , for x ∈ Ω.
Therefore, the weak solution (10) exists for almost any t ∈ (0, +∞).

Proof. First observe that if u 1 is positive then since u 0 2 is positive the function u 2 is positive almost everywhere. Actually multiplying (1b) by u - 2 = max(0, -u 2 ) and integrating by parts implies

1 2 ∂ t u - 2 (t, •) 2 L 2 ( Ω) = -d 2 ∇u - 2 2 L 2 ( Ω) -λ Ω u 1 u - 2 + λu 1 (u - 2 ) 2 ≤ λM (u - 2 ) 2 L 2 ( Ω ,
hence u - 2 equal zero by Gronwall's lemma. Similarly, if u 1 is positive, then since u 2 is therefore also positive we infer that v is positive by multiplying (1c) by v -and integrating by parts.

Prove now that u - 1 = max(0, -u 1 ) vanishes too. Multiply (1a) by u - 1 and integrate by parts to obtain for almost any t ∈ [0, T ]:

1 2 ∂ t u - 1 (t, •) 2 L 2 (Ω) ≤ -d 1 ∇u - 1 (t, •) 2 L 2 (Ω) + λ(1 + M ) u - 1 (t, •) 2 L 2 ( Ω) + χ 0 Ω u - 1 v 1 + |v| (1 -u - 1 )∇v∇u - 1 dx .
Moreover applying estimate [START_REF] Eymard | Finite Volume Methods, Handbook of Numerical Analysis[END_REF] to v implies

Ω u - 1 (1 -(u 1 + u 2 ))∇v∇u - 1 dx ≤ C Ω t 1/4 γ 1 u 1 (t, •) + γ 2 1 Ω u 2 (t•) L ∞ (Ω) × Ω u - 1 (t, x)|1 -u - 1 |(t, x)∇u - 1 (t, x) dx ≤ C Ω T 1/4 M (γ 1 + γ 2 )(1 + 2M ) Ω u - 1 (t, x)∇u - 1 (t, x) dx ≤ CT 1/4 1 4α ∇u - 1 2 L 2 (Ω) + α u - 1 2 L 2 (Ω)

RR n°7998

by Cauchy-Schwarz estimates and the well-known Peetre's inequality with α > 0 large enough. Thus since |s|/(1 + |s|) ≤ 1 for any s ∈ R we infer

∂ t u - 1 (t, •) 2 L 2 (Ω) ≤ α CT 1/4 u - 1 (t, •) 2 L 2 (Ω) . Gronwall's lemma implies therefore that u - 1 (t, •) 2 L 2 (Ω) = 0, since u - 1 (0, •) equals zero. Prove now that u 2 ≤ 1. Let U 2 = u 2 -1: ∂ t U 2 = ∆U 2 -λu 1 U 2 ,
hence multiplying by U + 2 = max(0, U 2 ) the above equation and since U + 2 (0, •) = 0 we infer that U 2 vanishes everywhere thanks to Gronwall's lemma, and therefore u 2 ≤ 1.

Similarly let

U 1 = u 1 -1.
Then, U 1 satisfies

∂ t U 1 = d 1 ∆U 1 -λ1 Ω (1 -u 2 ) (U 1 + 1) + χ 0 ∇. (U 1 + 1) v 1 + |v| U 1 ∇v , in Ω, (12)
Once again multiply ( 12) by U + 1 = max(U 1 , 0) and integrate by parts to obtain

1 2 ∂ t U + 1 2 L 2 (Ω) = -d 1 ∇U + 1 2 L 2 (Ω) -λ Ω (1 -u 2 )U + 1 (U + 1 + 1) dx -χ 0 Ω (U + 1 + 1) v 1 + |v| U + 1 ∇v • ∇U + 1 dx.
Since 1 -u 2 is positive and using Cauchy-Schwarz estimate and Peetre inequality for α large enough (as used above to prove that u 1 ≥ 0) implies that

∂ t U + 1 2 L 2 (Ω) ≤ α C U + 1 2 L 2 (Ω) .
Therefore, Gronwall lemma implies that U + 1 vanishes almost everywhere in (0, T ) × Ω hence u 1 ≤ 1.

To obtain the positivity of v, first multiply (1c) by v -and integrate by part to infer, since u 1 and u 2 are positive that:

∂ t v -2 L 2 (Ω) ≤ 0. Then, the function V = v -η -1 (γ 1 + γ 2 ) satsifies ∂ t V = ∆V -ηV + γ 1 (u 1 -1) + γ 2 (u 2 -1).
Since γ 1 (u 1 -1) + γ 2 (u 2 -1) ≤ 0, we infer that V + identically vanishes after multiplication and integration by parts, hence

0 ≤ v ≤ η -1 (γ 1 + γ 2 ).
From the implicit representation integral of u 1 and u 2 we deduce easily that if T M is the maximal time of existence, then there exists a sequence (t n ) n∈N tending to T M , with

t n < T M such that lim n→+∞ u 1 (t n , •) L ∞ = +∞,
hence u 1 and u 2 exists for almost any t ∈ (0, +∞) by contraposition. We obtain a good agreement with the expected evolution. Indeed, the density of the cells on the adhesive areas increases in time, whereas outside it becomes very small. Cells are stucked on the strips and stop moving once they are over them. As a consequence the density of the attractant on the strips also increases.

Behavior on large strips

We now consider a domain composed by two large strips. The geometry is presented in Fig. 8(b) In Fig. 9, we present the total density of endothelial cells (u = u 1 + u 2 ) at times 0.3 (in Fig. 9(a)) and 1.0 (in Fig. 9(b)) obtained for the choice of parameters :

d 1 = d 2 = χ 0 = γ 2 = 1, γ 1 = 0.5, λ = 100.
In Fig. 10, we present the behavior of v for the same choice of parameters at time t = 0.3 and t = 0.6.

As previously we observe a behavior in good agreement with the experiments. When consid- ering two large adhesive areas the velocity of the migration is smaller than for a large number of thin strips. Indeed, at the time step t = 1.0 we observe that with thin strips the migration seems to be more advanced than in the case with large strips. This could be explained by the fact this last case some cells are far away from a strip and their migration toward the strips take more time.

Influence of the number of strips on the migration

We want to study the influence of the geometry on the migration. We set the surface of the adhesive domain, and let the number of strips, N s , vary. The average of u 2 in term of the time for N s = 1, 2, and 4 is presented in Fig. 11 . We observe that when considering four strips the migration is quicker. Moreover the mean density reached is higher, which corroborates the experiments.

Inria 5 Conclusion

In this paper, a macroscopic model for endothelial cells migration is presented. Its major biological assumption is that cells are not actively attracted by the adhesion patterns but just adhere on it and try to gather.

Mathematically, mass conservation and global existence is shown. Numerically, the model behaves in good agreement with the biological knowledge. Despite the lack of active attraction by the adhesion patterns, the non-washed out endothelial cells end up on the patterns. We have observed two facts that have been reported by the experiments: 1) For a given surface of active principle the process of cell migration is more efficient with a large number of thin strips than with a small number of large strips.

2) There exists a minimum value of the initial density of endothelial cells to impose in order to have an optimal cell migration towards the active principle.

We therefore believe that this model is a first step towards better understanding of cell migration on micropatterns, the long-term goal being optimal designing of patterns in order to build biological tissues.

Figure 1 :

 1 Figure 1: Endothelial cells alignment onto micropatterned polymer (PET) (10µm (A) and 300 µm (B) stripes of SVVYGLR peptides)[START_REF] Lei | Micropatterning of polyethylene terephthalate with RGD peptides for induction of endothelial cell alignment and morphogenesis[END_REF]. The distance between bioactive patterns is 100 µm.

Definition 3 . 1 .

 31 The kernels B, G and G of the respective operators ∂ t -∆ + η, ∂ t -d 1 ∆ on Ω, and ∂ t -d 2 ∆ on Ω, all with homogeneous Neumann conditions, are respectively defined by

Figure 5 :

 5 Figure 5: Initial setup: endothelial cells (left) and adhesion substrat (right).

Figure 6 :

 6 Figure 6: The density of endothelial cells u at two different time steps for d 1 = d 2 = χ 0 = γ 2 = 1, γ 1 = 0.5, λ = 100.

Figure 7 :

 7 Figure 7: The density v at two different time steps for d 1 = d 2 = χ 0 = γ 2 = 1, γ 1 = 0.5, λ = 100.

Figure 8 :

 8 Figure 8: Initial setup: endothelial cells (left) and patterns (right).

Figure 9 : 6 Figure 10 :Figure 11 :

 961011 Figure 9: The density of endothelial cells u at two different time steps for d 1 = d 2 = χ 0 = γ 2 = 1, γ 1 = 0.5, λ = 100

RR n°7998

Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

Theorem 2.1 is an easy consequence of the above results. The following result is a straightforward consequence of proposition 4. It ensures that the mass of the cells tends to concentrate on the micropatterns.

Corollary 2. Let (u 0 1 , u 0 2 ) ∈ L ∞ (Ω) × L ∞ ( Ω) such that 0 ≤ u 0 1 ≤ 1, 0 ≤ u 0 2 ≤ 1, and let (u 1 , u 2 ) the weak solution to problem [START_REF] Anderson | Continuous and discrete mathematical models of tumorinduced angiogenesis[END_REF]. Then

Numerical results

We provide now numerical schemes that are used to compute problem [START_REF] Anderson | Continuous and discrete mathematical models of tumorinduced angiogenesis[END_REF], and then we show the simulations that corroborate the experimental results.

Approximation of the problem

We consider a cartesian mesh (composed by quadrilaterals). We discretize the model using finite volume method [START_REF] Eymard | Finite Volume Methods, Handbook of Numerical Analysis[END_REF] and we use an implicit Crank-Nicolson scheme for the time discretization. We solve the model using a decoupled approach [START_REF] Hundsdorfer | Numerical solution of time-dependent advectiondiffusion-reaction equations[END_REF]. In particular, the first equation is split into advection and diffusion parts. Let us recall the expression of this equation :

To simplify the notations we define A and B as:

, and B(u

Let us denote the time step by ∆t, set t n = n∆t and let (u n 1 , u n 2 , v n ) be the solution at the time t n . At each time step we first solve the diffusive part :

. For all the diffusive terms, the spatial discretization is handled by a centered finite volume scheme, all the species being computed at the centre of each element of the mesh. We then solve the advection part :

The high order WENO 5 (Weighted Essentially Non-Oscillatory) finite difference scheme introduced in [START_REF] Liu | Weighted essentially non-oscillatory schemes[END_REF] and improved in [START_REF] Jiang | Efficient implementation of weighted ENO schemes[END_REF] has been used to approach the convective term. These solvers are implemented in the academic library eLYSe 1 .

For the numerical computations, we consider a bounded domain Ω = [0, L] × [0, L], with L > 0. Initially, endothelial cells are randomly distributed on the computational domain.

In the following, the initial conditions write

where u 0 is a function of x ∈ Ω. Hence at the initial time the support of u 1 and u 2 are disjoint.

1 http://www.math.u-bordeaux1.fr/~saut/ RR n°7998

Mathematical behavior of the model

In this paragraph, we present numerical results illustrating the good mathematical behavior of our model. We want to check the properties of the model, when the maximal density on the adhesive area is reached. We consider a domain Ω = [0, 1] × [0, 1] and meshes composed by 100 × 100 quadrilaterals. The domain is composed by a unique adhesive area (cf Fig. 2).

Figure 2: Geometry of the micropattern.

At the initial time the cells are uniformly distributed meaning u 0 of ( 14) is constant. We consider two different values of u 0 :

We plot the results along the axis {y = 0.5} in order to have a profile of the distribution of u 1 and u 2 in the domain. The densities u 1 along the axis at different time steps are given by Fig. 3 and the densities u 2 at the same time steps are represented in Fig. 4. When considering u 0 = 0.08, the maximal density on the adhesive area is never reached. We observe that u 1 is decreasing, while u 2 is increasing with respect to the time. In the second case, u 0 = 0.25, the maximal density for u 2 is reached for t = 0.3. As expected the migration stops and no more movement are observed. These simulations show that a minimum amount of endothelial cells is needed at initial time in order to reach an optimal concentration on the strips at the end of the experiment. If this initial concentration is too small the final density of endothelial cells is suboptimal.

Behavior on realistic benchmarks

We now provide simulations in realistic confirgurations. Therefore, throughout this subsection the function u 0 of ( 14) is a random spatial distribution between 0 and 1 using a normal distribution.

Behavior on thin strips

We first consider a micropattern composed by six thin adhesive areas (in red on Fig. 5(b)).

The simulation Fig. 6, represents the total density of endothelial cells (u = u 1 + u 2 ) at times 0.


Mounira Bachir
email: mounira.bachir@lip6.fr

Albert Cohen
email: albert.cohen@inria.fr

Sid-Ahmed-Ali Touati
email: sid.touati@inria.fr

On the Effectiveness of Register Moves to Minimise Post-Pass Unrolling in Software Pipelined Loops

Keywords: Instruction Level Parallelism, Compiler, Register Allocation, Software Pipelining, Loop Unrolling, Code Optimisation

Software pipelining is a powerful technique to expose fine-grain parallelism, but it results in variables staying alive across more than one kernel iteration. It requires periodic register allocation and is challenging for code generation: the lack of a reliable solution currently restricts the applicability of software pipelining. The classical software solution that does not alter the computation throughput consists in unrolling the loop a posteriori [12], [11]. However, the resulting unrolling degree is often unacceptable and may reach absurd levels. Alternatively, loop unrolling can be avoided thanks to software register renaming. This is achieved through the insertion of move operations, but this may increase the initiation interval (II) which nullifies the benefits of software pipelining. This article aims at tightly controling the post-pass loop unrolling necessary to generate code. We study the potential of live range splitting to reduce kernel loop unrolling, introducing additional move instructions without inscreasing the II. We provide a complete formalisation of the problem, an algorithm, and extensive experiments. Our algorithm yields low unrolling degrees in most cases -with no increase of the II.

I. INTRODUCTION

Our focus is on the exploitation of instruction-level parallelism (ILP) in embedded VLIW processors [START_REF] Faraboschi | Embedded Computing: a VLIW Approach to Architecture, Compilers and Tools[END_REF]. Increased ILP translates into higher register pressure and stresses the register allocation phase(s) and the design of the register files. In the case of software-pipelined loops, variables can stay alive across more than one kernel iteration, which is challenging for code generation and generally addressed through: (1) hardware support -rotating register files -deemed too expensive for almost embedded processors, (2) insertion of register moves with a high risk of reducing the computation throughputinitiation interval (II) -of software pipelining, and (3) postpass loop unrolling that does not compromise throughput but often leads to unpractical code growth.

We investigate ways to keep the size of the generated code compatible with embedded system constraints without compromising the throughput benefits of software pipelining. Having a minimal unroll factor reduces code size, which is an important performance measure for embedded systems because they have a limited memory size. Regarding high performance computing (desktop and supercomputers), loop code size may not be important for memory size, but may be so for I-cache performance. In addition to the minimal unroll factors, it is necessary that the code generation scheme for periodic register allocation does not generate additional spill; the number of required registers must not exceed MAXLIVE [START_REF] Huff | Lifetime-sensitive modulo scheduling[END_REF] (the maximum number of values simultaneously alive).

When the instruction schedule is fixed then the circular lifetime intervals (CLI) and MAXLIVE are known. In this case, known methods exist for computing unroll factors: [START_REF] Bachir | Loop Unrolling Minimisation for Periodic Register Allocation[END_REF] Modulo Variable Expansion (MVE) [START_REF] Faraboschi | Embedded Computing: a VLIW Approach to Architecture, Compilers and Tools[END_REF], [START_REF] Lam | Software pipelining: an effective scheduling technique for vliw machines[END_REF] which computes a minimal unroll factor but may introduce spill (since MVE may need more than MAXLIVE registers without proving an appropriate upper-bound); (2) Hendren's heuristic [START_REF] Hendren | A register allocation framework based on hierarchical cyclic interval graphs[END_REF] which computes a sufficient unroll factor without introducing spill, but with no guarantee in terms of minimal register usage or unroll factor; (3) the Meeting Graph framework [START_REF] De Werra | Circular-arc graph coloring: On chords and circuits in the meeting graph[END_REF] whose unroll factor is formally proven to minimize register usage (MAXLIVE), but which does not guarantee a minimal unroll factor.

In [START_REF] Bachir | Post-pass periodic register allocation to minimise loop unrolling degree[END_REF], [START_REF] Bachir | Using the meeting graph framework to minimise kernel loop unrolling for scheduled loops[END_REF], the authors claim that the loop unrolling minimisation (LUM) using extra remaining registers after a periodic register allocation is an efficient method to bring loop unrolling as low as possible -with no increase of the initiation interval (II). However, some kernel loops may still require high unrolling degrees. These occasional high unrolling degrees suggest it may be worthwhile to consider combining the insertion of move operations with kernel loop unrolling.

This article proposes a technique to minimise the unroll factor after periodic register allocation by inserting move instructions without compromising the throughput benefits of software pipelining 1 . We study all the possible cases of adding move operations without altering the initiation interval (II) for already software pipelined loops. This is done by splitting a chosen variable into two new variables at a given clock cycle if a free functional unit exists and can execute the added move operation.

The paper is structured as follows. Section II discusses the limitations of using the interference graph for register allocation and briefly recalls the use of circular lifetime intervals (CLI) for periodic register allocation. Section III presents the closest related work on code generation. Section IV illustrates our idea intuitively on simple examples. In Section V, we explain in more detail how it is possible to minimise loop unrolling through the insertion of extra moves. Section VI formally describes our exhaustive algorithm. In Section VII, we demonstrate the effectiveness of our approach in a standalone tool by hightlighting the different results. Finally, we present our conclusion in Section VIII.

II. PERIODIC REGISTER ALLOCATION: INTERFERENCE GRAPH VS CIRCULAR LIFETIME INTERVALS

Register allocation in many compilers performs the classical graph coloring method originally proposed by Chaitin [START_REF] Chaitin | Register allocation and spilling via graph coloring[END_REF]. The different nodes in the graph represent the different variables that are candidates for register allocation. Two nodes are connected by an arc if the corresponding live ranges variables are simultaneously alive. This graph is commonly called interference graph [START_REF] Chaitin | Register allocation and spilling via graph coloring[END_REF]. The register allocation attempts to find a k coloring of the interference graph. That is, an assignement of k colors to the nodes such as two connected nodes have different colors. If such coloring is found then it can map the k colors into registers. However, if the graph is not k-colorable, some nodes are removed from the interference graph until this latter becomes k-colorable. The removed nodes are spilled into the memory because they are not assigned to registers.

The interference graph representation of register allocation has some weaknesses. In fact, given a set of variables live ranges, the interference graph contains only the overlapping information of any two variables live ranges. By the way, it does not reveal the overlapping information for more than two variables and it does not provide the times of variables overlapping. The lack of this information between variables ovelapping may lead to some limitations.

For instance, Chaintin's register allocator may fail to find the minimum coloring even for simple loops. Figure 1(a) illustrates a simple loop program where MAXLIVE = 2 and the traditional Chaitin's [START_REF] Chaitin | Register allocation and spilling via graph coloring[END_REF] register allocator fails to find a register allocation with 2 registers. Some state-of-art work [START_REF] De Werra | Circular-arc graph coloring: On chords and circuits in the meeting graph[END_REF], [START_REF] Hendren | A register allocation framework based on hierarchical cyclic interval graphs[END_REF] studied the circular lifetime intervals as an alternative representation for register allocation.

Let us consider the interference graph given in Figure 1(b) which is similar to the one presented in [START_REF] Briggs | Coloring heuristics for register allocation[END_REF]. Even though this graph is clearly 2-colorable, Chaintin's heuristic fails because there is no node with degree < 2, and thus spill code has to be introduced.

To address this problem, we represent in Figure 1(c) an alternative representation called circular lifetime intervals (CLI) [START_REF] De Werra | Circular-arc graph coloring: On chords and circuits in the meeting graph[END_REF], [START_REF] Lelait | A new fast algorithm for optimal register allocation in modulo scheduled loops[END_REF], [START_REF] Hendren | A register allocation framework based on hierarchical cyclic interval graphs[END_REF]. The x-axis represents the different clock cycles and the y-axis represents the live ranges of the different variables. The barrier in the diagram represents the definition point of the variable while the arrow represents their last use. For example, variable X1 is defined at clock cycle 1 and finishes at clock cycle 3. Note that the lifetime of the variable can span more than one iteration, as shown with variable X4. Since the relation between the intervals is periodic, we characterise this relation within one iteration.

FOR i=0, i<N, i=i+1 a[i+2] = b[i]+1 ; b[i+2] = c[i]+2 ; c[i+2] = a[i]+3 ; ENDFOR
As in any general graph coloring problem, finding the minimum coloring of a circular lifetime intervals is NP-hard [START_REF] Garey | The complexity of coloring circular arcs and chords[END_REF]. Many heuristics, exploiting the information provided in such representation are proposed in the litterature [START_REF] De Werra | Circular-arc graph coloring: On chords and circuits in the meeting graph[END_REF], [START_REF] Hendren | A register allocation framework based on hierarchical cyclic interval graphs[END_REF]. The number of minimum registers needed for a circular interval is related to the width of the circular lifetime intervals. The main caracteristics related to circular lifetime intervals [START_REF] Hendren | A register allocation framework based on hierarchical cyclic interval graphs[END_REF] are as follows:

• Two intervals I1, I2 are overlapping if there exists a clock cycle t where both of the intervals I1 and I2 are simultaneously alive. • The width of a given circular lifetime intervals F at time t, written as W idth(F, t), is the number of intervals which are simultaneously alive at this time t. • The maximum width of a given circular lifetime intervals F , written as W (F ), is the maximum W idth(F, t), for all clock cycle t ∈ [0, II[. The maximum width is called also MAXLIVE, the maximum number of variables which are simultaneously alive. The following section reviews the main state-of-art periodic register allocation.

III. CODE GENERATION: BACKGROUND AND CHALLENGES

We review the main issues and approaches to code generation for periodic register allocation using the loop example described in Listing 1.

i=0, i<N, i=i+1 R3 = b[i]+1 ; b[i+2] = c[i]+2 ; c[i+2] = R1+3 ; R1 = R2 ; R2 = R3
; ENDFOR There are two ways to deal with periodic register allocation: using special architecture support such as rotating register files, or without using such support. This latter may require the insertion of move operations or loop unrolling

A. Rotating Register File

A rotating register file (RRF) [START_REF] Dehnert | Overlapped loop support in the cydra 5[END_REF] is a hardware mechanism to prevent successive lifetime intervals from being assigned to the same physical registers.

In Listing 1, variable a[i] spans three iterations (defined in iteration i -2 and used in iteration i). Hence, at least 3 physical registers are needed to carry simultaneously a[i], a[i + 1] and a[i + 2]. A rotating register file R automatically performs the move operation at each iteration. R acts as a FIFO buffer. The major advantage is that instructions in the generated code see all live values of a given variable through a single operand, avoiding explicit register copying. Below R[k] denotes a register with offset k from R.

Iteration i Iteration i+2 R=b

[i]+1 ; R[+2]=b[i]+1 ; b[i+2]=c[i]+2 ; b[i+2]=c[i]+2 ; c[i+2]=R[-2]+3 ; c[i+2]=R+3 ;
Using a RRF avoids increasing code size due to loop unrolling, or to decrease the computation throughput due to the insertion of move operations.

B. Move Operations

This method is also called register renaming. Considering the example of Listing 1 for allocating a[i], we use 3 registers and perform move operations at the end of each iteration [START_REF] Llosa | Reduced code size modulo scheduling in the absence of hardware support[END_REF], [START_REF] Nicolau | Register allocation, renaming and their impact on fine-grain parallelism[END_REF], [START_REF] Lelait | A new fast algorithm for optimal register allocation in modulo scheduled loops[END_REF]:

a[i] in register R1, a[i + 1] in register R2 and a[i + 2] in register R3.
Then we use move operations to shift registers across the register file at every iteration as shown in Listing 2. However, it is easy to see that if variable v spans d iterations, we have to insert d -1 extra move operations at each iteration. In addition, this may increase the II and may require rescheduling the code if these move operations do not fit into the kernel. This is generally unacceptable as it negates most of the benefits of software pipelining.

C. Loop Unrolling

Another method, loop unrolling, is more suitable to maintain II without requiring hardware support such as RRF. The resulted loop body itself is bigger but no extra operations are executed in comparison with the original code. Here different

Listing 3. Example of Loop Unrolling FOR i=0, i<N, i=i+3 R1 = b[i]+1 ; b[i+2] = c[i]+2 ; c[i+2] = R2+3 ; R2 = b[i+1]+1 ; b[i+3] = c[i+1]+2 ; c[i+3] = R3+3 ; R3 = b[i+2]+1 ; b[i+4] = c[i+2]+2 ; c[i+4] = R1+3 ; ENDFOR
registers are used for different instances of the variable a of Listing 1. In Listing 3, the loop is unrolled three times.

a[i+2] is stored in R1, a[i + 3] in R2, a[i + 4] in R3, a[i + 5] in R1,
and so on.

By unrolling the loop, we avoid inserting extra move operations. The drawback is that the code size will be multiplied by 3 in this case, and by the unrolling degree in the general case. This can have a dramatic impact by causing unnessary instruction cache misses when the code size of the loop happens to be larger than the size of the instruction cache. For simplicity, we did not expand the code to assign registers for b and c. In addition, brute force searching for the best solution using loop unrolling has a prohibitive cost, existing solutions may either sacrifice the register optimality [START_REF] Hendren | A register allocation framework based on hierarchical cyclic interval graphs[END_REF], [START_REF] Lam | Software pipelining: an effective scheduling technique for vliw machines[END_REF], [START_REF] Rau | Register allocation for software pipelined loops[END_REF] or incur large unrolling overhead [START_REF] De Werra | Circular-arc graph coloring: On chords and circuits in the meeting graph[END_REF], [START_REF] Touati | Early Periodic Register Allocation on ILP Processors[END_REF].

1) Modulo Variable Expansion: Lam designed a general loop unrolling scheme called modulo variable expansion (MVE) [START_REF] Lam | Software pipelining: an effective scheduling technique for vliw machines[END_REF]. In fact, the major criterion of this method is to minimize the loop unrolling degree because the memory size of the i-WARP processor is low [START_REF] Lam | Software pipelining: an effective scheduling technique for vliw machines[END_REF]. The MVE method defines a minimal unrolling degree to enable code generation after a given periodic register allocation. This unrolling degree is obtained by dividing the length of the longest live range (max v LT v ) by the number of cycles of the kernel α = ⌈ maxv LTv II ⌉. Once the loop is unrolled, MVE uses the interference graph for allocation.

Having MAXLIVE the maximum number of values simultaneously alive, the problem with MVE is that it does not guarantee a register allocation with minimal number of register equal to MAXLIVE [START_REF] De Werra | Circular-arc graph coloring: On chords and circuits in the meeting graph[END_REF], [START_REF] Hendren | A register allocation framework based on hierarchical cyclic interval graphs[END_REF], and in general it may lead to unnecessary spills breaking the benefits of software pipelining.

A concrete examples of this limitation can be found in [START_REF] Bachir | Loop Unrolling Minimisation for Periodic Register Allocation[END_REF]. In Listing 1, the longest live range lasts 8 cycles and the number of cycles of the loop is 3 cycles, so α = ⌈ 8 3 ⌉, and we should unroll the loop 3 times. Then we can assign to each variable a number of registers equal to the least integer greater than the span of the variable that divides II. In Listing 4, each variable a, b, c is assigned 3 registers using MVE: R1, R2, R3 for a, R4, R5, R6 for b, R7, R8, R9 for c, and the loop is unrolled 3 times.

One can verify that it is not possible to allocate the different variables on less than 9 registers when unrolling the loop 3 times. But MVE does not ensure a register allocation with a minimal number of registers, and hence it is not optimal. As we will see in the next section, we need 8 registers to In MVE, the round up to the nearest integer for choosing the unrolling degree may miss opportunities for achieving an optimal register allocation.

2) Meeting Graphs: The algorithm of Eisenbeis et al. [START_REF] De Werra | Circular-arc graph coloring: On chords and circuits in the meeting graph[END_REF] can generate a periodic register allocation using a minimal number of registers equal to MAXLIVE if the kernel is unrolled, thanks to a dedicated graph representation called the meeting graph. It is a more accurate graph than the usual interference graph, as it has information on the number of clock cycles of each variable lifetime and on the succession of the lifetimes all along the loop. It is based on circular lifetime intervals (CLI). A preliminary remark is that without loss of generality, we can consider that the width of the interval representation is constant at each clock cycle. If not, it is always possible to add unit-time intervals in each clock cycle where the width is less than MAXLIVE [START_REF] De Werra | Circular-arc graph coloring: On chords and circuits in the meeting graph[END_REF].

The formal definition of the meeting graph is as follows.

Definition 1 (Meeting Graph). Let F be a set of circular lifetime intervals graph with constant width r. The meeting graph related to F is the directed weighted graph G = (V, E).

V is the set of circular intervals. Each edge e ∈ E represents the relation of meeting. In fact, there is an edge between two nodes v i and v j iff the interval v i ends when the interval v j begins. Each v ∈ V is weighted by its lifetime length in terms of processor clock cycles.

The meeting graph (MG) allows us to compute an unrolling degree which enables an allocation of the loops with RC=MAXLIVE registers. It can have several connected components of weights µ 1 , . . . , µ k (if there is only one connected component, its weight is µ 1 = RC), this leads to the upper bound of unrolling α = lcm(µ 1 , ..., µ k ) (RC if there is only one connected component). Moreover a possible lower bound of loop unrolling is computed by decomposing the graph into as many circuits as possible and then computing the least common multiple (lcm) of their weights [START_REF] De Werra | Circular-arc graph coloring: On chords and circuits in the meeting graph[END_REF]. The circuits are then used to compute the final allocation. This method can handle variables that are alive during several iterations. This allocation always finds an allocation with a minimal number of registers (MAXLIVE). The main drawback of the meeting graph method is that the loop unrolling degree can be high in practice although the number of registers used is minimal. That may cause spurious instruction cache misses or even be impracticable due to the memory constraints, like in embedded processors. In order to minimise loop unrolling, Bachir et al. [START_REF] Bachir | Using the meeting graph framework to minimise kernel loop unrolling for scheduled loops[END_REF], [START_REF] Bachir | Post-pass periodic register allocation to minimise loop unrolling degree[END_REF] propose a method called loop unrolling minimisation (LUM method) that minimises the loop unrolling degree by using extra remaining registers for a given periodic register allocation. The LUM method brings loop unrolling as low as possible -with no increase of the II. However, some kernel loops may still require high unrolling degrees. This occasional cases suggest that it may be worthwhile to consider combining the insertion of move operations with kernel loop unrolling.

IV. SPLITTING VARIABLES LIFETIMES: MOTIVATION

As illustration in this section, we first present the simple example described in Figure 3 to illustrate the weaknesses of the loop unrolling minimisation (LUM method) described in [START_REF] Bachir | Using the meeting graph framework to minimise kernel loop unrolling for scheduled loops[END_REF], [START_REF] Bachir | Post-pass periodic register allocation to minimise loop unrolling degree[END_REF]. We then outline how these problems can be handled by adding move operations.

Let us consider the two lifetime intervals a and b as described in Figure 3. Let us assume that the variable a is resulted from a multiplication and b is resulted from an addition. From these lifetimes, the corresponding meeting If we want to not unroll the loop, we notice that the loop unrolling minimisation does not provide better loop unrolling degree bacause the minimal loop unrolling degree α * ≥ max i µ(C i ). In our example α * = µ(C 1 ) = 2.

One possibility to reduce loop unrolling degree is called register renaming [START_REF] Nicolau | Register allocation, renaming and their impact on fine-grain parallelism[END_REF]. This latter proposes to allocate the variable a with 2 registers and perform move operation to carry the value from the first register to the second one. By the way, if we split the variable a into two variables (a ′ , a ′′ ) where each spans less than one kernel iteration, then the problem of loop unrolling can be resolved.

In Figure 3, we assume that the adder unit is free at the clock cycle II. If we decide to split the variable a into two variables a ′ and a ′′ at the clock cycle II then we add only one move operation (a ′′ ← a ′ ) which is performed by using the free adder unit. The resulting allocation can be done with two registers as shown in Figure 3 without unrolling the loop (α * = lcm(µ(C 1 ), µ(C 2 )) = 1). Another running example of splitting variables is represented in Figure 4. In fact, Figure 4 shows how to minimise the loop unrolling degree by splitting the variables of the circular lifetime intervals of Figure 2. In Figure 2(a), variable lifetimes are depicted by circular lifetime intervals. Each variable is alive during 8 cycles, this means 3 iterations as II = 3. From these lifetimes, the corresponding meeting graph is drawn in Figure 2(b). The interval family has a width equal to µ(C) = 8. So hence, we need 8 registers to carry successfully three variables a, b, c if we unroll the loop 8 times. In Figure 4, we assume that adding one move operation at each clock cycle does not alter II. Each variable is splitted into one 6 cycles interval and one 2 cycles interval. In Figure 4(a), variable a is cut after 6 cycles and hence it results two The following section describes in details the exhaustive solution to minimise the loop unrolling degree by splitting variables without increasing the computation throughput.

V. SPLITTING VARIABLES TO MINIMISE LOOP UNROLLING DEGREE

The impact of a loop register allocation scheme can be measured by 3 parameters. The first one is the number RC of allocated registers. A lower bound for RC is the maximal number of simultaneously alive variables, denoted as MAXLIVE. The second one is the unrolling degree α. Large unrolling degrees are induced by the fact that a least common multiple (lcm) is computed. Roughly speaking, if you must unroll p times for one set of variables and q times for another one, then you have to unroll at least lcm(p, q) times. When we use the loop unrolling minimisation (LUM method) by using the remaining registers [START_REF] Bachir | Using the meeting graph framework to minimise kernel loop unrolling for scheduled loops[END_REF], [START_REF] Bachir | Post-pass periodic register allocation to minimise loop unrolling degree[END_REF], the minimal loop unrolling degree can still be high. This is due to the formulation of the LUM method which proposes to update the different circuits weights µ 1 , . . . , µ k generated by the meeting graph by adding the remaining registers. The resulting final loop unrolling degree α * is greater or equal to the maximum circuit weight (α * >= max(µ 1 , . . . , µ k )).

This paper proposes to find a periodic register allocation with a small circuit weights which leads to a small lcm. Our exploratory method consists in splitting the different variables by adding m extra move instructions. The impact of this third parameter is hard to measure because sometimes the added move instructions can be executed in parallel with the other operations. Based on these parameters our study follows these different phases:

1) Consider a software pipelined loop 2) Find the number of free functional units at each clock cycle which execute move instruction while preserving the Initiation Interval (II) 3) Extract the circular lifetimes intervals (CLI) from the software pipelined loop devoted to be allocated by the meeting graph 4) Compute all the possible combinations of splitting the different variables of CLI following the number of possible move at each clock cycle 5) For each configuration of splitting variables, compute the periodic register allocation and then choose the one with the minimal loop unrolling degree Step 2 aims to collect information about the number of move operations at each clock cycle that can be executed without altering II. Once this information available, step 4 computes all possible splitting variables following the number of possible move operations at each clock cycle. Each splitting variable generates a new CLI. Finally, having all the different CLIs, we perform for each one the periodic register allocation using the meeting graph technique and then we choose the best CLI with the minimal loop unrolling degree.

Figure 5 portrays all the different steps. An arrow between two phases means that the target phase needs the information of the source phase. We note CLI for Circular Lifetimes Intervals and PRA for Periodic Register Allocation (PRA i means the meeting graph of the i th CLI).

VI. EXHAUSTIVE SEARCH ALGORITHM

In this section, we explain our exhaustive algorithm to find a periodic register allocation with a minimal loop unrolling size. First of all, we assume that we have the initial circular lifetimes intervals (CLI) and we provide the number of possible added move instructions at each clock cycle t = 0, II -1. These move instructions, if added, can be executed in parallel with other operations without altering II. In our study, the number of possible added move instructions at each clock cycle t is denoted m t and computed as m t = max(0, MAXMOVEdef t )

Where MAXMOVE is the number of functional units performing move operations and def t is the number of variables which are defined at clock cycle t.

In order to compute all the possible cases of splitting variables for a given CLI, we look for the number of variables which can be splitted at each clock cycle t. In fact, a variable v can be splitted into two variables v ′ and v ′′ at clock cycle t as explained in the following definition. The variable v is splitted into two variables v ′ , v ′′ at clock time t iff and only iff S(v) < t < E(v). The resulting variables v ′ and v ′′ have the following lifetimes:

• S(v ′ ) = S(v) and E(v ′ ) = t • S(v ′′ ) = (t mod II) and E(v ′′ ) = E(v) -t
When the variable v is splitted, new circular lifetime intervals called Result is created such as:

Result = {l ∈ CLI | l <> v} ∪ {v ′ , v ′′ }
From Definition 2, the number w t of variables which can be splitted at clock cycle t is the number of variables which are alive but not defined at this clock cycle. Let us assume that W idth(CLI, t) is width of CLI at clock cycle t and def t is the number of variables which are defined at this clock cycle then the number of variables which can be splitted at time t is as w t = W idth(CLI, t)def t Consequently, for a given CLI and a given time t, we can split at most m t variables among w t alives one. That is, we can decide either to not split any variable (saving the original CLI) or to split In order to find the set G of all possible splitting variables configurations, we need to store the results of splitting variables at each clock cycle t, which are then reused later at cycle = t + 1, II -1. In fact, we need this new CLIs to compute further splitting configurations. Formally, the cardinal of the set G of all possible splitting variables configurations is as follows:

|G| = II-1 t=0 mt i=0 C i wt
Algorithm 1 implements our exhaustive search for splitting variables following possible added move intructions. In this algorithm, we require the initial CLI, the initiation interval II and the information m t about available units at each clock cycle t. This algorithm returns all the possible splitting variables configurations (CLIs) which are stored in the set G by computing the different combinations. In fact, at time t, splitting i variables among w t variables (the generation of C i wt CLIs) consists in splitting one variable for the different CLIs generated by splitting (i -1) variables among w t variables (the already generated C i-1 wt CLIs). That is, at time t, we use two other sets G c which contains the CLIs resulting from a given combination (C i wt ) and the set G t which contains the CLIs resulting from the previous combination C i -1 wt or which initially contains the generated CLIs at clock time t -1.

As we can see in Algorithm 1, the principal algorithm calls the sub-algorithm SPLIT-VARIABLE with the following parameters: the given CLI, the variable v which we want to split into two variables, the cycle window where exactly we want to split the variable v (variable v can span more than one kernel iteration). This sub-algorithm returns the generated CLI called Result and a boolean equal to TRUE if the variable v is actually splitted at clock cycle window, FALSE otherwise. Algorithm 1 delivers all the possible splitting configurations in finite time because it is dominated by the number of generated CLIs:

II-1 t=0 mt i=0 C i wt VII. EXPERIMENTAL RESULTS
To study the efficiency of our empirical study, we developed a standalone tool to generate all the possible new CLIs following the number of possible move operations at each clock cycle and then perform for each generated CLI the meeting graph framework in order to find the periodic register allocation with a minimal loop unrolling.

We did extensive experiments on 1935 CLIs extracted from various known benchmarks, namely Livermore loops, Linpack and Nas. In our experiments, the processor has one functional unit which can execute move operations (adder unit).

In order to have an idea about the complexity of the exhaustive solution, we draw in Table I the number of experimented CLIs for each benchmarks, the number of CLIs where it is possible to add move operation and the number of generated CLIs needed to find the best CLI with a minimal number of registers.

Table I shows that from 1935 there are only 343 CLIs where it is possible to add move operations without increasing II. The number of the generated CLIs is 108934. On average, the running time of our exhaustive solution on a dual core 2 GHz Linux PC is 272 milli-seconds. The maximal observed runtime of Algorithm 1 is about 84 seconds.

In order to highlight the improvements of splitting variables to minimise loop unrolling, we performed a periodic register allocation for each generated CLI and then compared the initial loop unrolling before splitting and the minimal loop unrolling after splitting. In addition, we wanted to know if the loop unrolling minimisation method (LUM method) described in [START_REF] Bachir | Using the meeting graph framework to minimise kernel loop unrolling for scheduled loops[END_REF], [START_REF] Bachir | Post-pass periodic register allocation to minimise loop unrolling degree[END_REF] provides lower loop unrolling.

In these experiments, we take into account only CLIs where it is possible to add move operations. We varied the number of architectural registers (R arch ) from 16 to 256 Table II displays for each method the following observations : the lower quartile (Q1 = 25%), the median (Q2 = 50%), the upper quartile (Q3 = 75%), and the largest loop unrolling degree (max). Furthermore, we use for comparison the initial loop unrolling, the final loop unrolling using the loop unrolling minimisation (LUM method), the final loop unrolling using splitting variables and final loop unrolling using splitting variables method followed by the LUM method.

From the different results, we remark that adding move operations decreases dramatically the initial loop unrolling degree. Actually, thanks to splitting method, we do not unroll 50% of loops and at most we unroll twice 75% of loops. We also see in Table II that combining the loop unrolling minimisation method with splitting method decreases greatly the maximum loop unrolling degree. For example, in a given machine with 64 registers, we minimise the loop unrolling from 330 to 8.

VIII. CONCLUSION

Minimising loop unrolling by using extra remaining registers is an efficient method to reduce unroll factors for periodic register allocation. However, some real-world loops still require high unrolling degrees [START_REF] Bachir | Post-pass periodic register allocation to minimise loop unrolling degree[END_REF], [START_REF] Bachir | Using the meeting graph framework to minimise kernel loop unrolling for scheduled loops[END_REF]. This occasional cases suggest that it may be worthwhile to combine the insertion of move operations with kernel loop unrolling.

R arch

Loop Unrolling Q1 Median Q3 MAX We study the potential for live range splitting to reduce kernel loop unrolling, describing how to introduce new moves without inscreasing II. Moreover, we show that an exhaustive solution for splitting variables can produce satisfactory results. For example, when we assume that we can add at most one move operations per cycle, then half of the loops do not need to be unrolled and 25% more need an unroll factor of 2. In addition, the maximun loop unrolling is highly reduced while combining splitting variables method with loop unrolling minimisation method.

In most cases, our method is a satisfactory solution to generate compact code after periodic register allocation. Nevertheless, we noticed that some loops still require high unrolling degrees despite our optimisation. This shows that optimal exploitation of remaining registers plus the insertion of IIpreserving moves may not always be sufficient to bring the unroll factor down to acceptable levels. This provides an interesting fundamental limit of the applicability of software pipelining without hardware support for rotating register files.

Fig. 1 .

 1 Fig. 1. Interference Graph vs Circular Lifetime Intervals

Listing 2 .

 2 Example of Register Renaming

FOR

Listing 4 .

 4 Example of MVE FOR i=0, i<N, i=i+3 R1 = R5+1 ; R4 = R8+2 ; R7 = R2+3 ; R2 = R6+1 ; R5 = R9+2 ; R8 = R3+3 ; R3 = R4+1 ; R6 = R7+2 ; R9 = R1+3 ; ENDFOR allocate the different variables.

Figure 2 (

 2 a) displays the circular lifetime intervals representing the different variables (a, b, c) of the loop example described in Listing 1, the maximum number of variables simultaneously alive MAXLIVE = 8. As shown in Figure 2(b), the meeting graph is able to use 8 registers to allocate the different variables instead of 9 with Modulo Variable Expansion by unrolling the loop 8 times. For the Listing 5. Example of Meeting Graph FOR i=0, i<N, i=i+8 R1 = R2+1 ; R4 = R5+2 ; R7 = R7+3 ; R2 = R3+1 ; R5 = R6+2 ; R8 = R8+3 ; R3 = R4+1 ; R6 = R7+2 ; R1 = R1+3 ; R4 = R5+1 ; R7 = R8+2 ; R2 = R2+3 ; R5 = R6+1 ; R8 = R1+2 ; R3 = R3+3 ; R6 = R7+1 ; R1 = R2+2 ; R4 = R4+3 ; R7 = R8+1 ; R2 = R3+2 ; R5 = R5+3 ; R8 = R1+1 ; R3 = R4+2 ; R6 = R6+3 ; ENDFOR loop described in Listing 1, the meeting graph generates the code shown in Listing 5.

Fig. 2 .

 2 Fig. 2. Meeting Graph of the Loop Example in Listing 1

2 bFig. 3 .

 23 Fig. 3. Splitting Variables Example

Fig. 4 . 2 Listing 6 .

 426 Fig.[START_REF] Briggs | Coloring heuristics for register allocation[END_REF]. Running Splitting Variables Method on Circular Lifetime Intervals described in Figure2

Fig. 5 .

 5 Fig. 5. Splitting Variables to Minimise Loop Unrolling

  1 variable among w t variables (generating C 1 wt possible CLIs) or to split 2 variables among w t variables (generating C 2 wt possible CLIs) or . . . or to split m t variables among w t ones (generating C mt wt possible CLIs). As results, the set of generated CLIs for a given CLI at time t denoted G t , is the sum of the different combinations C i wt where i is the number of variables which we want to split and w t the number of maximum splitted variables at time t. That is, the cardinal of the set G t is as follows:

	mt
	|G t | =	C i wt
	i=0
	Furthermore, if we have k different CLIs at time t then we
	generate for each CLI, number of generated CLIs at time t is k × mt i=0 C i wt CLIs. So hence, the total mt i=0 C i wt .

Algorithm 1

 1 Splitting Variables Algorithm Require: CLIinit , II and at each clock cycle t number of possible added move mt without atering II Ensure: The set G of all the possible generated CLIs G ← {CLIinit} {G contains the initial CLI} for t = 0 to II -1 do Gt ← G {Gt contains all generated CLIs at time t -1} while mt <> 0 do

	Gc ← ∅ {initialisation of Gc which contains the generated combination C i w t of CLIs (i = 0, mt)}
	for each CLI ∈ Gt do
	for each v ∈ CLI do
	window ← t {variables spanning more than one itera-
	tion}
	while window < E(v) do
	if	SPLIT-VARIABLE(CLI, v, window, Result)
	then
	Gc ← Gc ∪ {Result}
	end if
	window ← window + II
	end while
	end for
	end for
	G ← G ∪ Gc ; Gt ← Gc ; mt ← mt -1
	end while
	end for
	RETURN(G)

TABLE I INITIAL

 I CLIS, THE POSSIBLE SPLITTING CLIS AND THE GENERATED SPLITTING CLIS

TABLE II STATISTICS

 II ON LOOP UNROLLING WHEN AT MOST ONE MOVE IS ADDED PER CLOCK CYCLE

In general, inserting move instructions during register allocation is paired with the split of live ranges to reduce the constraints on register assignement.


O Bénichou

N Meunier

S Redner

R Voituriez

Non-Gaussianity and Dynamical Trapping in Locally Activated Random Walks

Keywords: numbers: 05, 40, Jc, 05, 40, Fb

We propose a minimal model of locally activated random walks, in which the di↵usion coe cient of a one-dimensional Brownian particle is modified in a prescribed way -either increased or decreased -upon each crossing of the origin. Importantly, the case of a local decrease of the motion ability is at work in the process of formation of the atherosclerotic plaque, when describing the dynamics of a macrophage cell that grows when accumulating localized lipid particles. We show in the general case that localized perturbations have remarkable consequences on the dynamics of the di↵usion process at all scales, such as the emergence of a non-Gaussian multi-peaked probability distribution and a dynamical transition to an absorbing state. In the context of atherosclerosis, this dynamical transition to an absorbing state can be viewed as a minimal mechanism leading to the segregation of macrophages in lipid enriched regions and therefore to the formation of the atherosclerosis plaque.

Many-particle systems that consume energy for selfpropulsion -active particle systems -have received growing attention in the last decade, both because of the new physical phenomena that they display and their wide range of applications. Examples include molecular motors, cell assemblies, and even larger organisms [1]. The intrinsic out-of-equilibrium nature of these systems leads to remarkable e↵ects such as non-Boltzmann distributions [2], long-range order even in low spatial dimensions [3] and spontaneous flows [4].

At the single-particle level, the active forcing of a Brownian particle leads to non-trivial statistics. For example, it has been recently shown [5,6] that a random walk which is reset to its starting point at a fixed rate has a non-equilibrium stationary state, as opposed to standard Brownian motion. Another example is given by selfpropelled Brownian particles [7], which can yield sharply peaked probability densities for the particle velocity.

In this letter, we consider a new class of problems in which the active forcing of a Brownian particle is localized in space. While the impact of localized perturbations on random walks has been investigated [START_REF] Hughes | Random Walks and Random Environments[END_REF], in part because of its relevance to a wide range of situations, such as localized sources and sinks [START_REF] Redner | [END_REF]10], trapping [START_REF] Weiss | Aspects and Applications of the Random Walk[END_REF][START_REF] Tejedor | [END_REF] or diffusion with forbidden [START_REF]Markov Chains with Stationary Transition Probabilities[END_REF], hop-over [START_REF] Zia | [END_REF] or defective [15] sites, the role of local activation on Brownian particles dynamics remains open. We present a minimal model of locally activated di↵usion, in which the di↵usivity of a Brownian particle is modified in a prescribed wayeither increased or decreased -at each crossing of the origin. A prototypical example is a bacterium in the presence of a localized patch of nutrients, which enhances the ability of the bacterium to move, or, alternatively, toxins that impair bacterial mobility.

Importantly, this type of localized decrease of the mo-tion ability is also at work when describing the dynamics of a cell (a macrophage e.g.) that grows when accumulating other smaller particles (lipid e.g.) that are localized at some point (Fig. . 1). As the cell size increases its ability to move decreases. Such a situation occurs in the process of formation of the atherosclerotic plaque [16]. The space localization is then due to the presence of lipids at some specific points of the arterial tree, which are determined, as is now well accepted, by the properties of the blood flow. Observations show that macrophages that have accumulated lipids are slowed down and eventually stopped in lipid enriched regions, resulting in the formation of the atherosclerotic plaque [17,18]. Here we propose a simple model accounting for this locally decreased process and address in particular the questions of (i) the potential trapping of cells on locally lipid enriched regions and (ii) the kinetics of the resulting segregation process when it exists.

Our formalism allows us to describe both situations of decreased and increased localized processes. We show that this localized type of perturbation has remarkable consequences on the di↵usion process at all scales. From the theoretical point of view, we stress that the di↵usion coe cient of the active particle at any time depends on the entire history of the trajectory. Thus the evolution of the particle position is intrinsically non-Markovian [START_REF] Kampen | Stochastic Processes in Physics and Chemistry[END_REF][START_REF] Majumdar | [END_REF][21][START_REF] Redner | A Guide to First-Passage Processes[END_REF].

Our main findings are: (i) The probability distribution of the position has a non-Gaussian tail. (ii) For local acceleration, a di↵using particle is repelled from the origin, so that the maximum in the probability distribution is at non-zero displacement. (iii) For local deceleration, a dynamical transition to an absorbing state occurs. For sufficiently strong deceleration, the particle can get trapped at the origin at a finite time. The exact time dependence

(3)

t x (1) (2)
(3)

a b
FIG. 1: a. Sketch of the di↵erent stages involved in the process of atherosclerosis plaque formation: (1) a "free" macrophage cell di↵uses "rapidly"; (2) upon crossing of a localized lipid enriched region, the macrophage accumulates lipids, which makes it grow and slow down; (3) after a sufficient number of crossings of the lipid enriched region, the macrophage eventually gets trapped, resulting in the formation of the atherosclerotic plaque. In this work, we discuss two questions: (i) Is the slow-down in the lipid enriched region efficient enough to eventually induce trapping of macrophage cells ? (ii) If so, what is the kinetics of the resulting segregation process ? b. Sketch of a one-dimensional particle trajectory of the model of locally decelerated random walk.

for the particle survival probability is determined explicitly. Conversely, if the deceleration process is su ciently weak, the particle never gets trapped. (iv) This dynamical transition to an absorbing state provides a minimal mechanism that could help understand the kinetics of formation of the atherosclerosis plaque.

The Model. A one-dimensional di↵using particle is accelerated or decelerated whenever it crosses the origin x = 0 according to the following Langevin equations:

ẋ = p 2D ⇠(t) , Ḋ = f (D) (x) , (1)
where ⇠ is a Gaussian white noise of intensity one, D the particle di↵usion coe cient, x the particle position, (x) the Dirac distribution, and f (D) an arbitrary prescribed function that accounts for the local activation. For simplicity, we assume that the particle is initially at x = 0 with D = D 0 > 0. Note that: (i) both the position x and the di↵usion coe cient D are random variables; (ii) as mentioned previously, the evolutions of x or D alone are non-Markovian; (iii) the function f (D) can be positive (local acceleration) or negative (local deceleration), but with f (0) = 0 so that D remains non negative. Following standard steps, the corresponding Fokker-Planck equation [START_REF] Gardiner | Handbook of Stochastic Methods for Physics[END_REF] for the joint distribution of position

x and di↵usion coe cient D at time t, P (x, D, t), is:

@P @t = D @ 2 P @x 2 (x) @[f (D)P ] @D (t) (x) (D), (2)
where the last term of the right side accounts for the absorbing state at (x = 0, D = 0). The explicit expression for (t) is determined demanding that P is normalized, from which we obtain

(t) = lim D!0 ⇥ f (D)P (0, D, t) ⇤ . (3)
When f (D) is positive, then D is always non zero. In this case, the particle is never trapped and (t) = 0 at all times. While intuitively obvious for local acceleration (f (D) > 0), we show below that (t) can equal zero for local deceleration processes.

Local Acceleration, f (D) > 0. Laplace transforming Eq. ( 2) gives

s b P + D @ 2 b P @x 2 = (x) h @(f b P ) @D (D D 0 ) i , (4)
where b P = b P (x, D, s) is the Laplace transform of the probability distribution. For x 6 = 0, the solution is

b P (x, D, s) = A(D, s) e |x| p s/D , (5)
where the coe cient A(D, s) is determined by integrating Eq. ( 4) across x = 0 to obtain the jump of the first derivative of b P with respect to x at this point:

D " @ b P @x x=0+ @ b P @x x=0 # = @[f b P (x = 0)] @D (D D 0 ) .
Using Eq. ( 5), we have

f @A @D + h f 0 + 2 p sD i A = (D D 0 ). ( 6
)
When f (D) is positive, then A(D, s) = 0 for D < D 0 , while for D > D 0 the solution to ( 6) is

A = B(s) f (D 0 ) f (D) e p 4s F , (7)
with

F (D) ⌘ Z D D0 p D 0 |f (D 0 )| dD 0 .
The unknown function B(s) is determined by the jump of A at D 0 :

A(D + 0 , s) A(D 0 , s) = 1 f (D 0 ) , which finally yields b P (x, D, s) = ⇥(D D 0 ) 1 f (D) e Z p s , ( 8
)
where ⇥ is the Heaviside step function and

Z ⌘ |x| p D + 2F (D) .
Laplace inverting this expression, we obtain the joint distribution

P (x, D, t) = ⇥(D D 0 ) Z f (D) p 4⇡t 3 e Z 2 /4t . (9)
The marginal distribution with respect to x, that is, the probability distribution of positions, is obtained by integrating Eq. ( 9) over all D in the range [D 0 , 1]. While it does not seem possible to evaluate this integral analytically, the large-x behavior can be obtained by the Laplace method. For the illustrative case where f (D) is a constant (that we define as a), this method gives

P (x, t) ⇠ 1 t r |x| 3a exp  8|x| 3/2 9 p a t x ! 1 , (10)
which we numerically checked is close to the exact value P (x, t) over a wide spatial range. We wish to emphasize two important features of this result for P (x, t) that are in marked contrast with the Gaussian propagator of the usual Brownian motion: (i) P (x, t) generally has a non-Gaussian tail; (ii) P (x, t) reaches its maximum at a non-zero displacement. Equation (10) shows that the location of this maximum asymptotically grows as t 2/3 when f (D) = a. Thus local acceleration pushes a di↵using particle away from the origin. From the general expression ( 9), the marginal distribution with respect to D can also be easily obtained by integration over x. We find

P (D, t) = ⇥(D D 0 ) 2 p D f (D) p ⇡t e F 2 /t . (11)
In the particular case of f (D) = a, Eq. [START_REF] Weiss | Aspects and Applications of the Random Walk[END_REF] shows that the di↵usion coe cient of the particle asymptotically grows as t 1/3 . As a byproduct, Eq. ( 11) also provides the distribution of the local time ⌧ (t) spent by the particle in the active zone (the origin for the present case) up to time t. Using the second of Eqs. (1), this basic observable in the theory of di↵usion [START_REF] Majumdar | [END_REF] is related to the di↵usion coe cient at time t by

⌧ (t) ⌘ Z t 0 (x(t 0 )) dt 0 = Z D D0 dD 0 f (D 0 ) . ( 12
)
Thus the distribution of the local time, defined as P(⌧, t), is given by P(⌧, t) = f (D)P (D, t), with P (D, t) given by Eq. ( 11) and D implicitly defined as a function of ⌧ in Eq. [START_REF] Tejedor | [END_REF]. For the illustrative case of f (D) = a, the distribution of the local time at time t therefore is

P(⌧, t)= 2 p a⌧ + D 0 p ⇡t exp ( 4 ⇥ (a⌧ +D 0 ) 3/2 D 3/2 0 ⇤ 2 9a 2 t ) . ( 13
)
This result strongly contrasts with the Gaussian distribution that arises in the case of Brownian motion, which can be recovered from Eq. ( 13) in the limit a ! 0:

P BM (⌧, t) = 2 p D 0 p ⇡t e D0⌧ 2 /t . (14)
Notice, in particular, the typical local time for an accelerated particle with f (D) = a grows as t 1/3 instead of t 1/2 in the case of Brownian motion.

It is worth noting an intriguing dichotomy with a discrete-time version of local acceleration -the "greedy" random walk [25]. In this discrete model, the step length `k after the k th return of a random walk to the origin is given by `k = k ↵ . To match with the continuous model with f (D) = D ↵ , one must choose the value ↵ = 1/2. With this choice, Eq. ( 13) of [25] gives, ignoring all multiplicative factors, P (x, t) / x 1/3 /t exp[ x 4/3 /t], which is di↵erent from (10). The source of this discrepancy is that the probability of being at the origin is not a↵ected by the enhancement mechanism of greedy walks [25], while this return probability is fundamentally modified in the case of locally-activated random walks, as seen explicitly from the distribution of the local time [START_REF]Markov Chains with Stationary Transition Probabilities[END_REF]. Thus our locally activated di↵usion model cannot be viewed as the continuous limit of the greedy random walk.

Local Deceleration, f (D) < 0. Following the same analysis as that used for local acceleration, the Laplace transform of the joint distribution is

b P = ⇥(D 0 D) |f (D)| e Z p s b (s) s (x) (D) , (15)
where b is the Laplace transform of (t) defined in Eqs.

(2) and (3). Using these defining relations for (t), Eq. ( 15) gives b (s) = lim

D!0 ⇥ f (D) b P (0, D, s) ⇤ = e p 4s F , (16)
where we define

F (D) ⌘ Z D0 0 p D 0 |f (D 0 )| dD 0 .
In this result for b , we have used (D)f (D) = 0, since f (0) = 0 by the definition of our model. The important feature of Eq. ( 16) is that b (s) = 0 as soon as F diverges.

Thus our final result is

Here

T (t) = erfc ⇣ 1 p t Z D0 0 p D 0 |f (D 0 )| dD 0 ⌘ (19)
is the trapping probability, namely, the probability that the particle becomes stuck at x = 0 by time t because the di↵usion coe cient has reached zero. As a corollary, the survival probability is given by S(t) = 1 T (t), and we have obtained this quantity for an explicitly non-Markovian process. We also mention that, as in the case of local acceleration, the joint distribution easily gives the marginal distributions of the position and the di↵usion coe cient, as well as the local time.

A fundamental consequence of the local deceleration of a Brownian particle is that two di↵erent dynamical regimes emerge. We illustrate these regimes for the particular case where f (D) = D as D ! 0. If the deceleration is su ciently strong, which occurs when < 3/2, there is a non-zero probability for the particle to get trapped at the origin. More precisely, the survival probability has the asymptotic behavior

S(t) ⇠ 4D 3/2 0 p ⇡t (3 2 ) ! 0 t ! 1 . (20)
Thus in this regime of strong deceleration, the survival probability has the same scaling with time as in the case of a usual Brownian particle in the presence of a perfect trap. In the opposite case of 3/2, then S(t) = 1 for all t > 0 and the particle never gets trapped at the origin. Thus a locally decelerated Brownian particle undergoes a dynamical transition to the absorbing state (x = 0, D = 0) as the deceleration strength increases. Mathematically, this transition occurs at the point where F is no longer divergent.

In conclusion, we introduced a minimal model of locally activated di↵usion, in which the di↵usion coe -cient of a Brownian particle is modified in a prescribed way at each crossing of the origin. In one dimension, a purely di↵using particle hits the origin of the order of p t times after a time t. Consequently, the local activation mechanism is repeatedly invoked during the trajectory of a Brownian particle. Thus the asymptotic dynamics of a Brownian particle is globally a↵ected, leading to markedly di↵erent behavior than that of pure di↵usion. Since the unusual properties of local activation rely on the recurrence of Brownian motion, we anticipate that qualitatively similar, but quantitatively distinct, behavior would arise in two dimensions.

Our model encompasses both the situations where the Brownian particle is locally accelerated or decelerated. For local acceleration, the probability distribution is non-Gaussian and multi-peaked, with maxima away from the origin no matter how weak the acceleration. For suciently weak local deceleration, a Brownian particle manages to avoid getting trapped at the origin in spite of its recurrence. However, for strong deceleration, there is a dynamical transition to an absorbing state in which the particle ultimately gets trapped at the origin.

In the context of atherosclerosis mentioned in introduction, this dynamical transition to an absorbing state can be viewed as a minimal mechanism leading to the segregation of macrophages in lipid enriched regions and therefore to the formation of the atherosclerosis plaque. Interestingly, our model suggests that there exists a critical intensity of the mobility decrease that could be related to the local lipid concentration for the cluster formation to occur. In addition, it could help understand the kinetics of the formation of the atherosclerosis plaque. SR gratefully acknowledges financial support from NSF Grant No. DMR-0906504 and the hospitality of UPMC where this work was initiated.