Gissell Estrada-Rodriguez
email: estradarodriguez@ljll.math.upmc.fr

Benoit Perthame
email: benoit.perthame@sorbonne-universite.fr

Motility switching and front-back synchronisation in polarized cells

The combination of protrusions and retractions in the movement of polarized cells leads to understand the effect of possible synchronisation between the two ends of the cells. This synchronisation, in turn, could lead to different dynamics such as normal and fractional diffusion. Departing from a stochastic single cell trajectory, where a "memory effect" induces persistent movement, we derive a kinetic-renewal system at the mesoscopic scale. We investigate various scenarios with different levels of complexity, where the two ends of the cell move either independently or with partial or full synchronisation. We study the relevant macroscopic limits where we obtain diffusion, drift-diffusion or fractional diffusion, depending on the initial system. This article clarifies the form of relevant macroscopic equations that describe the possible effects of synchronised movement in cells, and sheds light on the switching between normal and fractional diffusion.

Introduction

Mathematical modelling of cell motility has been largely studied, specially through the use of partial differential equations (PDEs). Depending on the biological context, cell trajectories can be described by a persistent random walk [START_REF] Othmer | Models of dispersal in biological systems[END_REF], where the individual tends to keep moving in the same direction as observed in [START_REF] Fricke | Persistence and adaptation in immunity: T cells balance the extent and thoroughness of search[END_REF][START_REF] Li | Persistent cell motion in the absence of external signals: a search strategy for eukaryotic cells[END_REF]. Often, directional persistence is also described as a correlated random walk where the direction of previous steps influences the direction of next ones.

This type of movement was studied in [START_REF] Huda | Lévy-like movement patterns of metastatic cancer cells revealed in microfabricated systems and implicated in vivo[END_REF] where, due to the synchronisation of protrusions and retractions in the front and back of metastatic cells, the authors observed a strong presence of long runs, interspersed by a sequence of short steps. These characteristics are analogous to Lévy walk trajectories, where the probability of a long run, i.e. a trajectory in the same direction for a long time, is non negligible. It was also observed that in non-metastatic cells the front and back movements are independent and then they follow a classical random walk. In contrast to a Brownian motion, where the distribution of the individuals' trajectories follows a Gaussian, Lévy walk trajectories asymptotically follow a power-law distribution [START_REF] Metzler | The random walk's guide to anomalous diffusion: a fractional dynamics approach[END_REF][START_REF] Zaburdaev | Lévy walks[END_REF]. Moreover, while for the Brownian case the mean square displacement x 2 grows linearly with respect to time ( x 2 ∼ t), for the Lévy walk case we have x 2 ∼ t ζ where ζ ∈ (1, 2). The exponent ζ = 2 corresponds to ballistic transport while ζ = 1 corresponds to normal diffusion. When ζ ∈ (1, 2) we are in the superdiffusive regime. For the ubiquitous appearance of Lévy walk models in biological systems we refer to [START_REF] Ariel | Swarming bacteria migrate by Lévy walk[END_REF][START_REF] Harris | Generalized Lévy walks and the role of chemokines in migration of effector cd8 t cells: Supplementary Information[END_REF][START_REF] Korobkova | From molecular noise to behavioural variability in a single bacterium[END_REF] at the cellular level, [START_REF] Focardi | Adaptive Lévy walks in foraging fallow deer[END_REF][START_REF] Reynolds | The weierstrassian movement patterns of snails[END_REF][START_REF] Reynolds | Displaced honey bees perform optimal scale-free search flights[END_REF][START_REF] Sims | Scaling laws of marine predator search behaviour[END_REF][START_REF] Viswanathan | Lévy flight search patterns of wandering albatrosses[END_REF] for animals and [START_REF] Raichlen | Evidence of Lévy walk foraging patterns in human hunter-gatherers[END_REF][START_REF] Reynolds | Lévy foraging patterns of rural humans[END_REF] for humans.

In this work, we start from the most general description: when front and back can make independent, non-synchronised steps to the right and to the left. In this setting, the model records the persistence time in each direction, thus leading to a complex system which can better be understood in terms of cell elongation and movement of the center of gravity.

To reduce the complexity, we assume that the cell length is fixed at the mesoscopic scale under investigation; the model is now amenable to multi-scale analysis and depending on the switching rates of both ends, we compute a diffusion coefficient at the macroscopic scale. Such a model with a given cell length can also be derived from the full system assuming fast switches in the dynamics of cell elongation compared to those for the movement. We also investigate the synchronisation between back and front, where we can observe a full range of behaviours. With low persistence, normal diffusion occurs, possibly with a drift and explicit formulas are computed. They show the possibility of a "backward diffusion" regime which can be interpreted as instability of the constant steady state at the kinetic (mesoscopic) scale. More interesting is when persistence is higher, then fractional diffusion occurs exhibiting the phenomena of long excursions which motivates our study.

When going back to individual cell (microscopic) dynamics, we perform numerical simulations and, as explained earlier, the long time behaviour quantifies the fractional diffusion in accordance to the developed theory.

To motivate our modelling, we are going to describe, based on biological evidence presented in [START_REF] Huda | Lévy-like movement patterns of metastatic cancer cells revealed in microfabricated systems and implicated in vivo[END_REF], two systems where none or full synchronisation of front and back leads to diffusive or superdiffusive dynamics, respectively.

Cell persistent movement through front-back synchronisation As previously introduced, the full synchronisation in space and time of the front and back movements of cells leads to Lévy walk dynamics, as observed in [START_REF] Huda | Lévy-like movement patterns of metastatic cancer cells revealed in microfabricated systems and implicated in vivo[END_REF] for the case of metastatic cancer cells. On the other hand, when the cells' front and back movements are independent, the trajectories follow a normal diffusion.

The synchronisation is translated into a persistence of the movement in a given direction. This means that at a given time step, the ends of the cell move forward or backwards simultaneously. The non-synchronisation is the result of independent movements in the front and back which might result in an intermittent cell movement pattern, where the cell length can vary. Hence, we assume that when the steps are not independent, and the cell keeps some "memory" of previous steps (non-Markovian process), the distribution of persistence length corresponds to a power-law and the movement is described by a Lévy walk. If the steps are independent, then the trajectories correspond to a diffusive movement. A more detailed description of this persistent and non-persistent movement is given in Section 1 (see also Supplementary Information in [START_REF] Huda | Lévy-like movement patterns of metastatic cancer cells revealed in microfabricated systems and implicated in vivo[END_REF]).

The aim of this article is to derive macroscopic equations that characterise these dynamics, starting from a kinetic description of the individual movement. To study this behaviour we propose the following setting. We consider that a cell is approximated by a one dimensional system of two identical point masses attached by an elastic cord. Each point mass is going to represent the front and back of the cell and the elastic cord represents the cell length. We aim to describe the trajectories of the whole system as in Figure 1. Each mass is considered as a material point that takes discrete (in space and time) infinitesimally small steps to the left or to the right with a certain probability. For the non-synchronisation case (Figure 1a) we consider the independent movement of the front (represented by y) and back (represented by x) and the cell length change, where c is the length at rest (see Subsection 2.1). While studying the synchronised movement leading to superdiffusion, since the front and back simultaneously move in the same direction at each time, we only consider the change in position of the centre of mass, represented by x in Figure 1b. Since the cell length is fixed in this case, we do not take it into account.

Outline of the paper

In Section 1 we present a detailed description of the individual movement as well as the main modelling assumptions. In Section 2 we introduce a general model where front and back movements are not synchronised. We also discuss some general notions as conservation of particles and realistic cell length for the model. Section 3 presents a simplified version of the previous model where we fix the cell length at the mesoscopic scale. We derive macroscopic equations depending on the switching rates at the ends of the cell. This model is also obtained by considering a fast switching dynamics in the full system as described in Subsection 3.1. In Section 4 we consider the system with synchronised front and back movement. We study the diffusive regime, when the persistence is low, and the superdiffusive regime in the opposite case (Section 5). Here we also show the possibility of a "backwards diffusion" regime. Finally, in Section 6 we perform numerical simulations at the individual level to compare with the developed theory.

Description of mathematical model

Front-back non-synchronisation In this case, the front and back movements are independent but with probabilities sampled from the same probability distribution. We consider that the front (y) gives a step k y of size δ at a given time τ , and similarly for the back x. The cell length is given by |y -x| where c denotes the equilibrium length. If the whole system is moving to the right initially, the cell front is allowed to reverse direction and move to y -δ with probability q ky = τ p ky or it can keep moving in the same direction with probability qky = 1 -τ p ky . Similarly, the back of the cell can move in the positive direction x + δ with probability qkx = 1 -τ p kx or reverse direction with probability q kx = τ p kx . The length of the cell varies in a specific range so that we preserve the physical properties, as described later in Section 2.1.

Front-back synchronisation For this case, since the front and back of the cell move simultaneously, i.e. at each time step they move to the left or to the right one step size with the same probability, the cell length is fixed at all times and we can consider the whole system as a point particle. The elastic cord can be considered as a solid road of fixed length and we study the movement of the center of mass only as in Figure 1b.

Analogous to the previous description, if we assume that the system is initially moving to the right, then it can reverse direction to x -δ with probability q k = τ p k or it can keep moving forward to x + δ with a probability qk = 1 -τ p k . Every time the cell changes direction we set k = 0 and we start counting again. The probability of keep moving without reversing will algebraically decrease with k. Note: From now on, we use the notation p k when k is discrete (see also Subsection 6.1) and p(k) when k is continuous.

Reversing direction probability

As described earlier, the rate at which the cell changes direction from left to right movement, given by p k , depends on the number of steps k given in that direction. The reversing rate is associated with a probability ψ(k), which is given by

ψ(k) = e -k 0 p(k * ) dk * . ( 1
)
This function ψ is often referred to as the survival probability, i.e., it gives the probability that the event of interest, in this case the reverse in direction, has not occurred for k steps. Equation [START_REF] Ariel | Swarming bacteria migrate by Lévy walk[END_REF] means that the probability of moving for k steps without changing is equal to the exponential of the cumulative reversing frequency. As indicated in [START_REF] Huda | Lévy-like movement patterns of metastatic cancer cells revealed in microfabricated systems and implicated in vivo[END_REF], for the case of metastatic cells this probability decays algebraically with k, therefore, here we consider that

ψ(k) = k 0 k 0 + k µ for µ ∈ (1, 3) . (2)
The reversing direction rate p(k) can also be expressed as the ratio

p(k) = φ(k) ψ(k) = -∂ k ψ(k) ψ(k) , (3)
where φ(k) is a probability density function. The above expression means that the reversing rate at step k equals the density of the event divided by the probability of keep moving in the same direction for k steps.

Non-synchronised movement description

We assume that the probabilities of the protrusions and retractions are independent from each other and therefore we have four different scenarios as in Figure 2. We denote by α(t, x, y, k x , k y ) and β(t, x, y, k x , k y ) the cells that are moving to the right and left, respectively, and by δ(t, x, y, k x , k y ) and γ(t, x, y, k x , k y ) the cells that only change their length by elongation and contraction. Moreover, the steps at the cell front, denoted by k y , are independent from the steps at the back, k x , and consequently, we have to take into account the rates p α x , p α y , p β x , p β y , p γ x , p γ y , p δ x and p δ y . These rates do not only depend on k x and k y but also on the distance |x -y| to preserve the cell physical size, as we discuss in Section 2.1. When the front and back of the cell change direction simultaneously, p x = p y , we denote the corresponding switching rates as p α , p β , p γ , p δ .

Note that cases (a) and (b) in Figure 2 are analogous to the synchronisation case later discussed in Section 4. In full generality, the number density of each cell population is described by the following systems coupled through the boundary terms,

         (∂ t + v∂ x + v∂ y + ∂ kx + ∂ ky )α = -(p α x + p α y )α-p α α , α(t, x, y, k x = 0, k y ) = ∞ 0 p γ x γ dk x , α(t, x, y, k x , k y = 0) = ∞ 0 p δ y δ dk y , α(t, x, y, k x = 0, k y = 0) = ∞ 0 ∞ 0 p β β dk x dk y , (4)
         (∂ t -v∂ x -v∂ y + ∂ kx + ∂ ky )β = -(p β x + p β y )β-p β β , β(t, x, y, k x = 0, k y ) = ∞ 0 p δ x δ dk x , β(t, x, y, k x , k y = 0) = ∞ 0 p γ y γ dk y , β(t, x, y, k x = 0, k y = 0) = ∞ 0 ∞ 0 p α α dk x dk y , (5)
         (∂ t -v∂ x + v∂ y + ∂ kx + ∂ ky )γ = -(p γ x + p γ y )γ-p γ γ , γ(t, x, y, k x = 0, k y ) = ∞ 0 p α x α dk x , γ(t, x, y, k x , k y = 0) = ∞ 0 p β y β dk y , γ(t, x, y, k x = 0, k y = 0) = ∞ 0 ∞ 0 p δ δ dk x dk y , (6)
         (∂ t + v∂ x -v∂ y + ∂ kx + ∂ ky )δ = -(p δ x + p δ y )δ-p δ δ , δ(t, x, y, k x = 0, k y ) = ∞ 0 p β x β dk x , δ(t, x, y, k x , k y = 0) = ∞ 0 p α y α dk y , δ(t, x, y, k x = 0, k y = 0) = ∞ 0 ∞ 0 p γ γ dk x dk y . (7)
The above systems of equations describe the different jumping combinations represented in Figure 2. For instance, in the case of population α, if x changes direction with certain rate p α

x , then we have a transition from population α to γ (Figure 2 (a)→(c)), which is given by γ(t, x, y, k x = 0, k y ).

Similarly, if y changes direction in population α we have a change from α to δ, and the individuals that leave population α appear in δ(t, x, y, k x , k y = 0). The reverse process also happens and transitions from population γ to α (cells appear at α(t, x, y, k x = 0, k y )) and from δ to α (cells appear at α(t, x, y, k x , k y = 0)) are also considered. When the jumping direction changes at x and y simultaneously then it always happens that populations switch directly from α ←→ β and γ ←→ δ.

Notation: For simplicity of notation, in the rest of the paper we use

(•) = ∞ 0 ∞ 0 (•) dk x dk y .
In the following we are going to check some physical properties that the systems (4)- [START_REF] Focardi | Adaptive Lévy walks in foraging fallow deer[END_REF] must satisfy.

Coordinates of the center of mass and cell elongation At a first stage, a desirable property is that cell polarization is preserved, that means x < y all along the movement assuming it is true initially. To examine the conditions which enforce this property it is easier to use the coordinates of the center of mass and the distance between back and front. For that we let x+y 2 = X and y -x = z, the cell length. From now on, and for simplicity in the notation, we keep the same functions α, β, γ, δ that will depend on the new variables (t, X, z, k x , k y ). From the system (4)- [START_REF] Focardi | Adaptive Lévy walks in foraging fallow deer[END_REF] we get

(∂ t + 2v∂ X + ∂ kx + ∂ ky )α = -(p α x + p α y )α -p α α , (8)
(∂ t -2v∂ X + ∂ kx + ∂ ky )β = -(p β x + p β y )β -p β β , (9)
(∂ t + v∂ z + ∂ kx + ∂ ky )γ = -(p γ x + p γ y )γ -p γ γ , (10)
(∂ t -v∂ z + ∂ kx + ∂ ky )δ = -(p δ x + p δ y )δ -p δ δ , (11)
with the boundary conditions in k

x , k y , α(•, k x = 0, k y ) = ∞ 0 p γ x γ dk x , α(•, k x , k y = 0) = ∞ 0 p δ y δ dk y , α(•, k x = 0, k y = 0) = p β β , (12)
β(•, k x = 0, k y ) = ∞ 0 p δ x δ dk x , β(•, k x , k y = 0) = ∞ 0 p γ y γ dk y , β(•, k x = 0, k y = 0) = p α α , (13)
γ(•, k x = 0, k y ) = ∞ 0 p α x α dk x , γ(•, k x , k y = 0) = ∞ 0 p β y β dk y , γ(•, k x = 0, k y = 0) = p δ δ , (14)
δ(•, k x = 0, k y ) = ∞ 0 p β x β dk x , δ(•, k x , k y = 0) = ∞ 0 p α y α dk y , δ(•, k x = 0, k y = 0) = p γ γ . (15)
Here (•) denotes the dependence on (t, X, z). In order to guarantee that z > 0 is preserved, we also need to ensure that δ(X, z = 0, k x , k y ) = 0 (similarly γ(X, z = 0, k x , k y ) = 0), which means that the jump rate p δ x + p δ y → ∞ as z → 0 and

• 0 (p δ x + p δ y )(z) dz = ∞.
Using the notation in [START_REF] Metzler | The random walk's guide to anomalous diffusion: a fractional dynamics approach[END_REF], this means that

• 0 µ δ (z) dz = ∞.
Conservation of particles Integrating with respect to k x and k y we define the macroscopic density

ᾱ(t, X, z) = α(•, k x , k y ) , (16)
and similarly for β, γ and δ. Moreover, integrating with respect to k x and k y equations ( 8)- [START_REF] Harris | Generalized Lévy walks and the role of chemokines in migration of effector cd8 t cells: Supplementary Information[END_REF] and adding them together we obtain the following macroscopic conservation equation

∂ t (u + w) + 2v∂ X j + v∂ z m = 0 . (17)
Here u(t, X, z) = ᾱ + β is the moving population, w(t, X, z) = γ + δ is the resting population, j(t, X, z) = ᾱ -β is the mean direction of motion and m(t, X, z) = γ -δ is the mean extension rate.

Biologically relevant switching probabilities

The switching rate is not going to depend only on the persisting steps k x and k y but also on the cell length z. Using ( 2) and ( 3) we can write the general expression

p(k) = µ 1 + k . ( 18
)
For the non-synchronised movement this rate is given by, for population α (and similarly for the rest),

p α x (k x ) = µ α (z) 1 + k x , p α y (k y ) = µ α (z) 1 + k y . ( 19
)
The dependence of µ on z guarantees that we keep a realistic cell length as we will discuss below.

For the synchronised case, since the cell does not change shape, we consider the switching rate given by [START_REF] Mellet | Fractional diffusion limit for collisional kinetic equations[END_REF]. As described in [START_REF] Huda | Lévy-like movement patterns of metastatic cancer cells revealed in microfabricated systems and implicated in vivo[END_REF] the length of a cell can only vary in a certain range. Considering that the resting length is c = |y -x|, we define L max = 3.5 c and L min = 0.5 c and the switching rates satisfy the following properties,

p α x = if z L max , µα(z) 1+kx if z L min , p α y = µα(z) 1+ky if z L max , if z L min , (20)
p β x = µ β (z) 1+kx if z L max , if z L min , p β y = if z L max , µ β (z) 1+ky if z L min , (21)
and finally,

p γ x = µγ (z) 1+kx if z L max , if z L min , p γ y = µγ (z) 1+ky if z L max , if z L min , (22)
p δ x = if z L max , µ δ (z) 1+kx if z L min , p δ y = if z L max , µ δ (z) 1+ky if z L min . (23)
Here is a small parameter. Let us take for instance the case when cells are moving right (population α). In the limit when z L max , since we want to preserve the physical length of the cell, the front has to change direction with a higher rate (p α y ) while the back should keep moving in the same direction (small p α x ). On the other hand, in the limit when z L min the back of the cell has to switch direction (p α x ) while the front should keep moving without changing (small p α y ). The opposite happens when the cell is moving to the left (population β).

For the case when the cell is at rest (γ and δ), if z L max for the case of the population γ, the switching rate has to be very high at both ends so that the cell recovers the resting length c . If z L min , then p γ x , p γ y are very small. The opposite happens in the case of the population δ.

Simplified system with resting population

The systems ( 8)-( 11) are very complex to analyse since they involve different dynamics such as left and right movement for four different populations, and additionally, the change in cell length. Therefore, in this section we consider a simplified model of three populations: cells moving left (β), cells moving right (α) and resting cells (γ 0 ). The population γ 0 represents the average of the populations γ and δ described before assuming the mean cell length is constant, thus ignoring the variable z. The dynamics are described now by

ε 2 ∂ t α ε + εv∂ x α ε + ∂ k α ε = -p(k)α ε , α ε (t, 0, x) = ζ β ∞ 0 p(k)β ε dk + ζ γ γ ε 0 (t, x, 0) , (24)
ε 2 ∂ t β ε -εv∂ x β ε + ∂ k β ε = -p(k)β ε , β ε (t, 0, x) = ζ α ∞ 0 p(k)α ε dk + (1 -ζ γ )γ ε 0 (t, x, 0) , (25)
ε 2 ∂ t γ ε 0 -∂ k γ ε 0 = (1 -ζ α )p(k)α ε + (1 -ζ β )p(k)β ε . (26)
The switching rate p(k) describes both the change in direction of the center of mass and transition to rest state γ 0 where memory is gradually lost. This replaces the movement of the front and the back of the cell as in Section 2 and is given by ( 18). Since we are not considering front and back movements then α, β and γ 0 depend only on (t, x, k). In the above system, ζ ∈ (0, 1) is a probability and note that we have introduced a diffusive scaling (t, x) → ( t/ε 2 , x/ε). The individuals from population α that switch direction with rate p(k) either start moving in the opposite direction with probability ζ α , represented by the first term in β ε (t, x, 0), or they go into a resting phase with probability 1-ζ α , given by the first term in the right hand side of [START_REF] Reynolds | Displaced honey bees perform optimal scale-free search flights[END_REF]. A similar dynamic is followed by individuals in population β. On the other hand, individuals that are at rest, population γ 0 , start to move to the right, with probability ζ γ , or to the left with probability 1 -ζ γ . As for the systems ( 4)-( 7), we can easily check that ( 24)-( 26) preserves the number density of individuals.

When the rate p(k) is large enough for large k, more precisely when ke -k 0 p(k * ) dk * is integrable, then the large scale dynamics is normal diffusion. To explain this, we define the survival probability ψ as

ψ(k) = Ze -k 0 p(k * ) dk * , Z -1 := ∞ 0 e -k 0 p(k * ) dk * dk , ∞ 0 p(k)ψ(k) dk = Z . (27)
For

p(k) = µ 1+k with µ > 2, then indeed ke -k 0 p(k * ) dk * = k (1+k) µ is integrable.
The first question is to determine under which conditions diffusion occurs in the small scale regime for ε. To do that, we compute the limiting α, β, γ 0 as ε → 0, using the solution of

                     ∂ k α = -p(k)α , α(t, x, 0) = ζ β ∞ 0 p(k)β dk + ζ γ γ 0 (t, x, 0) , ∂ k β = -p(k)β , β(t, x, 0) = ζ α ∞ 0 p(k)α dk + (1 -ζ γ )γ 0 (t, x, 0) , -∂ k γ 0 = (1 -ζ α )p(k)α + (1 -ζ β )p(k)β . (28)
Therefore, we obtain the limits

α(t, x, k) = α(t, x, 0) ψ(k) Z = ᾱ(t, x)ψ(k) , β(t, x, k) = β(t, x, 0) ψ(k) Z = β(t, x)ψ(k) , (29)
γ 0 (t, x, k) = (1 -ζ α ) ᾱ(t, x)ψ(k) + (1 -ζ β ) β(t, x)ψ(k) , where ᾱ(t, x) = ∞ 0 α dk and β(t, x) = ∞ 0 β dk. The γ 0 (t, x, 0) is given by γ 0 (t, x, 0) = (1 -ζ α ) ᾱ(t, x)ψ(0) + (1 -ζ β ) β(t, x)ψ(0) . (30)
With these expressions, we can compute a relation between ᾱ(t, x) and β(t, x) starting from

ᾱ(t, x) = α(t, x, 0) Z = ζ β Z ∞ 0 p(k)β(t, x, k) dk + ζ γ Z γ 0 (t, x, 0) = ζ β β(t, x) + ζ γ [(1 -ζ α )ᾱ(t, x) + (1 -ζ β ) β(t, x)] ,
that finally gives ᾱ(t, x)

= ζ β +ζ γ (1-ζ β ) 1-ζ γ (1-ζ α ) β(t, x)
. Therefore, the condition for a diffusive limit, i.e., ᾱ(t, x) = β(t, x) turns out to be

ζ γ = 1 -ζ β 2 -ζ α -ζ β . (31)
Then, the second question is to compute the diffusion coefficient. The macroscopic conservation equation is obtained using

∂ t ( ᾱε + βε + γε 0 ) + v∂ x J ε = 0 , where J ε := ᾱε -βε ε .
The difficulty here is to compute the flux J ε when p(k) is not constant. In the following we are going to compute the diffusive limit under the condition (31), i.e., such that, as ε → 0

ᾱ(t, x) = β(t, x) , for ζ γ = 1 -ζ β 2 -ζ α -ζ β , and J ε → -Dv∂ x β(t, x) .
We start by computing the Taylor expansion of α ε using the equation

ε 2 ∂ t α ε (t, k, x) ψ(k) + εv∂ x α ε (t, k, x) ψ(k) + ∂ k α ε (t, k, x) ψ(k) = 0 .
Integrating with respect to k we find,

α ε (t, x, k) ψ(k) = α ε (t, x, 0) ψ(0) -εv∂ x k 0 α ε (t, x, k * ) ψ(k * ) dk * + O(ε 2 ) = α ε (t, x, 0) ψ(0) -εkv∂ x α ε (t, x, 0) ψ(0) + O(ε 2 ) = α ε (t, x, 0) ψ(0) -εkv∂ x ᾱε + O(ε 2 ) ,
since, integrating in k after multiplying by ψ(k), we find ᾱε = αε(t,x,0) ψ(0)

+ O(ε). For β ε , we obtain

β ε (t, x, k) ψ(k) = β ε (t, x, 0) ψ(0) + εkv∂ x βε + O(ε 2 ) . ( 32
)
From the above equations the aim is to compute

J ε = ᾱε -βε ε = ∞ 0 ψ(k) α ε (t, x, 0) -β ε (t, x, 0) εψ(0) dk -v∂ x ∞ 0 kψ(k)(ᾱ ε + βε ) dk + O(ε) = α ε (t, x, 0) -β ε (t, x, 0) εψ(0) -v∂ x (ᾱ ε + βε ) ∞ 0 kψ(k) dk + O(ε) . (33)
In order to compute the term α ε (t, x, 0) -β ε (t, x, 0), we first multiply by p(k)ψ(k) in (32) and integrate in k. This gives, using that

∞ 0 kp(k)ψ(k) dk = 1, ∞ 0 p(k)β ε (t, x, k) dk = β ε (t, x, 0) + εv∂ x βε ∞ 0 kp(k)ψ(k) dk + O(ε 2 ) = β ε (t, x, 0) + εv∂ x βε + O(ε 2 ) . ( 34
)
The boundary condition α ε (t, x, 0) in ( 24) can be written as, after using ( 30) and (34),

α ε (t, x, 0) = (ζ β + ζ γ (1 -ζ β )) ∞ 0 p(k)β ε (t, x, k) dk + ζ γ (1 -ζ α ) ∞ 0 p(k)α ε (t, x, k) dk = (ζ β + ζ γ (1 -ζ β ))[β ε (t, x, 0) + εv∂ x βε ] + ζ γ (1 -ζ α )[α ε (t, x, 0) -εv∂ x ᾱε ] + O(ε 2 ) ,
and we obtain, for c :=

1 1-ζ γ (1-ζ α ) = 1 ζ β +ζ γ (1-ζ β ) thanks to (31), α ε (t, x, 0) = β ε (t, x, 0) + εv∂ x βε -cζ γ (1 -ζ α )εv∂ x ᾱε + O(ε 2 ) = β ε (t, x, 0) + εv∂ x βε + (1 -c)εv∂ x ᾱε + O(ε 2 ) .
From here we obtain, with ᾱ = β the limits of ᾱε and βε ,

lim ε→0 α ε (t, x, 0) -β ε (t, x, 0) ε = (2 -c)v∂ x ᾱ .
Finally, we write (33) and compute the diffusion coefficient

J = -v∂ x ᾱ(t, x) 2 ∞ 0 ψ(k)k dk - 2 -c Z , D = 2 ∞ 0 ψ(k)k dk - 2 -c Z . ( 35
)
Note that, by opposition to the formalism developed in [START_REF] Frank | On a generalized Boltzmann equation for non-classical particle transport[END_REF], this diffusion coefficient is not always positive. This is because the authors also rescale k in such a way to give more weights to the large values of k. When D is negative, some instability arises for the kinetic model, which is analysed in Section 4.2.

Partial synchronisation limit

We may also assume that moving forward is more effective than elongating and shortening. To represent that, we may derive a partial synchronisation limit starting from ( 8)-( 11). This limit consists on introducing fast transition rates for conformations γ and δ so that the whole system approximately converges to the two moving populations α and β (as in ( 50) and (51) in Section 4). We start by introducing the following scaling

(∂ t + v∂ z + ∂ kx + ∂ ky )γ = - 1 ε (p γ x + p γ y )γ -p γ γ , (36)
(∂ t -v∂ z + ∂ kx + ∂ ky )δ = - 1 ε (p δ x + p δ y )δ -p δ δ , (37)
and we change accordingly the boundary conditions ( 12) and [START_REF] Jara | Limit theorems for additive functionals of a Markov chain[END_REF]. With this scaling, we find γ → 0 and δ → 0 as ε → 0. The difficulty is to compute the limiting contribution to the boundary terms

α(•, k x = 0, k y ) = 1 ε ∞ 0 p γ x γ(•, k x , k y ) dk x .
To do so, we use the method of characteristics in [START_REF] Fricke | Persistence and adaptation in immunity: T cells balance the extent and thoroughness of search[END_REF] and [START_REF] Harris | Generalized Lévy walks and the role of chemokines in migration of effector cd8 t cells: Supplementary Information[END_REF] and neglect the initial contribution which is immediately absorbed due to our scaling. We find, respectively,

γ(•, k x , k y ) =        γ(t -k x , X, z -vk x , k x = 0, k y -k x )e -kx 0 p γ x +p γ y ε +p γ dk * x for k x < k y , γ(t -k y , X, z -vk y , k x -k y , k y = 0)e - ky 0 p γ x +p γ y ε +p γ dk * y for k y < k x , (38)
δ(•, k x , k y ) =        δ(t -k x , X, z + vk x , k x = 0, k y -k x )e -kx 0 p δ x +p δ y ε +p δ dk * x for k x < k y , δ(t -k y , X, z + vk y , k x -k y , k y = 0)e - ky 0 p δ x +p δ y ε +p δ dk * y for k y < k x . ( 39
)
We may estimate these quantities thanks to the Laplace approximation in the regime where ε → 0,

1 ε ∞ 0 p γ x γ(•, k x = 0, k y )e -kx 0 p γ x +p γ y ε +p γ dk * x dk x → p γ x (k x = 0) p γ x (k x = 0) + p γ y γ(•, k x = 0, k y ) , 1 ε ∞ 0 p γ x γ(•, k x , k y = 0)e - ky 0 p γ x +p γ y ε +p γ dk * y dk x → δ 0 (k y ) ∞ 0 p γ x p γ x + p γ y (k y = 0) γ(•, k x , k y = 0) dk x . (40)
Note that δ 0 (k y ) represents a Dirac delta function in k y while δ(•, k x , k y ) is the density of individuals moving as in Figure 2. Substituting these expressions into the first equation in [START_REF] Huda | Lévy-like movement patterns of metastatic cancer cells revealed in microfabricated systems and implicated in vivo[END_REF] and using [START_REF] Korobkova | From molecular noise to behavioural variability in a single bacterium[END_REF], we find the boundary condition

α(•, k x = 0, k y ) → p γ x (k x = 0) p γ x (k x = 0) + p γ y γ(•, k x = 0, k y ) + δ 0 (k y ) ∞ 0 p γ x p γ x + p γ y (k y = 0) γ(•, k x , k y = 0) dk x = p γ x (k x = 0) p γ x (k x = 0) + p γ y ∞ 0 p α x α dk x + δ 0 (k y ) ∞ 0 p γ x p γ x + p γ y (k y = 0) ∞ 0 p β y β dk y dk x . (41)
Similarly, we obtain from the second equation in ( 12)

α(•, k x , k y = 0) = δ 0 (k x ) ∞ 0 p δ y p δ x (k x = 0) + p δ y ∞ 0 p β x β dk x dk y + p δ y (k y = 0) p δ x + p δ y (k y = 0) ∞ 0 p α y α dk y .
(42) We can simplify the above expressions by using the following notation

P γ x = p γ x p γ x + p γ y (k y = 0) , P γ y = p γ x (k x = 0) p γ x (k x = 0) + p γ y , P δ x = p δ y (k y = 0) p δ y (k y = 0) + p δ y , P δ y = p δ y p δ x (k x = 0) + p δ x . (43)
For the population β, we follow the same steps, starting from

β(•, k x = 0, k y ) = 1 ε ∞ 0 p δ x δ dk x . Using (40) we have β(•, k x = 0, k y ) = p δ x (k x = 0) p δ x (k x = 0) + p δ y ∞ 0 p β x β dk x + δ 0 (k y ) ∞ 0 p δ x p δ x + p δ y (k y = 0) ∞ 0 p α y α dk y dk x , (44)
β(•, k x , k y = 0) = δ 0 (k x ) ∞ 0 p γ y p γ x (k x = 0) + p γ y ∞ 0 p α x α dk x dk y + p γ y (k y = 0) p γ y (k y = 0) + p γ x ∞ 0 p β y β dk y . (45)
Here, we similarly define

Pγ x = p γ y (k y = 0) p γ x + p γ y (k y = 0) , Pγ y = p γ y p γ x (k x = 0) + p γ y , Pδ x = p δ x p δ y (k y = 0) + p δ y , Pδ y = p δ x (k x = 0) p δ x (k x = 0) + p δ x . (46)
Using (41), (42), ( 44) and (45) and integrating with respect to k x , k y in ( 8) and ( 9) we obtain

(∂ t + 2v∂ X )ᾱ = - P 1 (k x , k y )α(•, k x , k y ) + P 2 (k x , k y )β(•, k x , k y ) , (47)
(∂ t -2v∂ X ) β = - P 3 (k x , k y )β(•, k x , k y ) + P 4 (k x , k y )α(•, k x , k y ) , (48)
where P 1 , P 2 , P 3 , P 4 are expressed in terms of ( 43) and (46). This system is analogous to (50)-(51) below that describes the left and right movement only, corresponding to the synchronisation case, with a modified jumping rate p(k) coming form the populations γ and δ.

To achieve conservation of particles we add (47) and (48) to obtain

∂ t u + 2v∂ X j = (-1 + P γ y + Pγ y )p α x α + (-1 + P δ x + Pδ x )p α y α + (-1 + P γ x + Pγ x )p β y β + (-1 + P δ y + Pδ y )p β x β = 0 , (49)
since the terms inside the brackets cancel out.

In conclusion, assuming fast transition in the "asynchronous states" γ and δ, we recover a simplified system where only the states α, β occur, that is the back and front are always synchronized. The new phenomena is the possibility of two fast transitions α → β → α (or symmetrically exchanging β and α) which modifies the boundary conditions compared to the model initially postulated.

Synchronised movement description

When the front and back movement of cells are synchronised, as described in Section 1, the system is much simpler but still can exhibit several remarkable features as oriented drift, instability or superdiffusive movement. As before, we denote by α(t, x, k) the probability that the cell moves to the right and by β(t, x, k) the probability that the cell moves to the left. The rate of changing the direction is denoted by p(k), defined in [START_REF] Mellet | Fractional diffusion limit for collisional kinetic equations[END_REF]. The system of equations that describes the synchronisation movement was derived in Appendix 6.1 from a discrete description and is given by

∂ t α(t, x, k) + v∂ x α(t, x, k) + ∂ k α(t, x, k) = -p(k)α(t, x, k) , α(t, x, 0) = ∞ 0 p(k)β(t, x, k) dk , (50)
∂ t β(t, x, k) -v∂ x β(t, x, k) + ∂ k β(t, x, k) = -p(k)β(t, x, k) , β(t, x, 0) = ∞ 0 p(k)α(t, x, k) dk . (51)

Normal diffusion limit of the synchronised system

We first study the scale in the memory term p(k) which leads to a usual diffusion equation, following the lines of Section 3. Because of the simplicity of the system, we may analyse the drift-diffusion behaviour in a more general context. To do so we re-scale (50)-(51) as follows,

         ε 2 ∂ t α ε (t, x, k) + ε v∂ x α ε (t, x, k) + ∂ k α ε (t, k, x) = -(p(k) + εp α (k))α ε (t, x, k) , ε 2 ∂ t β ε (t, x, k) -ε v∂ x β ε (t, x, k) + ∂ k β ε (t, x, k) = -(p(k) + εp β (k))β ε (t, x, k) , α ε (t, x, 0) = ∞ 0 (p(k) + εp β (k))β ε dk, β ε (t, x, 0) = ∞ 0 (p(k) + εp α (k))α ε dk .
(52)

The drift-diffusion limit is obtained using again the identity

∂ t ( ᾱε + βε ) + ∂ x J ε = 0 , J ε := ᾱε -βε ε , (53)
where we need to compute the x-flux J ε . We are going to compute the constants V and D such that, as ε → 0, ᾱ(t, x) = β(t, x) , and

J ε → V ᾱ(t, x) -Dv∂ x ᾱ(t, x) . (54)
We complete the system (52) with initial data such that α ε (0, x, k) = β ε (0, x, k), this is because the definition of the flux J ε requires a bounded quantity ᾱε-βε ε . As ε vanishes, we find limits that we denote by α, β and that satisfy

   ∂ k α(t, x, k) = -p(k)α(t, x, k) , α(t, x, 0) = ∞ 0 p(k)β(t, x, k) dk , ∂ k β(t, x, k) = -p(k)β(t, x, k) , β(t, x, 0) = ∞ 0 p(k)α(t, x, k) dk ,
which means that the limits are given by ( 27) and [START_REF] Viswanathan | Lévy flight search patterns of wandering albatrosses[END_REF]. Consequently, we deduce that

ᾱ(t, x) = α(t, x, 0) ψ(0) = 1 ψ(0) ∞ 0 p(k)β(t, x, k) dk = β(t, x) .
To compute the flux J ε we follow the steps in Section 3 and we write

α ε (t, x, k) ψ(k) = α ε (t, x, 0) ψ(0) -ε k 0 p α (k * ) α ε (t, x, k * ) ψ(k * ) dk * -ε v∂ x k 0 α ε (t, x, k * ) ψ(k * ) dk * + O(ε 2 ) ,
which can be re-arranged as

α ε (t, x, k) = ψ(k) α ε (t, x, 0) ψ(0) -εψ(k) k 0 p α (k * ) dk * α ε (t, x, 0) ψ(0) -ε kψ(k)v∂ x α ε (t, x, 0) ψ(0) + O(ε 2 ) . ( 55
)
Arguing in a similar way for β ε , and using the expression for α ε (t, x, 0) in ( 52) we find

α ε (t, x, 0) = ∞ 0 (p(k) + εp β (k))β ε dk = β ε (t, x, 0) 1 + ε Z β ψ(0) -ε β ε (t, x, 0) ψ(0) Y β + εv∂ x β ε (t, x, 0) ψ(0) ∞ 0 kp(k)ψ(k) dk + O(ε 2 ) . ( 56
)
Here we have used the normalization of ψ(k) introduced in ( 27) and

Z α := ∞ 0 p α (k)ψ(k) dk, Z β := ∞ 0 p β (k)ψ(k) dk, Y α,β := ∞ 0 p(k)ψ(k) k 0 p α,β (k * ) dk * dk.
We can further notice that ∞ 0 kp(k)ψ(k) dk = 1 as before. Writing an analogous expression of (56) but for β ε (t, x, 0) we can compute

α ε (t, x, 0)-β ε (t, x, 0) = β ε (t, x, 0) -α ε (t, x, 0) + ε Z β ψ(0) β ε (t, x, 0) -ε Z α ψ(0) α ε (t, x, 0) -ε β ε (t, x, 0) ψ(0) Y β - α ε (t, x, 0) ψ(0) Y α + εv∂ x β ε (t, x, 0) ψ(0) + α ε (t, x, 0) ψ(0) + O(ε 2 ) .
This shows that in the limit α(t, x, 0) = β(t, x, 0), ᾱ(t, x) = β(t, x) and thus α(t, x, k) = β(t, x, k), giving

2 lim ε→0 α ε (t, x, 0) -β ε (t, x, 0) ε = (Z β -Z α )ᾱ(t, x) -ᾱ(t, x)(Y β -Y α ) + 2v∂ x ᾱ(t, x) = 2v∂ x ᾱ(t, x) ,
since (and the same argument shows that

Y β = Z β ), Y α = - ∞ 0 dψ(k) dk k 0 p α (k * ) dk * = ∞ 0 ψ(k)p α (k) dk = Z α .
Back to (55), using [START_REF] Viswanathan | Lévy flight search patterns of wandering albatrosses[END_REF], we find, with

V = ∞ 0 ψ(k) k 0 (p α (k * ) -p β (k * )) dk * dk, lim ε→0 J ε = lim ε→0 α ε (t, x, 0) -β ε (t, x, 0) εψ(0) -ᾱ(t, x) ∞ 0 ψ(k) k 0 (p α (k * ) -p β (k * )) dk * dk -2v∂ x ᾱ(t, x) ∞ 0 kψ(k) dk = V ᾱ(t, x) -2v∂ x ᾱ(t, x) ∞ 0 kψ(k) 1 - p(k) 2Z dk .
We finally obtain the transport coefficients V and D

lim ε→0 J ε = V ᾱ -v D∂ x ᾱ(t, x), D = 2 ∞ 0 kψ(k) dk - 1 Z .
We recover the result of Section 3 when p α,β = 0, and ζ α,β = 1, then we find V = 0, c = 1 and D = D. The same comment on the positivity of D applies here. Remark 4.1. When ψ(k) is given by ( 2), then for k 0 = 1 we have that

∞ 0 kψ(k) dk = 1 (µ-2)(µ-1) > 0 for µ > 2. Moreover, it holds that 1 (µ -2)(µ -1) > 1 2Z for µ ∈ (2, 3)
and the diffusion coefficient D > 0.

Stability analysis

Consider, for simplicity, the system (52). When the above diffusion coefficient is negative, D < 0, we expect instability for the kinetic system when ε small enough. This phenomena has been already observed for chemotaxis and semilinear parabolic equations in [START_REF] Moussa | Backward parabolicity, cross-diffusion and Turing instability[END_REF][START_REF] Perthame | Stiff-response-induced instability for chemotactic bacteria and flux-limited Keller-Segel equation[END_REF] and we explain it in the present context.

To study the stable/unstable modes, we consider a simple Fourier mode α(t, x, k) = ᾱ(t, x) + a n (k)e λt e inx that we substitute in the equation for α in (52). We get

ε 2 λa n (k) + εvina n (k) + ∂ k a n (k) + p(k)a n (k) = 0 , which we can solve to obtain    a n (k) = a n (0)e -k 0 (p(k * )+ε 2 λ+εivn) dk * , a n (0) = ∞ 0 p(k)b n (k) dk .
Hence we have

a n (k) = ∞ 0 p(k)b n (k) dk e -k 0 (p(k * )+ε 2 λ+εivn) dk * , b n (k) = ∞ 0 p(k)a n (k) dk e -k 0 (p(k * )+ε 2 λ-εivn) dk * .
Substituting b n (k) into a n (k), multiplying by p(k) and integrating we obtain the dispersion relation

1 = ∞ 0 p(k)e -k 0 (p(k * )+ε 2 λ+εivn) dk * dk ∞ 0 p(k)e -k 0 (p(k * )+ε 2 λ-εivn) dk * dk . ( 57
)
For ε very small, we use a Taylor expansion and re-write (57) as

1 = ∞ 0 p(k)e -k 0 p(k * ) dk * (1 -εivnk -ε 2 λk - ε 2 2 n 2 k 2 v 2 ) dk × ∞ 0 p(k)e -k 0 p(k * ) dk * (1 + εivnk -ε 2 λk - ε 2 2 n 2 k 2 v 2 ) dk .
As before, see [START_REF] Sims | Scaling laws of marine predator search behaviour[END_REF], we may use that

∞ 0 p(k)e -k 0 p(k * ) dk * dk = 1, ∞ 0 kp(k)e -k 0 p(k * ) dk * dk = 1 Z .
Thus the terms of order 1 and ε cancel and the second order terms give

2 Z λ = -n 2 v 2 ∞ 0 p(k)e -k 0 p(k * ) dk * k 2 dk + 1 Z 2 v 2 n 2 . Computing ∞ 0 p(k)ψ(k)k 2 dk = 2 ∞ 0 ψ(k)k dk, we obtain λ = n 2 v 2 2 1 Z -2 ∞ 0 ψ(k)k dk = -D . ( 58
)
This condition shows that when D < 0, for ε small, the kinetic model is Turing unstable as in [START_REF] Moussa | Backward parabolicity, cross-diffusion and Turing instability[END_REF][START_REF] Perthame | Stiff-response-induced instability for chemotactic bacteria and flux-limited Keller-Segel equation[END_REF].

Fractional equation for the synchronised movement

In the context of system (50)-(51), using an appropriate scaling, we obtain a macroscopic fractional diffusion equation describing the persistent movement of the total population when the front and back of the cell are synchronised. When ψ(k) has a "fat tail", meaning µ ∈ (1, 2) for p(k) = µ 1+k , fractional diffusion occurs as already pointed out in [START_REF] Frank | Fractional diffusion limits of non-classical transport equations[END_REF]. We also recall that superdiffusion regimes are well established in different contexts of purely kinetic theory since the seminal works [START_REF] Jara | Limit theorems for additive functionals of a Markov chain[END_REF][START_REF] Mellet | Fractional diffusion limit for collisional kinetic equations[END_REF].

Kinetic system

We start by integrating (50) and (51) with respect to k, taking into account the boundary condition at k = 0,

∂ t ᾱ(t, x) + v∂ x ᾱ(t, x) = - t 0 p(k)α(t, x, k) dk + t 0 p(k)β(t, x, k) dk , (59)
∂ t β(t, x) -v∂ x β(t, x) = - t 0 p(k)β(t, x, k) dk + t 0 p(k)α(t, x, k) dk , (60)
where ᾱ(t, x) = t 0 α(t, x, k) dk and β(t, x) = t 0 β(t, x, k) dk . We also consider initial conditions β 0 (0, x, k) = β0 (x)δ(k) and α 0 (0, x, k) = ᾱ0 (x)δ(k). Now the aim is to write the right hand side of ( 59)-(60) in terms of the macroscopic densities ᾱ(t, x) and β(t, x). For that purpose we follow some steps from [START_REF] Estrada-Rodriguez | Fractional Patlak-Keller-Segel equations for chemotactic superdiffusion[END_REF] and [START_REF] Fedotov | Persistent random walk of cells involving anomalous effects and random death[END_REF]. Using the method of characteristics we find the solution of ( 50) and (51) for k < t where we neglect the initial data:

α(t, x, k) = α(t -k, x -vk, 0)e -k 0 p(k * ) dk * , (61) β(t, x, k)
= β(t -k, x + vk, 0)e -k 0 p(k * ) dk * . (62)
Next, from (59) let us define the escape and arrival rates of individuals at position x at time t as

j α (t, x) = t 0 p(k)α(t, x, k) dk , j β (t, x) = t 0 p(k)β(t, x, k) dk . (63)
Recalling the definitions (1) and ( 3) and following the steps in Appendix A-I we write

j α (t, x) = t 0 φ(t -s)e -v(t-s)∂x α(s, x, 0) ds + α 0 (x -vk)φ(k) , (64)
j β (t, x) = t 0 φ(t -s)e v(t-s)∂x β(s, x, 0) ds + β 0 (x + vk)φ(k) . (65)
Using the Laplace transform L[f ](t) = f (λ) where λ is the Laplace variable, we have

ĵα (λ, x) = φ(λ + v∂ x )α(λ, x, 0) + α 0 φ(λ + v∂ x ) . (66)
Moreover, using the Laplace transform of the characteristic solution (61) and the definition of ᾱ(t, x) we write

α(λ, x) = α(λ, x, 0) ψ(λ + v∂ x ) + α 0 ψ(λ + v∂ x ) . (67)
Substituting α(λ, x, 0) from (67) into (66) we finally get

ĵα (λ, x) = φ(λ + v∂ x ) ψ(λ + v∂ x ) α(λ, x) = Q(λ + v∂ x ) α(λ, x) . (68)
The operator Q(λ + v∂ x ) can be explicitly computed in the Laplace space. Transforming back to the (t, x)-space we have, for j α and j β ,

j α (t, x) = t 0 Q(t -s)ᾱ(s, x -v(t -s)) ds , (69)
j β (t, x) = t 0 Q(t -s) β(s, x + v(t -s)) ds . (70)
Using the expressions (69) and (70) we write the system (59)-(60) in term of the macroscopic quantities ᾱ and β. In the following we obtain explicit expressions for j β and j α by using the distribution of persistence steps k given in (2).

Left and right persistent movement

Using the results from the previous section we write the kinetic system as follows

∂ t ᾱ + v∂ x ᾱ = -j α + j β , (71)
∂ t β -v∂ x β = -j β + j α , (72)
where j α and j β are given by ( 69) and (70), respectively. The quantities in the right hand side of ( 71) and ( 72) are best expressed in the Fourier-Laplace space, where the Fourier-Laplace transform is defined as

FL[f ](t, x) = f (λ, ξ) = R ∞ 0 e iξx-λt f (t, x) dt dx .
Transforming the system (71)-(72) we write

λ α + α 0 + viξ α = -jα + jβ , λ β + β 0 -viξ β = -jβ + jα , (73)
where jα = Q(λ + viξ) α(λ, ξ) and jβ = Q(λ -viξ) β(λ, ξ). To obtain Q(λ ± viξ) we first have to compute the quantities φ(λ ± viξ) and ψ(λ ± viξ), previously defined in (2) and (3). Letting λ ± = λ ± viξ, φ± = φ(λ ± ) and ψ± = ψ(λ ± ) we write,

ψ± = k µ 0 λ µ+1 ± e k0λ± Γ(-µ + 1, k 0 λ ± ) , φ± = µ(k 0 λ ± ) µ Γ(-µ, k 0 λ ± )e k0λ± .
Using an asymptotic expansion of the Gamma function [2] and following the steps in [START_REF] Estrada-Rodriguez | Fractional Patlak-Keller-Segel equations for chemotactic superdiffusion[END_REF] we get

ψ± = - k 0 1 -µ - k 2 0 λ ± (1 -µ)(2 -µ) + k µ 0 λ µ-1 ± Γ(-µ + 1) + O(k 3 0 λ 2 ± ) , φ± = 1 + k 0 λ ± 1 -µ + k µ 0 λ µ ± + O(k µ+1 0 λ µ+1 ± ) . (74)
Using (74) we can write

Q(λ ± viξ) µ -1 k 0 - λ ± viξ 2 -µ -k µ-2 0 (λ ± viξ) µ-1 (µ -1)Γ(-µ + 1) . (75)
Hence, system (73) is now written in the (t, x)-space, for b = k µ-1 0 (µ -1)Γ(-µ + 1),

∂ t ᾱ + v∂ x ᾱ = - µ -1 k 0 (ᾱ -β) + ∂ t + v∂ x 2 -µ ᾱ - ∂ t -v∂ x 2 -µ β + b (∂ t + v∂ x ) µ-1 ᾱ -(∂ t -v∂ x ) µ-1 β , ∂ t β -v∂ x β = - µ -1 k 0 ( β -ᾱ) + ∂ t -v∂ x 2 -µ β - ∂ t + v∂ x 2 -µ ᾱ + b (∂ t -v∂ x ) µ-1 β -(∂ t + v∂ x ) µ-1 ᾱ .
Here we have used the fact that

FL ∂ t ± v∂ x µ-1 f = (λ ± viξ) µ-1 f .
Remark 5.1. The tempered fractional material derivative [START_REF] Fedotov | Persistent random walk of cells involving anomalous effects and random death[END_REF][START_REF] Sokolov | Towards deterministic equations for Lévy walks: The fractional material derivative[END_REF], defined as

∂ t ± v∂ x µ-1 f (t, x) = 0 D µ-1 t f (t, x ± vt) ,
generalises the standard material derivative d dt f (t, x ± vt) = (∂ t ± v∂ x )f for µ = 2 through the introduction of the Riemann-Liouville operator [START_REF] Sokolov | Towards deterministic equations for Lévy walks: The fractional material derivative[END_REF].

Macroscopic PDE for the total population

We may now write a macroscopic equation for the total density ρ(t, x) = ᾱ(t, x) + β(t, x). From the definitions of j α and j β in (63) we know that j α (t, x) = β(t, x, 0) and j β (t, x) = α(t, x, 0) ,

and therefore, using ( 61) and (62) we have

α(t, x, k) = j β (t -k, x -vk)ψ(k) , β(t, x, k) = j α (t -k, x + vk)ψ(k) . (77)
Hence, from (77) we can write

j α (t, x) = t 0 p(k)α(t, x, k) dk = t 0 φ(k)j β (t -k, x -vk) dk + α 0 (x -vk)φ(k) , (78)
j β (t, x) = t 0 p(k)β(t, x, k) dk = t 0 φ(k)j α (t -k, x + vk) dk + β 0 (x + vk)φ(k) . (79)
On the other hand we have

ᾱ(t, x) = t 0 α(t, x, k) dk = t 0 j β (t -k, x -vk)ψ(k) dk + α 0 (x -vk)ψ(k) , (80)
β(t, x) = t 0 β(t, x, k) dk = t 0 j α (t -k, x + vk)ψ(k) dk + β 0 (x + vk)ψ(k) . (81)
Re-arranging the above expressions and using the notation introduced in Section 5.2 for ψ ± and φ ± we get

α(λ, ξ) = φ- ψ- β(λ, ξ) + α0 ψ+ , β(λ, ξ) = φ+ ψ+ α(λ, ξ) + β0 ψ-, (84)
or equivalently,

α(λ, ξ) = φ-β0 ψ+ + α0 ψ+ 1 -φ+ φ- , β(λ, ξ) = φ+ α0 ψ-+ β0 ψ- 1 -φ-φ+ . ( 85
)
Note that if we substitute the jα from the expression for β in (83) into jβ = φjα + β0 φwe obtain the relation (68) in Section 5.

Fractional scaling We consider the following scaling

(t n , k n , x n ) → (t/ε θ , k/ε κ , x/ε ν ) , (86)
where θ, κ, ν > 0. We introduce the scaling in the expressions (2) and (3)

ψ ε (k) = ε κ k 0 ε κ k 0 + k µ , φ ε (k) = µ(ε κ k 0 ) µ (ε κ k 0 + k) µ+1 , p ε (k) = µε κ ε κ k 0 + k , (87)
and from now on, we take a = ε κ k 0 . Now consider the case when the cell starts to move to the right at t = 0 from the point x = 0, then α 0 = ε v δ(x) where v > 0 is a constant and β 0 (x) = 0. Since ρ = α + β we have, in the Fourier-Laplace space, using (84),

ψ+ ε ψ- ε (ρ -ε v ) = ψ+ ε φ-β ψ+ ε + ψ- ε φ+ ε α ψ- ε . (88)
Using the expansions (74) and following the steps in Appendix A-II the above expression can be written as

ε v + aλ - β 1 -µ + α 2 -µ + aλ + β 2 -µ + α 1 -µ = a µ-1 Γ(-µ + 1)(1 -µ)(λ µ-1 + β + λ µ-1 -α) . ( 89
)
Replacing λ ± = λ ± ivξ and including the scaling we have

ε v + k 0 (ε θ+κ λ -ε ν+κ ivξ) β 1 -µ + α 2 -µ + k 0 (ε θ+κ λ + ε ν+κ ivξ) β 2 -µ + α 1 -µ = k µ-1 0 ε (κ+ν)(µ-1) Γ(-µ + 1)(1 -µ) (ivξ) µ-1 β + (-ivξ) µ-1 α . ( 90
)
Note that on the right hand side we have used a quasi-static approximation (ε θ λ ± ε ν ivξ) µ-1 ε ν(µ-1) (±ivξ) µ-1 , assuming θ > ν. Grouping terms and using the definitions for the macroscopic density and the local flux ρ = α + β , Jε = α- β ε respectively, we obtain

ε θ+κ λρ -c µ ε v + 1 (3 -2µ) ε ν+κ+1 ivξ Jε = d µ ε (κ+ν)(µ-1) (ivξ) µ-1 β + (-ivξ) µ-1 α , (91)
where

c µ = - (1 -µ)(2 -µ) k 0 (3 -2µ) > 0 , d µ = k µ-2 0 (1 -µ) 2 (2 -µ) (3 -2µ) > 0 for 1 < µ < 3/2 . (92)
Choosing θ = κ(µ -2) + ν(µ -1), v = (κ + ν)(µ -1) and noting that ν + κ + 1 > (κ + ν)(µ -1) (which means that the normal diffusion is of lower order) for µ < 2 + 1 κ+ν we get

λρ -c µ = d µ (ivξ) µ-1 β + (-ivξ) µ-1 α . ( 93
)
Using FL[∂ t ρ(t, x)] = λρ(λ, ξ) -ρ 0 (0) , where we assume that c µ = ρ 0 (0) and the following relations for s ∈ (0, 1) [START_REF] Ferrari | Some nonlocal operators in the first Heisenberg group[END_REF][START_REF] Ferrari | Weyl and Marchaud derivatives: A forgotten history[END_REF]]

F[D s -f ] = (iξ) s f , F[D s + f ] = (-iξ) s f , we have ∂ t ρ = d µ (D µ-1 - β + D µ-1 + ᾱ) for µ ∈ (1, 3/2) .
Here D s -and D s + are Riemann-Liouville fractional derivatives defined as

D s -f = -1 Γ(1 -s) ∂ ∂x ∞ x f (y) (y -x) s dy , D s + f = 1 Γ(1 -s) ∂ ∂x x -∞ f (y) (x -y) s dy .
Following the steps in Appendix A-III we finally write the macroscopic equation as

∂ t ρ(t, x) = C - d 2 dx 2 µ-1 2 ρ(t, x) , (94)
where

C = d µ µ-1 2Γ(2-µ) 1 c(1, µ-1 2 ) > 0 .

Numerical results

We present some numerical results for the discrete synchronised system where we show the diffusive and superdiffusive regimes, in agreement with the results in Sections 4 and 5.

We start with the discrete description of the synchronised movement, which leads, in the limit, to (50)-(51).

Discrete description of the fully synchronised cell movement

The system (50)-(51) can be derived from a point particle when the probability of moving depends on previous steps taken in the same direction. We only treat the full synchronisation case, this derivation can be extended to the non-synchronised system.

As before, we denote by α(N, x, k) the probability that the cell moves to the right. Here N is the total number of steps, k are the number of steps given by the cell in the same direction, and x is the position. Analogously, we denote by β(N, x, k) the probability that the cell moves to the left. We recall that the probability of changing the direction is denoted by q k = τ p k where τ is a small time step. Therefore the probability of keep moving in the same direction is qk = 1 -τ p k . Since the cell has "memory" of the direction of the previous steps, we assume that the probability of changing direction decreases with the number of steps k according to a power-law. This models the directional persistence observed in experiments in [START_REF] Huda | Lévy-like movement patterns of metastatic cancer cells revealed in microfabricated systems and implicated in vivo[END_REF].

At each time a particle makes a step to the left or to the right according to its status, and then decides to keep moving in the same direction or reverse direction.

Discrete jumping

We first consider a cell moving to the right, after N steps, where it gave k steps in this direction. In the previous step N -1, this cell had done k -1 steps to the right and thus the probability of keep moving to the right is

α(N, x, k) = (1 -τ p k )α(N -1, x -δ, k -1) . (95)
We also have to consider the events when the cell was moving to the left at step N -1, described by β(N -1, x, k) and reverses direction with probability τ p k . Since the particle changed direction, it is set at k = 0 moving to the right, and thus we have

α(N, x, 0) = τ N -1 k=1 p k β(N -1, x -δ, k) . (96)
From (95) we can write, after dividing by τ α

(N, x, k) -α(N -1, x -δ, k -1) τ = -p k α(N -1, x-δ, k -1) . (97)
In the limit, for τ, δ → 0 and v = δ/τ we get,

∂ t α(t, x, k) + ∂ k α(t, x, k) + v∂ x α(t, x, k) = -p(k)α(t, x, k) , α(t, x, 0) = ∞ 0 p(k)β(t, x, k) dk . (98)
The second relation in (50) is obtained from (96), in the limit. Following the same steps for the left movement of the particle we start from

β(N, x, k) = (1 -τ p k )β(N -1, x + δ, k -1) , β(N, x, 0) = τ N -1 k=1 p k α(N -1, x+δ, k) , (99)
and in the limit we obtain

∂ t β(t, x, k) + ∂ k β(t, x, k) -v∂ x β(t, x, k) = -p(k)β(t, x, k) , β(t, x, 0) = ∞ 0 p(k)α(t, x, k) dk . (100)

Numerical set up and main numerical results

We consider a discrete velocity jump model which describes the left and right movement as in Section 6.1, in an infinite one dimensional domain. We assume that the speed of the cell is constant given by v = ±1 and the probability of changing direction from left to right is governed by (2).

To decide whether the cell changes direction or not, we use the rejection method. We randomly generate a number between (0, 1), if that number is bigger than a probability P = ψ(k) ψ(k-1)

1 , then the cell changes direction, otherwise it keeps moving without changing. The steps k are updated in each iteration and therefore P , where we always start with k = 1. The cell updates its position according to x(t i+1 ) = x(t i ) + v. This same description can be extended for the non-synchronisation case, where the movement of the front (y) and the back (x) are independent. Every time the cell changes direction we count the number of steps k given in the same direction. For the non-synchronisation case we take into account the biologically relevant switching probabilities given in Section 2.1 to preserve the realistic cell length.

With this toy example we are able to compute the mean square displacement (MSD) x 2 of the cells. As stated in the Introduction, normal diffusion processes are characterised by x 2 ∼ t, while for the case of superdiffusion x 2 ∼ t ζ for ζ ∈ (1, 2), where ζ = µ/2.

In Figure 3 we have the average of the MSD where this average is taken over 10 000 runs and the trajectories of the cell follows the discrete velocity jump process described before. As obtained in (92), the superdiffusion movement for the synchronised case is observed when µ ∈ (1, 3/2) which agrees with the results in Figure 3a. On the other hand, we consider the normal diffusion limit of the synchronised system derived in Section 4.1 where we observed normal diffusion for µ ∈ (2, 3). From Figure 3b, we see that the slope of the MSD is approximately 1, corresponding to the normal diffusion case.

Moreover, these findings are in agreement with [START_REF] Huda | Lévy-like movement patterns of metastatic cancer cells revealed in microfabricated systems and implicated in vivo[END_REF], where the authors observed superdiffusion for Lévy exponents µ = 1.39, 1.58, 1.50 and normal diffusion for µ = 2.17, 2.36, 3.57 (see Table 1 in [START_REF] Huda | Lévy-like movement patterns of metastatic cancer cells revealed in microfabricated systems and implicated in vivo[END_REF]).

Finally, for completeness we also present the numerical results for the non-synchronised case in Figure 4. Here we observe a similar behaviour as for the synchronised with the difference that now the superdiffusion is "weaker" in the sense that even for very small values of µ the slope of the MSD is close to one.

Conclusion and perspectives

We developed a formalism allowing to take into account how eukariotic cells move by protrusions (front of the cell) and retractions (back of the cell), keeping the simplicity of one space dimension for motion. Full generality, assuming that back and front are independent leads to a mathematical model hardly amenable to analysis, but various synchronisation levels lead to simpler models for which macroscopic effects can be observed. Among them we found normal drift-diffusion but more interestingly, instability can occur and, in the fully synchronised case, fractional diffusion characterised by long jumps. This is in accordance with experimental observations in [START_REF] Huda | Lévy-like movement patterns of metastatic cancer cells revealed in microfabricated systems and implicated in vivo[END_REF] where the trajectories of metastatic cells, which move in a synchronised way, followed a power-law distribution, characteristic of a superdiffusion process.

1 ψ(k) is the probability of a run of length at least k. I would like to achieve this distribution by independent decisions whether to turn or not (based on rand(1)). The probability to continue the run after the first time step is P (1), the probability to continue after the second time step is P (2), etc. The probability that the cell has not turned within the first k time steps is P (1)P (2)...P (k). This is in fact equal to ψ(k). The formula ψ(k) = P (1)P (2)...P (k) ∀k has a unique solution for the probabilities P : P (j) = ψ(j)/ψ(j -1). From a modelling and analytical point of view, several questions are left open. For instance, a better understanding of the full model and of possible model reduction. Also, the introduction of more biological details, for example, in the switching direction probability (2). We could tailor this function to a specific system by knowing the internal mechanisms that leads to synchronisation in cells. Moreover, we could extend our model to several dimensions and connect it to models of cell polarisation such as [START_REF] Loy | Kinetic models with non-local sensing determining cell polarization and speed according to independent cues[END_REF][START_REF] Loy | Modelling physical limits of migration by a kinetic model with non-local sensing[END_REF]. Finally, it would be interesting to look at the effect of the interactions with the environment and collective effects.

A Miscellaneous

(I) We compute the escape and arrival rates introduced in (63) by using the characteristic solutions (61) and (62). We start from j β (t, x) = t 0 p(k)β(t, x, k) dx which, by using (3) and (62), can be re-written as j β (t, x) = t 0 φ(k)β(t -k, x + vk, 0) dk + β0 (x + vk)φ(k) = t 0 φ(t -s)e v(t-s)∂x β(s, x, 0) ds + β0 (x + vk)φ(k) .

(101) The last equality is obtained using the change of variables k = t -s along with the following Taylor expansion

e v(t-s)∂x f (x) = ∞ m=0 (v(t -s)∂ x ) m m! f (x) = ∞ m=0 1 m! (v(t -s)) m ∂ m x f (x) = f (x + v(t -s)) .
Analogously we can obtain (64) for j α (t, x).

(II) Now we aim to derive the expression (88). From (84) we write, after multiplying both sides by ψ+ ψ-ψ+ ψρ = ψ+ φ-β + α0 ψ+ ψ-+ ψ-φ+ α ψ-+ β0 ψ-ψ+ .

Using the initial conditions α 0 ε = ε z δ(x) and β 0 (x) = 0 we obtain (88). Now, we introduce the scaling to (74) and we write

ψ± = - a 1 -µ - a 2 λ ± (1 -µ)(2 -µ)
+ a α λ µ-1 ± Γ(-µ + 1) + O(a 3 λ 2 ± ) , φ± = 1 + aλ ± 1 -µ + a µ λ µ ± + O(λ µ+1 ) .

Hence from here we compute

ψ+ φ-= - a 1 -µ - a 2 λ + (1 -µ)(2 -µ)
+ a µ λ µ-1 + Γ(-µ + 1) -

a 2 λ - (1 -µ) 2 + O(a µ+1 ) , ψ-φ+ = - a 1 -µ - a 2 λ - (1 -µ)(2 -µ)
+ a µ λ µ-1 -Γ(-µ + 1) -

a 2 λ + (1 -µ) 2 + O(a µ+1 ) , ψ+ ψ-= a 2 (1 -µ) 2 -a µ+1 Γ(µ + 1) 1 -µ (λ µ-1 - + λ µ-1 + ) + O(a 3 ) .
Substituting these three quantities in (102) we arrive at (89).

Figure 1 :

 1 Figure 1: Schematics of the front-back cell movement.

Figure 2 :

 2 Figure 2: Non-synchronised movement where the cell moves to the right (a), to the left (b) or changes its length by stretching (c) and shrinking (d).

  Next we apply the Fourier-Laplace transform to (78)-(81) to obtainjα (λ, ξ) = φ(λ + ivξ) jβ (λ, ξ) + α0 (ξ) , jβ (λ, ξ) = φ(λ -ivξ) jα (λ, ξ) + β0 (ξ) ,(82)α(λ, ξ) = jβ (λ, ξ) + α0 (ξ) ψ(λ + ivξ) , β(λ, ξ) = jα (λ, ξ) + β0 (ξ) ψ(λ -ivξ) .

Figure 3 :

 3 Figure 3: Average of the MSD taken from 10000 individual trajectories when we let them run for t = 10000. (a) describes the superdiffusive regime while (b) gives the normal diffusion case.

Figure 4 :

 4 Figure 4: MSD for the non-synchronised case. The blue and yellow lines are for µ = 1.3 and 1.5 respectively, and the orange and green are for the cases µ = 2.3 and 2.5.

1 B.P. has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 740623). G.E.R. acknowledges the support of the Fondation Sciences Mathématiques de Paris (FSMP) for the postdoctoral fellowship.

(III) Finally, we are going to work only with the fractional operators. Following [START_REF] Ferrari | Some nonlocal operators in the first Heisenberg group[END_REF][START_REF] Ferrari | Weyl and Marchaud derivatives: A forgotten history[END_REF] we have

The above relation is true if ᾱ, β ∈ C 1 (R) and ᾱ, β = o(|x| µ-2-), x → +∞ for > 0 (equivalence between Marchaud derivative and Riemann-Liuoville derivative). Now we are going to use the fact that the sum D µ-1 -f + D µ-1 + f gives the fractional Laplace operator in one dimension, also known as the Riesz derivative,

where c(1, µ-1 2 ) is a normalization constant.