Christophe Deroulers
email: deroulers@imnc.in2p3.fr

Marine Aubert
email: aubert@imnc.in2p3.fr

Mathilde Badoual
email: badoual@imnc.in2p3.fr

Basil Grammaticos
email: grammaticos@univ-paris-diderot.fr

Modeling tumor cell migration: From microscopic to macroscopic models

Keywords: numbers: 87.18.Gh, 87.10.Ed, 05.10.-a, 87.10.Hk, 87.19.xj, 87.19.lk

published or not. The documents may come

Modeling tumor cell migration: from microscopic to macroscopic

I. INTRODUCTION

The migration of tumor cells plays a principal role in tumor malignancy, particularly for brain tumors. The fact that glioma cells migrate fast contributes greatly to glioblastoma lethality and thwarts therapy strategies based on tumor resection [1]. Understanding the mechanisms of migration can be of utmost importance when it comes to devising efficient treatments [2,3]. In vitro studies combined with mathematical models constitute a promising exploration path, a fact supported by several high-quality results [4,5].

How does one go about modeling cell migration? There are essentially two approaches: on one hand, "macroscopic" models that take the mathematical form of one or several partial differential equations (PDE) [6,[START_REF] Murray | Mathematical biology. II: Spatial models and biomedical applications[END_REF], and, on the other hand, "microscopic" models where individual cells are represented and evolve according to some (stochastic) rules: the so-called individual based models, agent based models, cellular automata, and so on [START_REF] Deutsch | Cellular Automaton Modeling of Biological Pattern Formation: Characterization, Applications, and Analysis[END_REF][START_REF] Hatzikirou | [END_REF][START_REF] Alber | Mathematical Systems Theory in Biology, Communication, and Finance[END_REF][START_REF]Single-cell-based Models in Biology and Medicine[END_REF].

In a macroscopic model, one eschews all reference to the elementary constituent, the cell, and introduces macroscopic quantities, like the average cell density, which evolve according to a PDE. One expects stochastic deviations from the average to be negligible at the scale of a large population of cells, provided this population is homogeneous (quite often, cell populations in cancer are made of cells having distinct phenotypes or genotypes, because e.g., of genetic mutations [START_REF] Nowell | [END_REF][START_REF] Weinberg | The Biology of Cancer[END_REF][START_REF] Bellomo | [END_REF], and this should be taken into account in a more general model of cancer development). As a matter of fact, the first brain tumor models (and the vast majority of models existing today, which are essentially refinements of these first ones) consider the dynamics of a small number of macroscopic quantities [15,16,17]. The tumor cell density is, of course, incontrovertible, the remaining quantities being usually concentrations of chemicals, which, one surmises, play a role in the migration.

Using a PDE has many advantages over a "microscopic" model. First, even with the ultrasimplification of a cellular automaton, one can treat a relatively small number of cells, at best a few thousand, a number very far from the millions one encounters in a tumor lump. Then, there are ready-to-use PDE solvers. Since the PDEs are deterministic, one simulation is enough where one should perform many simulations of the cellular automaton to average over the stochastic noise. The number of discretization steps of the space is fixed only by the mesoscopic geometry of the medium in which the cells diffuse; there is no need to have as many steps as cells. Finally, it is easier to get analytical results from a PDE (even approximate) than from a cellular automaton.

In this spirit, the simplest approach to cell migration is to assimilate the process to linear diffusion. Then a diffusion equation may be used to find how the cell density evolves with time,

∂ρ ∂t = D ∇ 2 ρ. (1)
This was, for instance, the approach of Burgess et al. [18], who modeled the evolution of glioblastomas with the help of a diffusion equation (plus an exponential growth representing cell proliferation). Refinements of the diffusion model can be, and have been, introduced. One can, for instance, consider a space-dependent diffusion coefficient that varies depending on the nature of the tissue in which migration takes place [19,20,21,22]. Terms can be added to Eq. ( 1) representing various chemotactic effects, either as arbitrary external source terms or as quantities that couple Eq. ( 1) to other evolution equations [23,24]. When the cells have a velocity or a persistence of the direction of motion, one has to study (e.g., at short time scales of the cell migration) a different type of equation, like a hyperbolic PDE or a degenerate parabolic PDE [25,26,27,28]. However, modeling of cell migration as simple, heatlike diffusion (noninteracting random walks [29,30]) neglects important phenomena, even in the absence of chemotaxis: cells cannot penetrate each other, and they often interact via contacts with other cells and/or with the extracellular matrix.

But what PDE should one use to take contact interactions into account in the migration process? It can be easier to devise realistic evolution and interaction rules in a cell-based model [START_REF] Hatzikirou | [END_REF][START_REF]Single-cell-based Models in Biology and Medicine[END_REF]31]. Even though there are situations in which it is not possible to discriminate between several macroscopic models, especially when only data for the evolution of the cell density are at hand [31,32,33], it is of little interest to use a PDE that can later prove to be wrong when more data become available and, in particular, when one is able to compare predictions for the trajectories of individual cells [5,31,32].

Sometimes, a microscopic model is derived from a PDE through a discretization procedure [34]. But there are several ways to discretize space that lead to very different behavior of the microscopic agents, some of which can be unrealistic. (For example, consider the discretized population pressure model of [33]: one can show that, for some discretization procedure, cells move toward the high-density region they are supposed to flee, but for some other procedure they move away from it -and both correspond to the very same PDE.) At least, one needs a biological criterion to choose the "right" discretized version of the PDE.

Therefore, to establish a macroscopic model for tumor cell migration with contact interactions via the so-called gap junctions between cells, we start from the microscopic model introduced and successfully compared to in vitro migration experiments in [35]. In this setting, cell proliferation, apoptosis (cell death), and mutations (changes of cell behavior) are negligible (hence the population is homogeneous). The number of cells is constant and cells simply move on a collagen substrate. The model is defined in Sec. II. (There we also discuss special cases in which our model is a kinetically constrained model, borrowing some knowledge from the theory of glasses, supercooled liquids, and jamming systems.)

The derivation of the macroscopic model, the so-called hydrodynamic limit, is carefully explained in Sec. III. The derivation of hydrodynamic limits is a well-known technique in mathematical physics [START_REF] Masi | Mathematical Methods for Hydrodynamic Limits[END_REF][START_REF] Deutsch | [END_REF][START_REF] Kipnis | Scaling Limits of Interacting Particle Systems[END_REF][START_REF] Lebowitz | [END_REF][START_REF] Spohn | Large Scale Dynamics of Interacting Particles[END_REF], and it has also been used in mathematical biology [6,25,26,27,28]. See [START_REF] Turner | [END_REF]42,43,44] for examples in the field of cell movement and cancer modeling. Since our derivation relies on a mean-field like approximation, we review in Sec. III the circumstances under which the approximation fails, and we provide thorough comparisons to simulation results.

The PDE we obtain is a nonlinear diffusion equation. Contrary to what is customary in the phenomenological models to which we alluded above, the terms of the diffusion equation depend crucially on the cell-cell interaction present in the microscopic model. As our analysis shows, these interactions introduce nonlinearities in the diffusion equation. These nonlinearities are essential in order to reproduce the density profiles obtained in in vitro migration experiments (Sec. IV).

II. THE STOCHASTIC CELLULAR AUTOMATON MODEL

In [35], a cellular automaton (lattice gas) was introduced to mimic the migration of cancer cells with contact interactions. The dynamics of migratory motion are described through the evolution of points (representing the cells) on a grid (or lattice), while the interactions are introduced through evolution rules. The impossibility of cells to penetrate each other is trivially taken into account (no two cells can share the same position on the grid), and the scale length of the lattice step is the typical size of one cell.

It is easy to extend such a cellular automaton to cells that occupy several sites on the lattice, using e.g., a cellular Potts model [START_REF]Single-cell-based Models in Biology and Medicine[END_REF][START_REF] Turner | [END_REF]43,44]. But then we would need to introduce many parameters (surface tension, binding energies, ...), whereas experimental data at hand were not sufficient to measure so many values and avoid dubious fits. Thus the model of [35] is an elementary cellular automaton with nonextended cells, involving a single parameter that takes contact interactions into account. This approach has made it possible to model with success the migration of glioma cells over substrates of collagen or of astrocytes and draw conclusions on the existence of homotype and heterotype interactions between the cells [35,45].

In order to define a cellular automaton model, one needs two things. First, one must choose the geometry of the lattice. One may work in one, two, or three dimensions and consider elementary cells of any shape and size. Second, one must specify the update rules of the automaton. With the adequate choice of these rules, one can mimic the interactions between the elementary entities (glioma cells in [35]) which are modeled by the automaton. In this paper, we shall follow the prescriptions we introduced in [35]. We recall them briefly below.

A. The geometry

The cellular automaton introduced in [35] for the description of the glioma cell migration is based upon a regular hexagonal tiling (or, equivalently, a triangular lattice). This choice is dictated by the fact that the hexagonal tiling is the less anisotropic among all regular tilings of the plane. Thus one can have a very simple algorithmic definition of the geometry while keeping the constraints due to the symmetry of the lattice at a minimum. Moreover, in such a lattice each cell is surrounded by six sites, a number sufficiently high to allow a certain freedom of motion at each update. A generalization to three space dimensions is obtained by considering the face-centered-cubic (fcc) lattice or the hexagonal compact (hc) lattice, where each site has 12 nearest)-neighbor sites. It is known that some difficulties arise in threedimensional lattices as compared to two-dimensional lattices when one aims at finding a cellular automaton, the continuous limit of which would be the Navier-Stokes equation [46,47], but here we have no such constraint: instead of starting from a macroscopic equation (such as the Navier-Stokes equation) and guessing what cellular automaton would mimick it, we start from a given cellular automaton and ask for (an approximation to) the corresponding macroscopic equation. The only constraint is that the macroscopic equation be a reasonable model of the behavior of real cells in experiments. Additionally, it turns out that both the fcc and hc lattices, although they do not have the same symmetry properties [79], lead to the same macroscopic equation. Therefore, we expect that refinements in the choice of the three-dimensional lattice would not significantly change our results.

In each case, we denote by a the distance between the centers of two nearest-neighbor sites (lattice step).

B. The stochastic evolution rules

Each hexagon can be occupied by at most one cell at a time. As a consequence, a cell can move only to a free site. The total number of cells is locally conserved (but cells may be introduced into or removed from the system on the boundaries). Thus our automaton is a kind of lattice gas or "box and ball" system [48]. Notice that our model is not a so-called lattice-gas cellular automaton (LGCA), where individual particles have both a position and a speed: here the cells lose the memory of their speed and previous position after each move. In the classification of [30], it is a "position jump process," not a "velocity jump process."

At each update of the automaton (discrete time step), the process of picking a lattice site at random and applying the stochastic evolution rules to this site is repeated a number of times equal to the total number of sites, so that, on average, each site is updated once at each time step [80]. If the site is empty, nothing happens. If it is occupied by one cell, the following rule is applied. It is inspired by what we believe is the biological mechanism of contact interaction. The interaction paremeter p is a fixed, constant number between 0 and 1. We pick at random one of the six neighboring sites, the target site. If the target site is occupied, nothing happens. Otherwise, the move from the cell to the target site (leaving its departure site empty) is done with probability p if at least one of the two common neighboring sites of the departure site and of the target site is occupied, and with probability 1 -p if the two common neighboring sites are empty. With these assumptions, p = 1 means that a cell can only move if it is in close contact with at least one other cell, and it can only move to a position where it will stay in contact with a least one of its former neighbors, while, for p = 0, a cell moves only in such a way that all contacts to its former neighbors are broken. If p = 1 2 , cells are indifferent to their environment (apart from mutual exclusion), and we expect that the cell density obeys a diffusion equation. Values of p > 1 2 model situations in which cells are reluctant to break established junctions with their neighbors, whereas p < 1 2 involves a kind of short-range repulsion but only between cells that were in close contact, unlike what is usually called repulsion between, e.g., electric charges. If these rules seem odd from a physical point of view (they cannot be related to some energy or free energy and convey a strange notion of repulsion, and in particular they do not reduce to a special case of the model of [44]), they incorporate the biologically reasonable assumption that a cell makes its decision regarding where it goes by considering only information it can get from its immediate vicinity, e.g., through established communication junctions with the cells with which it is in contact.

The preceding rules can readily be extended to the three-dimensional hexagonal compact or fcc lattices; the only change is that a pair of a departure and a target site has four neighbors in common instead of two.

The setting considered in [35] consists in a central part of the lattice occupied by the equivalent of the glioma cell spheroid, which, we assumed, acts as a cell source. Once a free position in the hexagons surrounding the center appears, it is immediately occupied by a cell "released from the center." In [35], we have explained how one can (and must) calibrate the model so as to allow quantitative comparisons with experiment.

C. Kinetically constrained models

Our cellular automaton contains as special cases (for p = 1 and 0) two kinetically constrained models (KCMs). These models were introduced in statistical physics to mimic the dynamics of glasses, supercooled liquids, and jamming systems [49,50].

In a KCM, a microscopic object (in our context, cells) can move to a free position only if some requirement about its neighborhood is fullfilled (e.g., only if it has at least another empty position around itself). As a con-sequence, at least for some concentrations of surrounding cells, the motion of an individual cell is not a Brownian motion, as it is in a diffusive model such as the simple symmetric exclusion processes (p close to 1 2 ). However, if one does a "zoom out" or coarse-graining of the trajectory of an individual cell, using characteristic length and time scales (which depend on the concentration of cells) as new unit length and time, the result looks Brownian [51]. Above these scales, a diffusion equation should hold.

The origin of this separation of scales and the existence of the characteristic length and time scales is easy to understand for the so-called noncooperative KCMs (modeling strong glass-former materials). There, cells are blocked unless they are close to some model-dependent defect or excitation that has a finite size. These excitations travel across the system in a kind of diffusive motion, and the configuration of the cells in a region will be changed or relaxed only after the region has been traversed by an excitation. Even though the excitations, and the individual cells while they are close to an excitation, diffuse "quickly," a cell goes "slowly" away from its initial position because its trajectory is an alternation of long nonmoving phases (when the cell is away from any excitation) and short diffusive phases (when the cell is close to an excitation) [51,52]. The characteristic time scale is the typical time between two stays of an excitation in the neighborhood of a cell, and the characteristic length is the typical size of the region explored by a cell during such a stay. Over times much greater than the characteristic time, the trajectory of a cell is made of many nearly independent non-moving/quickly diffusive patterns and, according to the central limit theorem, the cell undergoes a Brownian motion with an effective diffusion constant. Fick's law is obeyed and the coarse-grained model is like a lattice gas without kinetic constraints [50].

In the case of the cooperative KCMs (modeling fragile glass-formers), no traveling finite-size excitation facilitating the motion of cells exists. Instead, the moves of the cells are organized hierarchically and form structures that can be infinitely large [53]: some changes of the system require the participation of an extensive fraction of the cells, hence the term "cooperative". A complete blocking at a finite concentration of cells can even be observed for some models [50]. Still, in the absence of such blocking, it is possible to define a characteristic time and a characteristic scale above which a diffusive behavior is recovered and Fick's law is obeyed [54].

In our cellular automaton, both the noncooperative and cooperative behaviors take place, for p = 1 and 0, respectively (and the behavior for p close to 0, on the one hand, and for p close to 1, on the other, is qualitatively the same as in these two limit situations).

For p = 1, our model is a noncooperative KCM: according to the rules of the automaton, an isolated cell cannot move, but a cluster of two neighboring cells has a finite probability to go anywhere in the lattice and it can be "used" to move any isolated cell (it plays the role of the aforementioned excitations).

For p = 0, our model is a cooperative KCM: this is the very same model as the "2-triangular lattice gas" or (2)-TLG studied in the context of the slow dynamics of glasses [53,54,55,56].

To conclude, for general p, our model interpolates between a cooperative KCM (p = 0), a simple symmetric exclusion process (p = 1 2 ) and a noncooperative KCM (p = 1).

III. THE HYDRODYNAMIC LIMIT

In this section, we explain in detail how we derive an approximate macroscopic, deterministic model (which takes the form of a partial differential equation) from the microscopic, stochastic cellular automaton. This kind of technique has been used in physics since the 1980s to take the so-called hydrodynamic limit of a stochastic model. We recall the basics for readers who are not familiar with it. Other readers may be interested in our discussion of rigorous results and of the quality and limits of our approximation.

The hydrodynamic limit exists for some model defined on a lattice if, as the size of the lattice goes to infinity, one can cut the lattice into parts such that (i) the parts are negligibly small with respect to the whole lattice, (ii) the values that macroscopic quantities (like the cell concentration) take for each part have vanishingly small fluctuations, and (iii) the averages of these values obey some partial differential equation(s). This process is a kind a coarse-graining. Similarly, water, though made of discrete molecules, looks like a fluid composed of homogeneous droplets, because each droplet contains a huge number of molecules that are typically spread homogeneously in the droplet rather than all packed in one-half of the droplet.

The establishment of a PDE for macroscopic quantities after coarse-graining has been extensively studied and can be made mathematically rigorous in a number of cases (see, e.g., [START_REF] Masi | Mathematical Methods for Hydrodynamic Limits[END_REF][START_REF] Kipnis | Scaling Limits of Interacting Particle Systems[END_REF][START_REF] Lebowitz | [END_REF][START_REF] Spohn | Large Scale Dynamics of Interacting Particles[END_REF]). We would like to show that, although one might think that this continuous time and continuous space approximation is useful only in the case of very large lattices, it yields results in excellent agreement with the results from our cellular automaton for lattices as small as 16×16. Indeed, even if in one simulation (or in one experiment) the microscopic concentration of tumor cells is discrete (it can take only two values on each lattice site, depending on whether the site is full or empty), the average of this concentration over independent simulations of the cellular automaton (or over independent experiments) is well predicted by the PDE. Furthermore, there are many ways to perform the approximation, some of which being rather involved, but we show that, here, an elementary procedure that can easily be applied to other cellular automata is sufficient.

A. On the existence and scale of the validity of the hydrodynamic limit

Several authors have proven rigorously that some cellular automata have a nonlinear diffusion equation as their macroscopic scaling limit [57,58,59,60]. In particular, the theorem of [60] proves that our model has a hydrodynamic limit for all values of p between 0 and 1 excluded (we have to exclude 0 and 1 because, for them, some transition rates are 0 and the theorem does not apply). However, it seems rather complicated to get an explicit formula for the corresponding nonlinear diffusion equation from this theorem. On the contrary, our computation provides approximate but explicit formulas that are very useful most of the time.

When p = 1 or 0, one needs to be more careful: then, we have seen in Sec. II that the cellular automaton is a KCM where Fick's law is obeyed only above some characteristic length and time scales.

Very recently, Gonçalves et al. constructed some KCM, the hydrodynamic limit of which would be the porous media equations, and they were able to prove this rigorously [61]. For p = 1, our model is only slightly different from theirs and we believe that generalizing their proof to ours would be possible. This is an indication that (nonlinear) diffusion should hold above some spatial scale, but in practice it is necessary to know how large this scale is. We found that, for p = 1, and more generally for all values of p not close to zero, this scale is of the order of the elementary lattice spacing -the agreement between simulations of the cellular automaton and solutions of the PDE is excellent already at this scale (see later).

For p = 0, we know from the previous section that our model is a cooperative KCM (the "2-triangular lattice gas" [55,56]). To the best of our knowledge, there is no rigorous proof of the hydrodynamic limit in this case or in a similar case, but it is highly probable that this limit exists because the length and time scales we have discussed above exist in particular for this model [54] (they diverge as the cell concentration goes to 1, but are otherwise finite). We found that our approximation breaks down as p gets close to 0, which is a sign that the kinetics cannot be understood only by looking at one cell and its nearest neighbors. Still, p needs to be as small as 0.05 for one to see a significant discrepancy between our approximation and the automaton.

B. Formalism and computation

We now proceed with our computation. On the basis of the previous considerations, this approximate computation should be valid only if the spatial scale above which a macroscopic, coarse-grained diffusive behaviour is of the order of the lattice size -fortunately, this is the case for most values of p.

Our procedure involves first taking averages over the stochastic noise, neglecting correlations, in order to yield deterministic equations on the lattice, then taking the continuous space limit to go from the microscopic to the macroscopic scale. Of course, averaging over the noise amounts to losing information about the true process.

Describing the behavior of a few cells on the lattice by the sole average of the cell concentration on each lattice site would be a crude approximation, of little interest. However, one can hope that this approximation is not so bad in the presence of a large number of cells, and that the typical fluctuations of interesting macroscopic quantities such as the cell concentration become very small as the system grows. Conversely, large fluctuations of the macroscopic quantities would become very rare as the system grows. As we shall see, this does happen for our model and the deterministic, continuous space aproximation is already excellent for rather small systems.

Averaging over the stochastic noise

To perform the averaging and neglect correlations, we use a technique that is similar to the Hubbard-Stratonovich transform and which was already used to study another cellular automaton, the contact process [62]. This technique may be used to build a perturbative expansion. Here, we will actually stop at the first order, since it yields results that are already in excellent agreement with our numerical simulations. Therefore, introducing the technique is not necessary and one can get the same results using other means. But we explain it because its formalism is convenient and compact, especially if one wants to use computer algebra to compute the macroscopic PDE on complicated or highdimensional lattices. In addition, it yields access to systematic perturbative corrections (we can, in principle, get closer and closer to the exact result by keeping more and more terms) on one hand, and to large deviations properties, on the other hand (see [62] for an example where large deviations properties are crucial).

Using the now standard "quantumlike" formalism for master equations, introduced a long time ago [63,64,65,66], we start by representing the probability distributions of the configurations of the cellular automaton as vectors in a simplex. The stochastic rules are encoded as a linear operator acting on this space.

For each site i of the lattice, we introduce two basis vectors |0 i and |1 i to encode the situations in which the site is empty and occupied by one cell, respectively. If the lattice had only one site, say i, only two configurations would be possible ("the site is occupied by a cell" or "the site is empty") and it would be enough to specify the occupation probability ρ of the (unique) site to characterize entirely the statistical distribution of the configurations under some fixed environment. In our formalism, we represent this distribution by the vector v i := ρ |1 i + (1 -ρ) |0 i ; the coefficient of each basis vector is the probability of observing the corresponding con-figuration (there is no need to take the square modulus as in quantum mechanics). If the lattice is composed of two sites, say i and j, there are four configurations and the most general probability distribution can be represented by the vector α |0 i , 0 j + β |0 i , 1 j + γ |1 i , 0 j + δ |1 i , 1 j with α + β + γ + δ = 1 (probabilities must sum up to 1). Here, the four basis vectors are the tensor products of the basis vectors for the site i and of the site j: |1 i , 0 j = |1 i ⊗ |0 j and so on. In the special case in which the sites i and j are statistically independent (in which case the probability that i is occupied and j is occupied is the product of the occupation probabilities of i and j, say ρ i and ρ j , respectively), this vector factorizes as [(1

-ρ i )|0 i + ρ i |1 i ] ⊗ [(1 -ρ j )|0 j + ρ j |1 j ].
More generally, if there are N sites, there are 2 N possible configurations and our vector space has dimension 2 N .

To explain how we encode the stochastic rules, let us take again the simple case in which the lattice has only one site. Let us address first the stochastic process of "radioactive decay": when the site is full, it has a constant probability per unit time α to become an empty site, and, when it is empty, it remains so. The occupation probability ρ i (t) of the site obeys the differential equation

dρ i /dt = -α ρ i (t). ( 2
)
This equation is called the master equation or the Chapman-Kolmogorov equation. The vector that encodes the probability distribution, v i (t), obeys

dv i (t)/dt = α(c i -c + i c i )v i (t), (3)
where c i and c + i are 2 × 2 matrices: in the basis |0 , |1 ,

c i = 0 1 0 0 and c + i = 0 0 1 0 . ( 4
)
The symbols for these matrices are borrowed from quantum mechanics, since c and c + behave like annihilation and creation matrices, respectively (for so-called hardcore bosons). The evolution operator c i -c + i c i of vector v i has two terms: the first one, c i , lets the probability of the configuration |0 increase proportionnally to the probability of the configuration |1 , and, without the second one, probability conservation would not be ensured because the probability of the configuration |1 would be unchanged. Let us give a second example in which two sites i and j are coupled: consider the stochastic process in which a particle is transferred from the site i (if site i is occupied) to the site j (if it is free) at rate κ, i becoming empty and j full. The evolution equation of the vector v encoding the probability distribution of the 4 configurations of the couple of sites reads

dv/dt = κ[c + j c i -(1 I -c + j c j )c + i c i ]v(t). ( 5
)
Here we extend the definition of the matrices. For instance, we understand as c i a 4 × 4 matrix, acting on the whole space of dimension four where v lives, equal to the tensor product of the former matrix c i (which acts on the two-dimensional space of v i ) and of the identity matrix of the two-dimensional space of v j . 1 I is the identity matrix acting on the whole space. The first term of the evolution operator, c + j c i , lets the probability of occupation of j increase proportionally to the probability that i was full and that j was empty. The second term ensures overall probability conservation.

It is easy to generalize this to any stochastic process. If the process has several rules that apply in parallel (e.g., "for any couple of neighboring sites i and j, a particle can jump from i to j if i is full and j empty, making i empty and j full"), the evolution matrix is the sum of the evolution matrices encoding each individual rule [in the example above, it would be the sum over the couples of neighboring sites i and j of the matrices c

+ j c i -(1 - c + j c j )c + i c i ].
It is equivalent to defining a process by its stochastic rules and by the expression of the evolution matrix of its configuration probability vector v.

Applying this technique, we find that our cell migration process has the following evolution matrix:

Ŵ = 1 z i,j n.n. [c + j c i -(1-c + j c j )c + i c i ] F {c + k c k } k n.
n. of i and j .

(6) "n.n." stands for nearest neighbors on the lattice, and z is their number for a given site (coordination number, equal to 6 on the hexagonal tiling and to 12 on the threedimensional fcc or hc lattices). F is a polynom of two or several variables; for each i and j, it is applied to the matrices c + k c k , c + l c l and so on where k, l, ... are the common neighbors of i and j on the lattice. The expression of F is chosen so that, for each possible configuration of the occupation numbers n k , n l , ... of the sites k, l, ..., F (n k , n l , . . .) is equal to the transition rate of one cell from site i to site j (provided that site i is full and site j is empty). It would be easy to study other rules, but, to conform to the model defined in the previous section, we took in two dimensions on the hexagonal tiling

F (n k , n l ) = p (n k +n l ) (3-n k -n l )/2+(1-p)(1-n k )(1-n l ) (7) which is equal to F (n k , n l ) = p (n k +n l -n k n l )+(1-p)(1-n k )(1-n l ) (8)
since n 2 k = n k and n 2 l = n l and, in three dimensions on the fcc or hc lattices,

F (n k , n l , n m , n n ) = p s(5 -s)(s 2 -5s + 10)/24+ (1 -p)(1 -n k )(1 -n l )(1 -n m )(1 -n n ), ( 9
)
where s is short for n k + n l + n m + n n . In both dimensions, we want F to be equal to p if at least one of the nearest neighbors of i and j is full, and to 1 -p if all are empty. The previous expressions are the simplest ones that achieve this. We use the technique of Ref. [62] to derive from the expression of the evolution operator Ŵ the mean-field (correlations are neglected) evolution equations for the occupation probability ρ i (t) of each site i. There are equivalent techniques [32,[START_REF] Turner | [END_REF]43] to get the bare meanfield evolution equations of ρ i (t); this one can additionally yield results for the large deviations and systematic perturbative corrections. The evolution equation of ρ i (t) in the typical (i.e., most probable) configurations reads

∂ t ρ i (t) = -∂ ψi(t) Ŵ , (10)
where ∂ t denotes time derivative, ∂ ψi(t) is the functional derivative with respect to ψ i at time t, taken at ψ i = 0, and Ŵ is the average of the evolution operator. This average is obtained by replacing each 1 I i with 1, each c + i c i with ρ i (t), each c i (t) with ρ i (t) exp[ψ i (t)], and con-

sistently each c + i (t) with [1-ρ i (t)] exp[-ψ i (t)].
Using the expression of Ŵ , Eq. ( 6), one gets the evolution equation of the typical site occupation probability on the discrete lattice,

∂ t ρ i (t) = 1 z j n.n. of i [ρ j (t)-ρ i (t)]F [{ρ k (t)} k n.n. of i and j ] . (11
)

Continuous space limit

So far, we have obtained N evolution equations for the N average occupation probabilities ρ i (t) of the sites from the initial 2 N evolution equations for the probabilities of the 2 N configurations of the lattice. We shall replace the N quantities ρ i (t) by a single field ρ( r, t) such that ρ( r i , t) = ρ i (t) at all times t, r i being the position of the center of site number i. If the lattice step a is finite, there are infinitely many such fields, but to the a → 0 limit there is only one that is regular (twice differentiable). In reality, a is of the order of a single cell's diameter and it is not vanishing (this is a natural length scale in the problem, and it remains relevant after our procedure): taking formally the a → 0 limit means actually that one is interested in phenomena that take place on a length scale of many individual cells, as is standard in statistical mechanics, by restricting to spatial variations of the average cell density ρ that are slow on the length scale a.

For a given site i, we replace each quantity ρ j (t) = ρ( r j , t) with the Taylor expansion of ρ( r j , t) in space around r i : introducing the unit vector from i to j, u i→j := ( r j -r i )/a, we can write

ρ( r j , t) = ρ( r i , t) + a u i→j • ∇ρ( r i , t) + (a 2 /2) u i→j • H( r i , t) u i→j + O(a 3 ) (12)
where the center dot denotes the usual scalar product, ∇ is the gradient operator, and H is the Hessian matrix of ρ (here taken at the position r i and at time t). After this substitution, all terms of order 0 in a vanish because there is locally conservation of the number of cells (no proliferation nor apoptosis), and all terms of order a vanish because the evolution rules are symmetric (induce no bias between left and right) and the lattice is reflectioninvariant. Also, the lattice is translation-invariant and the sites are coupled all in the same way (to their nearest neighbors), thus each evolution equation yields the same relation between the field ρ and its space and time derivatives. Since we are interested in the continuous space (a → 0) limit, we disregard terms of order 3 and above in a, and after some algebra we find the evolution equation of the site occupation probability,

∂ t ρ( r, t) = a 2 div[D (ρ( r, t)) ∇ρ( r, t)], (13)
where

D(ρ) = (1 -p)/4 + (2p -1)ρ(1 -ρ/2)/2 (14)
on the hexagonal tiling in two dimensions (2D) and

D(ρ) = (1 -p)/6 + (2p -1)ρ(4 -6ρ + 4ρ 2 -ρ 3 )/6 ( 15
)
on the fcc lattice in 3D. For the hc lattice, we find the very same expression of D(ρ) as for the fcc lattice. The expression of D(ρ) is the value of the function F [Eqs. ( 8) and ( 9)] when all its arguments are equal to ρ, up to a geometrical factor. In the a → 0 limit, the righthand side of Eq. ( 13) vanishes, leaving the trivial equation ∂ t ρ( r, t) = 0: this is a sign that nothing interesting happens, in the limit of continuous space, at the time scale of the updates of individual sites, t. The only time scale leading to a nontrivial equation has a unit time of the order of 1/a 2 [81]. This scale is characteristic of diffusion (in the case of an advective phenomenon, we would have chosen the time scale 1/a [START_REF] Lebowitz | [END_REF]). Equation ( 13) is a diffusion equation with a concentration-dependent diffusion coefficient D(ρ). This equation belongs to the family of the porous media equations [START_REF] Vázquez | The Porous Medium Equation[END_REF]. A striking difference with the heat equation (or with Fick diffusion) comes from the possibility that D vanishes at some site occupation probability ρ, which may create abrupt fronts in the diffusion profile instead of large, Gaussian tails.

To completely free ourselves from the lattice, which in most practical applications does not exist, but was convenient to define the model, we write down the equation for the spatial concentration of cells c instead of the site occupation probability ρ. This is done be dividing ρ by the volume of a single site, √ 3a 2 /2 for the hexagonal tiling, and a 3 /(3 √ 6) for the fcc lattice. Because ρ and c are proportional, we can discuss either one in the rest of this Section; when we discuss experimental results, we shall use only c.

Before proceeding to a comparison of this result to numerical simulations, let us discuss the vector "density of cell current," . We expect it to satisfy Fick's law, with a nonlinear diffusivity. To compute , we first compute the net cell current along the lattice link i → j, i.e., the current from site i to site j,

c + j c i -c + i c j = a z F ({ρ( r, t), . . .}) u i→j • ∇ρ( r, t) + O(a 3 ), (16)
where r := ( r i + r j )/2 and all arguments of F are equal to ρ( r, t). Then we add the contributions of all links (on the hexagonal tiling, there are three types of contributions since the links can have three directions), and we go to the density of current by multiplying each contribution by the density of links with that direction on a hyperplane orthogonal to that direction. On the hexagonal tiling, if some links are, say, vertical, then we count the density of intersections of a horizontal line with links. These intersections have a periodic pattern and can be grouped four by four; a group has the length √ 3a and contains two vertical links (each of these contributes one), and two links making an angle π/3 with the vertical direction [each of these contributes cos(π/3)], hence the factor √ 3/a. We find on the hexagonal tiling

( r, t) = - 2 √ 3 D[ρ( r, t)] ∇ρ( r, t). ( 17
)
Putting this expression into the conservation equation for the number of cells (there is neither apoptosis nor proliferation in our model),

∂ t c( r, t) = -div ( r, t), ( 18
)
and using that c = 2ρ/( √ 3a 2 ) (on the hexagonal tiling), we find again Eq. ( 13).

C. Comparison to steady-state simulations of the automaton

In order to test the analytical results ( 13) and ( 14), we did some simulation of the cellular automaton in a simple geometry. More realistic simulations are the topic of the next Section. Since, as will be discussed later on, the equilibrium state of the cellular automaton is trivial (all configurations of the tumor cells are equally probable), we chose a setting where a nonequilibrium steady state of the tumor cells can exist: a cylinder made of n L + 1 rings of cells (its total length is (n L + 1)a √ 3/2) and base circumference H = n H a connected to a reservoir full of cells (ρ = 1) at the left end and to an empty reservoir (ρ = 0) at the right end, as shown in Fig. 1 (n L and n H are two integers). The first and the last ring of cells of the cylinder belong to the full and to the empty reservoirs respectively; hence the effective length over which diffusion takes place is L := n L a √ 3/2. We chose a cylinder rather than a simple rectangle in order to make the boundary effects in the direction perpendicular to the flow small, if not vanishing. We introduce two coordinates on the cylinder: x is measured in the direction parallel to the axis and y along the shortest circles. To simulate the reservoirs, on the boundaries of the system, cells behave under special evolution rules: each time one of the sites on the left border of the lattice is left by a cell, it becomes occupied again by a new cell with no delay, and each time a cell arrives on a site of the right border of the lattice, it is removed at once. In this geometry, if we let the system evolve starting from a random configuration, it relaxes to a steady state with a permanent current where gains and losses of cells from and to the two reservoirs compensate in average (but the instantaneous total number of cells still fluctuates as time advances). On an infinite lattice, the steady state might be reached only after an infinite duration. In practice, we simulate in parallel two families of independent systems, one where the site occupation probability is initially close to one and one where it is initially close to zero, and we stop our simulation when the values of the macroscopic observables (averaged over the systems of each family) are undistinguishable up to the error bars. The number of systems in each family is chosen so that the error bars are shorter than the precision that we request in advance.

In the steady state, the concentration profile can be computed in the macroscopic limit thanks to Fick's law (17) or to the PDE [START_REF] Weinberg | The Biology of Cancer[END_REF]. In the stationary regime (where no macroscopic quantity depends on time), any initial asymmetry has been "forgotten" by the system and ρ( r) does not depend on y because the system is translation-invariant in the y direction. Furthermore, ρ( r) is such that the current of cells is uniform: , the density of cell current, is independent of x and y. Then Fick's law (17) reads

(2/ √ 3)D[ρ(x)] ∂ x ρ(x) = -j. ( 19
)
With the boundary conditions ρ(0) = 1 and ρ(L) = 0 FIG. 2: (Color online) Dependence on the interaction parameter p of the stationary profiles ρ(x/L) of the site occupation probability along the cylinder between the full reservoir (left, x = 0, site occupation probability ρ = 1) and the empty reservoir (right, x = L, ρ = 0) -see Fig. (1). Since the cell concentration c is proportional to the site occupation probability ρ, the profiles of c would be the same. For each value of p, we plot the simulation results for the cellular automaton (circles) and our prediction from the analytic approximation Eq. ( 21) (solid lines). From bottom to top, p = 0.01, 0.05, 0.2, 0.4, 0.7 and 1. The error bars are smaller than the circles' size. The cylinder is made of 64 × 64 sites except for p = 0.01 and 0.05 where there are 512 × 512 sites.

and the expression for D(ρ) above, we find in 2D

j = (1 + p)/(6 √ 3L) (20)
and

x(ρ) = {1 + ρ[(2p -1)ρ(ρ -3) -3(1 -p)]/(p + 1)}L. ( 21
)
Then one can get the curve for ρ(x) by plotting the parametric curve (x(ρ), ρ). In Fig. 2, we plotted the stationary profile of the site occupation probability cell ρ on the cylinder, between the two reservoirs, as obtained from simulations of the cellular automaton and from the macroscopic PDE. Notice the effect of nonlinear diffusion: for p = 1/2 (not shown in Fig. 2), the profile is linear as predicted by the usual Fourier or Fick law, but, for other values of p, interactions between cells make it more complicated. For most values of p, there is an excellent agreement between the simulation results and the analytical approximation for all lattice sizes we tried (starting from L = 16). However, we can distinguish two situations in which the microscopic and the macroscopic models disagree. The first one is related to a boundary effect: at the vicinity of the empty reservoir when p = 1 (respectively at the vicinity of the full reservoir when p = 0), a clear discrepancy of FIG. 3: (Color online) Finite-size effects for the stationary profile of the cell concentration between the two reservoirs. Main curve: dots: profile ρ(x/L) of the site occupation probability from simulation results for p = 0.01 and system sizes nL = nH -1 = 15, 31, 63, 127, 255, and 511 (from bottom to top). It shows a strong finite-size dependence (for p = 0.05, not shown, a weaker finite-size dependence is found). The solid line is the analytic approximation, clearly in disagreement with the simulation results. Bottom inset: for p=0.01, the extrapolation of the occupation probability ρ at position x = L/4 to an infinite lattice shows that the discrepancy between the simulation data (dots) and the analytical approximation (dashed line) is not due to finite-size effects. Indeed, there is convergence toward a value of the site occupation probability that is clearly below the analytic approximation. Top inset: profile ρ(x/L) for p = 1, with numerical data for the sizes nL = 15, 31, and 63 only (dots) and analytic prediction (solid line). The finite-size effects are very weak and there is a good agreement between simulation data for large systems and the analytic prediction. This is true more generally for all values of p away from 0 (say, larger than 0.1).

the average occupation probability ρ of the last-but-one lattice site (respectively of the first few lattice sites) is seen between the cellular automaton results and the corresponding value from our analytic formula. The second situation is what happens when p gets close to 0; see, e.g., p = 0.01 and 0.05 in Fig. 2. Here, the concentration profiles of the finite cellular automata show a strong finite-size dependence (main plot of Fig. 3, compared to the top inset of Fig. 3), and they seem to converge to some shape that is clearly different from the profile predicted from the PDE. To exclude the possibility that the simulation results agree with the analytic approximation at very large sizes, but not at the sizes we have simulated, we plot in the bottom inset of Fig. 3 an extrapolation to infinite systems: it shows that the discrepancy should persist even for huge systems. These phenomena are confirmed in Fig. 4, where the permanent current of cells FIG. 4: (Color online) Stationary current of cell density between the full and the empty reservoirs, across the cylinder (geometry defined in Fig. 1), as a function of the interaction parameter p. To be able to compare different system sizes (bottom to top : nH -1 = nL = 15, 63 and 255), we actually plot the current density j multiplied by the system length L (dots). For nL = 255, only the points that differ significantly from the other sizes are represented. The error bars are smaller than the symbol's sizes, except for nL = 255, where they are drawn. Solid line: analytic approximation, Eq. ( 20) times L. Like for the profile of the cell concentration, there is a large disagreement between the analytic approximation and the simulation results when p is close to 0 and a small disagreement when p is close to 1.

across the cylinder is plotted as a function of p: one sees a good agreement for most values of p, with a slight deviation and practically no finite-size effects around p = 1, and some more important deviation, together with strong finite-size effects, around p = 0. This will be understood in the next subsection.

D. Discussion of our approximations: The occurrence of correlations

Our approximate computation of the hydrodynamic limit fails if the correlations between the occupation of neighboring sites are not negligible or if the continuous space limit has singularities (e.g., an infinite gradient) that cannot exist on the lattice.

Singularities in continuous space

Singularities are easy to predict and to address; they could happen for our model since the diffusion coefficient vanishes (for p = 1 and ρ = 0 on one hand and for p = 0 and ρ = 1 on the other). Indeed, Fick's law reads

( r, t) = -D[ρ( r, t)] ∇ρ( r, t), (22)
and at a point where D vanishes while  is finite (which happens for p = 1 or 0 between the two reservoirs of the previous subsection), the concentration gradient is infinite. To get the correct, physical behavior if necessary, it may suffice to reintroduce two or a few lattice sites around the point where the singularity appears and to solve separately the discrete equations (obtained after averaging but before the Taylor expansion) on these sites, and the PDE on the rest of the lattice.

In our case, this simple procedure does not enable one to quantitatively (or even qualitatively) understand the discrepancies between the automaton and the approximate results close to the reservoirs. This is because these discrepancies have to do with correlations.

Correlations

If i and j are two sites, let us consider the so-called connected correlation function [START_REF]This name comes from the fact that these correlation functions have connected diagrams in a diagrammatic expansion[END_REF] of i and j,

ρ i (t)ρ j (t) c := [ρ i (t) -ρ i (t) ][ρ j (t) -ρ j (t) ] = ρ i (t)ρ j (t) -ρ i (t) ρ j (t) . ( 23
)
This quantity is zero if the statistical distribution of the occupation numbers of the sites i and j is the same as if the sites i and j would be filled independently at random with given mean occupation probabilities. More precisely, in the case of a finite lattice of N sites with a fixed total number of cells, this definition must be replaced with

ρ i (t)ρ j (t) c := [ρ i (t)-ρ i (t) ][ρ j (t)-N ρ j (t)-1 /(N -1)] (24
) if one wants to detect correlations due to the interactions between cells that could be "hidden" behind the contribution -ρ(1 -ρ)/(N -1) due to the constraint that the total number of cells is fixed.

Our approximate computation of the hydrodynamic limit makes use of an assumption of statistical independence ("all connected correlations functions are zero"), and we found in our simulations that the correlations are indeed small (they vanish as L → +∞) in almost all cases.

A situation in which we can compute exactly the correlations is the equilibrium state of a translation invariant lattice. Since, in our stochastic automaton, any movement of a cell may be reversed and since, if the lattice is translation invariant, any transition from one configuration of the cells on the lattice to another has the same rate, 1/(6N ), the detailed balance condition [49] is fulfilled, the equilibrium state exists, and the probability distribution of the configurations is uniform. Then the connected correlation functions (with the definition FIG. 5: (Color online) Schematic explanation of the correlations that appear for p close to 0 near the full reservoir (left) and for p close to 1 near the empty reservoir (right). Left: when two cells sit next to one another close to the full reservoir and if p is close to 0, they prevent one another from moving in one of the two directions that would be allowed if they were alone (black arrows represent permitted moves). Thus accompanied cells stay longer close to the reservoir than single ones, and positive correlations occur (see text). Right: when one cell sits close to the empty reservoir and if p is close to 1, it cannot move unless it has a neighboring cell. Thus single cells stay longer close to the empty reservoir than accompanied ones, and negative correlations occur (see text).

above, relevant to the case of conserved total cell number) between any two sites are zero. We simulated the cellular automaton in this situation and we checked that our numerical data are compatible with vanishing correlations (up to the error bars) for several values of p and of the cell concentration; this is a check of correctness of the simulation program. In this case, our analytical approximation is exact, but the result is of course trivial (ρ is uniform).

Let us come back to the nonequilibrium steady state between the two reservoirs. In such a case, it is well known that correlations do exist, can be quantitatively important, and, above all, are generically longrange [START_REF] Garrido | [END_REF]69]. For our automaton, we found numerically that, for p close to 1 (0) there exists a significant negative (positive) correlation of the occupation probability of adjacent sites in the last (first) row of sites, parallel to the reservoir. This is why the macroscopic limit disagrees, in these situations, from the simulation. In other situations, they decrease at least as fast as 1/L as the system length L goes to +∞. Since studying these correlations in more detail would go beyond the scope of the present article, we limit ourselves here to explaining qualitatively what happens close to the boundaries. The correlation there can be interpreted using the microscopic rules (see Fig. 5). If p = 1, when two cells are close to the empty, absorbing reservoir (Fig. 5 right), sooner or later one of the two cells will fall into the reservoir, leaving the other one alone. This other cell will not move until another cell comes on a neighboring site. Therefore, isolated cells stay longer close to the reservoir than accompanied cells, and there is more chance to find an empty site close to a cell than one would expect if cells were distributed at random, hence the negative correlations. Far away from the reservoir, the cell concentration is higher and there is always some cell to "unblock" an isolated cell, so the same phenomenon cannot be observed. Conversely, if p = 0, when two cells are close to the full reservoir (Fig. 5, left), they are in some sense in each other's way and they may progress away from the reservoir only in one of the two directions they could use if they were alone, the other direction being forbidden by the kinetic rule of the automaton since they would stay in contact with a cell they were already in contact with. If one could consider only these cells on the nearest line of sites from the reservoir and forget about the cells that lie on the next-to-nearest line of sites, one could conclude that accompanied cells stay longer close to the reservoir than single cells, that there is more chance to find another cell close to a cell than one would expect if cells were distributed at random, and that the correlation is positive. This is not so simple because the next-to-nearest line of sites has many cells that can compensate this effect, and we do observe that correlations extend over a few lines of sites away from the full reservoir, but the qualitative idea is correct. Far away from the full reservoir, the site occupation probability is not close to 1 and there is always some hole to "unblock" a cluster of cells, hence no such correlation is observed. We expect that taking into account these correlations into the analytic computation would improve the agremeent between simulations and analytic formulas. Now let us come to the disagreement of the simulated and predicted concentrations of cells in the whole system when p comes close to zero. Our interpretation is the following. We know from a preceding discussion that, at p = 0, the cellular automaton is a cooperative KCM for which Fick's law is obeyed only above some spatial scale that diverges with the cell concentration. Since, between the two reservoirs, all concentrations can be observed, Fick's law might never hold even in arbitrarily large systems, but this is beyond the scope of the present paper. If p is close to zero but strictly positive, we expect that the threshold spatial scale for effectively Brownian diffusion should be large (diverging with p) but finite at all values of the cell concentration. Therefore, and contrary to what happens for larger values of p where the threshold scale must be comparable to the lattice spacing, there is a whole range of system sizes where the Brownian regime is not yet reached and the concentration profile depends strongly on the system size. Above this scale, diffusion sets in but, as a consequence of the cooperative nature of the kinetics, correlations always play a significant role and our evaluation of the diffusion constant, which neglects them, is irrelevant, hence the disagreement observed even at large system sizes between the predicted profile and the actual one.

IV. RESULTS FOR CELLS DIFFUSING OUT OF A SPHEROID

In this section, we come back to the radial geometry, where cells diffuse away from a center (cf. Sec. II B), in order to compare the results from the macroscopic model to experimental results. We use the expression of the diffusion coefficient obtained in Sec. III B 2, cf. Eq. ( 14), for a 2D system on a hexagonal lattice. We calculate numerically the solution of the nonlinear diffusion equation in a radial geometry by using a fully implicit method where the spatial derivatives in the diffusion equation are evaluated at time step t+ 1 (for the sake of stability) [START_REF] Press | Numerical Recipes: The Art of Scientific Computing[END_REF]. The site occupation probability ρ is constrained to be equal to 1 in the center at all times.

A comparison between the density profiles from the automaton and from the numerical solution of the diffusion equation at the same time is presented in Fig. 6. As in the case of diffusion between a full and an empty reservoir, the agreement is very good except for p = 0.05, close to the full reservoir, where the numerical solution does not agree with the density profile of the automaton as nicely as in the cases of larger p. This effect is due to correlations between close neighbors which are not negligible when p is close to 0 (cf. Sec. III D 2). In general, the agreement between both profiles is nice. This is a supplementary check of our derivation of a nonlinear equation from the automaton. This good agreement allows us to make a direct comparison between the numerical solution of the nonlinear diffusion equation and experimental cell densities. The experiments consisted in placing a spheroid of glioma cells (GL15) on a collagen substrate and following the evolution of the migration pattern during a few days. A detailed protocol has been described in [35]. Experiments have been performed by Christov at the Henri Mondor Hospital. Figure 7 is a photograph of an experimental spheroid, 36 h after cells started migrating. In order to compare the experimental density profiles at different times of evolution (cf. the full lines in Fig. 8) and the numerical solution of the diffusion equation, we rewrite Eq. ( 14) as (25) where D * is the value of D when p = 1 2 (i.e., without interactions between cells). Note that D * = 1 8 when the time step of the numerical solution coincides with that of the automaton (cf. Sec. III B 2). For comparison with experimental data, the numerical time step is chosen equal to 1 h (and thus D * can be different from ). The spatial length step is set to 35.2 µm as in the automaton, which corresponds to a characteristic size of cells. The diffusion equation was integrated for p = 0.95. The choice p = 1 induces sharp edges at low densities that do not correspond to our experimental profiles.

D(ρ) = D * [2(1 -p) + 4(2p -1)ρ -2(2p -1)ρ 2 ],
We find that, for the coefficient D * = 1240 ± 100 µm 2 h -1 , numerical solutions at 12 and 36 h agree well with the experimental data at the same times, cf. Fig. 8. One can notice that this value falls within the range of values of the diffusion coefficient values obtained in vivo [START_REF] Swanson | [END_REF] (between 4 and 88 mm 2 /d = 450 -10000 µm 2 h -1 ), although this comparison must be taken with a grain of salt due to the different nature of the systems. For this type of in vitro experiments, where cells are followed individually, the automaton remains the best way to model the system: time is counted as the number of cells that have left the spheroid, the lattice step corresponds to the mean size of a cell, and interaction between cells is represented by rules defined between neighboring cells. The only parameter to be fixed is p. On the contrary, for real tumors where the number of cells is of the order of 10 7 , the automaton becomes unmanageable, and in this case, the approach through a diffusion equation is more efficient.

V. CONCLUSION AND OUTLOOK

In this paper we set out to establish a link between a microscopic description of tumor cell migration with contact interactions (which we developed in [35]) and a macroscopic ones, based on diffusion-type equations used in most phenomenological models of tumor growth.

A microscopic description has obvious advantages. First, it allows a straightforward visualisation of the results: one "sees" the cells migrating. Moreover it can be directly compared to experiments and thus reveal special features of the migration like, for instance, cell attraction which we postulated in [35] in order to interpret the experimental cell distribution. It allows predictions for the trajectories of individual cells that can be observed in time-lapse microscopy experiments [32,72,73]. However the microscopic approach is not without its drawbacks. In all implementations we presented we have limited ourselves to a two-dimensional geometry. While a three-dimensional extension is, in principle, feasible it would lead to a substantial complication of the model. Moreover, even in a 2-D setting, it is not computationally feasible to accomodate the millions of cells present in a tumor. Our typical migration simulations involve at best a few thousand cells.

A macroscopic treatment on the other hand deals with mean quantities, like densities. Thus limitations such as the number of elementary entities (cells, in the present work) simply disappear. However the trace of microscopic interactions equally disappears unless it is explicitly coded into the macroscopic equations. This is a very delicate matter: the present work is dedicated to bridging just this gap between the microscopic and macroscopic approaches.

Our starting point was the microscopic description of cell migration which was validated through a direct comparison to experiment. Starting from a cellular automaton (with a given geometry and update rules) we have derived a nonlinear diffusion equation by taking the adequate continuous limits. Our main finding was that the interactions between cells modify the diffusion process inducing nonlinearities.

Such nonlinearities have been observed in models of ecosystems where individuals tend to migrate faster in overcrowded regions [START_REF] Murray | Mathematical biology: I. An Introduction[END_REF][START_REF] Okubo | Diffusion and Ecological Problems[END_REF]. This mechanism (fleeing overcrowded regions) has also been the subject of speculation in the context of cell migration [START_REF] Witelski | [END_REF], but, in our experimental situation, it is irrelevant: on one hand, cells in overcrowded regions cannot migrate fast, even if they want to, because they are fully packed (thus blocked); on the other hand it has been shown experimentally that the nonlinearity is due to attractive interactions between cells (when interactions are suppressed, the nonlinearity diminishes) [45].

While the present model is most interesting, it is fair to point out the precise assumptions that entered its derivation. We were based on a specific cellular automaton: many more do exist, the limit being only one's imagination. Of course for the problem at hand it was essential to have a geometry allowing every site to possess a sufficiently large number of neighbours so as not to artificially hinder the cell motion. Moreover the update rules were chosen in fine through a comparison to experiment. But it would be interesting to see how robust the nonlinear diffusion coefficient is with respect to changes of the geometry: introducing a random deformation of the lattice or mimicking the migration of cancer cells on a susbtrate of astocytes in the spirit of [45]. Moreover, if one forgets about the application to tumor cell migration, there may exist several interesting cases of cellular automaton geometry and kinetic rules leading to nonlinear diffusion equations worth studying, in addition to the kinetically constrained models we have seen and that have their own interest for the physics of glasses. We intend to come back to these questions in some future publications.

Finally, an incontrovertible generalization of the present work would be that of three-dimensional models. While the extension of microscopic models to three dimensions may be computationally overwhelming, the solution of three-dimensional diffusion equations, be they nonlinear, does not present particular difficulties. Thus once the automaton rules (mimicking the cell-cell interaction) have been used for the derivation of a diffusion equation like Eq. ( 13) with coefficient (15), we can proceed to the study of the macroscopic model, coupling diffusion to proliferation (plus, perhaps, other macroscopically described effects) for the modeling of real-life tumors (in particular, glioblastomas). This is a path we intend to explore in some future work.

FIG. 1 :

 1 FIG.1:(Color online) 3D view of the first lattice geometry we used to test our results: a cylinder (lattice with periodic boundary conditions in one direction) connected to a full reservoir of cells (a source) on its left end and to an empty reservoir of cells (a sink) on its right end. Here, the system is made of nL + 1 = 16 rings of sites (including the two rings in the reservoirs), and each ring is made of nH = 16 sites. After some transient regime, a steady current of cells establishes along the cylinder. A typical configuration of the cells in the steady state for interaction parameter p = 0.8 is shown.

FIG. 6 :

 6 FIG. 6: Density profiles from the automaton (dotted line) and from the numerical solution of the diffusion equation (full line) for p = 1, 0.5, and 0.05. Error bars for the automaton density profile are smaller than the markers and thus are not represented.

FIG. 7 :

 7 FIG. 7: Experimental pattern of migration after 36 h. The real size of the image is 1.90 × 1.77 mm.

Acknowledgments

We acknowledge financial support from the "Comité de financement des théoriciens" of the IN2P3 (CNRS).


Alexandre Gaffet
email: alexandre.gaffet@vitesco.com

Nathalie Barbosa Roa
email: nathalie.barbosa.roa@vitesco.com

Pauline Ribot
email: pauline.ribot@laas.fr

Elodie Chanthery
email: elodie.chanthery@laas.fr

Christophe Merle
email: christophe.merle@vitesco.com

A Hierarchical XGBoost Early Detection Method for Quality and Productivity Improvement of Electronics Manufacturing Systems

This paper presents XGBoost classifier-based methods to solve three tasks proposed by the European Prognostics and Health Management Society (PHME) 2022 conference. These tasks are based on real data from a Surface Mount Technologies line. Each of these tasks aims to improve the efficiency of the Printed Circuit Board (PCB) manufacturing process, facilitate the operator's work and minimize the cases of manual intervention. Due to the structured nature of the problems proposed for each task, an XGBoost method based on encoding and feature engineering is proposed. The proposed methods utilise the fusion of test values and system characteristics extracted from two different testing equipment of the Surface Mount Technologies lines. This work also explores the problems of generalising prediction at the system level using information from the subsystem data. For this particular industrial case: the challenges with the changes in the number of subsystems. For Industry 4.0, the need for interpretability is very important. This is why the results of the models are analysed using Shapley values. With the proposed method, our team took the first place, capable of successfully detecting at an early stage the defective components for tasks 2 and 3.

INTRODUCTION

The 2022 PHME Data challenge encourages participants to solve multiple classification problems for a real production line from Bitron Spa. The dataset includes data from Solder Paste Inspection (SPI) and Automatic Optical Inspection (AOI) equipment of a real industrial production line equipped with automated, integrated and fully connected machines (Industry 4.0). A detailed description of the dataset is given in 3. The challenge is to design an algorithm to predict test labels for the components. Specifically, the goal is to develop a hierarchical classification predicting: 1. whether the AOI classifies the component as defective; 2. in the case of a defect, the label applied by the operator; 3. in the case of confirmation of the defect by the operator, the repair label.

To tackle this challenge, we pursue the following steps: data exploration and domain knowledge extraction, data cleaning, data preparation (normalization and encoding), data modelling (model training and validation) and results in analysis. The four last steps were made recursively while trying different approaches, as shown in Section 4. Data exploration allowed us to identify three main issues within the given dataset: missing information, highly imbalanced classes (for all tasks) and high cardinality on the categorical features. The latter is not necessarily an issue but implies that a special treatment needs to be done to these features a priori. We will elaborate on the issues in Section 4.1 and on the categorical encoding in Section 4.2.

To solve each task, we take different information units formed by feature tuples, corresponding to different levels on the data hierarchy. At the same time, different features are kept as relevant for each task and followed by specific normalization or encoding. The specific tools used for each task are described in detail in Section 4.3. Finally, after presenting the experimental setup used to tune the model hyperparameters (Section 4.4), the achieved results are discussed in Section 5.

RELATED WORK

Several scientific articles already present machine learning applications for Surface Mount Technology production lines. [START_REF] Richter | On the development of intelligent optical inspections[END_REF] proposes a convolution neural networks deep learning application working on the AOI system to automatically detect defects. In [START_REF] Tavakolizadeh | Industry 4.0: mining physical defects in production of surface-mount devices[END_REF], some binary classifiers are tested to detect defects inside products using simulated production data. These data are simulated from several SMT lines and give a good classification score. In [START_REF] Parviziomran | Data-driven prediction model of components shift during reflow process in surface mount technology[END_REF], a component shift prediction method is proposed to predict the shift of the pad during the reflow process. In [START_REF] Park | D 3 pointnet: Dual-level defect detection pointnet for solder paste printer in surface mount technology[END_REF]) SPI data are used to predict at an early stage the potential defects of the prediction. This work is based on a dual-level defect detection method. In [START_REF] Jabbar | Supervised learning approach for surface-mount device production[END_REF]) some treebased machine learning methods are used to predict the defects found in AOI using SPI data. [START_REF] Gaffet | Data-driven capability-based health monitoring method for automative manufacturing[END_REF] proposes an unsupervised univariate method for monitoring the In-Circuit Testing machine (located at the end of the Surface Mount Technology lines) and components. Another large topic of interest for this study is prognosis and health management at different levels: systemlevel or sub-system level. In our case, we have to use information from a pin level to retrieve the health at a system level which is the product component. This topic is linked to the decentralized diagnosis approach. [START_REF] Zhang | Decentralized fault detection for a class of large-scale nonlinear uncertain systems[END_REF] proposes a decentralized model-based approach with a simulation example of automated highway systems. [START_REF] Ferdowsi | A decentralized fault detection and prediction scheme for nonlinear interconnected continuous-time systems[END_REF] proposes a decentralized fault diagnosis and prognosis methodology for large-scale systems adapted for aircraft, trains, automobiles, power plants and chemical plants applications. [START_REF] Tamssaouet | Online joint estimation and prediction for system-level prognostics under component interactions and mission profile effects[END_REF] proposes a component interaction-based method to provide the prognosis of multi sub-system model.

DATA DESCRIPTION

The PHME provides the dataset used in this article as part of the 2022 conference data challenge (PHM Society, 2022). The dataset includes measurement information from two different steps of the PCB production (see Figure 1). The first step is the SPI, in which each solder pad is checked to verify its compliance and, accordingly, a sanction is generated depending on the quality of the solder (evaluated on several physical aspects). The second step is the AOI. In addition to checking the finished solder pads, their position, shape, etc., this process also inspects the component itself, looking for defects like missing or misaligned components. The AOI inspection has two types of sanctions: the automatic sanction provided by the machine itself and the one given by an operator that verifies the first one in case of spotted defects. For the AOI, each line corresponds to a unique entry of the set {PanelID, FigureID, MachineID, ComponentID, PinNumber}, where the PinNumber can be fill as NaN. On such cases, we believe the AOI does not inspect the solder paste but the component itself.

As introduced before, the challenge is divided into three tasks.

Task 1 is to predict whether or not a component will be classified as defective by the AOI using only the inputs provided by the SPI, i.e. the pad measurements. Task 2 is about predicting the operator's label using the SPI test results and the automatic defect classification provided by the AOI (AOIlabel). Finally, task 3 concerns the reparation operation. Again, the objective is to predict whether the component detected as a defect by the operator can be repaired or not. For this task, the information used for the prediction is the SPI test result, the AOIlabel and also the OperatorLabel.

METHODOLOGIES

Challenges

The exploratory analysis of the training data has revealed several issues that need to be tackled to solve correctly the different tasks:

1. Missing values: for three different PanelID, the proposed data missed some information in the SPI dataset. We choose to exclude the lines with no information. 2. Class imbalance: the number of components detected as defects by the AOI is much lower than the number of components from the SPI dataset. Similarly, the number of components really classified as a fault by the operator is much lower than the number of components classified as correct.

High cardinality of the categorical features: the categorical features

XGBoost algorithm and categorical features encoding

For tabular data applications, Gradient Boosting Decision Trees (GBDT) are widely used, being XGBoost [START_REF] Chen | Xgboost: extreme gradient boosting[END_REF],

CatBoost [START_REF] Prokhorenkova | Catboost: unbiased boosting with categorical features[END_REF] and LightGBM [START_REF] Ke | Lightgbm: A highly efficient gradient boosting decision tree[END_REF] the algorithms with the best results. Among these algorithms, we decide to use XGBoost, which is a scalable, paralleled and distributed implementation of the original gradient boosting tree algorithm.

GBDT is an ensemble model algorithm, i.e. it combines several decision trees to perform a better prediction than a single model. XGBoost uses, in particular, the idea of boosting: it uses a collection of weak models to generate a strong model. In practice, for XGBoost, the idea is to use a gradient descent algorithm over a cost function to iteratively generate and improve weak models. On each iteration, a new weak decision tree is generated based on the error residual of the previous weak model. The final prediction is a weighted sum of all the iterated weak trees. Among ensemble methods, boosting can minimize the model's bias. We propose one model for each task. In this section, these models as well as the used features are presented.

XGBoost is a very performing algorithm, although some caution has to be taken when categorical features are used. This is true for all tree-based or boosted tree methods. In particular, one hot encoding can lead to very poor results when the categorical features have many levels. Indeed, a large number of levels leads to sparsity as for each level, a new variable is created. These new variables have only a small fraction of data points that shall have the value 1 and the other the value 0, which is a problem for tree-based methods because tree split searching for the purest nodes. Indeed, a hot encoded variable is not very likely to lead to the purest nodes if it is very sparse. That is why, the tree split will not be done using this one hot encoded variable. Even if the original categorical feature has a lot of importance for the prediction. In our cases, other encoding techniques should be used.

First, a common technique is the Hash encoding. It is already present in our dataset with a numerical value for each PadID level. The PadID level values depend on the FigureID of the pad. Here the Hash encoding is done with only one feature but in general, it could be encoded into more features. One of the most used hashing methods is described in (Yong-Xia & Ge, 2010).

The next approach is the frequency-based encoding method: it is based on the frequency of the levels as the label value. If the frequency is linked to the target, it will help the prediction of the variable. For instance, in tasks 2 and 3, the frequency of one component is probably linked to the issues that can exist for each component. This encoding can be useful in that case.

Finally, the last type of encoding is label-based. The idea of label encoding is to replace each categorical value with the conditional probability of the class to be predicted by knowing the categorical features. This can be done by several methods such as Leave One Out Encoding, and CatBoost en-coding [START_REF] Prokhorenkova | Catboost: unbiased boosting with categorical features[END_REF]. The main issue with this method is to learn the conditional probability without overfitting. It can be realized by not taking the observation into account in the learning of the probability for each observation, as in the Leave One Out Encoding. This can also be done more efficiently using CatBoost encoding. For this case, we found that the CatBoost encoding performs best.

The frequency-based encoding and Hash encoding have less success.

Solving the challenges

Task 1

The first task is the most challenging of all. The main difficulty arises from the fact that some defects are related to the pin, and others to the component itself. In our opinion, the most important question for each task is "Should we predict (model) by component or by pin ?". For this first task, we decide to go for a per-pin prediction. This is mainly guided by the difficulty to generate coherent labels and features at the component level for this task. Indeed, only one pin can have an issue. It does not seem right to affect the same label to a component with only one pin detected as a defect by AOI equipment and another component with multiple pins with defects. Moreover, the number of pins varies too much depending on the studied component. As a result, any aggregation of continuous variables will probably hide important information if only one pin has a defect. For instance, the aggregation with the mean of the Volume will not contain a lot of information if there are many pins, and only one defective pin for the considered component.

For each tuple, we want to predict if the tuple is detected as a defect by the AOI equipment or not. Accordingly, the training target column is 1 if the tuple appears in the AOI dataset and 0 otherwise. It is worth noting that we are not considering as defective tuples for which only PanelID, FigureID, ComponentID appear with PinNumber = NaN. Both, categorical and continuous features, are used as input. We use the following continuous variables: Volume(%),Area(%),OffsetX(%), OffsetY(%),Shape(um),PosX(mm), PosY(mm),SizeX, SizeY that we simply call "numerical features" in the following of the article. As a categorical value, we only keep ComponentID that we encode using a CatBoost encoding method [START_REF] Prokhorenkova | Catboost: unbiased boosting with categorical features[END_REF].

Task 2

For this second task, we first split the AOI dataset into two parts according to whether or not PinNumber = NaN. In the case where PinNumber is specified, we can join the AOI and the SPI easily using the columns PanelID, FigureID, Compon-entID, PinNumber in each dataset. From the joint dataset, we use the "numerical features" from the SPI test results as de-fined in task 1. For the categorical features, we keep three features: AOILabel, ComponentID and FigureID ComponentID where the later is the string concatenation of the variables with the same name. For these features, we use a CatBoost encoding based on the categorical encoders python library. We believe a better result can be achieved with deeper work on the optimization of the encoder hyper-parameters. Finally, we also create two new meta-features (not encoded) Count Pin Component and Count Pin Figure . These two variables are counting the number of pins detected as a defect by the AOI respectively for the component and figure of the tuple PanelID, FigureID, ComponentID, PinNumber. The XGBoost classifier algorithm classes each tuple into the class "Bad" or "Good" of the OperatorLabel target column.

For the AOI defect without any PinNumber associated, we propose to also use an ensemble of categories and continuous variables. We use the same categorical and meta features, while for the numerical features we use the mean values per component of the following variables: Volume(%), Area(%), OffsetX(%), OffsetY(%) to keep the information only on the PanelID, FigureID, ComponentID tuple level.

Finally, we also use an XGBoost classifier algorithm to class each tuple PanelID, FigureID, ComponentID, PinNumber, with PinNumber referenced as NaN as described before.

For the final sanction (that must be given at the PanelID, FigureID, ComponentID three-tuple level), we use the following rule: if one of the four-tuple PanelID, FigureID ComponentID, PinNumber entry is predicted as "Bad", then the associated three-tuple will also be considered as "Bad". If not, the label "Good" is assigned to the three-tuple.

Task 3

Task 3 is about the prediction of one categorical value presented in the AOI dataset RepairLabel. This label can take two values, FalseScrap or NotPossibleToRepair. For this task, we tried the same approach as in task 2 predicting each four-tuple PanelID, FigureID, ComponentID, PinNumber and merging the results per component, but the approach was not successful. To improve the result, we choose to do the prediction for the PanelID, FigureID, ComponentID threetuple directly. As input, we use the same method as for task 2, grouping the SPI values per three-tuple using the mean as an aggregation method for the "numerical features". The categorical variables used are ComponentID, FigureID Component ID. As before, these variables were encoded using the Cat-Boost encoding method.

The meta-features generated are Count Pin Component, Count Pin Figure and Count Pin Panel, respectively, the number of pins detected as defects by the AOI per component, figure and panel. We also created one-hot encoded features from the AOILabel. For each labelled type of error, the asso-ciated feature has the value of 1 if the three-tuple is detected as having this error by the AOI machine (in at least one pin) and 0 otherwise. To predict the class of each tuple, we use the XGBoost classifier.

Hyper-parameter tuning

The XGBoost model has been tuned using the Optuna python library [START_REF] Akiba | Optuna: A next-generation hyperparameter optimization framework[END_REF]. This package is an optimization framework that searches for the best hyper-parameters of the space defined by the user. It uses some distributed computation and early stopping to improve the speed of the solution. It is implemented with various optimization algorithms. In our case, we use the sampling algorithm Tree Parzen Estimator Sampler (TPES) [START_REF] Bergstra | Algorithms for hyper-parameter optimization[END_REF]. It is a sequential modelbased optimization. As a bayesian optimization algorithm, it computes a probability model of the optimization function and selects the best hyper-parameter according to this probability model and the real cost result. For each step, the tuning is done using the mean of the F1-score of a 4 cross-validation method. The dataset is split into four parts, each of these parts is recursively the testing dataset while the other is used for training the model. For the sake of reproducibility, hyperoptimization techniques and training algorithms are available in [START_REF] Gaffet | Phme data contest 2022[END_REF]. Our prediction can probably be easily improved with more iterations of the tuning phase, as we did not spend much time on it.

RESULTS AND DISCUSSION

The results obtained for the three tasks are detailed in this section. Table 2 shows the F1-score for each task of the challenge. The training and testing set results are close, showing good generalization capability. It seems that our models avoid overfitting issues that are difficult to handle with this dataset. Indeed, the imbalance issue and the fact that some variables such as ComponentID have a lot of importance can lead to large bias. ComponentID could be thought of as a domainspecific key factor, if, for example, depending on the type of component, the reparation is possible / allowed or not. For instance, if there is an issue with a microchip, this is far more difficult and costly to solve than an issue with a simple resistance. The following most important features are the metafeatures being in the order of importance more critical the number of defective pins per figure, then component and final panel. This also seems correct as the number of issues increases the potential damage to the product. The "Missing" AOILabel (one-hot encoded) also has a considerable impact. Maybe it is not possible to replace a missing component due to oven operation.

Score

Feature importance

SHAP values

The interpretation of the machine learning model described as a black block model is a really important topic in the industry. Indeed, for the acceptance of a model, it is mandatory to explain the model decision to the process experts. The interpretation allows validating the model by comparing the model and experts' explanation of a phenomenon. It also allows recommending some repair actions to the experts. SHAP (SHapley Additive exPlanation) interpretation [START_REF] Lundberg | A unified approach to interpreting model predictions[END_REF] is based on the game theory [START_REF] Štrumbelj | Explaining prediction models and individual predictions with feature contributions[END_REF]. The idea is to compute an interpretation value called SHAP value and denoted ϕ j for all variables and each sample of the dataset. The output of the model is described by the sum of the SHAP values ϕ j .

ϕ j = S⊆J\{j} |S|!(M -|S| -1)!(f (S ∪ {j}) -f (S)) M ! (1)
where M is the number of variables, J is the ensemble of variables, f the model output and j the variable index. SHAP is an additive method. The output of a classifier can be described as the cumulative sum of the impact of all variables. This result was expected as the more the AOI equipment found defective pins, the more likely a real defect in the component is. The impact of continuous variables is very low for this task. For Task 3, the output value 1 is associated with the "Not Possible to Repair" label and the output value 0 with "FalseScrap". The SHAP values indicate that the more the number of defective pins per component per figure, the more likely the product is impossible to repair. If a component or a pin is missing, then the product seems also more likely to be not possible to repair. The same result is found for "Unsoldered" AOILabel level. For the other level of AOILabel, the results are unclear.

Issues with Task 1

In this subsection, the issues encountered in the prediction of task 1 are detailed. As previously reported, the results of the prediction of task 1 are highly biased by one component "BC1". Table 3 and Table 4 present the confusion matrix results of the training of task 1. It can be seen that most of the correctly detected defects by the classifier model of task 1 is coming from the results of the component "BC1". Only five other tuples are correctly classified for the defects. Obviously, this is a concern for our method. Another observation is that the F1-score obtained with the results coming from Table 4 is 0.76 whereas the F1-score obtained considering all the tuples as defects is 0.75. Both scores are very close and depending on the cross-validation random selection, the simple rule considering all the tuples of the the"BC1" component as defects can be better than a more advanced classifier. To solve this issue, more advanced analyses are required. From our experience, it seems easier to predict the defects associated with the "Unsoldered" AOILabel than the other. The first task results are probably too biased to be used in production. However, the results of task 2 and task 3 are promising for a production application.

In future work, it could be interesting to add some information to improve the prediction. The time between operations for instance is missing (probably to not share cycle time information). A product with a larger cycle can lead to an issue with dust deposit, humidity change... This information could improve the performance of all tasks. Some extra information like the temperature curves of the oven can also help to better predict future issues. Finally, for task 1, unsupervised monitoring methods may be a better solution due to the class imbalance issue.

Figure 1 .

 1 Figure 1. Surface Mount Technologies Production Line. (PHM Society, 2022)

  PadID have more than 1000 modes. Without any sort of variable encoding, classifiers are very difficult to use for such variables. PadID is already encoded in the sort because the values are ordered by FigureID. In a way, the variable is encoded by component area. 4. High bias in continuous features: the continuous features such as volume, area, height... are highly correlated to the categorical feature PadID. 5. Level of prediction: the prediction has to be done at a component level, whereas the available data are given at a PadID. This leads to a lot of issues in creating the target training and prediction. Indeed, for instance, it is unclear if the target training has to be created by component or by pad. 6. Different numbers of pins: the number of pins depends on the component. It is difficult to use all the pin results as input of a classifier for each component, as the number of pins and therefore features varies. The generalisation of the training depending on the component is very difficult.

Figure 2

 2 Figure 2 presents the feature importance of the classification model for task 1. The most important features are the component's position on the panel and the encoded ComponentID variable. Because almost 30 per cent of the defects found in the AOI dataset are coming from one component ("BC1"),

Figure 3 .

 3 Figure 3. Task 2: feature importance of the XGBoost classifier for the four-tuple (AOI defect is linked to a pin)

Figure 3 and

 3 Figure 3 and Figure 4 present the feature importance for task 2 on the four-tuple and three-tuple classifiers respectively. The position and the ComponentID are also important for this task, but the continuous variables Volume(%), OffsetX(%), OffsetY(%) and Shape(um) seem to have also a great impact on the prediction. Even more, the generated meta-features Count Pin Component, representing the number of pins with a defect per component, and Count Pin Figure, representing the number of defective pins per figure, have also a large impact. Actually, these two features allow improving the F1score for this task by almost 0.2.

Figure 5

 5 Figure5shows the feature importance of the XGBoost clas-

Figure 6 ,

 6 Figure 6, Figure 7, Figure 8, Figure 9 present the computed impact of the features on the model output. For task 1, a

Table 1 .

 1 Summary characteristics of the datasets

	Feature	SPI	AOI
	Number of lines	5 985 382	31 617
	Number of panels	1924	1 924
	Number of components	129	102
	Figures/panel	{1,8}	{1,. . . ,8}
	Components/panel	{128,129}	{2,. . . ,27}
	Lines/panel	{3112,. . . } {2,. . . ,203}
	Simple data exploration allows to discover anomalous entries
	and clean the datasets. A summary of the information found
	in each dataset is shown in Table 1. The cleaned SPI dataset
	is composed of 1921 panels, each with 3112 entries. A panel
	is an ensemble of 8 PCBs, also called figures, grouped. Each
	figure is composed of 129 components. The component refer-
	ence is found on the feature ComponentID. These components
	have several pins that can be identified using the PinNumber
	feature. The SPI test results provide the volume, area, height,
	size, shape and offset of each PadID as well as a final re-
	sult flag. A PadID corresponds to a unique combination of
	{FigureID, ComponentID, PinNumber}. Each panel provided
	by the competition has at least one component detected as a
	defect by AOI automatic sanction. Each line of the datasets
	describes a PadID of one electronic board.

Table 2 .

 2 F1-score for the three tasks with training and testing data set

	Dataset	Task 1 Task 2 Task 3 Score
	Training 0.43	0.66	0.90	0.66
	Testing	0.41	0.67	0.77	0.62

Table 3 .

 3 Confusion matrix for all training sets

	Predicted Label
	Defects Not Defects

Table 4 .

 4 Confusion matrix for the training set observation of the component "BC1" larger value of the CatBoost encoded variable ComponentID logically results in larger model output (the encoding is done by replacing the level value with its conditional probability of target). The variable SizeX does not seem to be connected to the model output. A more interesting larger value of the continuous variable of Volume and

Offset seems to be representative of a defect for task 1. For task 2, a low number of defective pins per component and per figure is more likely to be a "Good" OperatorLabel.

Table 5 .

 5 Task 1 score results for different AOILabel levels6. CONCLUSION AND FUTURE WORKTo solve the different challenges in the Printed Circuit Board production, we proposed three different methods. These methods use the XGBoost classifier and are based on variable encoding and feature engineering. We have shown that generating meta-features representing the circuit state at different levels (component, figure and panel) improved the generalization of the classifiers. The use of different encoding techniques for the categorical features has shown CatBoost as the most promising one, due to the high cardinality issues.

	AOILabel / Score Score Task 1
	Coplanarity	0.00
	Translated	0.00
	Soldered	0.00
	Unsoldered	0.55
	Size	0.00
	LeanSoldering	0.27
	Misaligned	0.15
	Missing	0.00
	Jumper	0.00