Soledad Berríos
email: sberrios@med.uchile.cl

Julio Lopez Fenner

Aude Maignan
email: aude.maignan@univ-grenoble-alpes.fr

Julio López Fenner

J López Fenner

Simulating aggregates of bivalents in 2n=40 mouse meiotic spermatocytes through inhomogeneous percolation processes

Keywords: MSC 2010 subject classifications: Primary 60K35; secondary 60C05, 62P10, 92B05, 92-08. Inhomogeneous percolation, Fullerene, Mus M. Domesticus 2n=40 meiotic spermatocytes

come

Simulating aggregates of bivalents in 2n = 40 mouse meiotic spermatocytes through inhomogeneous site percolation processes

Introduction

The question about randomness in biological processes constitutes a source of many interesting cross-disciplinary discussions among biologists and mathematicians. More often than not, uniform probability distributions seem to constitute the very paradigm of randomness, while the appearance of other probability distributions are often considered to be the result of dynamical rules applied to uniform random variables.

On other side, it has been frequently discussed whether biological models at the scales considered in this work, in the range of 1 to 10 µm, should be better approached from the continuous rather than the discrete side of modeling.

The challenge of producing a mathematical model of randomness for explaining the (apparent) bivalent's association clusters observed in the Meiosis's Pachytene for the house mouse 2n = 40 -Mus m. domesticus spermatocytes was initiated in [START_REF] Berríos | Model of chromosome associations in Mus domesticus spermatocytes[END_REF] and continued in [START_REF] López-Fenner | Bivalent associations in Mus domesticus 2n= 40 spermatocytes. Are they random?[END_REF], in an eminently discrete approach, by considering first surface hexagonal tilings as approximation of the nuclear envelope and next by a discrete regular graph with a finite number of vertices. The reason of this sort of approach is simple: Objects being modeled are distinct entities and the sought probability distributions referred to the (discrete) partitions of a number. Hence, more than absolute positions of points in a sphere, the mathematical/biological constraints for the stochastic distribution of a finite number of indistinguishable objects, related to each other by a neighborhooddependent condition, lead us naturally to a description using graphs, with satisfactory results [START_REF] López-Fenner | Bivalent associations in Mus domesticus 2n= 40 spermatocytes. Are they random?[END_REF].

In this work we propose an enhanced discrete model along the lines of the earlier works, but in which the nuclear envelope is now modelled by the dual C 1200 of a fullerene, [START_REF] Weisstein | Fullerene MathWorld-A Wolfram Web Resource[END_REF], which is a particular planar graph embedded in the sphere. On this graph the objects of interest are percolation domains obtained by a suitable 'Inhomogeneous Percolation Process', which we call P-percolation (see [START_REF] Grimmett | Inhomogeneous bond percolation on square, triangular and hexagonal lattices[END_REF]), that are put into correspondence with the units of biological interest (bivalents, or, more precisely, CTC's, see below). Induced clustering generated by overlapping domains can be now explored and analyzed using the model.

Hence, the article is structured as follows: Section 2 contains the biology context, in which a minimal set of notions and working definitions (necessary for understanding the underlying biological problem) are presented. It includes some new derived statistics from the original data set, that were omitted in the previous works, which we use as test frames for our model. The interested reader and the biological audience will find the complete theory with details in the references.

In section 3 we argue against using the theory of random graphs G(n, p), as introduced by Erdös [START_REF] Erdös | On random graphs[END_REF], for explaining the clustering probabilities obtained in the dataset. We show that there can be no value for the probability p that adequately fit or reproduces the data.

Section 4 deals with the actual construction of the graph used to model the nuclear envelope. It will turn out to be the dual of fullerene C 1200 (subsection 4.1). This graph appears naturally by imposing size and connectivity constraints upon the number of vertices and planarity of the graph, that correlates well with both, the size and form of the biological objects being modeled. The construction of random domains of vertices via P-percolation, (called 'ranches') that represent the CTC's (subsections 4.2 and 4.3), coupled with a simple Laplace dynamics for the parallel evolution of the positions of vertices with attributes (subsection 4.4), produce a structure with which the probability distribution of the association clusters is correctly reproduced. Finally, subsection 4.5 presents the results from the simulation process.

Section 5 hands the word back to the biologist by presenting a biological discussion and interpretation of the mathematical and simulation results obtained.

The Biology Behind

Meiosis is an extraordinary process that produces haploids and genetically diverse gametes in organisms of sexual reproduction and differs significantly to the more widely known process of Mitosis. It is characterized by a plethora of stages, of which one of the conspicuous, remarkable ones is called the pachytene in early prophase I, which we are interested in. At pachytene, the homologous chromosomes synapse along (by or through) a proteinacious structure called synaptonemal complex (SC), which enables recombination between them, a process that is known to produce genetic variation. Chromosomes in synapse are called bivalents and describe arcs 'floating' inside the nucleous, bound to the nuclear envelope by both their extremes, called telomeres.

Unions or associations among bivalents, mediated by overlapping domains of their own constitutive heterochromatin, are frequently observed during this stage.

Enter the house mouse: 2n = 40 Mus m. domesticus. This species, which is widely represented in Europe, exhibits all their chromosomes (save for the X and Y), to be of similar morphology and size, hence making them particularly well suited for establishing a reference biological subject with which it should be possible to contrast mathematical models. In our case, we will only consider the 19 bivalents of 2n = 40 Mus m. domesticus (hence leaving the X and Y out) and will consider them to be indistinguishable units.

The centromere -telomere -complex (CTC) is a structure of the short arm of each one of the 19 autosomal bivalents and -under appropriate staining techniques, as immunocytochemical staining -several structures of them can be brought to light. So, for example the Synaptonemal Complexes (SC's) and also the constitutive pericentromeric heterochromatin domains (CPCH's).

Spreads are obtained as the result of removing the nuclear envelope and projecting the spermatocyte nucleus's to a flat 2D surface. A typical spread is presented in Figure 1. a), in which the pericentromeric heterochromatin domains (CPCH's) are stained red, while the Synaptonemal Complexes (SC's) are stained green. Notice the apparent clustering induced by the overlapping domains of CPCH. These clusters will be interpreted as domains in which the corresponding CTC's are in asso-ciation and, in order to describe the frequencies of their appearance in the data, they will be associated 1:1 to a partition of the number 19 (see [START_REF] Berríos | Model of chromosome associations in Mus domesticus spermatocytes[END_REF]):

The 1 st -Class (or simply 'Class') of a spread is a characterization of the induced partition given by the number of bivalents present in the biggest association cluster observed in it. The 2 nd -Class is the size of the second biggest subcluster observed, and so forth.

So, for example, if a given spread exhibits one association cluster of size 4, three of size 2, and all others in singletons, then the associated partition of 19 is written as 19 = 4 + 2 + 2 + 2 + 1 + • • • + 1 and we say that the spread belongs to 1 st -Class 4, 2 nd -Class 2, etc.

In order to make this article as self-contained as possible, we reprocessed the original data from [START_REF] Berríos | Model of chromosome associations in Mus domesticus spermatocytes[END_REF] and reproduced the empiric 1 st -Class and 2 nd -Class distribution, which we completed by adding the intrinsic variation in each class, computed by simple bootstrap.

The data is given in Table 1, results are reported in Figure 1,b), where the shaded areas represent two times (the square root of) the corresponding bootstrap variances. For each spread in the data set, the corresponding partition was determined and the frequencies of their appearance computed. Figure 2. a) depicts the full set of frequencies, where the x-axis denotes the set of all partitions ordered canonically, i.e., the first partition is 19 = 1+1+1+• • •+1, the second one is 19 = 2+1+1+• • •+1, until partition number 489, which is 19 = 18 + 1 and the last one, 19 = 19, the 490th partition. One cannot help but notice a suggestive fractal-like form for the distribution showing up. Remark: The reader will notice a small abuse in notation while using the letter n: Each time we refer to the mouse under study, we shall write 2n = 40, since this is standard notation in biology: 2n = 40 Mus M. Domesticus. In contrast, whenever we talk about placing the 19 autosomal bivalents (not considering the X and Y) upon the nuclear envelope and determining the corresponding association clusters, we shall put n = 19, since it is the partitions of 19 we will be referring to.

Random Graphs

In this section we show that the observed partitions do not follow the distribution obtained when applying the theory of random graphs G(n, p).

Indeed, the exact probability of G(n, p) being connected, which we denote by f (n), can be computed recursively by means of [START_REF] Gilbert | Random Graphs[END_REF]

f (1) = 1 f (n) = 1 - n-1 i=1 f (i) n -1 i -1 (1 -p) i(n-i) , n > 1. (3.1)
The probability P(r) for a random graph composed of n nodes to be split into k components of size r = (r 1 , r 2 , ...r k ) is known to be [START_REF] Allard | Exact solution of bond percolation on small arbitrary graphs[END_REF]:

P(r) = n! Π m d m !Π i r i ! Π l f (r l )Π i< j (1 -p) r i r j ,
in which d m is the number of containers of size m, which essentially is a combinatoric result obtained in the determination of the so called 'Bell's numbers' (see [START_REF] Berend | Improved bounds on Bell numbers and on moments of sums of random variables[END_REF], [START_REF] Rota | The number of partitions of a set[END_REF]). So we ask ourselves whether there is a value p, for which the observed partition frequencies corresponds to the theoretical frequencies of the cluster distributions of G(19, p).

To answer this we take, for example, 1st-Class 4, of which there are 93 observations in the dataset, corresponding to 28 different partitions.

Under this model, the probability of partition r

= (4, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1) in G(19, p) is given by P(r) = a -a (1 -p) 3 -b p(1 -p) 4 -b 1 -(1 -p) 2 -2 p(1 -p) 2 (1 -p) 3 p 3 (1 -p) 162 ,
in which the coefficients are, respectively, a = 290990700, b = 872972100.

Solving for p such that r∈Cl(4) (P(r)) = 93/400, the observed frequency for 1st Class 4, yields three possible values: p = -.121, p = -0.0125, p = 1.807. Hence we conclude that no p will ever furnish a correct approximation for the observed distribution using G(n = 19, p).

For this reason, we now take the nuclear architecture of the spermatocytes into consideration and propose a random partition model that suits better the experimental data.

Model Construction

We take first into account the scales of the biological units considered in our study (as reported in [START_REF] Berríos | Model of chromosome associations in Mus domesticus spermatocytes[END_REF]): It is known that approximately 1µm 2 portion of the CTC's surface adheres to the internal face of the nuclear envelope, whose size -in turn -range between 10µm -15µm in diameter. This means that, for example, in a nucleus of 14µm diameter, the total surface of the nuclear envelope measures approximately 615µm 2 , hence accepting a tiling containing 615 cells of 1µm 2 each, while in a nucleus of 13µm diameter, a similar tiling would contain 530 cells of 1µm 2 .

We propose a discrete structure (a connected undirected graph G = (V, E), with V a set of vertices (or nodes) and E a set of edges, E ⊆ V × V) for modeling the nuclear envelope. The difference with our previous model (reported in [START_REF] López-Fenner | Bivalent associations in Mus domesticus 2n= 40 spermatocytes. Are they random?[END_REF]) is that now the bivalent's structures attached to the envelope are spatially extended structures which will be represented as attributes to the vertices, and we replace the assumption of a rigid CTC by a random structure representing two of the highlighted components in the spreads: the Constitutive Pericentromeric Heterochromatine (CPCH) and the Synaptonemal Complex (SC).

To capture the biological constraints posed above, we search for a graph G satisfying (at least) the conditions: C1. Regular and maximally connected. C2. Embedded on a sphere. C3. |V| ∼ 600.

Condition C1 represents a biological constraint, since the SC's are spacially extended objects placed upon a surface, which imposes naturally a six-degree of regularity upon the nuclear architecture, considered as a graph. To respect the observed sizes of the nuclei, we require C3 and, while a six-regular embedding of a surface to a sphere is certainly not possible, the constraint of being maximally connected yields a candidate: The dual fullerene.

Dual Fullerene

Embedding a graph in the plane is equivalent to embedding it on the sphere, hence we use planar graph theory [START_REF] Trudeau | Introduction to Graph Theory[END_REF]. We need the following results, which we state without proof: A graph satisfying conditions C1 and C2 above with exactly N -12 vertices of degree 6 and 12 vertices of degree 5 can be constructed in a straightforward manner: Indeed, from the previous theorem and corollary, e ≤ 3N -6. If we consider N -12 vertices of degree 6 and 12 vertices of degree 5, then e = 1 2 (6(N -12) + 5 × 12), which is exactly 3N -6.

A Fullerene is a cubic planar graph having all faces 5-or 6-cycles, i.e., a 3-regular planar graph such that all the face's sizes are pentagons or hexagons [START_REF] Weisstein | Fullerene MathWorld-A Wolfram Web Resource[END_REF]. Hence the graph we consider is a fullerene dual and is constructed as follows: Start with a node, attach to it a ring made of 6 connected vertices, then a second ring of 12 vertices, etc . . ., until attaching the i th ring with 6i vertices, see Figure 3. After having built this structure with i rings, which we call a semi-sphere, the graph has exactly N i = 1+3i(i+1) vertices, all of them of degree 6, save those that are located at the exterior ring, which are of degree 3 or 4.

Construct now a second half-sphere with i -1 rings and connect with the first one as to leave only 12 vertices of degree 5, all other being of degree 6.

The total number of vertices is therefore N i +N i-1 = 1+3i(i+1)+1+3i(i-1) = 2+6i 2 and all vertices of degree 5 are located on the i th ring of the first semi-sphere.

For i = 10 we finally obtain a graph of 602 vertices, which matches our expectations for the size of the nucleus (see Figure 3). Of them, 590 vertices are of degree 6 and the rest (twelve) are of degree 5. It is C 1200 , the Dual of the Fullerene C 1200 of 1200 vertices. The presence of these twelve degree 5 vertices will prove to be non significant for the computations/simulations, in the sense that without loss of accuracy, the planar graph can be assumed to be plain 6-regular.

Random Chromatin Neighborhoods

Now we turn to the CTC, which will be characterized by two of their main constituents: The SC and the CPCH. It is a 3-d structure that attaches to the nuclear envelope. Hence the contact zone will be described by the position of the SC, i.e., a vertex in C 1200 , and the positions of the surrounding CPCH, described by vertices in the neighborhood of the SC.

Let G = (V, E) be a graph, v an arbitrary but otherwise fixed vertex. Let N G (v) be the set of neighbors of v in G, i.e., the set of vertices that are connected to v by an edge in E. For v, w ∈ V, let d G (v, w) the distance between them to be defined as the length of the shortest path connecting them.

For each vertex v we assign a state: state(v), which belongs to the set {¬Ex, S , ¬S }, in which ¬Ex stands for "non explored", S for "selected" and ¬S for "not selected".

Let P = (P 1 , P 2 , . . . , P κ ) ∈ [0, 1] κ be a finite length vector of probabilities 0

≤ P i ≤ 1, 1 ≤ i ≤ κ.
In analogy to [START_REF] Grimmett | Inhomogeneous bond percolation on square, triangular and hexagonal lattices[END_REF], we construct an Inhomogeneous Bernoulli Percolation Process upon C 1200 , also called P-percolation process, as a site selection process according to the following algorithm: P-percolation

• Input: G ← C 1200 , v ∈ V, P = (P 1 , P 2 , . . . P κ ).
• Put vertices in V \ {v} to be in state ¬Ex, assign state(v) ← S and set

V P 0 (v) = {v}. • for j from 1 to κ do: ∀k ∈ V P j-1 (v), ∀l ∈ N G (k) : * if state(l) == ¬Ex then:
• Put X ∼ B(1, P j ), the Bernoulli random variable with probability P j ; • If X == 1 then add l to V P j (v) and assign state(l) ← S . * else state(l) ← ¬ S .

• Ouput: ∪ κ j=0 V P j (v). As a result, ranch P (v) ← ∪ κ j=0 V P j (v) ⊆ V is obtained layer by layer by adding successive rings V P j (v) of vertices centered at v in which the selection follows independent Bernoulli trials, hence building a random neighborhood that represents a natural probabilistic structure with which bivalents (CTC's) will be modeled. Clearly, ranch P (v) ≡ V(v) iff all the values of P are 1, otherwise it is a subset of nodes of V(v).

We assume that the CTC will be attached to the nuclear envelope at the points given by the vertices in ranch P (v), v ∈ V. v denotes the position of the SC, and the set ranch P (v) \ {v} the CPCH, see Figure 4.

represents the position of a S C and corresponds to the intersection between the chromatine and the nucleus surface discretization.

From the fact that -by construction -the sets V P j (v), j = 1, 2; . . . are mutually disjoint, the number of elements of the ranch is given by

|ranch P (v)| = κ j=0 |V P j (v)|.
We restrict ourselves to the graph G = (V, E) with |V| = N = 602 nodes, as introduced in section 4.1, see also Figure 3, and the vector P a given κ-dimensional vector of probabilities.

Let P = (P 1 , P 2 , . . . P κ ). We shall derive an expression for the first two moments of the probability distribution of the size of the ranch P (v) of a given vertex v ∈ V: E(|ranch P (v)|) and σ 2 (|ranch P (v)|). For this, we assign to the central position of the local planar graph in which v is contained the (local) coordinate (0,0). In this coordinate system, vertices located at distance 1 will be given the positions (1, 0), (1, 1), . . . [START_REF] Allard | Exact solution of bond percolation on small arbitrary graphs[END_REF][START_REF] Gilbert | Random Graphs[END_REF] counting clockwise, as in figure 5, and so forth for vertices located at greater distances: for i ≥ 1, j ∈ {0, 1, . . . 6i -1}, the vertex (i, j) will be the j-th vertex counting clockwise from the central vertex located at distance i from the center.

Denote by Z(i, j) the Bernoulli random variable Z(i, j) ∼ B(1, p i, j ), 0 ≤ p i, j ≤ 1 that takes the value 1 with probability p i, j if the vertex (i, j) belongs to the ranch P , 0 otherwise. This means that (i, j) is reached from the central vertex by the P-percolation described above, i.e., by a percolation path of length i.

Notice that at each level i = 1, 2, . . . , κ, once the preceding level has been determined by the P percolation algorithm, the set of reachable vertices depends strongly upon the previous level. In terms of the random variables Z(i, j), this means that Z(i, j) depends always upon some Z(i -1, l) that lies below, for some l. The dependency/independency of the Z variables is illustrated in Table 2. Due to the particular construction and symmetry of the local architecture of the graph, being 6-regular, variables Z(i, 0) and Z(i, i • k), for k ∈ {1, 2, . . . 5}, are independent, identically distributed Bernoulli random variables with parameter π i := i j=1 P j . As for the probability p i, j of Z(i, j) ∼ B(1, p i, j ) for j i • k, k ∈ {1, 2, . . . 5}, it can be readily seen that it satisfies the following finite difference equation

p i, j = P i p i-1, j-1 + p i-1, j -p i-1, j-1 • p i-1, j , 1 ≤ i ≤ κ. (4.1)
This equation follows from the fact that, by the graph architecture, nodes at distance i of the central node v can be reached either by a unique path of length i, which is the case for nodes at (i, 0), (i, i), . . . , (i, 5i) or by the two sole adjacent paths reaching (i, j): one that pass through the node situated at (i -1, j), the other one through (i -1, j -1), so that Eqn. (4.1) is nothing else than the probability of the union. Theorem 4.2. Let v = (0, 0) be a given node of the graph G, P = (P 1 , P 2 ) Then the expected size and variance of the size of ranch P 1 ,P 2 (v) is given by

µ(|ranch P (v)|) = 1 + 6P 1 + 6P 1 P 2 (3 -P 1 ) (4.2) σ 2 (|ranch P (v)|) = 6P 1 (1 -P 1 ) + 6P 1 P 2 9 -11P 1 + 4P 2 1 + 6P 1 P 2 2 (6 -19P 1 + 14P 2 1 -3P 3 1 ). (4.3)
Proof. The mean size and variance of |ranch P (v)| will be given by the number of vertices in the neighborhood of v = (0, 0) which are selected by the procedure. With our notation, (see figure 5)

|ranch P (v)| = 1 + 2 i=1 6i-1 j=0 Z(i, j),
from which

E(|ranch P (v)|) = 1 + 2 i=1 6i-1 j=0 p i, j = 1 + 6P 1 + 6P 1 P 2 (3 -P 1 )
because all of them are Bernoulli random variables, which is the easy part and uses equation (4.1). The variance, though, is more involved, as the random variables are not mutually independent, save those of the first ring Z(1, 0), . . . Z [START_REF] Allard | Exact solution of bond percolation on small arbitrary graphs[END_REF][START_REF] Gilbert | Random Graphs[END_REF].

Variables Z(1, i) ∼ B(1, P 1 ), i = 0, 1, . . . 5, are independent identically distributed Bernoulli random variables with probability P 1 . It is immediate to recognize that Z(2, 0), Z(2, 2), . . . Z(2, 10) are all i.i.d ∼ B(1, P 1 P 2 ). From equation (4.1), variables Z(2, 1), Z(2, 3), . . . Z(2, 11) are all B(1, P 1 P 2 (2 -P 1 )). But they are not all independent among them, nor among the variables directly below following any path of length 2 to the central vertex (0, 0):

Depends on Independent of Z(2, 0) Z(2, 11), Z(2, 1) Z(2, 2) --Z(2, 10) Z(2, 1) Z(2, 11), Z(2, 0), Z(2, 2), Z(2, 3) Z(2, 4) --Z(2, 10) . . . . . . . . .

Table 2 Z dependencies at level 2

This is relevant for the computation of the variance, since then the joint probability distribution must be determined. Since

σ 2 (ranch P 1 ,P 2 (v)) = E |ranch P 1 ,P 2 (v)| 2 -E 2 (|ranch P 1 ,P 2 (v)|),
we expand the terms involved and get, for example

        1 + 2 i=1 6i-1 j=0 Z(i, j)         2 = 1 + 2 5 j=0 Z(1, j) + 2 11 j=0 Z(2, j) + 2 5 k=0 11 l=0 Z(1, k) • Z(2, l) + 5 k=0 5 l=0 Z(1, k) • Z(1, l) + 11 k=0 11 l=0 Z(2, k) • Z(2, l).
Also,

E 2 (|ranch P 1 ,P 2 (v)|) = E 2         1 + 5 j=0 Z(1, j) + 11 j=0 Z(2, j)         = (1 + 6P 1 + 6P 1 P 2 (3 -P 1 )) 2 , so that σ 2 (|ranch P 1 ,P 2 |) = 1 + 2 • 6P 1 + 2(6P 1 P 2 + 6P 1 P 2 (2 -P 1 )) + (6P 1 + 30P 2 1 ) + 2 5 k=0 11 l=0 E (Z(1, k) • Z(2, l)) + 11 k=0 11 l=0 E (Z(2, k) • Z(2, l)) -(1 + 6P 1 + 6P 1 P 2 (3 -P 1 )) 2 .
It remains only to compute the terms 5

k=0 11 l=0 E (Z(1, k) • Z(2, l))
for the interactions between level 1 and 2, and 11 k=0 11 l=0 E (Z(2, k) • Z(2, l)) for level 2, which is done by using the symmetry of the local graph structure (6-regularity) and the fact that the product of Bernoulli random variables is again Bernoulli.

Thus, for example, by conditioning upon the random variables at the level below, we determine that Z(1, 0) • Z(2, 0) ∼ B(1, P 1 P 2 ) and it has the same distribution as Z(1, 0) • Z(2, 1) or Z(1, 0)•Z [START_REF] Berend | Improved bounds on Bell numbers and on moments of sums of random variables[END_REF][START_REF] Core | R: A Language and Environment for Statistical Computing R Foundation for Statistical Computing[END_REF] and that Z(1, 0) is independent of Z(2, 2), Z(2, 3), . . . Z [START_REF] Berend | Improved bounds on Bell numbers and on moments of sums of random variables[END_REF][START_REF] Scherthan | A bouquet makes ends meet[END_REF].

Hence the products yield either a B(1, P

E (Z(1, 0) • Z(2, j)) = 3P 1 P 2 + 5P 2 1 P 2 + 4P 2 1 P 2 (2 -P 1 ).
Since there are six terms of this kind, we conclude that

5 k=0 11 l=0 E (Z(1, k) • Z(2, l)) = 6(3P 1 P 2 + 5P 2 1 P 2 + 4P 2 1 P 2 (2 -P 1 )) = 6P 1 P 2 (3 + 13P 1 -4P 2 1
).

X 0 = Z(1, 0) and X 1 = Z(1, 1). Then

P (Y 0 = 1 ∧ Y 1 = 1) = i, j∈{0,1} P (Y 0 = 1 ∧ Y 1 = 1| X 0 = i ∧ X 1 = j) • P (X 0 = i ∧ X 1 = j) = P 2 2 • P (X 0 = 1 ∧ X 1 = 1) + P 2 2 • P (X 0 = 1 ∧ X 1 = 0) + 0 • P (X 0 = 0 ∧ X 1 = 1) + 0 • P (X 0 = 0 ∧ X 1 = 0) = P 2 2 P 2 1 + P 1 (1 -P 1 ) = P 1 P 2 2 .
If we denote by M = (M i, j ) the matrix with entries m i, j = E (Z(2, i -1) • Z(2, j -1)), we obtain

M =                                                        
                                                       
in which we only wrote the upper diagonal part of the matrix, since it is symmetric. Hence the second sum is found to be

E         11 k=0 11 l=0 Z(2, k) • Z(2, l)         = 6π 2 + 6 π2 + 24a + 30b + 48c + 12d + 18e = 6P 1 P 2 (3 -P 1 ) + 6P 1 P 2 2 (4 + 23P 1 -16P 2 1 + 3P 3 1 ). (4.4)
Adding up all terms, we conclude that

σ 2 (|ranch P 1 ,P 2 (v)|) = 6P 1 (1 -P 1 ) + 6P 1 P 2 (9 -11P 1 + 4P 2 1 ) + 6P 1 P 2 2 (6 -19P 1 + 14P 2 1 -3P 3 1 ).
As a small check of consistency, we can point out that if P 1 = P 2 = 1, the percolation process reach every vertex of the neighborhood at distance 2 with probability one, hence the variance is zero: 0 + 6(9 -11 + 4) + 6(6 -19 + 14 -3) = 0; if P 1 = 1, all vertices at distance 1 are reached by the percolation process, rendering the varianceat level 2 -of the sum of 12 independent Bernoulli random variables B(1, P 2 ), which yields 6P 2 (9 -11 + 4) + 6P 2 2 (6 -19 + 14 -3) = 12P 2 (1 -P 2 ), or, if P 2 = 0, P 1 0, 1, only the first layer of vertices is reached, which produces the variance of the sum of 6 independent Bernoulli's B(1, P 1 ), i.e., 6P 1 (1 -P 1 ).

Corollary 4.2. The homogeneous (P 1 , P 1 ) percolation process produces ranches with average size 1 + 6P 1 + 12P 2 1 + 6P 2 1 (1 -P 1 ) and variance

6P 1 + 48P 2 1 -30P 3 1 -90P 4 1 + 84P 5 1 -18P 6 1 .
The maximum variance is obtained with P 1 = 0.59511, at which the expected ranch size is 9.68 and variance 8.43. In contrast, the nonhomogeneous Bernoulli percolation process with P = (P 1 , P 2 ) = (1, 0.47) produce ranches with average size 12.64 and variance 2.98, respectively. Theorem 4.3. Let v = (0, 0) be a given node of the graph G, P = (P 1 ,

P 2 , • • • , P κ ). Let π i := i j=1 P j , i = 1, 2, • • • κ, and
p i, j = π i if j = k • i, k ∈ {0, 1, . . . 5} P i p i-1, j-1 + p i-1, j -p i-1, j-1 • p i-1, j j k • i, k ∈ {0, 1, . . . 5}, (4.5)
Define ∂ κ := 6κ-1 j=0 Z(κ, j), so that

|ranch P 1 ,••• ,P κ (v)| = |ranch P 1 ,••• ,P κ-1 (v)| + ∂ κ . (4.6)
Then the expected size and variance of |ranch

P 1 ,P 2 ,••• ,P κ (v)| is given by µ(|ranch P 1 ,••• ,P κ (v)|) = 1 + κ i=1 6i-1 j=0 p i, j = µ(|ranch P 1 ,••• ,P κ-1 (v)|) + µ(∂ κ ), (4.7)
σ 2 (|ranch P 1 ,••• ,P κ (v)|) = σ 2 (|ranch P 1 ,••• ,P κ-1 (v)|) + σ 2 (∂ κ ) + 2COV |ranch P 1 ,••• ,P κ-1 (v)|, ∂ κ . (4.8)
In light of the computations made so far, the first result (expected value) is by now straightforward. The second result is just the variance of the sum of two random variables, the detailed computation requires computing the joint distribution of all variables that are dependent with each other, so for example Z(i, j) with Z(l, j) whenever 1 ≤ l ≤ i, j = k • i or Z(i, j) and Z(i, l) for 0 ≤ j ≤ l ≤ i -1, etc. We omit this computation.

Farms

Once ranches have been appropriately defined into our graph structure, we turn over to the associations between them, which will be given by overlapping clouds of CPCH and can hence undertake the problem of determining the association distribution: Definition 4.2. Let P = (P 1 , P 2 , . . . , P κ ), G = (V, E) as before.

1. Two ranches r 1 and r 2 are said to overlap if and only if r 1 ∩ r 2 ∅.

2. Two ranches r 1 and r 2 are said to be connected or in association if and only if there exists r i 1 , r i 2 , . . . r i k such that r 1 overlaps r i 1 , r i j overlaps r i j+1 , 1 ≤ j ≤ k -1 and r i k overlaps r 2 . 3. A Farm is a partition induced upon any set of ranches by the equivalence relation of "being connected" defined between ranches.

Remark: From the biological point of view, an association cluster according to this definition consists of a set of vertices whose overlapping domains represent actual sites in which CPCH from different bivalents come into contact, hence providing a scaffold in which genetic material can be paired, compared, copied, etc., (we do not fix ourselves or pursue this topic further). The observed partitions are the sizes of the equivalence classes determined upon the constructed farm. From first principles, any random distribution of ranches will provide a different partition distribution, it remains to check whether we can reproduce the observed distribution with sufficient accuracy.

For example, in Figure 6 we depict the results obtained for the simulation of 19 SC's distributed randomly upon the set of 602 vertices, along with their corresponding ranches and for the probability vectors a) P = (2/3, 1/2, 2/3) and b) P = (1, 0.5). Boundaries of chromatin have been drawn in order to emphasize the similarity with the observations.

SC Dynamics

Notice that in the field of dynamics upon a fixed regular net, cellular automata and agent based models have been widely used, see for example [START_REF] Leng | An extended floor field model based on regular hexagonal cells for pedestrian simulation[END_REF], and also to a lesser extent some dynamics have been defined upon planar tilings made of hexagonal and/or pentagonal cells, as in [START_REF] Porzycki | Agent based approach and cellular automata: A promising perspective in crowd dynamics modeling?[END_REF]. From this point of view, we can consider that the SC's can be thought of as being agents interacting upon the surface of the nucleus (a finite connected field constituted by 590 hexagonal cells and 12 pentagonal cells).

It is known that in an earlier stage of the prophase, previous to pachytene, the size of the nuclear envelope is a fraction of the one used in our setting, and that at this

SC-Dynamics

• Input G ← C 1200 , n.IT ← number of iterations. • Initialize: α 0 i ← 1, for i ∈ {1, 2, . . . , 19}, otherwise α 0 i ← 0, i ∈ {20, . . . , 602}.

• Iterate for t ∈ {0, . . . , n.IT -1}:

• Step 1 : for all 1 ≤ i ≤ 602 do:

-If α t i == 1 then choose v ∈ N G (i) such that α t v == 0 if any; put (i, v) in a non ordered list L.
• Step 2 : α t+1 ← α t ; while L ∅ do:

-extract randomly from L an element (i, v) -If (α t+1 v == 0) then α t+1 i ← 0 and α t+1 v ← 1 • Output: MC ← {i ∈ V : α n.IT i == 1}
At each time unit the complete set of SC's is allowed to simultaneously choose a free neighbor and move to it, or to remain at the present position if either there are no free neighbors or the chosen one was already selected marked for movement. In this dynamics, the sum 602 i=1 α t i = 19 is a conserved quantity. It is a conflict free parallel discrete time upgrading dynamics.

Figure 8 illustrates one possible step from time t to time t + 1 for a parallel evolution of three nodes in a subgraph of G. Simulations produce convergence to a uniform distribution for the positions of the SC's upon the nuclear envelope, already after ∼ 10 4 iterations. Hence, we can implement now an algorithm for assessing the percolation process after a uniform sample of nodes has been reached.

Simulation results

We simulate the random bivalent's associations via the following Farm-Partition algorithm:

Farm-Partition

• Input: G ← C 1200 , P = (P 1 , P 2 , . . . P κ ).

-MC <-SC-Dynamics (G, n.IT ← 10 4 ).

-For v ∈ MC: ranch P (v) ←P-Percolation(G, v, P).

-Farm ← {F 1 , . . . , F r } vertex sets of the connected components of ∪ v∈MC ranch P (v) in G.

• Output: Part ← sort( |F j |, j = 1, . . . , r, decreasing = TRUE).

This algorithm was implemented using R software [START_REF] Core | R: A Language and Environment for Statistical Computing R Foundation for Statistical Computing[END_REF] and iterated ad libitum in order to produce a simulated set of nuclei with which available data could be contrasted. Notice that in our notation, the biggest subcluster found in the partition corresponds to Part [START_REF] Allard | Exact solution of bond percolation on small arbitrary graphs[END_REF], i.e., 1 st -Class, the second one in Part [START_REF] Berend | Improved bounds on Bell numbers and on moments of sums of random variables[END_REF], i.e., 2 nd -Class, etc. Best approximation to the observed distribution for the 1 st -and 2 nd -Class was obtained for different values of P = (P 1 , P 2 ), where a bootstrapping technique was used for assessing the intrinsic variation in the data/simulations in both cases, the homogeneous and the non homogeneous one.

Figure 9 shows, for example, the results obtained for the homogeneous percolation with P = (0.68, 0.68) after simulating 10 4 spermatocytes (i.e., by running Farm-Partition 10 4 times with G ← C 1200 , P = (0.68, 0.68)). Best approximation for 1 st -Class and 2 nd -Class distributions for the homogeneous case P = (0.68, 0.68). Similar approximations are obtained in the non homogeneous case for various values of the P-vector, see Table 3 3 A sample of (P 1 , P 2 ) values for which the associated P percolation process gives the best approximation for the 1 st -and 2 nd -Class probability distribution, together with theoretical and simulated ranch sizes and their variances.

Other values with good agreement with respect to the data are the non homogeneous case P = (1, 0.47), or P = (0.8, 0.6), etc., see Table 3. Included are the theoretical values for the mean value and variance of the ranch's sizes in the 6-regular approximation (equations (4.2) and (4.3) from Theorem 4.2), as well as those obtained by simulation in the dual fullerene C 1200 . The homogeneous percolation with minimal variance P = (0.595, 0.595) does not furnish a satisfactory 1 st -Class probability distribution, so it is not included in the table. Figure 10 depicts the zone in (P 1 , P 2 )-space with the best fitting to the observed frequencies in both, 1 st -Class and 2 nd -Class. The criteria used for determining the lower and upper curves was that the entire simulated mean distribution does not fall outside the variation zone around the observed mean values (frequencies) in the data set. The red curve depicts the best fit in which both, the simulated and the observed mean distribution lie close together. Depicted here is also the homogeneous percolation region P = (P 1 , P 1 ), which -intersected with the region in between the lower and upper curves -yield the best approximation range for the distribution of 1 st -Class and 2 nd -Class in the homogeneous case.

Comments and Discussion

Chromosomal bivalent's associations through intersecting CPCH create rich dynamic and diverse scenarios via the participating elements. They are triggered by the corresponding associations of CPCH located at the short arms of the bivalents, but also by the resulting convergence of the rest of the constituent chromatin along them. This supra chromosomal bonds allows, for example, joint expression of genes coming from different bivalents.

Undoubtedly, the better understanding of the general principles behind bivalent's associations in prophase meiotic nuclei, as well as precising the type of randomness being at play at this stage could bring us also a step closer to a better understanding of the different chromosome combinations present in the gametes. Since these associations and combinations persist until the meiotic divisions, the chromosomal associations as described here necessarily leaves some imprint in the chromosomal sets passed on to gametes and hence their importance to evolution. While available data provides no distinctive insights upon prophase progression from the bouquet state to the observed associations, whose distributions are studied here, it is by simulating different scenarios that we obtain a simple conflict free parallel dynamics for randomly selecting free neighbors, with which these distributions, as measured by 1 st -and 2 nd -Classes, can be correctly reproduced. For this, a suitable choice of vector P has to be made, which is a parameter needed for the ranch determination and therefore is an architectural parameter worth considering.

The approach pursued in this work possess clear advantages over earlier modeling attempts, since now we have a tool that not only reproduces available data, but can also drive further biological research in topics that have not been explored yet, for example expressing or assessing individual properties attached to the individual chromosomes via different choices of their ranches, i.e., via assigning individual vector P's to each CTC. To our knowledge this has not been reported in the literature and will be undertaken next. We conjecture a correlation between our ranch-sizes, expressed in terms of the number of neighboring vertices with chromatin attributes, and the actual sizes of the chromosomes (respectively their chromatin content). These conditions should allow us to propose new experiments for assessing individual chromatin content (for example as a ratio volume/surface) or to characterize it by means of appropriate choices of vectors P. This approach would allow us to effectively relax our previous hypothesis of indistinguishable CTC's, or at least, to put it to a test.

From the biological point of view, the introduction of chromatin neighborhoods with (random) variable sizes reflects more properly the true nature of the objects being modeled and, perhaps more significantly, it is the very intersection of chromatin domains, which appears here naturally for describing the interactions of CTC's via its induced clustering relation, that represents better the biological importance of the described process, because, while in hetero-chromatin there is no actual genic expression, it certainly contributes to bringing together genes located at other chromosomal domains, so that the associations themselves could be conceived as a form of dynamical organization that contributes to the functionality of the joint genic expression, which -ultimatelymay be looking at us hidden behind the available data.

Fig 1 .

 1 Fig 1. a) A representative pachytene spermatocyte spread of 2n = 40 Mus domesticus mice treated by immunocytochemical techniques. In red the pericentromeric heterochromatin domains, in green the synaptonemal complexes. The (sex) bivalent XY is indicated. b) Observed 1 st -Class and 2 nd -Class frequencies with bootstrap variation

Figure 2 ,

 2 b) shows the observed distribution for 1st-Class 4.

Fig 2 .

 2 Fig 2. a) The full set of 400 observations, with their frequencies. b) The subset of 1st-Class 4

Theorem 4 . 1 .

 41 If G is a planar graph of N ≥ 3 vertices and e edges, then e ≤ 3N -6.

Corollary 4 . 1 .

 41 Every planar graph of N ≥ 3 vertices contains a vertex of degree at most 5.

Fig 3 .

 3 Fig 3. A local planar view of the (almost) 6-regular sphere and the resulting graph C 1200 .

Definition 4 . 1 .

 41 Let P = (P 1 , P 2 , . . . , P κ ) ∈ [0, 1] κ . Let v ∈ V be a given (fixed) vertex of G. The random chromatin neighborhood of v, denoted by ranch P (v) ⊆ V, is the set of vertices in state S resulting from a P-percolation processes around v.

Fig 4 .

 4 Fig 4. (a) Chromatin neigborhood created with P = (1, 0.5). (b) Chromatin neigborhood created with P = (0.8, 0.6).

π 2 aπ2 a d c e π 2 a b c π2 a d π 2 a π2

 2 b c b c b c b c b a π2 a d c e c e c e c d π 2 a b c b c b c b c π2 a d c e c e c e π 2 a b c b c b c π2 a d c e c e π 2 a b c b c

Fig 6 .

 6 Fig 6. Artist's rendering of simulated spermatocytes (a) P = (1, 1/2), partition 19 = 6+4+3+2+1+1+1+1. (b) P = (2/3, 1/3, 7/24), partition 19 = 7 + 3 + 2 + 2 + 1 + 1 + 1 + 1 + 1. (Lengths are not to scale).

Fig 8 .

 8 Fig 8. Description of 1 step of (parallel) evolution in a subgraph of G. Red lines indicate direction of displacement.

Fig 9 .

 9 Fig 9.Best approximation for 1 st -Class and 2 nd -Class distributions for the homogeneous case P = (0.68, 0.68). Similar approximations are obtained in the non homogeneous case for various values of the P-vector, see Table3anf Figure10.

Fig 10 .

 10 Fig 10. Best approximation regions in P = (P 1 , P 2 )space for 1 st -Class and 2 nd -Class distributions obtained with P-percolation.

  anf Figure10.

	P 1	P 2	µ	σ 2	μ	σ2
	1.0	0.47	12.64	2.99	12.55	3.00
	0.8	0.60	12.14	5.84	12.05	5.97
	0.68	0.68	11.52	8.03	11.40	8.03
	0.6	0.78	11.34	10.01	11.22	10.06
	0.52	1.0	11.86	13.50	11.82	13.86
				Table

* Partially supported by Chilean MINEDUC Grant MECE-SUP 2016-2017. † Partially supported by the LabEx PERSYVAL-Lab (ANR-11-LABX-0025-01) funded by the French program Investissement d'avenir.

‡ "This is a pre-print of an article published in Journal Math. Biology.

In the similar lines of reasoning, the products Z(2, i) • Z(2, j) are determined to be Bernoulli, with probabilities π 2 := P 1 P 2 , π2 := π 2 (2 -P 1 ), a := P 1 P 2 2 = π 2 • P 2 , b := P 2 1 P 2 2 = π 2 2 , c := P 2 1 P 2 2 (2 -P 1 ) = π 2 • π2 , d := P 1 P 2 2 (1 + P 1 -P 2 1 ) = a • (1 + P 1 -P 2 1 ), and e := P 2 1 P 2 2 (2 -P 1 ) 2 = π2 2 . We explain a few of these expressions, the rest follows from symmetry considerations:

, which is deduced by conditioning upon Z(1, 0) and Z(1, 1), see below. In order to see why -for example -Z(2, 0)

2 ), it suffices to recall that once a vertex, say w, at a given level i has been selected by the percolation process, all vertices at level i + 1 that can be reached from it are available for the next Bernoulli selection B(1, P i+1 ).

Hence, the joint probability P (Z(2, 0) = 1 ∧ Z(2, 1) = 1) is determined upon conditioning with respect to variables Z(1, 0) and Z(1, 1): Let Y 0 = Z(2, 0), Y 1 = Z(2, 1), stage bivalents appear to be grouped into a single "bouquet" (Figure 7), see [START_REF] Scherthan | A bouquet makes ends meet[END_REF], it is during the pachytene stage that the bivalents (and more specifically the SC's) seem to be randomly (i.e;, uniformly) distributed.

While the available data does not provide information about the wandering processes of the SC's gliding upon the internal surface of the nucleus, we nevertheless propose and explore a simple dynamics for the displacement of the SC's upon the internal surface of the fullerene dual. It will turn out that it actually converges to the uniform distribution: Let G = (V, E) be the graph representing the nuclear envelope of a spermatocyte. A node v ∈ V corresponds to a position upon the nuclear envelope that can accommodate either the SC or the CPCH of a given bivalent.

To decide whether at time t this position is occupied by an SC or not, we add an attribute to each vertex: Definition 4.3. The attribute at time t of a vertex v ∈ V, denoted by α t v ∈ {0, 1} is defined by: α t v = 1 if and only if a SC is located at (the position of) vertex v ∈ V at time t.

At time t = 0, we assume that all the SC are grouped and occupy the 19 first nodes, see Figure 7. (The assumption is not restrictive, as any other compact set of vertices forming a disc can also be chosen.) Thus the vector α 0 = (α 0 1 , . . . , α 0 602 ), which describes the "bouquet" configuration, is defined as α 0 i = 1, i ∈ {1, 2, . . . , 19}, and α 0 i = 1 for all i ∈ {20, 21, . . . , 602}. The only constraint is that two SC cannot share the same position at the same time. Thus the random trajectory of an SC is only affected by the position of the other SC's.

At time t > 0, with vector α t already been computed, the positions of the SC's will be given by the following algorithm:


Salma Jiddi

Philippe Robert

Eric Marchand

Technicolor

ESTIMATION OF POSITION AND INTENSITY OF DYNAMIC LIGHT SOURCES USING CAST SHADOWS ON TEXTURED REAL SURFACES

Keywords: Illumination, Reflectance, Textures, Shadows, Augmented reality

In this paper, we consider the problem of estimating the 3D position and intensity of multiple light sources without using any light probe or user interaction. The proposed approach is twofold and relies on RGB-D data acquired with a low cost 3D sensor. First, we separate albedo/texture and illumination using lightness ratios between pairs of points with the same reflectance property but subject to different lighting conditions. Our selection algorithm is robust in presence of challenging textured surfaces. Then, estimated illumination ratios are integrated, at each frame, within an iterative process to recover position and intensity of light sources responsible of cast shadows. Estimated lighting characteristics are finally used to achieve realistic Augmented Reality (AR).

INTRODUCTION

The image brightness of a 3D scene is a function of three components: scene geometry, surface reflectance and illumination distribution. Provided that these three components are accurately estimated, one is able to photo-realistically render a virtual scene. In this paper, we focus on estimating illumination characteristics. Specifically, we aim at recovering the 3D position and intensity of light sources in indoor scenes. Existing solutions often use additional devices such as light probes [START_REF] Debevec | Rendering synthetic objects into real scenes: Bridging traditional and image-based graphics with global illumination and high dynamic range photography[END_REF] [START_REF] Nowrouzezahrai | Light factorization for mixedfrequency shadows in augmented reality[END_REF] or fisheye cameras [START_REF] Rohmer | Interactive near-field illumination for photorealistic augmented reality on mobile devices[END_REF] [START_REF] Knecht | Interactive brdf estimation for mixed-reality applications[END_REF]. Also, when only RGB-D data is considered, lighting is often assumed to be distant [START_REF] Sato | Illumination distribution from brightness in shadows: Adaptive estimation of illumination distribution with unknown reflectance properties in shadow regions[END_REF] or reduced to a single point light [START_REF] Arief | Realtime estimation of illumination direction for augmented reality on mobile devices[END_REF] [START_REF] Boom | Point light source estimation based on scenes recorded by a rgb-d camera[END_REF].

The proposed framework handles both static and dynamic lighting. It takes three inputs: (a) coarse 3D model of the scene acquired by an active sensor (Intel R200). (b) nearambient reference image to which we will refer as reference image. The latter is acquired by simulating an ambient lighting which, in theory, does not generate shading or shadowing. In practice, it can be easily produced by considering a fairly uniform indirect lighting. (c) color images of the scene from which lighting will be recovered.

The assumptions we made for the proposed approach are: (i) scene geometry is assumed to be static and contains a main planar surface on which shadows are cast. (ii) scene reflectance is described by the Lambertian reflection model. In this work, the viewpoint is fixed within the analysis part to avoid camera pose registration errors but can vary in AR.

The main contributions of our method are as follows: Illumination ratio maps. we propose a novel approach to factor albedo/texture and illumination in arbitrary textured surfaces. Our algorithm relies on a fine selection of pairs of points with the same albedo but subject to different illumination conditions. Lightness ratios of selected pairs provide an illumination ratio map which is texture-free. This map contains mainly shadowing information and is recovered for each incoming frame as the lighting can be static and/or dynamic. Position and intensity of multiple light sources. Illumination distribution is initially represented by a set of equally distributed point lights. Within an iterative process, we extract the 3D position of lights whose shadow maps correlate with the illumination ratio map. Finally, a robust estimator allows to recover their respective intensities. Near real-time processing. On average, the implementation of our method is able to perform at approximately 4fps. Recovered illumination characteristics are further integrated in a real-time rendering pipeline for AR with virtual shadows that are visually coherent in terms of shape and attenuation.

RELATED WORK

Early work on illumination estimation using shadowing information was presented by Sato et al. [START_REF] Sato | Illumination distribution from brightness in shadows: Adaptive estimation of illumination distribution with unknown reflectance properties in shadow regions[END_REF]. They used cast shadows to recover the direction of multiple distant lighting. In [START_REF] Sato | Illumination distribution from brightness in shadows: Adaptive estimation of illumination distribution with unknown reflectance properties in shadow regions[END_REF], cast shadows are detected by manually removing occluding objects and capturing a shadow-free image. Panagopoulos et al. [START_REF] Panagopoulos | Illumination estimation and cast shadow detection through a higherorder graphical model[END_REF][9] used a graphical model and a voting procedure to detect shadows cast on planar textured surfaces under the assumption of distant lighting. Their framework takes about 3 minutes per image, which is critical for AR scenarios. In [START_REF] Arief | Realtime estimation of illumination direction for augmented reality on mobile devices[END_REF], Arief et al. analyzed the shadow cast by a cube to estimate the position of a single light source. The surface on which the shadow is cast must be uniform. Our methodology differs from previous shadow-based estimation work in three main points: first, our approach does not require any user intervention or reference objects to recover illumination. On the contrary, it targets middle scale scenes including various 3D objects with arbitrary geometry. Secondly, we do not make the assumption of distant lighting and instead, estimate both 3D position and intensity of multiple light sources. Finally, our approach handles uniform and textured surfaces.

OUR APPROACH TO ILLUMINATION ESTIMATION

In the following, we choose Phong reflection model [START_REF] Phong | Illumination for computer generated pictures[END_REF] to describe the way a diffuse point p in the scene reflects light:

I p = k p d (L a + M i=1 (N p • ω p i )L i O p i ) (1)
where I p is the color intensity of a 3D point p, k p d is its albedo and N p is its normal vector. L a and L i are respectively the intensities of ambient lighting and light source i, ω p i is the incoming light direction vector of light source i, and M is the number of light sources present in the scene. O p i is a binary visibility term, equal to 1 if point light i is visible from the 3D point corresponding to pixel p and equal to 0 if occluded.

Estimation of illumination ratio maps

Shadows are caused by the occlusion of incoming light, and thus contain various pieces of information about the illumination of the scene. For this reason, shadow detection has long been considered within the image processing research community [START_REF] Finlayson | Entropy minimization for shadow removal[END_REF][12][13][14][8] [START_REF] Panagopoulos | Robust shadow and illumination estimation using a mixture model[END_REF]. If we consider only local texture features, we are not capable of accurately separating texture/albedo and illumination. Up to now, it remains an extremely challenging problem and no generic solution exists. In this section, our goal is to implicitly detect cast shadows for a given frame. Our approach is twofold:

The first step is achieved only once since geometry is static: we detect the main planar surface on which shadows are cast. Specifically, we compute surface normals N and apply a region growing algorithm to cluster similarly oriented surfaces together (a deviation value of 3 degrees is typically allowed between normals). Then, we use a RANSAC estimator to fit each cluster to a planar surface model. Finally, the cluster including the largest number of inliers is considered as the main plane. Points above the detected plane are further grouped as belonging to the occluding objects.

Secondly, we aim at separating texture/albedo and illumination in the current frame. The analysis is limited within the 2D projection of the previously detected plane. Our factorization is achieved through a voting scheme using pairs of points (p, p) with the same reflectance k d but subject to different lighting conditions. We propose a fine selection based on two features: (i) L 2 norm in the CIELAB color space within the reference image R as it provides accurate similarity measures compared to using only the current frame (Fig. 1). In fact, the use of the reference image R makes our algorithm robust in presence of challenging textures, poor lighitng conditions and/or sensor noise. (ii) the lightness channel L of the CIELAB color space within the current image enables us to compare pairs illumination-wise. Hence, points with lower lightness values are prone to belong to shadowed regions. Using these two features, the voting scheme is applied as follows:

v(p) = +1, if R p -R p ≤ R and L p ≥ L p + L 0, otherwise (2)
where R and L are respectively thresholds with respect to color vectors distance in R and lightness difference in L. A value of 2.5 is chosen for R as it corresponds to an almost noticeable difference [START_REF] Sharma | Digital color imaging handbook[END_REF]. L is kept at a low value (typically 10) to handle both weak and strong shadows. Finally, pixels holding a significant voting value v(p) are further considered to estimate their respective illumination ratio δ(p):

δ(p) = L p L with: L = p L p v(p) ( 3
)
where L is the mean lightness value of matched points p with respect to p. Furthermore, pixels for which similar-reflectance pairs are found but received a low voting value hold an illumination ratio value equal to 1. In fact, these points are considered to be fully lit. Last but not least, pixels for which no match is found are discarded as we can not reliably estimate their illumination ratio (green pixels in Fig. 1).

Estimation of light sources 3D position

Since surface materials are assumed to be opaque, illumination distribution is approximated by a set of point lights equally distributed in the 3D space above the detected plane.

In this section, we aim at recovering the 3D position of light sources responsible of cast shadows in the scene. To achieve this, we extract a subset (S) of point lights whose shadow maps correlate with the estimated illumination ratio map (Section 3.1) for each incoming frame. The identification of subset (S) is carried within an iterative process as follows:

1. We initially compute correlation values by matching the illumination ratio map with the shadow maps of the sampled light candidates. The light source whose shadow map has the best correlation value is selected.

2. For each iteration, previously selected light sources are discarded. Also, the matching operation is carried out by discarding previously matched pixels.

3. The process ends either when the currently selected shadow map has a significantly low matching value or if the number of selected lights is higher than N .

The chosen correlation corresponds to Pearson's correlation coefficient ranging between 0 and 1. Perfectly matching maps have a coefficient equal to 1.

Estimation of light sources intensity

As we consider small to middle scale scenes, we assume that ((N p • ω p ) = cos θ) is equal across selected pairs, and use L i to represent L i cos θ. Subsequently, equation ( 1) becomes:

I p = k p d (L a + M i=1 L i O p i ) (4)
Our goal is to estimate the intensity of recovered light sources (Section 3.2). By considering the ratio of color intensities of a selected pair (p, p) as described in Section 3.1, we have:

I p I p = k d (L a + M i=1 (L i O p i )) k d (L a + M i=1 L i ) = δ(p) (5)
In fact, the albedo terms in the numerator and denominator cancel out. Subsequently, we obtain the illumination ratio value. As mentioned in Section 3.1, fully lit points hold a δ(p) equal to 1. Thus, we set L a + M i=1 L i = 1. Finally, we integrate illumination ratio values within a linear system:

AL = δ (6)
where:

A =        1 O p1 1 • • • O p1 M 1 O p2 1 • • • O p2 M . . . . . . . . . . . . 1 O p N 1 • • • O p N M        ;L =      L a L 1 . . . L M      ;δ =      δ(p 1 ) δ(p 2 ) . . . δ(p N )      (7)
The linear system (6) is solved using an iterative Least Squares with bounds and equality constraints:

L = min L ( 1 2 ωAL -ωδ 2
) such as:

0 ≤ L i ≤ 1 and 0 ≤ L a ≤ 1 L a + M i=1 L i = 1 (8)
The weights ω are computed using Tukey s bisquare loss function. Small weights are discarded throughout iterations.

EXPERIMENTAL RESULTS

The proposed approach has been tested on various real scenes with varying lighting conditions. In Fig. 2, we illustrate our results within a selection of four scenes (S1, S2, S3 and S4). Further results are shown in the video available at https : //youtu.be/kj dY N epRDg. In Fig. 2, scenes are grouped row-wise. In the first column, we overlay the contours of shadow maps corresponding to recovered lights on the input color frame. For instance, red and green contours are used respectively for the first and second detected lights. In the second column, we demonstrate estimated illumination ratio maps where background/noise are represented by red pixels and occluding objects by blue pixels. On the detected planar surface, white pixels represent fully lit 3D points (their illumination ratio is equal to 1), green pixels represent 3D points for which no pairs are found. Finally, grayscale pixels are the illumination ratio value. These pixels correspond to 3D points partially or fully occluded with regard to lighting.

Illumination ratio maps

As illustrated in Fig. 2, our algorithm delivers accurate illumination ratio maps where texture/albedo is accurately separated from illumination. The proposed framework handles both uniform surfaces (S1) and challenging textured surfaces (S2, S3 and S4). Also, we demonstrate robustness in presence of poor geometry, especially when the scene contains specular objects (S1 and S3). In fact, the Lambertian assumption must mainly hold on the principal plane. Finally, the use of the reference image enables us to be robust even when lighting and/or image quality lack (S4).

Recovered 3D position of light sources

In Fig. 2, we overlay the contours of the selected shadow maps on the current color frame. For S2, the illumination ratio map contains shadow information but also illumination ratio values corresponding to the narrow spot light effect. Nonetheless, our algorithm succeeds in recovering an accurate light position. Furthermore, our approach demonstrates good results in the presence of overlapping shadows (S3 and S4).

In order to evaluate the precision of recovered lighting positions, we used a special setup for several experimental scenes. First, we choose a world frame on the main plane and measure distance to real lights using a telemeter. Our algorithm is tested on various scenes and recovers light sources position with an average error of 17cm for a mean distance of 2.55m to the light source and a standard deviation of 3.5cm.

As Augmented Reality is our target application, temporal stability with regard to recovered light sources is of paramount importance. When the lighting is static, recovered lighting properties (position and intensity) must be temporally stable, otherwise virtual shadows might suffer from apparent flickering. In Fig. 3, we can notice that under the same lighting conditions, the selected light source (ID98)

Estimated light sources intensity

In this section, we demonstrate the effectiveness of our approach when estimating the intensity of light sources. First, we show the temporal stability of the intensity estimates throughout a sampled sequence for scenes S2 and S4 under varying lighting (Fig. 4). Furthermore, recovered lighting in- tensities are used to render virtual shadows. In Fig. 5, the first row represents reconstructed shading using scene geometry, lighting 3D position and intensity for scenes S3 and S4. A virtual sphere is introduced as well to demonstrate its interaction with the real scene. Furthermore, the second row in Fig. 5 shows augmented scenes where virtual shadows are visually coherent with real shadows in terms of shape and intensity.

CONCLUSION AND FUTURE WORK

We presented a new approach to estimate both position and intensity of multiple light sources. The algorithm succeeds in accurately factoring texture/albedo and illumination on challenging textured surfaces. For future work, we are interested in relaxing the Lambertian assumption. Furthermore, achieving a full GPU implementation would allow us to tackle online analysis and rendering.

Fig. 1 .

 1 Fig. 1. (a): reference image R of the scene. (b): current scene capture. (c)(d): recovered illumination ratio map using respectively (b) and (a). (d) demonstrate a better separation of texture and lighting. Note the shadow of the front book in (d) compared to (c) as well as discarded pixels (green).

Fig. 2 .

 2 Fig. 2. Col.1: Overlay of selected shadow maps contours on the input color frame. Col.2: Estimated illumination ratio maps for uniform (S1) and textured surfaces (S2, S3 and S4).

Fig. 3 .

 3 Fig. 3. First and second best correlation coefficients for a scene under static lighting.

Fig. 4 .

 4 Fig. 4. Recovered ambient intensity L a and point light sources intensity L 1 and L 2 for scenes S2 and S4.

Fig. 5 .

 5 Fig. 5. Row 1: Reconstructed shading using scene geometry and recovered lighting properties (position and intensity) for S3 (left) and S4 (right). Row 2: AR scenarios with visually coherent virtual shadows (red capsule in S1 (left) and brown cube in S2 (right)).


Joan Van Baaren
email: joan.van-baaren@univ-rennes1.fr

Ulrika Candolin
email: ulrika.candolin@helsinki.fi

Plasticity in a changing world: Behavioural responses to human perturbations

Most insect species are affected by Human Induced Rapid Environmental Changes (HIREC). Multiple responses to HIREC are observed in insects, such as modifications of their morphology, physiology, behavioural strategies or phenology. Most of the responses involve phenotypic plasticity rather than genetic evolution. Here, we review the involvement of behavioural plasticity in foraging, reproduction, habitat choice and dispersal; and how behavioural plasticity modifies social behavior and inter-specific interactions. Although important, behavioural plasticity is rarely sufficient to cope with HIREC. An increasing number of studies find species to respond maladaptively or insufficiently to various anthropogenic disturbances, and less often is large degree of plasticity linked to success.

Introduction

Most insect species are affected by Human Induced Rapid Environmental Changes (HIREC, defined by [START_REF] Sih | Evolution and behavioural responses to human-induced rapid environmental change[END_REF] [START_REF] Sih | Evolution and behavioural responses to human-induced rapid environmental change[END_REF], which include variable threats like climate change, habitat fragmentation, habitat loss, human harvesting, and pollution. The species can respond to HIREC, as to any other environmental changes, through plasticity, genetic evolution or dispersal. HIREC generally impact negatively the species. However some species cope well with HIREC and their populations increase to a point where they become themselves a threat to other species. Multiple responses to HIREC are observed in insects, such as modifications of their morphology (eg, size, wing area), physiology (eg, immune response, metabolic rate), behavioural strategies or phenology [START_REF] Sulmon | Abiotic stressors and stress responses: What commonalities appear between species across biological organization levels?[END_REF]. [START_REF] Hendry | Human influences on rates of phenotypic change in wild animal populations[END_REF] [START_REF] Hendry | Human influences on rates of phenotypic change in wild animal populations[END_REF] used meta-analysis dedicated to animals and found most of the responses to involve phenotypic plasticity rather than genetic evolution.

Moreover, when genetic changes or shifts in demography, distribution or phenology occurred, these were generally preceded by a modification involving phenotypic plasticity [START_REF] Beever | Behavioral flexibility as a mechanism for coping with climate change[END_REF]. Thus, behavioural plasticity appears important in explaining variation in the success of species to resist HIREC [START_REF] Sih | Understanding variation in behavioural responses to human-induced rapid environmental change: A conceptual overview[END_REF] (Figure 1). The behavioural responses can be either maladaptive, such as the incapacity to detect new predators that can precede species' decline [START_REF] Sih | Understanding variation in behavioural responses to human-induced rapid environmental change: A conceptual overview[END_REF], or be adaptive and improve fitness, such as the finding of new host plant species by Drosophila suzukii when invading new areas [START_REF] Poyet | The wide potential trophic niche of the Asiatic fruit fly Drosophila suzukii: The key of its invasion success in temperate Europe?[END_REF]. In this chapter, we examine how behavioural plasticity is involved in insect responses to HIREC during foraging, reproduction, habitat choice and dispersal; and how behavioural plasticity modifies social behavior and inter-specific interactions. We also examine if this plasticity is sufficient to respond adaptively to HIREC or not.

Foraging behavior

The Optimal Foraging Theory (OFT) predicts that generalist species should turn to specialist strategies when resources become rare, such as after habitat changes due to HIREC. Species that show plasticity in their temporal or spatial resource use, such as aphid parasitoids [START_REF] Eoche-Bosy | Seasonal variations of host resources influence foraging strategy in parasitoids[END_REF], may have a higher ability to resist HIREC. The main plastic response of foraging insects to changes in their habitats is to modify or enlarge their diet choice. Evans & Moustakas (2017) [START_REF] Evans | Plasticity in foraging behavior as a possible response to climate change[END_REF] showed with a model that if predators (or species at high trophic level), shift among prey species when the preferred prey becomes unavailable because of climate change, they survive longer. This modification of diet has been observed in phytophagous insects, and sometimes improves resistance to climate change, as shown by [START_REF] Raffa | Host selection by bark beetles: playing the odds in a high stakes game. Pine bark beetles[END_REF] [START_REF] Raffa | Host selection by bark beetles: playing the odds in a high stakes game. Pine bark beetles[END_REF] for bark beetles, and sometimes not, as shown by Pol et al. (2017) [START_REF] Pol | Behavioural flexibility does not prevent numerical declines of harvester ants under intense livestock grazing[END_REF] for ants. In some species, like the invasive Drosophila suzukii that shows extreme plasticity in diet choice, including more than 30 plant species, diet breath is probably responsible for their success [START_REF] Poyet | The wide potential trophic niche of the Asiatic fruit fly Drosophila suzukii: The key of its invasion success in temperate Europe?[END_REF]. The capacity of this species to use transgenerational medication (preference A c c e p t e d M a n u s c r i p t of oviposition in fruit containing an entomotoxic substance) contribute to its success [START_REF] Poyet | The invasive pest Drosophila suzukii uses trans-generational medication to resist parasitoid attack[END_REF]. However, D. suzukii also oviposits on plants that prevent the development of their larvae. Such traps could represent a cost of extreme plasticity [START_REF] Poyet | The wide potential trophic niche of the Asiatic fruit fly Drosophila suzukii: The key of its invasion success in temperate Europe?[END_REF].

Another solution to climatic stress is to selectively feed on particular types of food. The tropical butterfly Bicyclus anynana, decreases its activity on hot days and increases its intake of polyphenols (antioxidant) in order to maintain elevated antioxidant levels, which may confer fitness benefits by up-regulating endogenous antioxidant defenses. In grasshoppers, individuals shift their nutrient intake between carbohydrates and proteins depending on temperature and predation stress: indeed, it was shown that a chronic risk of predation (i.e. a chronic stress), induces an elevated metabolism and then a change in nutrient demand and resource consumption and that a higher temperature stands to exacerbate this stress and this diet modification. Growth or reproduction are protein-demanding whereas to cope with stress, a diet balanced towards carbohydrates is more appropriate [START_REF] Schmitz | Temperature dependence of predation stress and the nutritional ecology of a generalist herbivore[END_REF].

Reproductive behaviour

Changes in the environment can influence the ability of individuals to attract and locate mates, as well as their mate choice behaviour. Such changes can alter the number and quality of offspring produced, as well as lead to a demographic Allee effect, and, hence, influence population dynamics [START_REF] Shaw | Mate finding, Allee effects and selection for sex-biased dispersal[END_REF]. Traffic noise, for instance, interferes with the ability of female field crickets Gryllus bimaculatus to locate singing males during mate searching, which can influence their reproductive success [START_REF] Schmidt | Sexy voices -no choices: Male song in noise fails to attract females[END_REF]. To counteract negative effects of noise on mate location, individuals may alter their mate attraction or mate location behaviour. For example, male grasshoppers Chortippus biguttulus elevate the frequency of their courtship song so that it can be distinguished from traffic noise [START_REF] Lampe | Staying tuned: Grasshoppers from noisy roadside habitats produce courtship signals with elevated frequency components[END_REF].

Another common environmental change influencing reproductive behaviours is the increased use of artificial light at night. Male glow-worms Lampyris noctiluca, for instance, are less able to locate glowing females under artificial light [START_REF] Bird | Low levels of light pollution may block the ability of male glow-worms (Lampyris noctiluca L.) to locate females[END_REF], while females of the Australian black field cricket Teleogryllus commodus become less selective in their mate choice [START_REF] Botha | Effects of lifetime exposure to artificial light at night on cricket (Teleogryllus commodus) courtship and mating behaviour[END_REF]. Increased light levels influence also sex pheromone production. For example, females of a noctuid moth Mamestra brassica reduce their sex pheromone production when light levels are high [START_REF] Van Geffen | Artificial night lighting disrupts sex pheromone in a noctuid moth[END_REF].

Changes in habitat structure are common causes of altered reproductive behaviour. An example is the degradation of tropical cloud forests, which has caused an Aftrotropical buttefly Salamis parhassus to shift its mate search behaviour A c c e p t e d M a n u s c r i p t from perching to patrolling [START_REF] Bonte | Mate-locating behaviour, habitat-use, and flight morphology relative to rainforest disturbance in an Afrotropical butterfly[END_REF]. Rising temperature and climate change influence in turn the production and reception of sex pheromones [START_REF] Boullis | Will climate change affect insect pheromonal communication?[END_REF], as well as the choice of oviposition site [START_REF] Sentis | The smell of change: warming affects species interactions mediated by chemical information[END_REF].

Choice of micro-habitat

The choice of micro-habitat in response to habitat and climate change has received much attention. For example, [START_REF] Burdick | Density and distribution of soybean aphid, Aphis glycines Matsumura (Hemiptera : Aphididae ) in response to UV radiation[END_REF] [START_REF] Burdick | Density and distribution of soybean aphid, Aphis glycines Matsumura (Hemiptera : Aphididae ) in response to UV radiation[END_REF] showed that aphids modify their feeding location under increased UV radiation by feeding on the underside of leaves, as this protects against radiation. [START_REF] Alford | Why is there no impact of the host species on the cold tolerance of a generalist parasitoid?[END_REF] [START_REF] Alford | Why is there no impact of the host species on the cold tolerance of a generalist parasitoid?[END_REF] showed that aphid parasitoids increase their thermal tolerance and, hence, resistance to climate stress by manipulating the settling place on parasitized aphids.

Activity

As mentioned previously, [START_REF] Beaulieu | Behavioural antioxidant strategies to cope with high temperatures: a study in a tropical butterfly[END_REF] [START_REF] Beaulieu | Behavioural antioxidant strategies to cope with high temperatures: a study in a tropical butterfly[END_REF] recorded a decrease in daily activity of the tropical butterfly on hot days. Physical activity elevates metabolic rate and the production of reactive oxygen species (ROS) and reduced activity consequently limits these effects. The ant Iridomyrmex purpureus (Formicidae) employs a range of strategies to overcome high temperature stress, such as adjusting time of foraging to the colder hours of the day, climbing grass blades to cool down, and foraging only within shaded areas [START_REF] Andrew | Can temperate insect take the heat? A case study of the physiological and behavioural responses in a common ant, Iridomyrmex purpureus (Formicidae), with potential climate change[END_REF].

The reaction to thermal stress varies among species. For example, [START_REF] Kruse | Temperature and prey capture: opposite relationships in two predator taxa[END_REF] [START_REF] Kruse | Temperature and prey capture: opposite relationships in two predator taxa[END_REF] showed that spiders increase their activity at high temperatures, whereas carabid beetles decrease their activity. Since both are predators on the same prey (flies), these opposite reactions modify the predation rate on the prey.

At the seasonal level, Tougeron et al. (2017) [START_REF] Tougeron | Are aphid parasitoids from mild winter climates losing their winter diapause?[END_REF] observed a change in overwintering strategy following winter warming in aphid parasitoid wasps, from diapause to active state. This plasticity allows individuals to increase their progeny production. The mechanism behind the change is a decrease in responsiveness to environmental signals, rather than a genetic loss of diapause. [START_REF] Poethke | Predator-induced dispersal and the evolution of conditional dispersal in correlated environments[END_REF] showed with a model that any persistent deterioration of the environment that decreases expected fitness can induce dispersal behaviour in philopatric organisms. The model was developed for aphids responding to declining population viability due to increased predation pressure, but the model can be applied to any environment deterioration that reduces population viability.

Dispersion

Poethke et al. (2010)

Interspecific interactions

A c c e p t e d M a n u s c r i p t

Interactions between species depend on local conditions and, hence, are sensitive to environmental change [START_REF] Hoover | Behavioural Responses to a Changing World[END_REF][START_REF] Valiente-Banuet | Beyond species loss: the extinction of ecological interactions in a changing world[END_REF].

Pesticides, for instance, shift the competitive relationship between two thrips species so that an intrinsically inferior species Frankliniella occidentalis displaces an intrinsically superior competitor Thrips tabaci [START_REF] Zhao | Pesticide-mediated interspecific competition between local and invasive thrips pests[END_REF]. Climate change and increases in temperature are similarly altering species interactions. Ground beetles, for instance, attack more mobile prey at higher temperatures [START_REF] Vucic-Pestic | Warming up the system: higher predator feeding rates but lower energetic efficiencies[END_REF], while the parasitoid Aphidius rhopalosiphi attack fewer aphid hosts [START_REF] Lann | Thermal change alters the outcome of behavioural interactions between antagonistic partners[END_REF].

The invasion of alien species can profoundly alter species interactions and thereby the dynamics of populations [START_REF] Grether | Causes and consequences of behavioral interference between species[END_REF].

When the mosquito Aedes albopictus -a vector of Zika, dengue, and chikungunya viruses -invaded North America, it displaced the resident Ae. aegypti through interspecific matings that sterilized the resident species [START_REF] Bargielowski | Satyrization and satyrization-resistance in competitive displacements of invasive mosquito species[END_REF]. However, invaders can as well have positive effects on ecosystems. When insect pollinators were introduced to the Galapágos Islands, these made more visits to plants than the native species, which improved the stability of the interaction network [START_REF] Traveset | Invaders of pollination networks in the Galápagos Islands: emergence of novel communities[END_REF].

The loss of species can similarly alter species interactions. When large herbivores were lost from the African savanna, Acacia trees invested less in food and shelter rewards for mutualistic ants, which defend the trees against large herbivores and insect pests. A non-mutualistic ant could then occupy the trees, which caused the trees to suffer increased attack by insect pests [START_REF] Palmer | Breakdown of an ant-plant mutualism follows the loss of large herbivores from an African savanna[END_REF].

Social behaviour

Social behaviour is an important determinant of fitness in social insects. Thus, alteration of social behaviour because of environmental change can have profound demographic consequences [START_REF] Blumstein | Behavioural responses to a changing world: mechanisms and consequences[END_REF][START_REF] Schwander | Nature versus nurture in social insect caste differentiation[END_REF]. For instance, higher temperature induces a switch from solitary to social lifestyle in sweet bees Halictus rubicundus, which increases the number of pollinators, which in turn could mitigate the current pollinator crisis [START_REF] Schurch | Consequences of a warming climate for social organisation in sweat bees[END_REF]. On the other hand, a rise in temperature can as well disrupt social interactions. For example, higher temperature increases variation in worker size in the ant Temnothorax nylanderi, probably through reduced ability of nestmate workers to regulate larval development [START_REF] Molet | Effect of temperature and social environment on worker size in the ant Temnothorax nylanderi[END_REF].

Changes in the distribution and abundance of resources are common causes of altered frequency or type of social interaction. Changed shelter availability, for instance, influences the formation of social aggregations in the maritime earwig Anisolabis maritima [START_REF] Hack | Big wigs and small wigs: Time, sex, size and shelter affect cohabitation in the maritime earwig (Anisolabis maritima)[END_REF]. In the pharaoh's ant Monomorium pharaonic, changed resource distribution influences pheromone deposition along trails to food sources, which allows ants to communicate and reach adaptive collective decision in changing environments [START_REF] Czaczkes | Ants adjust their pheromone deposition to a changing environment and their probability of making errors[END_REF].

Is behavioural plasticity sufficient to cope with human perturbation?

A c c e p t e d M a n u s c r i p t

In this review, we have shown that behaviourally plastic response to HIREC are common. However, whether the plastic responses are sufficient to prevent population decline and extinction is poorly known. The research field is still in its infancy, and examples are scarce for insects. Some species are able to adjust their behaviour in an adaptive manner to changes, such as bark beetles that use multiple, integrated sensory modalities to adjust their choice of host trees to local conditions [START_REF] Raffa | Host selection by bark beetles: playing the odds in a high stakes game. Pine bark beetles[END_REF]. However, an increasing number of studies find species to respond maladaptively or insufficiently to various anthropogenic disturbances. For instance, plasticity in feeding behaviour of the harvester ant Pogonomyrmex mendozanus is insufficient for copying with deteriorating food quality [START_REF] Raffa | Host selection by bark beetles: playing the odds in a high stakes game. Pine bark beetles[END_REF]. Similarly, most species have limited potential to adjust behaviourally to extreme temperatures and many species are moving towards the poles to escape climate change [START_REF] Andrew | Can temperate insect take the heat? A case study of the physiological and behavioural responses in a common ant, Iridomyrmex purpureus (Formicidae), with potential climate change[END_REF][START_REF] Gunderson | Plasticity in thermal tolerance has limited potential to buffer ectotherms from global warming[END_REF].

Ecological traps can worsen survival prospects in changing environments. For instance, jewel beetles (Julodimorpha bakewelli) copulate with beer bottles because the texture of the bottles resembles that of females [START_REF] Gwynne | Beetles on the bottle -Male Buprestids mistake stubbies for females (Coleoptera)[END_REF], while mayflies oviposite on asphalt because the surface reflects polarized light in a similar manner as water bodies [START_REF] Kriska | Why do mayflies lay their eggs en masse on dry asphalt roads? Waterimitating polarized light reflected from asphalt attracts Ephemeroptera[END_REF].

Behavioural responses that are insufficient but still improve survival could facilitate evolutionary adaptation by providing more time for genetic changes [START_REF] Tuomainen | Behavioural responses to human-induced environmental change[END_REF]. This is especially likely in species with short generation time, as suggested for many pest and disease vectors [START_REF] Hoffman | Rapid adaptation of invertebrate pests to climatic stress?[END_REF].

Conclusions and perspectives

Behavioural plasticity is often not sufficient to cope with HIREC. One explanation is proposed by Sih (2013) [START_REF] Sih | Understanding variation in behavioural responses to human-induced rapid environmental change: A conceptual overview[END_REF]: if novel items have been poor options in the evolutionary history of a species, this could explain neophobic behaviours, while if novel items have been beneficial in the past, this could explain neophilic behaviours, which can help species to respond to HIREC. Because HIREC often results in novel conditions, not experienced in the recent evolutionary past, the responses may often be maladaptive and result in population decline. In rare cases, human induced changes can have a positive effect on an ecosystem's stability. For instance, in the Galapagos Islands, an introduced pollinators visited more plants than their native or endemics counterparts, which increased ecosystem stability [START_REF] Traveset | Invaders of pollination networks in the Galápagos Islands: emergence of novel communities[END_REF][START_REF] Wong | Behavioral responses to changing environments[END_REF].

Two fields of research emerge from this review that need more attention: (1) from a fundamental point of view, more research is needed on the effects of transgenerational plasticity on behavioural responses to environmental change, ie, when environments experienced by the parents influence the responses of offspring. Transgenerational plasticity is increasingly found to be involved in adjustments to variable environments, but the degree to which insects adjust to HIREC through transgenerational plasticity in behaviour is a largely unexplored topic. (2) From an applied point of A c c e p t e d M a n u s c r i p t view, the field of Conservation Behavior should receive more attention. This is a young discipline that investigates how the knowledge of animal behaviour can be taken in consideration in actions for preventing the loss of biodiversity [START_REF] Buchholz | Behavioural biology: an effective and relevant conservation tool[END_REF].

Figure legend Figure 1 (

 legend1 Figure legend

Acknowledgements

Thanks to Valérie Briand for bibliographic assistance. We would like to thank the section editors Eric Wajnberg and Emmanuel Desouhant for the invitation to contribute to the Behavioural ecology section 2018. This work was supported by the CNRS (Centre National de Recherche Scientifique) and by the University of Rennes 1 for Joan van Baaren, and by the Academy of Finland grant 277667 to Ulrika Candolin.