Jérémy Jover

Sarra Aissani

Laouès Guendouz

André Thomas

Daniel Canet
email: daniel.canet@univ-lorraine.fr

NQR DETECTION OF SODIUM NITRITE RECRYSTALLIZED IN WOOD

Cylindrical pieces of wood which can be placed in 10mm o.d. NMR sample tubes have been impregnated by an aqueous solution of sodium nitrite (NaNO 2 ). They were subsequently dried and examined by 14 N quadrupole resonance. To our surprise, we were able to observe a signal at the right frequency (4.64 MHz, the highest frequency NaNO 2 line). This implies that the material has been properly embedded but also that sodium nitrite was properly recrystallized. This latter point was more than speculative. Anyway, we could observe that it works for spruce, beech, ash, maritime pine, but not for oak (which is known to be difficult to impregnate). Instrumental parameters have been optimized for reducing the duration of the experiment. If the first measurements required around four hours, we obtain at the present time acceptable results in about four minutes or less. In addition, we have designed a simple NLLS (Non Linear Least Squares) algorithm by which the spectral parameters of the NQR signal can be retrieved even if the peak is hardly visible in the frequency domain.

Extensive measurements have been performed on NaNO 2 recrystallized in spruce. They show, among other things, that the totality of NaNO 2 which has penetrated the material has actually recrystallized and seemingly for ever. Various assays have been carried out as a a function of impregnation and drying conditions. The evolution of the width of the NaNO 2 line reflects defects in the crystal lattice of NaNO 2 due probably to alteration of the wood structure.

Introduction

This work represents an attempt for detecting by NQR a molecule which has been imbedded in wood, the far-off objective being material labelling. In fact, our objective is rather to demonstrate the feasibility of this approach. Moreover, because we are familiar with 14 N quadrupole resonance, we have focused on this quadrupolar nucleus in spite of its low sensitivity. Indeed, we are facing several challenges: how to insert a substance within the material? If wood is impregnated, what is the amount of the substance which has recrystallized (if any)? Will it be possible to use NQR as a sensor? This paper is devoted to these questions. For simplicity, we have chosen a molecule often used as a model in Nitrogen-14 Quadrupole Resonance, namely sodium nitrite NaNO 2 . It has the further advantage of being soluble in water and can be considered as bio-compatible in certain circumstances (food preservation).

Experimental

Impregnation

Small cylindrical pieces of wood of 70 mm height and 8 mm diameter were machined so that they could fit into 10 mm o.d. tubes (employed for NMR measurements). They were impregnated with a saturated aqueous solution of sodium nitrite with the help of a specially constructed device (figure 1). Impregnation can be possibly performed under high pressure (up to 15 bar) [START_REF] Malkov | Towards complete impregnation of wood chips with aqueous solutions[END_REF] , but also under partial vacuum. We mean by "partial vacuum" soft pumping of the vessel containing the wood rod immersed in the aqueous solution of sodium nitrite. This has the virtue of removing air thus facilitating the impregnation process. In the case of impregnation "under high pressure", pumping is first applied to the wood rod (again for removing air) and then the rod is immersed in the pressurised aqueous solution. Only raw material (without any preparation prior to impregnation) led to significant results (see below). open (ambient temperature and moisture) or in a dry kiln at 103°C (highest limit before wood degradation 2 ) for a stabilization period (controlled by weighing at regular intervals) which generally did not exceed 48 hours. and at ambient temperature with a home-made spectrometer described previously [START_REF] Hiblot | A fully home-made N-14 quadrupole resonance spectrometer[END_REF] . For the most recent measurements, a new probe was constructed: still with a solenoidal coil surrounding the sample but with a cylindrical shield at a distance of 8 cm from the solenoid center. This prevents electromagnetic reflections thus parasitic capacitors which exist in the previous probes. Also, tuning and matching of this probe involve equilibrated circuits. All measurements reported here have been performed with a simple low power pulse of 100 μs duration (corresponding to the first maximum of the nutation curve). A dead time (between the end of the pulse and the beginning of acquisition) of 300 μs along with low power pulses avoids baseline distortions. A repetition delay of 0.5 s was allowed between two consecutive acquisitions. Typically, with this new probe, 4096 transients were accumulated but, as shown later, much less accumulation can be envisioned.

Results

It was almost a surprise to observe a NQR signal in most impregnated wood samples. This implies that, within this material, sodium nitrite has recrystallized upon drying. Moreover, relatively sharp lines indicate few (if any) defects in the lattice. Such experiments proved to be successful for almost all species that we have tried. Spectra are reported in figure 2. The first obvious observation concerns the ability of the different species to be impregnated by aqueous solutions of sodium nitrite. This can be quantified by the peak area in either of the two series, thus by the amount of recrystallized NaNO 2 . This impregnation ability is in the order spruce>Sylvester pine>beech>ash>oak (for the latter, we just observe the spike at zero frequency and it is indeed well known that oak is difficult to impregnate 4 ). Two other properties concern the drying mode: less NaNO 2 is recrystallized in the case "out in the open" while line-widths are smaller than in the case "dry kiln". The interpretation of the first feature is immediate: upon heating at a higher temperature, more water is eliminated. The second feature can be explained by the fact that taking the material at a relatively elevated temperature modifies its structure entailing some defects in the crystal lattice of NaNO It can be seen that the amount of recrystallized NaNO 2 is hardly larger when the sample has been subjected to the dry kiln. This was already recognized in figure 2 that the line-width is even weaker than in the case of impregnation under partial vacuum. This means that the interactions of recrystallized NaNO 2 with the material are weakened due to the impregnation mode. Conversely, as all water has been eliminated by the drying process in the kiln (see above), and owing to the large line-width (right spectrum of figure 4), we can infer that the spruce sample was impregnated to the core and that high pressure has contributed to alter the material and, as a consequence, to enhance the interactions undergone by sodium nitrite. Indeed, alteration of the material is especially deleterious when, prior to impregnation, the piece of wood is dried in the kiln. This procedure was attempted, again with the hope to increase the amount of recystallized sodium nitrite [START_REF] Matsumura | Impregnation of radiate pine wood by vacuum treatment: effect of pre-steaming on wood structure and resin content[END_REF] . It turns out that, in these conditions, we totally failed to observe any NQR signal (with impregnation under partial vacuum or under high pressure). This could mean that water, which is naturally present in wood, is necessary for keeping a structure prone to recrystallization of embedded sodium nitrite.

Finally, we have checked that, when sodium nitrite recrystallizes in wood, its NQR spectrum remains the same for months and probably longer.

Data analysis

All results, in terms of peak area (A) or line-width at half-height (LW), which have been given

until now have been obtained through the procedure described below.

As seen above, obtaining a spectrum with a suitable signal-to-noise ratio (S/N) requires the accumulation of a large number of transients. Let us recall that the spectrum (frequency domain) results from the Fourier transform of accumulated free induction decays (fids; time domain). As a consequence, when the NQR line is relatively broad, the S/N is evidently lowered. Hence, a better strategy is to look at the time domain data since the parameter of interest is the amplitude A of the fid (which is also the peak area in the frequency domain) while the fid damping factor is related to the line-width (LW) in the frequency domain. It is also important to realize that, in our case, the spectrum involves a single line, this point facilitating any fitting procedure. However, for future applications, this measuring time might still be too long. With one minute measuring time, we can guess the presence of a peak but this is no longer true for a measuring time of 30 seconds. Now, we are going to show that a simple non linear least squares (NLLS) fitting procedure applied to the time domain data (fid) is able not only to reveal the presence of a peak but also to provide reliably the parameters of interest, namely A and LW. Such least-squares procedures are relatively common in multiline NMR spectra [START_REF] Barkhuijsen | retrieval of frequencies, amplitudes, damping factors and phases from time domain signals using a linear least-squares procedure[END_REF][START_REF] Montigny | Analysis of time domain NMR data by standard non linear least-squares[END_REF] and make use of the following mathematical model for the fid (in the case of a single line and assuming that the fid is devoid

of any dc component) ] ) ( exp[ ) exp( ) 2 exp( ) ( 0 t LW i t i A t S π ϕ πν - = (

Conclusion

When this project started, it was considered as very challenging and rather utopian. The first result was looked at with some mistrust because it required extensive time averaging (several hours) and involved some artefacts. In fact, with some tenacity, several instrumental improvements and with appropriate data processing, we were able to obtain very consistent and reproducible results in a very reasonable measuring time (lower limit: 30 seconds).

Extensive studies of spruce impregnated by sodium nitrite as a function of impregnation and drying conditions have shown that heating always modify the wood structure, as reflected by the NaNO 2 signal line-width. This latter property is indicative of interactions between the NaNO 2 crystal lattice and the material. In particular, it was noticed that drying wood at elevated temperature prior to impregnation totally precludes the formation of NaNO 2 crystals.

Otherwise, it was shown that in soft conditions (without going to elevated temperatures), a large amount of NaNO 2 could be recrystallized in spruce wood. This opens the way of material labelling by nitrogen containing compounds and its further detection by nitrogen-14 Quadrupole Resonance.

Figure 1 .

 1 Figure 1. Block diagram of the impregnation device.

SpruceBeechFigure 2 . 4 .

 24 Figure 2. 4.64 MHz spectra of recrystallized NaNO 2 in different wood species (impregnation under partial vacuum) with the two drying modes described in the text. Left column: out in the open; right column: dry kiln. Continuous curves (in red) correspond to recalculated spectra (see below). A stands for the peak area normalized to 4096 scans whereas LW is the line-width at half-height.

  2 5 (due to possible interactions of NaNO 2 with the material). Notice that these early experiments have Ash 65536 scans. A = 8534. LW = 2823 Hz 65536 scans. A = 14202. LW = 2109 Hz Oak 131072 scans been performed with the old probe. In fact, the above trends have been confirmed by accurate measurements carried out with the new probe on spruce samples. Relevant spectra are shown in figure 3.

Figure 3 .

 3 Figure 3. 4.64 MHz spectra of NaNO 2 obtained with the new probe. Top: reference spectrum of a NaNO 2 powder sample of the same size as the spruce samples (bottom). Impregnation was carried out under partial vacuum. Left: out in the open; right: dry kiln. Continuous curves (in red) correspond to recalculated spectra (see below). A stands for the peak area normalized to 4096 scans whereas LW is the line-width at half height.

  but seems to be specific of spruce (for other species, less NaNO 2 is recrystallised out in the open). This amount, as determined by peak areas (the quantities A), is five times larger in the powder than in wood. From the weight of NaNO 2 in the powder sample, we arrive at the absolute weight of NaNO 2 in the spruce sample (650 mg). It turns out that this value coincides almost perfectly with weighing data (before impregnation and after drying). The next issue is evidently the line-width of the recrystallized NaNO 2 NQR signal. In the case of spruce, there is an order of magnitude when going from crystallized powder to NaNO 2 recrystallized in wood (with 4096 scans. A = 456477. LW = 1232 Hz 4096 scans. A = 461949. LW = 2101 Hz 512 scans. A = 2319450. LW = 478 Hz drying "out in the open").As already mentioned, this can be reasonably attributed to defects in the crystalline lattice which arise from interactions with the material. These defects are still more important (line-width increases by a factor of two) when drying is carried out in a kiln.As assumed above, this is very likely due to structural modifications within the material and thus to modifications of the interactions with NaNO 2 .

Figure 4 .

 4 Figure 4. 4.64 MHz spectra of NaNO 2 in spruce. Impregnation was carried out under a pressure of 15 bars. Left: drying out in the open; right: drying in the kiln. Continuous curves (in red) correspond to recalculated spectra (see below). A stands for the peak area normalized to 4096 scans whereas LW is the line-width at half height.

Figure 5 .

 5 Figure 5. 4.64 MHz spectra of NaNO 2 in spruce (impregnation under partial vacuum; drying out in the open. Left column: raw spectra. Right column; the same raw spectra with the reconstructed spectra (continuous curves in red) with parameters derived from time domain data fitting.

1 )

 1 The parameters A and LW have been defined earlier; 0 ν is the resonance frequency with respect to the carrier frequency and ϕ the phase factor. Alternatively the damping function by a Gaussian function. A dedicated software has been developed under MATLAB. It requires starting values for these four parameters and provides refined values along with reconstructed spectrum in the frequency domain. The strating values can be derived from the experimental spectrum if it is sufficiently legible. Otherwise, as the spectral parameters of the considered NQR signal can be relatively easily guessed, they can be used manually as input parameters. As shown in the right column of figure 5, reconstructed spectra ere perfectly consistent with spectra obtained from Fourier transform of the experimental fids. It was even possible to detect, thanks to the fitting of time domain data, a signal which is invisible in the frequency domain (bottom spectra of figure 5). The A values are, as expected, approximately proportional to the number of scans while LW values remain of the same order of magnitude.


Vincent Calvez
email: vincent.calvez@umpa.ens-lyon.fr

Nicolas Meunier
email: nicolas.meunier@parisdescartes.fr

Nicolas Muller
email: nicolas.muller@parisdescartes.fr

Raphael Voituriez
email: voiturie@lptmc.jussieu.fr

Numerical simulation of the dynamics of molecular markers involved in cell polarisation

Keywords: Cell polarisation, global existence, blow-up, numerical simulations, Keller-Segel system

, in this work, we study mathematical models for cell polarisation. These models are based on nonlinear convection-diffusion equations. The nonlinearity in the transport term expresses the positive loop between the level of protein concentration localised in a small area of the cell membrane and the number of new proteins that will be convected to the same area. We perform numerical simulations and we illustrate that these models are rich enough to describe the apparition of a polarisome.

Introduction

Cell polarisation is a major step involved in several important cellular processes such as directional migration, growth, oriented secretion, cell division, mating or morphogenesis. When a cell is not polarised molecular markers (proteins CDC42) are uniformly distributed on the membrane while polarisation is characterized by the concentration of molecular markers in a small 1 area of the cell membrane. In [START_REF] Wedlich-Soldner | Spontaneous cell polarization through actomyosin-based delivery of the cdc42 gtpase[END_REF], it has been observed that if the external pheromone concentration is above a critical concentration, polarisation can occur spontaneously. It has also been observed that cell asymmetry can be driven by an external asymmetric stimulus.

Cell polarisation in yeast cells has been intensively studied during the past decade. Recently, many models describing cell polarisation have been developed. The majority of these models are based on reaction-diffusion systems where polarisation appears as a type of Turing instability [START_REF] Iglesias | Navigating through models of chemotaxis[END_REF], [START_REF] Onsum | A mathematical model for neutrophil gradient sensing and polarisation[END_REF], [START_REF] Levine | Directional sensing in eukaryotic chemotaxis: A balanced inactivation model[END_REF], or due to stochastic fluctuations [START_REF] Altschuler | On the spontaneous emergence of cell polarity[END_REF], other models include cytoskeleton proteins as a regulatory factor [START_REF] Eugenio | Principles for the dynamic maintenance of cortical polarity[END_REF], [START_REF] Wedlich-Soldner | Spontaneous cell polarization through actomyosin-based delivery of the cdc42 gtpase[END_REF]. Many biological studies have shown that the cytoskeleton plays an important role in polarisation. It has been suggested that the cytoskeleton has a positive feedback on molecular markers density. Indeed, disruption of transport along the cytoskeleton greatly reduces the stability of polar cap [START_REF] Wedlich-Soldner | Spontaneous cell polarization through actomyosin-based delivery of the cdc42 gtpase[END_REF]. The cell cytoskeleton is a network of long semi-flexible filaments made up of protein subunits [START_REF] Phillips | Physical biology of the cell[END_REF]. These filaments (mainly actin or microtubules) act as roads along which motor proteins are able to perform a biased ballistic motion and carry various molecules. Molecular markers play a key role in the formation of these filaments.

Following [START_REF] Hawkins | Rebuilding cytoskeleton roads: active transport induced polarisation of cells[END_REF], [START_REF] Calvez | A one-dimensional Keller-Segel equation with a drift issued from the boundary[END_REF] and [START_REF] Calvez | Analysis of a nonlocal model for spontaneous cell polarization[END_REF], in this work we study models that describe the dynamics of cell polarisation. In these models, molecular markers, such as proteins, diffuse in the cytoplasm and are actively transported along the cytoskeleton. The resulting motion is a biased diffusion regulated by the markers themselves. Using numerical simulations and mathematical heuristics, we observe that the coupling on the velocity field achieves an inhomogeneous distribution of molecular markers without any external asymmetric field. Such an inhomogeneous distribution is only due to interaction between molecular markers.

Throughout this paper, the density of molecular markers (resp. advection field) is denoted by ρ(t, x) (resp. u(t, x)). The advection is obtained though a coupling with the membrane concentration of markers. The cell is figured by the domain Ω ⊂ R n with n = 1, 2 and a part of the boundary of the domain will be the active membrane denoted by Γ. The time evolution of the molecular markers satisfies the following advection-diffusion equation, see [START_REF] Hawkins | Rebuilding cytoskeleton roads: active transport induced polarisation of cells[END_REF] and [START_REF] Calvez | Analysis of a nonlocal model for spontaneous cell polarization[END_REF]:

∂ t ρ(t, x) = D ∆ρ(t, x) -χ ∇. (ρ(t, x) u(t, x)) , t > 0, x ∈ Ω, ρ(0, x) = ρ 0 (x). ( 1
)
There is no creation nor degradation of molecular markers in the cell, so the quantity of molecular markers remains constant in time:

M = x∈Ω ρ 0 (x)dx = x∈Ω ρ(t, x)dx. ( 2
)
This condition is ensured by a zero flux boundary condition on the boundary.

A first simplified step is to assume that the cell is essentially bidimensional and to neglect curvature effects. The membrane boundary is then a 1D line along the y-axis and the cytoplasm is parametrized by

x = (x, y) ∈ R + × R.
The plan of this work is the following. First, we recall the main mathematical results of the simplified model in 1D for Ω = (0, ∞) and Γ = {x = 0}, see [START_REF] Calvez | A one-dimensional Keller-Segel equation with a drift issued from the boundary[END_REF], [START_REF] Calvez | Analysis of a nonlocal model for spontaneous cell polarization[END_REF] for more details. Then we study a more realistic model, that includes dynamical exchange of markers on the boundary for a general Ω. This model was introduced in [START_REF] Hawkins | Rebuilding cytoskeleton roads: active transport induced polarisation of cells[END_REF] and studied in [START_REF] Calvez | Analysis of a nonlocal model for spontaneous cell polarization[END_REF] in the one dimensional case. Here, we will perform a first numerical analysis of this model in the two dimensional case, for periodic (in one direction) and bounded (in the other direction) domain. Finally, we provide a methodology for parameter estimation by using mathematical heuristics and biological literature.

One dimensional case

In this section, we study the one dimensional case on the half line for Ω = (0, ∞). The membrane is then the point Γ = {x = 0}. For the first model, the advection field towards the membrane is equal to the density of molecular markers on the boundary ρ(t, 0). Then we improve this model by considering that only the trapped molecular markers on the membrane contribute to the advection field.

Simplified model set on the half line

In [START_REF] Calvez | Analysis of a nonlocal model for spontaneous cell polarization[END_REF] a first mathematical studies has been done on this model. We define an advection field u(t, x) for (1) u(t, x) = -ρ(t, 0), in such a case (1) reads as (with D = 1 and χ = 1):

∂ t ρ(t, x) = ∂ xx ρ(t, x) + ρ(t, 0) ∂ x ρ(t, x), t > 0, x > 0, (3)
with the following zero flux condition on the boundary {x = 0}, that ensures the mass conversation (2),

∂ x ρ(t, 0) + ρ(t, 0) 2 = 0. (4)
In [START_REF] Calvez | Analysis of a nonlocal model for spontaneous cell polarization[END_REF], it has been proved that solutions of (3) blow-up in finite time if their masses are above a certain critical mass, M > 1, and exist globally in time if M ≤ 1. Let us first recall the definition of weak solutions of (3).

Definition 1.1 We say that ρ(t, x) is a weak solution of (3) on (0, T ) if it satisfies: andρ(t, x) is a solution of (3) in the sense of distributions in D (R + ).

ρ ∈ L ∞ (0, T ; L 1 + (R + )) , ∂ x ρ ∈ L 1 ((0, T ) × R + ) ,
Let us now recall the main results for weak solutions of (3). Theorem 1.2 (Global existence: M ≤ 1) Assume that the initial data ρ 0 satisfies both ρ 0 ∈ L 1 ((1 + x) dx) and x>0 ρ 0 (x)(log ρ 0 (x)) + dx < +∞. Assume in addition that M ≤ 1, then there exists a global weak solution of equation (3).

Theorem 1.3 (Blow-up: M > 1) Assume M > 1.
Any weak solution of equation (3) with non-increasing initial data ρ 0 blows-up in finite time.

Remark 1 It would tempted to interpret blow-up of solutions of the one dimensional model as cell polarisation. But it is to be noticed that concentration of markers on the boundary doesn't mean polarisation. Indeed, consider a radially symmetric 2D cell case. Equation then reduces to the one dimensional one. Above a threshold on the total mass, the convection wins and markers concentrate on the boundary. In some situations, these markers may be homogeneously distributed on the boundary and in such a case there is no symmetry breaking.

The model with dynamical exchange of markers at the boundary

Such a direct activation of transport on the boundary seems to be unrealistic. Indeed possible occurrence of blow-up in finite time suggests this claim. We improve the previous model by distinguishing between cytoplasmic content ρ(t, x) and the concentration of trapped molecules on the boundary, that will be denoted by µ(t). The dynamical exchange of markers at the boundary is done with an attachment rate k on and a detachment rate k of f , hence the time evolution of µ(t) is

d dt µ(t) = k on ρ(t, 0) -k of f µ(t). (5)
The advection field u(t, x) in ( 1) is now defined by

u(t, x) = -µ(t), hence (1)
(with D = 1 and χ = 1) reads as:

∂ t ρ(t, x) = ∂ xx ρ(t, x) + µ(t) ∂ x ρ(t, x), t > 0, x > 0, (6)
with a modified boundary condition

∂ x ρ(t, 0) + ρ(t, 0) µ(t) = d dt µ(t). (7)
This ensures the following mass conservation shared among ρ(t, x) and µ(t):

M = R + ρ 0 (x)dx + µ 0 = R + ρ(t, x)dx + µ(t). (8)
With equation ( 5), the self-activation of transport by ρ(t, 0) is then delayed in time. Since the transport speed is bounded µ(t) ≤ M , the solution of the model with dynamical exchange on the boundary exists globally in time.

More precisely it is possible, see [START_REF] Calvez | Analysis of a nonlocal model for spontaneous cell polarization[END_REF], to prove that it converges towards a non trivial stationary state.

Two dimensional case : the model with dynamical exchange of markers at the boundary

Let Ω ⊂ R 2 be the cytoplasm domain, as in the one dimensional case (5) we consider dynamical exchange of markers at the boundary, so for x ∈ Γ we have the evolution in time of µ(t, x)

∂ t µ(t, x) = k on ρ(t, x) -k of f µ(t, x). (9)
with a modified boundary condition for ρ(t, x)

at point x ∈ Γ (D ∇ρ(t, x) -χ ρ(t, x) u(t, x)). n x = -∂ t µ(t, x). ( 10
)
where n x is the outward normal to Γ. This ensures the following mass conservation sharing by ρ(t, x) and µ(t, x):

M = Ω ρ 0 (x)dx + Γ µ 0 (x)dx = Ω ρ(t, x)dx + Γ µ(t, x)dx. ( 11
)
We consider the following advection field deriving from a harmonic potential modeling the transport by actin filaments (cytoskeleton):

u(t, x) = ∇c(t, x), where -∆c(t, x) = 0, if x ∈ Ω, ∇c(t, x). n x = S(x)µ(t, x), if x ∈ Γ. ( 12
)
This advection field orientation is due to the actin networks.

Actin filaments are attached on the membrane and randomly distributed, there orientations are mixed up. We also add the external pheromone concentration at x ∈ Γ which acts by the mating-pheromone MAPK cascade on the actin transport.

In dimension 2, we have global existence for the model without exchange on the boundary (replacing equation [START_REF] Onsum | A mathematical model for neutrophil gradient sensing and polarisation[END_REF] by ∇c(t, x). n x = S(x)ρ(t, x) if x ∈ Γ) with Ω = (0, +∞) × R and Γ = {0} × R. For clarity, we recall this result, see [START_REF] Calvez | Analysis of a nonlocal model for spontaneous cell polarization[END_REF] for more details.

Theorem 1.4 (Global existence in dimension 2) Assume that the advection field satisfies the two following conditions: ∇ • u ≥ 0 and u(t, 0, y) • e e = ρ(t, 0, y). Assume that the initial data ρ 0 satisfies both ρ 0 ∈ L 1 ((1 + |x| 2 ) dx) and ρ 0 L 2 is smaller than some constant c. Then there exists a global weak solution to equations (1)-( 2).

In the two dimensional case, for the model with exchange on the boundary, blow-up or global existence have not been proved yet. In this work, we make a first step in this direction by using a mathematical heuristic and numerical simulations.

Heuristics

The mathematical analysis performed in [START_REF] Calvez | Analysis of a nonlocal model for spontaneous cell polarization[END_REF] has demonstrated that a class of models exhibit pattern formation (either blow-up or convergence towards a non homogeneous steady state) under some conditions. However the main question still remains unanswered: do these models describe cell polarisation or not? Thus in order to provide a first answer to this question, we will perform numerical simulations. Our aim is to see if, under some conditions, the model leads to a concentration of markers, not only on the boundary, but on a small region of the boundary. In such a case polarisation occurs. In order to obtain more information on the critical value distinguishing the polarised case and the stable case, in the two dimensional case we will use a mathematical heuristics that we describe now. Let x = (x, y) be in Ω = R + × R and let Γ = {0} × R be the boundary, we have that

u(t, x) = ∇c(t, x), where -∆c(t, x) = 0, if x ∈ R + × R, -∂ x c(t, 0, y) = S(y)µ(t, y), if y ∈ R, (13)
hence, see [START_REF] Evans | Partial differential equations[END_REF] e.g., it is well known that c(x, y) = -1 π y ∈R log( (y -y ) 2 + x 2 )(Sµ)(y )dy .

The tangential component at the boundary is then given by u(t, 0, y) • e y = -H(Sµ)(y) , y ∈ R,

where H denotes the one-dimensional Hilbert transform that we recall now, see [START_REF] Calvez | Modified keller-segel system and critical mass for the log interaction kernel[END_REF] e.g., with respect to the y variable:

H(f µ)(y) = 1 π p.v. R 1 y -y f (y ) dy .
Integrating the main equation ( 1) with respect to x with zero flux condition on Γ = {x = 0}, we obtain:

∂ t x>0 ρ(t, x, y) dx = D ∂ yy x>0 ρ(t, x, y) dx -χ ∂ y x>0 ρ(t, x, y)(u(t, x, y) • e y )dx .
In the super-critical case, numerical simulations, see [START_REF] Muller | Mathematical and numerical studies of nonlinear and nonlocal models involved in biology[END_REF], suggest that the density ρ(t, x) concentrates on the boundary {x = 0}. Assuming ρ(t, x, y) = ν(t, y)δ(x = 0), we can formally write the dynamics of ν(t, y) as follows:

∂ t ν(t, y) = D ∂ yy ν(t, y) + χ ∂ y (ν(t, y)H(Sµ)(y)) .
Assuming S constant on R and µ(t, y) = kon k of f ν(t, y) for y ∈ R, it reads as

∂ t ν(t, y) = D ∂ yy ν(t, y) + χS k on k of f ∂ y (ν(t, y)H(ν)(y)) .
Hilbert transform has a critical singularity to offset the diffusion on this equation [START_REF] Calvez | Modified keller-segel system and critical mass for the log interaction kernel[END_REF]. We have a blow-up if R ν(t, y) dy = M is above

2πDk of f
Sχkon . This is a first step to observe a critical mass phenomenon and this may lead to blow-up if the mass is large enough. In this way, we define an order of magnitude for some parameters.

It is to be noticed that this latter criterion is valid for an infinite domain, namely y ∈ R. In the case of a cell, the domain will be finite and the existence of such a dichotomy has not been proved yet. In order to see if such a dichotomy holds true we will perform numerical simulations. This is the object of the following section.

Numerical analysis

We first give a discretization of the convection-diffusion model set on a 1D periodic domain. This first step allows us introducing the discretization of this model on a 2D domain which is periodic in one direction and bounded on the other direction.

One dimensional case

Let u(t, x) be a given function. We consider the following advection-diffusion equation on the periodic domain Ω = R/Z

∂ t ρ(t, x) = ∂ x (∂ x ρ(t, x) -u(t, x) ρ(t, x)), t > 0, x ∈ Ω. ( 14
)
Let t n = n dt be the time discretization and {x j = j dx, j ∈ {1, ..., N x }} be the space discretization of the periodic interval R/Z. Since the equations of the model are written in a conservative form, the natural framework to be used for the spatial discretization is the finite volume framework. We hence introduce the control volume defined for j ∈ {1, ..., N x }

V j = (x j-1 2 , x j+ 1 2 ). ( 15
) Let ρ n j (resp. u n j+ 1
2

) be the approximated value of the exact solution ρ(t n , x j ) (resp. u(t n , x j+ 1 2

)), the classical upwind scheme for ( 14) reads as

ρ n+1 j -ρ n j dt = F j+ 1 2 -F j-1 2 dx , j ∈ {1, ..., N x }, (16)
where the numerical flux F j+ 1 2 and F j-1 2 are defined by

F j+ 1 2 = ρ n+1 j+1 -ρ n+1 j dx -A up (u n j+ 1 2 , ρ n j , ρ n j+1 ), F j-1 2 = ρ n+1 j -ρ n+1 j-1 dx -A up (u n j-1 2 , ρ n j-1 , ρ n j ),
with the advection numerical flux is given by

A up (u, x -, x + ) = u x -, si u > 0, u x + , si u < 0. ( 17
)
The periodic flux condition on boundary reads as

F 1 2 = F Nx+ 1 2
and we set the value u n

1 2 = u n Nx+ 1 2
. The diffusion part is treated implicitly and it is then unconditionally stable, while the advection term is treated explicitly. The CFL condition of the scheme is

u n j+ 1 2 j∈{1,...,Nx} ∞ < dx dt .
We define the column vector

ρ n = ρ n 1 ρ n 2 . . . ρ n Nx T .
As usual, see e.g.

[1], the discrete heat matrix A ∈ M Nx (R) with periodic flux condition on the boundary is defined as

A =          2 + dx 2 dt -1 -1 -1 2 + dx 2 dt . . . . . . . . . . . . . . . 2 + dx 2 dt -1 -1 -1 2 + dx 2 dt          . ( 18
)
Periodic flux condition adds the top right term and the bottom left term. Next, in order to use A up defined by equation ( 17), we define (u) + = max(u, 0) and (u) -= min(u, 0). The discrete advection matrix B ∈ M Nx (R) with periodic flux condition on the boundary is then defined as in [START_REF] Allaire | Numerical analysis and optimization. An introduction to mathematical modelling and numerical simulation[END_REF]

B = dx 2 dt I Nx -dx             u n 3 2 + u n 3 2 - . . . . . . u n j+ 1 2 + (u n j+ 1 2 - . . . u n Nx-1 2 - u n Nx+ 1 2 - u n Nx+ 1 2 +             (19) +dx 
           u n 1 2 - u n 1 2 + u n 3 2 + . . . u n j-1 2 + u n j-1 2 - . . . . . . u n Nx-1 2 + u n Nx-1 2 -            
. We use a standard numerical method to invert the symmetric positive definite matrix A. Finally, at each time step we resolve

ρ n+1 = A -1 B ρ n .

Two dimensional case

We perform numerical simulations on the model with dynamical exchange of markers at the boundary. In this work, we assume that the cell occupies a disk of radius r > 0. Furthermore for simplicity, we consider a boundedperiodic domain Ω = [0, r]×R/2πrZ with Γ = {r}×R/2πrZ. This simplifies our numerical approach by using finite difference schemes on Cartesian grid. We start with the numerical study of the equation on ρ by assuming that the advection field u(t, x) = ∇c(t, x) is known. Then we perform the discretization of c. In this section, for simplicity we fix all the parameters values to 1 except M . Let us first recall the model with dynamical exchange of markers at the boundary on Ω = [0, r] × R/2πrZ:

∂ t ρ = ∇. (∇ρ -ρ ∇c) , in (0, r) × R/2πrZ, (20)
∂ x ρ -ρ ∂ x c = -∂ t µ, on {r} × R/2πrZ, (21)
∂ x ρ -ρ ∂ x c = 0, on {0} × R/2πrZ. (22)
Dynamical exchange markers on active boundary {r} × R/2πrZ is given by

∂ t µ = ρ -µ, on {r} × R/2πrZ. (23)
Laplace equation on c with non appropriate Neumann conditions on a bounded domain is ill-posed, see [START_REF] Allaire | Numerical analysis and optimization. An introduction to mathematical modelling and numerical simulation[END_REF] e.g. In order to handle this problem, we add a degradation term:

-∆c + α c = 0, in (0, r) × R/2πrZ, (24)
-∂ x c = µ, on {r} × R/2πrZ, (25)
-∂ x c = 0, on {0} × R/2πrZ. (26)
We take random initial conditions c, µ 0 and ρ 0 satisfying the following mass conservation

Ω ρ 0 + Γ µ 0 = M. ( 27
)
Let t n = n dt be the time discretization and {x j = j dx, j ∈ {1, ..., N x }} (resp. {y k = k dy, k ∈ {1, ..., N y }}) be the space discretization of the bounded interval [0, r) (resp. periodic interval R/2πrZ). We introduce the control volume

W (j,k) ⊂ R 2 W (j,k) = (x j-1 2 , x j+ 1 2 ) × (y k-1 2 , y k+ 1 2
).

Let P n (j,k) (resp. µ n k ) be the approximated value of the exact solution ρ(t n , x j , y k ) of equations ( 20),(21),( 22) and (27) (resp. µ(t n , y k ) of equations ( 23) and ( 27)). Let c (j,k) be the approximated value of the exact solution c(x j , y k ) of equations ( 24),(25) and (26).

Equation on µ

We can resolve at each time step for k ∈ {1, ..., N y }

µ n+1 k = µ n k + dt (ρ n k -µ n k ).

Equation on c

For simplicity, we call F the numerical flux as in the 1D case, we can write the following scheme: for (j, k) ∈ {1, ..., N x } × {1, ..., N y }

F (j+ 1 2 ,k) -F (j-1 2 ,k) dx + F (j,k+ 1 2 ) -F (j,k-1 2 ) dy -αc (j,k) = 0.
with numerical flux defined by

F (j+ 1 2 ,k) = c (j+1,k) -c (j,k) dx , F (j-1 2 ,k) = c (j,k) -c (j-1,k) dx .
The zero flux boundary condition (26) impose that F ( 1 2 ,k) = 0 and the boundary condition (25) F (Nx+ 1 2 ,k) = -µ n k for k ∈ {1, ..., N y }. Similarly, the periodic conditions impose F (j,Ny+ 1 2 ) = F (j, 1 2 ) for j ∈ {1, ..., N x }. We define the column vector C by C(k + (j -1)N y ) = C (j,k) with (j, k) ∈ {1, ..., N x } × {1, ..., N y }. As previously the rigidity matrix A 2D,α is defined by

A 2D,α =              A α + Id -Id -Id A α + 2 Id -Id -Id A α + 2 Id . . . . . . . . . . . . . . . A α + 2 Id -Id -Id A α + 2 Id -Id -Id A α + Id             
, where discrete Poisson matrix A α ∈ M Ny (R) in 1D with periodic flux condition on boundary is defined by

A α =          2 + α dx 2 -1 -1 -1 2 + α dx 2 . . . . . . . . . . . . . . . 2 + α dx 2 -1 -1 -1 2 + α dx 2         
.

The flux boundary condition {r} × R/2πrZ imposes this right hand side column vector of length N x N y :

R c = -dx (µ n k ) k 0 . . . 0 .
We use a standard numerical method to invert the symmetric positive definite matrix A 2D,α and then resolve at each time step

C = A -1 2D,α R c .

Equation on ρ

For simplicity, we call F the numerical flux as in the previous cases, we can write the upwind scheme as follows:

P n+1 (j,k) -P n (j,k) dt = F (j+ 1 2 ,k) -F (j-1 2 ,k) dx + F (j,k+ 1 2 ) -F (j,k-1 2 )
dy , and periodic flux boundary condition in the y-axis direction is defined by

B 2D =        B B . . . B B        -dx           U + 3 2 U - 3 2
. . . . . .

U + j+ 1 2 U - j+ 1 2 . . . U - Nx-1 2 U + Nx+ 1 2           + dx           U - 1 2
U + We use a standard numerical method to inverse the symmetric positive definite matrix A 2D and then resolve at each time step

P n+1 = A -1 2D (B 2D P n + R ρ ).

Graphics

With the previous numerical analysis, we implement all numerical simulations using MATLAB. We test different values of M :

Conclusion

In this work we have provided a first answer to the following question: do the nonlinear convection-diffusion models given in [START_REF] Hawkins | Rebuilding cytoskeleton roads: active transport induced polarisation of cells[END_REF] and [START_REF] Calvez | Analysis of a nonlocal model for spontaneous cell polarization[END_REF] describe cell polarisation or not? To do so we have used both a mathematical heuristic and numerical simulations. Numerical simulations were necessary because the heuristic is only valid for an infinite geometry while the cell is obviously finite. The numerical simulations ensure that solutions develop symmetry breaking over a critical value M * given us a first justification of the mathematical heuristic. In a further work, we will estimate an approximate value of this critical mass.

.

  The flux boundary condition {r} × R/2πrZ imposes this right hand side column vector of length N x N y : R ρ = -dx (

where u n

and the numerical flux is defined by

The zero flux boundary conditions (22

2 ) for j ∈ {1, ..., N x }. We define the column vector P n by P n (k + (j -1)N y ) = P n (j,k) with (j, k) ∈ {1, ..., N x } × {1, ..., N y }. In what follows for simplicity we consider that dx = dy. We define the rigidity matrix A 2D ∈ M NxNy (R) with A ∈ M Ny (R) defined by (18):


Vincent Calvez
email: vincent.calvez@umpa.ens-lyon.fr

Nicolas Meunier
email: nicolas.meunier@parisdescartes.fr

Nicolas Muller
email: nicolas.muller@parisdescartes.fr

Raphael Voituriez
email: voiturie@lptmc.jussieu.fr

Numerical simulation on a cell polarisation model: the polar case

Keywords: Cell polarisation, global existence, blow-up, numerical simulations, Keller-Segel system

When it is polarised, a cell develops an asymmetric distribution of specific molecular markers, cytoskeleton and cell membrane shape. Polarisation can occur spontaneously or be triggered by external signals, like gradients of signalling molecules... In this work, we use the models of cell polarisation introduced in [9] and [4] and we set a numerical analysis for these models. They are based on nonlinear convection-diffusion equations and the nonlinearity in the transport term expresses the positive loop between the level of protein concentration localised in a small area of the cell membrane and the number of new proteins that will be convected to the same area. We perform numerical simulations and we illustrate that these models are rich enough to describe the apparition of a polarisome.

Introduction

Cell polarisation is a symmetry-breaking event that occurs in cell division, mating or morphogenesis. Molecular markers play a central role in establishing this phenomenon. Indeed, there are two different behaviours: a non-polarised cell has its markers radially homogeneously distributed while markers are located in a small area of the cell membrane for a polarised 1 cell. Yeast cells are dynamically polarised in response to the extracellular gradients of signals (chemokynes). However, it has been observed in [START_REF] Wedlich-Soldner | Spontaneous cell polarization through actomyosin-based delivery of the cdc42 gtpase[END_REF] that polarisation can occur spontaneously without any external asymmetric stimulus.

During the past decade, many models describing cell polarisation have been developed. The majority of these models are based on reaction-diffusion systems where polarisation appears as a type of Turing instability [START_REF] Iglesias | Navigating through models of chemotaxis[END_REF], [START_REF] Onsum | A mathematical model for neutrophil gradient sensing and polarisation[END_REF], [START_REF] Levine | Directional sensing in eukaryotic chemotaxis: A balanced inactivation model[END_REF], or due to stochastic fluctuations [START_REF] Altschuler | On the spontaneous emergence of cell polarity[END_REF], other models include cytoskeleton proteins as a regulatory factor [START_REF] Eugenio | Principles for the dynamic maintenance of cortical polarity[END_REF], [START_REF] Wedlich-Soldner | Spontaneous cell polarization through actomyosin-based delivery of the cdc42 gtpase[END_REF]. Many biological studies have shown that the cytoskeleton plays an important role in polarisation. It has been suggested that there is a positive feedback on molecular markers density. Indeed, disruption of transport along the cytoskeleton greatly reduces the stability of polar cap [START_REF] Wedlich-Soldner | Spontaneous cell polarization through actomyosin-based delivery of the cdc42 gtpase[END_REF]. The cell cytoskeleton is a network of long semiflexible filaments made up of protein subunits [START_REF] Phillips | Physical biology of the cell[END_REF]. These filaments (mainly actin or microtubules) act as roads along which motor proteins are able to perform a biased ballistic motion and carry various molecules. Molecular markers play a key role in the formation of these filaments.

Following [START_REF] Hawkins | Rebuilding cytoskeleton roads: active transport induced polarisation of cells[END_REF], [START_REF] Calvez | A one-dimensional Keller-Segel equation with a drift issued from the boundary[END_REF] and [START_REF] Calvez | Analysis of a nonlocal model for spontaneous cell polarization[END_REF], in this work we study models that describe the dynamics of cell polarisation. In these models, molecular markers, such as proteins, diffuse in the cytoplasm and are actively transported along the cytoskeleton. The resulting motion is a biased diffusion regulated by the markers themselves. Using numerical simulations and mathematical heuristics, we observe that the coupling on the velocity field achieves an inhomogeneous distribution of molecular markers without any external asymmetric field. Such an inhomogeneous distribution is only due to interaction between molecular markers.

Throughout this paper, the density of molecular markers (resp. advection field) is denoted by ρ(t, x) (resp. u(t, x)). The advection is obtained through a coupling with the membrane concentration of markers. The cell is figured by the domain Ω ⊂ R n with n = 1, 2 and a part of the boundary of the domain will be the active membrane denoted by Γ. The time evolution of the molecular markers satisfies the following advection-diffusion equation, see [START_REF] Hawkins | Rebuilding cytoskeleton roads: active transport induced polarisation of cells[END_REF] and [START_REF] Calvez | Analysis of a nonlocal model for spontaneous cell polarization[END_REF]:

∂ t ρ(t, x) = D ∆ρ(t, x) -χ ∇. (ρ(t, x) u(t, x)) , t > 0, x ∈ Ω, ρ(0, x) = ρ 0 (x). ( 1
)
There is no creation nor degradation of molecular markers in the cell, so the quantity of molecular markers remains constant in time:

M = x∈Ω ρ 0 (x)dx = x∈Ω ρ(t, x)dx. ( 2
)
This condition is ensured by a zero flux boundary condition on the boundary.

A first simplified step is to assume that the cell is essentially bidimensional and to neglect curvature effects. The membrane boundary is then a 1D line along the y-axis and the cytoplasm is parametrized by

x = (x, y) ∈ R + × R.
The plan of this work is the following. First, we present the models, that are based on different expressions for u, and we recall the main mathematical results of the simplified model in 1D for Ω = (0, ∞) and Γ = {x = 0}, see [START_REF] Calvez | A one-dimensional Keller-Segel equation with a drift issued from the boundary[END_REF], [START_REF] Calvez | Analysis of a nonlocal model for spontaneous cell polarization[END_REF] for more details. Then we study a more realistic model, that includes dynamical exchange of markers on the boundary. This model was introduced in [START_REF] Hawkins | Rebuilding cytoskeleton roads: active transport induced polarisation of cells[END_REF] and studied in [START_REF] Calvez | Analysis of a nonlocal model for spontaneous cell polarization[END_REF]. We provide a methodology for parameter estimation and qualitative description of cell polarisation by using mathematical heuristics. Then, we perform a numerical analysis of this model. We introduce the numerical part by the one dimensional case. Finally, we give tools to study the numerical implementation of the model on an annulus domain.

2 Presentation of the models and mathematical results

One dimensional case

In this section, we study the one dimensional case on the half line for Ω = (0, ∞). The membrane is then the point Γ = {x = 0}. For the first model, the advection field towards the membrane is equal to the density of molecular markers on the boundary ρ(t, 0). Then we improve this model by considering that only the trapped molecular markers on the membrane contribute to the advection field.

Simplified model set on the half line

In [START_REF] Calvez | Analysis of a nonlocal model for spontaneous cell polarization[END_REF] a first mathematical studies has been done on this model. We define an advection field u(t, x) for ( 1)

u(t, x) = -ρ(t, 0),
in such a case (1) reads as (with D = 1 and χ = 1):

∂ t ρ(t, x) = ∂ xx ρ(t, x) + ρ(t, 0) ∂ x ρ(t, x), t > 0, x > 0, (3)
with the following zero flux condition on the boundary Γ = {x = 0}, that ensures the mass conversation (2),

∂ x ρ(t, 0) + ρ(t, 0) 2 = 0. ( 4
)
Interestingly enough in [START_REF] Calvez | Analysis of a nonlocal model for spontaneous cell polarization[END_REF], it has been proved that solutions of (3) blow-up in finite time if their masses are above a certain critical mass, M > 1, and exist globally in time if M ≤ 1. Let us first recall the definition of weak solutions of (3).

Definition 2.1 We say that ρ(t, x) is a weak solution of (3) on (0, T ) if it satisfies: andρ(t, x) is a solution of (3) in the sense of distributions in D (R + ).

ρ ∈ L ∞ (0, T ; L 1 + (R + )) , ∂ x ρ ∈ L 1 ((0, T ) × R + ) ,
Let us now recall the main results for weak solutions of (3). Theorem 2.2 (Global existence: M ≤ 1) Assume that the initial data ρ 0 satisfies both ρ 0 ∈ L 1 ((1 + x) dx) and x>0 ρ 0 (x)(log ρ 0 (x)) + dx < +∞. Assume in addition that M ≤ 1, then there exists a global weak solution of equation (3). Remark 1 It would tempted to interpret blow-up of solutions of the one dimensional model as cell polarisation but concentration of markers on the boundary doesn't automatically mean polarisation. Indeed, consider a radially symmetric 2D cell case then equation (1) reduces to the one dimensional one. Above a threshold on the total mass, the convection wins and markers concentrate on the boundary. In some situations, these markers may be homogeneously distributed on the boundary and in such a case there is no symmetry breaking.

The model with dynamical exchange of markers at the boundary

Such a direct activation of transport on the boundary seems to be unrealistic. Indeed possible occurrence of blow-up in finite time suggests this claim. We improve the previous model by distinguishing between cytoplasmic content ρ(t, x) and the concentration of trapped molecules on the boundary, that will be denoted by µ(t). The dynamical exchange of markers at the boundary is done with an attachment rate k on and a detachment rate k of f , hence the time evolution of µ(t) is

d dt µ(t) = k on ρ(t, 0) -k of f µ(t). (5)
The advection field u(t, x) in ( 1) is now defined by

u(t, x) = -µ(t), hence (1)
(with D = 1 and χ = 1) reads as:

∂ t ρ(t, x) = ∂ xx ρ(t, x) + µ(t) ∂ x ρ(t, x), t > 0, x > 0, (6)
with a modified boundary condition

∂ x ρ(t, 0) + ρ(t, 0) µ(t) = d dt µ(t). (7)
This ensures the following mass conservation shared among ρ(t, x) and µ(t):

M = R + ρ 0 (x)dx + µ 0 = R + ρ(t, x)dx + µ(t). (8)
With equation ( 5), the self-activation of transport by ρ(t, 0) is then delayed in time. Since the transport speed is bounded µ(t) ≤ M , the solution of the model with dynamical exchange on the boundary exists globally in time.

More precisely it is possible, see [START_REF] Calvez | Analysis of a nonlocal model for spontaneous cell polarization[END_REF], to prove that it converges towards a non trivial stationary state.

Theorem 2.4 (Global existence: dynamical exchange case) Assume that the initial data ρ 0 satisfies both

ρ 0 ∈ L 1 ((1+x) dx) and x>0 ρ 0 (x)(log ρ 0 (x)) + dx < +∞. Assume the mass is super-critical M > 1.
The partial mass m(t) = x>0 ρ(t, x) dx converges to 1 and the density ρ(t, x) strongly converges in L 1 towards the exponential profile (M -1) exp(-(M -1)x).

Two dimensional case

In this section, we study the two dimensional case on Ω ⊂ R 2 . As in the one dimensional case, the advection field towards the membrane depends on the density of molecular markers on the boundary, two different situations (actin and microtubules) will be described, then we improve this model by considering an exchange of markers at the membrane and only the trapped molecular markers contribute to the advection field.

Simplified model set on the half plan

We study the model on the half plan for x = (x, y) ∈ Ω = R + × R. The membrane is then the line Γ = {x = 0}×R. We have the following boundary condition for

ρ(t, x) at point x ∈ Γ (D ∇ρ(t, x) -χ ρ(t, x) u(t, x)). n x = 0, (9)
where n x is the outward normal to Γ. This ensures the following mass conservation:

M = Ω ρ 0 (x)dx = Ω ρ(t, x)dx. ( 10
)
First we consider the transversal case, the field u is normal to the boundary

u(t, x) = -S(y)ρ(t, 0, y) e x . ( 11
)
The microtubules of the cytoskeleton are normally oriented to the cell membrane and their growth depends on the density of molecular markers on the boundary. In the potential case, we consider the following advection field deriving from a harmonic potential modelling the transport by actin filaments:

u(t, x) = ∇c(t, x), where -∆c(t, x) = 0, if x ∈ Ω, ∇c(t, x). n x = S(x)ρ(t, x), if x ∈ Γ. ( 12
)
This advection field orientation is due to the actin networks.

Actin filaments are attached on the membrane and randomly distributed, there orientations are mixed up. We also add the external pheromone concentration at x ∈ Γ which acts by the mating-pheromone MAPK cascade on the actin transport. We have global existence and blow-up theorems for the simplified model with x = (x, y) ∈ Ω = R + × R. For clarity, we recall this result, see [START_REF] Calvez | Analysis of a nonlocal model for spontaneous cell polarization[END_REF] for more details.

Theorem 2.5 (Global existence in dimension 2) Assume that the advection field satisfies the two following conditions: ∇ • u ≥ 0 and u(t, 0, y) • e x = -ρ(t, 0, y). Assume that the initial data ρ 0 satisfies both ρ 0 ∈ L 1 ((1 + |x| 2 ) dx) and ρ 0 L 2 is smaller than some constant c. Then there exists a global weak solution to equations (1)-(2).

Theorem 2.6 (Blow-up in dimension 2) Assume that ρ(t, x) is a strong solution to (1) which verifies:

• ∂ x ρ(t, x) ≤ 0 for all x ∈ Ω and t > 0 when the advective field is given by (11).

• ∂ x ρ(t, x) ≤ 0 for all x ∈ Ω and t > 0 the matrix A(t, x) = x ⊗ ∂ x ∂ y log ρ(t, x) satisfies A T + A ≥ 0 (in the matrix sense) when the advective field is given by [START_REF] Onsum | A mathematical model for neutrophil gradient sensing and polarisation[END_REF].

Assume in addition that the second momentum is initially small enough: there exists a constant C such that x∈Ω |x| 2 ρ 0 (x)dx < CM 3 . Then the maximal time of existence of the solution is finite.

The model with dynamical exchange of markers at the boundary

Let Ω ⊂ R 2 be the cytoplasm domain, as in the one dimensional case (5) we consider dynamical exchange of markers at the boundary, so for x ∈ Γ we have the evolution in time of µ(t, x)

∂ t µ(t, x) = k on ρ(t, x) -k of f µ(t, x), (13)
with a modified boundary condition for ρ(t, x)

at point x ∈ Γ (D ∇ρ(t, x) -χ ρ(t, x) u(t, x)). n x = -∂ t µ(t, x), (14)
where n x is the outward normal to Γ. This ensures the following mass conservation sharing by ρ(t, x) and µ(t, x):

M = Ω ρ 0 (x)dx + Γ µ 0 (x)dx = Ω ρ(t, x)dx + Γ µ(t, x)dx. ( 15
)
As before, we consider the following advection field (cytoskeleton):

u(t, x) = ∇c(t, x), where -∆c(t, x) = 0, if x ∈ Ω, ∇c(t, x). n x = S(x)µ(t, x), if x ∈ Γ. ( 16
)
We can also consider the transversal case. For example for Ω = R + × R, we take similarly as before

u(t, x) = -S(y)µ(t, y) e x . (17)
For the model with exchange on the boundary, blow-up or global existence have not been proved yet. In this work, we make a first step in this direction by using a mathematical heuristic and numerical simulations. There is also one open question: does advection field (transversal or potential) create a break of symmetry ?

Heuristic

The mathematical analysis performed in [START_REF] Calvez | Analysis of a nonlocal model for spontaneous cell polarization[END_REF] has demonstrated that a class of models exhibit pattern formation (either blow-up or convergence towards a non homogeneous steady state) under some conditions. However the main question still remains unanswered: do these models describe cell polarisation or not? Thus in order to provide a first answer to this question, in the next section, we will perform numerical simulations. Our aim will be to see if, under some conditions, the model leads to a concentration of markers, not only on the boundary, but on a small region of the boundary.

Let us now describe the mathematical heuristic which is done on the half plan Ω = R + × R. We give formal arguments to motivate the differences arising in the dynamics leading by the two possible drifts u T (transversal case) and u P (potential case). Let u P = ∇c be the solution of [START_REF] Onsum | A mathematical model for neutrophil gradient sensing and polarisation[END_REF] we have (see [START_REF] Evans | Partial differential equations[END_REF]) c(x, y) = -1 π y ∈R log( (y -y ) 2 + x 2 )S(y )ρ(t, 0, y )dy .

We notice that the two possible drifts u T and u P share common features: they are both divergence free and their normal components at the boundary coincide. u T • e x = u P • e x = -S(y)ρ(t, 0, y) .

On the other hand, a key difference holds when looking at the tangential component at the boundary:

u T • e y = 0 , u P • e y = -πH(S(•)ρ(t, 0, •)) ,
where H denotes the one-dimensional Hilbert transform:

H(f )(y) = 1 π p.v. R 1 y -x f (x) dx .
We expect the solution to concentrate on the boundary in the super-critical case for both choices of u, numerical simulations suggest it (see section 4.3).

Postulating the ansatz ρ(t, x, y) = ν(t, y)δ(x = 0), we can formally write the dynamics of ν(t, y) for the two cases. Integrating the main equation ( 1) with respect to x with zero flux condition on Γ = {x = 0}, we obtain:

∂ t R + ρ(t, x, y)dx = D∂ yy R + ρ(t, x, y)dx -χ∂ y R + ρ(t, x, y)(u(t, x, y) • e y )dx .
• u T • e y = 0: diffusion equation for ν:

∂ t ν(t, y) = D ∂ yy ν(t, y) .

• u P • e y = -πH(S(•)ρ(t, 0, •)) and assuming S constant on R, it reads as:

∂ t ν(t, y) = D ∂ yy ν(t, y) + χS ∂ y (ν(t, y)H(ν)(y)) .
Transversal dynamics are very different: in transversal case boundary diffusion dominates while in potential case Hilbert transform has a critical singularity to offset the diffusion on this equation. This latter equation exhibits blow-up if R ν(t, y) dy = M is above the critical mass 2πD Sχ , see [START_REF] Calvez | Modified keller-segel system and critical mass for the log interaction kernel[END_REF] done in a peculiar variant of Keller-Segel equation for more informations.

Remark 2 This is a first step to observe a critical mass phenomenon and this may lead to blow-up if the mass is large enough. In this way, we also define an order of magnitude for some parameters. It is to be noticed that this latter criterion is valid for an infinite domain, namely y ∈ R. In the case of a cell, the domain is finite and the existence of such a dichotomy has not been proved yet. In order to obtain more information on the critical value distinguishing the polarised case and the stable case we will perform numerical simulations.

This heuristic shows us that potential case is able to break symmetry more readily than the transversal case. This difference between the two models is also discussed in the following section.

Numerical analysis

We introduce this numerical section by the discretization of the convectiondiffusion model set on a 1D periodic domain. This first step allows us introducing the discretization of this model on a 2D domain. We model the cell as an annulus, molecular markers cannot pass inside the nucleus. In this section, for simplicity we fix all the parameters values to 1 except M .

First step with the one dimensional case

Let u(t, θ) be a periodic function on the periodic domain Ω = R/2πZ. We consider the following advection-diffusion equation

∂ t ρ(t, θ) = ∂ θ (∂ θ ρ(t, θ) -ρ(t, θ)u(t, θ)) , t > 0, θ ∈ R/2πZ. ( 18
)
Let t n = n ∆t be the time discretization and {θ k = k ∆θ, k ∈ {1, ..., N θ }} be the space discretization of the periodic interval R/2πZ. Since the equations of the model are written in a conservative form, the natural framework to be used for the spatial discretization is the finite volume framework. We hence introduce the control volume defined for k ∈ {1, ..., N θ }

V k = (θ k-1 2 , θ k+ 1 2 ). ( 19
) Let ρ n k (resp. u n k+ 1
2

) be the approximated value of the exact solution ρ(t n , θ k ) (resp. u(t n , θ k+ 1 2

)), the classical upwind scheme for (18) reads as

ρ n+1 k -ρ n k ∆t = F k+ 1 2 -F k-1 2 ∆θ , k ∈ {1, ..., N θ }, ( 20
)
where the numerical flux F k+ 1 2 and F k-1 2 are defined by

F k+ 1 2 = ρ n+1 k+1 -ρ n+1 k ∆θ -A up (u n+1 k+ 1 2 , ρ n+1 k , ρ n+1 k+1 ), F k-1 2 = ρ n+1 k -ρ n+1 k-1 ∆θ -A up (u n+1 k-1 2 , ρ n+1 k-1 , ρ n+1 k ),
with the advection numerical flux is given by

A up (u, x -, x + ) = u x -, if u > 0, u x + , if u < 0. ( 21
)
The periodic flux condition on boundary reads as

F 1 2 = F N θ + 1 2 = ρ n+1 1 -ρ n+1 N θ ∆θ -A up (u n+1 1 2 , ρ n+1 N θ , ρ n+1 1 )
and we recall that u is periodic so we set the value u n

1 2 = u n N θ + 1 2
. The diffusion and advection terms are both treated implicitly, the scheme is then unconditionally stable. We define the column vector

ρ n = ρ n 1 ρ n 2 . . . ρ n N θ T .
As usual, see e.g. [START_REF] Allaire | Numerical analysis and optimization. An introduction to mathematical modelling and numerical simulation[END_REF], the discrete heat matrix A ∈ M N θ (R) with periodic flux condition on the boundary is defined as

A =          2 -1 -1 -1 2 . . . . . . . . . . . . . . . 2 -1 -1 -1 2          . ( 22
B n =             u n 3 2 + u n 3 2 - . . . . . . u n j+ 1 2 + u n j+ 1 2 - . . . u n N θ -1 2 - u n N θ + 1 2 - u n N θ + 1 2 +             (23) -             u n 1 2 - u n 1 2 + u n 3 2 + . . . u n j-1 2 + u n j-1 2 - . . . . . . u n N θ -1 2 + u n N θ -1 2 -             .
At each time step we have

ρ n+1 -ρ n ∆t = - 1 ∆θ 2 Aρ n+1 - 1 ∆θ B n+1 ρ n+1 .
We use a standard numerical method to invert the matrix A + ∆θ B n+1 + ∆θ 2 ∆t I N θ . Finally, at each time step we resolve

ρ n+1 = A + ∆θ B n+1 + ∆θ 2 ∆t I N θ -1 ∆θ 2 ∆t ρ n .

Two dimensional case: polar case

Let Ω ⊂ R 2 be the cytoplasm domain, the domain is obviously an annulus where the molecular markers cannot enter in the nucleus of the cell. Let us first recall the model on the annulus Ω = B(0, R max ) \ B(0, R min ) with n x the unit normal vector to Ω at point x ∈ ∂Ω (we note C(0, R) the circle of center (0, 0) and radius R):

∂ t ρ(t, x) = ∇. (∇ρ(t, x) -ρ(t, x) u(t, x)) , in Ω, (24) (∇ρ(t, x) -ρ(t, x) u(t, x))
.

n x = 0, on C(0, R max ), (25) (∇ρ(t, x) -ρ(t, x) u(t, x))
.

n x = 0, on C(0, R min ). (26)
Since the domain Ω is assumed to be an annulus, it is appropriate to introduce polar coordinates r and θ. Let x = (r cos(θ), r sin(θ)) ∈ Ω, we have the following equations on 1 r ρ(t, r, θ) = ρ(t, x) with (r, θ)

∈ [R min , R max ] × R/2πZ: ∂ t ρ(t, r, θ) = ∂ r r∂ r ρ(t, r, θ) r -ρ(t, r, θ)u r (t, r, θ) + ∂ θ 1 r 2 (∂ θ ρ(t, r, θ) -ρ(t, r, θ)u θ (t, r, θ)) , in Ω, (27) 0
= r∂ r ρ(t, r, θ) r -ρ(t, r, θ)u r (t, r, θ), on C(0, R max ),(28) 0 = r∂ r ρ(t, r, θ) r -ρ(t, r, θ)u r (t, r, θ), on C(0, R min ).( 29
)
We recall the potential case u(t, x) = ∇c(t, x). Laplace equation on c with non appropriate Neumann conditions on a bounded domain is ill-posed, see [START_REF] Allaire | Numerical analysis and optimization. An introduction to mathematical modelling and numerical simulation[END_REF] e.g. In order to handle this problem, we add a degradation term:

-∆c(x) + α c(x) = 0, in Ω, ( 30
) ∇c(x). n x = ρ(t, x), on C(0, R max ), (31) ∇c(x)
.

n x = 0, on C(0, R min ). ( 32
)
We also give the previous equations in polar coordinates by

1 r c(r, θ) = c(x) with (r, θ) ∈ [R min , R max ] × R/2πZ: -∂ r r∂ r c(r,θ) r -1 r 2 ∂ θθ c(r, θ) + α c(r, θ) = 0, in Ω, (33)
∂ r c(r,θ) r = ρ(Rmax,θ) Rmax , on C(0, R max ), (34)
∂ r c(r,θ) r = 0, on C(0, R min ). (35)
If we consider dynamical exchange of markers at the active boundary, for x ∈ Γ = C(0, R max ) we have the evolution in time of µ(t, x).

∂ t µ(t, x) = ρ(t, x) -µ(t, x), on C(0, R max ). (36)
We replace then (28) by

-∂ t µ(t, θ) = r∂ r ρ(t, r, θ) r -ρ(t, r, θ)u r (t, r, θ), on C(0, R max ), (37) and (34) by
∂ r c(r, θ) r = µ(t, θ), on C(0, R max ). (38)
The transversal case is then

u(t, r, θ) = ρ(t, R max , θ) R max e r ,
and we can also write the transversal case for the dynamical exchange model u(t, r, θ) = µ(t, θ)e r , Let t n = n ∆t be the time discretization and {r j = R min + j ∆r, j ∈ {1, ..., N r }} (resp. {θ k = k ∆θ, k ∈ {1, ..., N θ }}) be the space discretization of the bounded interval [R min , R max ] (resp. periodic interval R/2πZ). We introduce the control volume

W (j,k) ⊂ R 2 W (j,k) = r j-1 2 , r j+ 1 2 × θ k-1 2 , θ k+ 1 2 .
Let P n (j,k) (resp. µ n k ) be the approximated value of the exact solution ρ(t n , r j , θ k ) (resp. µ(t n , θ k )). Let c(j,k) be the approximated value of the exact solution c(r j , θ k ).

Equation on µ

In the dynamical exchange model, we can resolve at each time step the discretization of equation (36) for k ∈ {1, ..., N y }

µ n+1 k = µ n k + ∆t (ρ n k -µ n k ).

Equation on c

For simplicity, we call F the numerical flux as in the 1D case, we can write the following scheme for equation (33) for (j, k) ∈ {1, ..., N r } × {1, ...,

N θ } - F (j+ 1 2 ,k) -F (j-1 2 ,k) ∆r + F (j,k+ 1 2 ) -F (j,k-1 2 ) ∆θ + α c(j,k) = 0.
In order to use finite volume we define the numerical flux by

F (j+ 1 2 ,k) = r j+ 1 2 c(j+1,k) r j+1 -c(j,k) r j ∆r , F (j-1 2 ,k) = r j-1 2 c(j,k) r j -c(j-1,k) r j-1
∆r ,

F (j,k+ 1 2 ) = 1 r 2 j c(j,k+1) -c(j,k) ∆θ , F (j,k-1 2 ) = 1 r 2 j c(j,k) -c(j,k-1) ∆θ .
The zero flux boundary condition (35) impose that F ( 1 2 ,k) = 0 and the boundary condition (34)

F (Nr+ 1 2 ,k) = r Nr+ 1 2
µ n k for k ∈ {1, ..., N θ }. Similarly, the periodic conditions impose for j ∈ {1, ..., N r }

F (j,N θ + 1 2 ) = F (j, 1 2 ) = 1 r 2 j c(j,1) -c(j,N θ ) ∆θ .
We

       +          1 r 2 1 A 1 r 2 2 A . . . 1 r 2 Nr -1 A 1 r 2 Nr A          . ( 39
)
The flux boundary condition C(0, R max ) imposes this right hand side column vector of length N r N θ :

R n c = r Nr+ 1 2       0 . . . 0 ρn Nr ,k r Nr k       or in the exchange case R n c = r Nr+ 1 2      0 . . . 0 (µ n k ) k      .
The zero flux boundary condition (29) impose that F ( 1 2 ,k) = 0. In the simplified model, the boundary condition (28) F (Nr+ 1 2 ,k) = 0 for k ∈ {1, ..., N θ } and in the model with exchange, we have

F (Nr+ 1 2 ,k) = - µ n+1 k -µ n k ∆t
for k ∈ {1, ..., N θ }. Similarly, the periodic conditions impose for j ∈ {1, ..., N r }

F (j,N θ + 1 2 ) = F (j, 1 2 ) = 1 r 2 j P n+1 (j,1) -P n+1 (j,N θ ) ∆θ -A up u n+1 (j, 1
2 ) , P n+1 (j,N θ ) , P n+1 (j,1)

.

We define the column vector P n by P n (k + (j -1)N θ ) = P n (j,k) with (j, k) ∈ {1, ..., N r } × {1, ..., N θ }:

P n = P n (1,1) . . . P n (1,N θ ) P n (2,1) . . . P n (2,N θ ) . . . P n (Nr,N θ ) T .
We define the following diagonal matrices for j ∈ {1, ..., N r }, U +

j+ 1 2 ∈ M N θ (R) and U - j+ 1 2 ∈ M N θ (R): U + j+ 1 2 =           . . . (u n (j+ 1 2 ,k-1) ) + (u n (j+ 1 2 ,k) ) + (u n (j+ 1 2 ,k+1) ) + . . .           , U - j+ 1 2 =           . . . (u n (j+ 1 2 ,k-1) ) - (u n (j+ 1 2 ,k) ) - (u n (j+ 1 2 ,k+1) ) - . . .           .
Thus we can define:

B n =           U + 3 2 U - 3 2
. . . . . .

U + j+ 1 2 U - j+ 1 2 P n+1 -P n ∆t = - 1 ∆r 2 AP n+1 - 1 ∆r B n+1 P n+1 .
In the exchange model, the flux boundary condition on C(0, R max ) imposes this right hand side column vector of length N r N θ :

R n exchange = -       0 . . . 0 ( µ n+1 k -µ n k ∆t ) k       .
We have at each time step

P n+1 -P n ∆t = - 1 ∆r 2 AP n+1 - 1 ∆r B n+1 P n+1 + 1 ∆r R n exchange .
We use a standard numerical method to invert the matrix A + ∆r B n+1 + ∆r 2 ∆t I NrN θ and then resolve at each time step

P n+1 = A + ∆r B n+1 + ∆r 2 ∆t I NrN θ -1 ∆r 2 ∆t P n + R n exchange .

Graphics

We use the numerical analysis done in this article to implement it with Matlab. Simulations have first been done in the transversal case u T . The following behaviour was obtained: The heuristic done in section 3 allowed us thinking that the potential case could break symmetry more readily than the transversal case. We have also done numerical simulation in potential case for small M and we found that the numerical behaviours are similar in the two possible drifts, see Fig 2 and4.

Figure 4: Numerical Simulations on Ω = [0.2, 2.5] × R/2πZ and all parameters equal to 1 with random initial conditions. Fo. For M = 0.01 small, steady state is isotropic.

We have assumed that the cell occupies a circle of radius R > 0. Furthermore for simplicity, we consider a bounded-periodic domain Ω = [0, R]×R/2πRZ.

With the numerical analysis done in [START_REF]Numerical simulation of the dynamics of molecular markers involved in cell polarisation[END_REF] for the potential case, we see that the behaviours are similar with the annulus case developed in this article.

Conclusion

Polar description of the cell has been described in this article, this improvement fitted the real cell shape and was a first step before establishing a model for membrane deformation. In this work we have provided a first answer to the following question: do the nonlinear convection-diffusion models given in [START_REF] Hawkins | Rebuilding cytoskeleton roads: active transport induced polarisation of cells[END_REF] and [START_REF] Calvez | Analysis of a nonlocal model for spontaneous cell polarization[END_REF] describe cell polarisation or not? To do so we have used both a mathematical heuristic and numerical simulations, which have ensured us that solutions develop symmetry breaking over a critical value M * . This has given us a first justification of the mathematical heuristic. On this point, the numerical behaviours are close to cell behaviours during biological experiences. In order to fit biological measurements, the choice of parameters is essential and we refer to biological literature. Several measurements on polarisation time and localisation of polar cap have been made, we will describe them in a further work.

Theorem 2 . 3 (

 23 Blow-up: M > 1) Assume M > 1. Any weak solution of equation (3) with non-increasing initial data ρ 0 blows-up in finite time.

)

  Periodic flux condition adds the top right term and the bottom left term. Next, in order to use A up defined by equation (21), we define (u) + = max(u, 0) and (u) -= min(u, 0) so A up (u, ρ k , ρ k+1 ) = (u) + ρ k + (u) -ρ k+1 . The discrete advection matrix B n ∈ M N θ (R) with periodic flux condition on the boundary is then defined as in[START_REF] Allaire | Numerical analysis and optimization. An introduction to mathematical modelling and numerical simulation[END_REF]

  define the column vector C by C(k + (j -1)N θ ) = c(j,k) with (j, k) ∈ {1, ..., N r } × {1, ..., N θ }: C = c(1,1) . . . c(1,N θ ) c(2,1) . . . c(2,N θ ) . . . c(Nr,N θ ) T For ∆r = ∆θ the rigidity matrix A is defined by

Figure 1 :

 1 Figure 1: Numerical Simulations on Ω = [0.2, 2.5] × R/2πZ and all parameters equal to 1 with random initial conditions. For M = 10 greater enough, no symmetry breaking appears, molecular markers are uniformly distributed on the membrane.

Figure 2 :Figure 3 :

 23 Figure 2: Numerical Simulations on Ω = [0.2, 2.5] × R/2πZ and all parameters equal to 1 with random initial conditions. Fo. For M = 0.01 small, steady state is isotropic.

Acknowledgement: this research has been supported by ANR program JCJC project MODPOL. N. Meunier and N. Muller want to thank M. de Buhan for helpful discussions.

We use a standard numerical method to invert the symmetric positive definite matrix 1 ∆r 2 A + αI NrN θ and then resolve at each time step

Equation on ρ

For simplicity, we call F the numerical flux as in the 1D case, we can write the following scheme for equation (27): for (j, k) ∈ {1, ..., N r } × {1, ..., N θ }

∆θ .

We define the numerical flux by

2 ) , P n+1 (j,k) , P n+1 (j,k+1)

,

2 ) , P n+1 (j,k-1) , P n+1 (j,k)

.

In the transversal case [START_REF] Levine | Directional sensing in eukaryotic chemotaxis: A balanced inactivation model[END_REF], we take u r = ρ(Rmax,θ) r and u θ = 0 we define at time

r Nr or in the exchange case u n

2 ) = 0,

In the potential case [START_REF] Onsum | A mathematical model for neutrophil gradient sensing and polarisation[END_REF], we have