Piotr Biler

Alexandre Boritchev

Lorenzo Brandolese

Sharp well-posedness and blowup results for parabolic systems of the Keller-Segel type

We study two toy models obtained after a slight modification of the nonlinearity of the usual doubly parabolic Keller-Segel system. For these toy models, both consisting of a system of two parabolic equations, we establish that for data which are, in a suitable sense, smaller than the diffusion parameter τ in the equation for the chemoattractant, we obtain global solutions, and for some data larger than τ , a finite time blowup. In this way, we check that our size condition for the global existence is sharp for large τ , up to a logarithmic factor.

Introduction and main results

This paper is concerned with parabolic systems (TM) and (TM') below, depending on a diffusion parameter τ > 0:

     u t = ∆u -u∆ϕ, τ ϕ t = ∆ϕ + u, u(0) = u 0 , ϕ(0) = ϕ 0 , x ∈ R d , t > 0, (TM) and      u t = ∆u + (∆ϕ) 2 , τ ϕ t = ∆ϕ + u, u(0) = u 0 , ϕ(0) = ϕ 0 , x ∈ R d , t > 0, (TM')
The above systems degenerate into the quadratic nonlinear heat equation (NLH) when τ = 0 and a compatibility condition, ∆ϕ 0 = -u 0 , is put on ϕ 0 :

u t = ∆u + u 2 , u(0) = u 0 , x ∈ R d , t > 0. (NLH)
These model systems are introduced in order to show the influence of the parameter τ in the second equation (a linear nonhomogeneous heat equation) on the size of admissible initial data leading to global-in-time solutions, and to finite time blowup, respectively. Our main motivation is to analyze those issues for the Keller-Segel system (PP), continuing the analysis started in [START_REF] Biler | Large global solutions of the parabolic-parabolic Keller-Segel system in higher dimensions[END_REF]. Here (PP) denotes the doubly parabolic Keller-Segel system describing chemotaxis, given by

     u t = ∆u -∇ • (u∇ϕ), τ ϕ t = ∆ϕ + u, u(0) = u 0 , ϕ(0) = ϕ 0 , x ∈ R d , t > 0, (PP)
where τ > 0, and (PE) is the parabolic-elliptic Keller-Segel system with τ = 0

     u t = ∆u -∇ • (u∇ϕ), ∆ϕ + u = 0, u(0) = u 0 , x ∈ R d , t > 0. (PE)
The existence of solutions for (PP) has been studied in the recent work [START_REF] Biler | Large global solutions of the parabolic-parabolic Keller-Segel system in higher dimensions[END_REF]. The issue of blowup is largely open for (PP), except for [START_REF] Winkler | Finite-time blow-up in the higher-dimensional parabolic-parabolic Keller-Segel system[END_REF][START_REF] Winkler | Blow-up profiles and life beyond blow-up in the fully parabolic Keller-Segel system[END_REF] where the radially symmetric problem is considered in a ball and in [START_REF] Winkler | Single-point blow-up in the Cauchy problem for the higher-dimensional Keller-Segel system[END_REF] in the whole space. Note that some concentration phenomena for (PP) have been shown in [START_REF] Calvez | Blow-up, concentration phenomenon and global existence for the Keller-Segel model in high dimension[END_REF]. There, a supplementary information on L 1 solutions is derived from entropy functionals and other specific properties of those drift-diffusion systems. Note that (TM) and (TM') have the same structure of steady states as (NLH), so they are related in a way similar to that as (PP) relates to (PE). However, (NLH) has a number of specific properties (such as variational structure, energy functional) that are not immediately extended to systems like (TM) and (TM').

All the considered systems feature cross-diffusion terms which make their analysis delicate. This justifies, in a sense, the application of Besov spaces in their analysis, cf. [START_REF] Iwabuchi | Global well-posedness for Keller-Segel system in Besov type spaces[END_REF][START_REF] Lemarié-Rieusset | Small data in an optimal Banach space for the parabolic-parabolic and parabolic-elliptic Keller-Segel equations in the whole space[END_REF] for well-and ill-posedness issues.

The analysis of the toy models (TM) and (TM') shed some light on the optimal spaces where it is reasonable to address the existence problem for (PP). Indeed, we prove in Sec. 2.2 existence for (TM) in the framework of Besov-type spaces, very close to optimal spaces with respect to the admissible initial data leading to local-in-time solutions. Next, in Sec 3.1, a couple of such existence and regularity results is shown for (TM), and for (TM') in Sec. 3.2, in the framework of pseudomeasures, similarly to the presentation in [START_REF] Biler | Large global solutions of the parabolic-parabolic Keller-Segel system in higher dimensions[END_REF]Sec. 2]. Section 4 is devoted to proofs of blowup, for (TM) in Sec. 4.1 and for (TM') in Sec. 4.2. The proofs are based on the idea going back to [START_REF] Montgomery-Smith | Finite time blow up for a Navier-Stokes like equations[END_REF] (and used for (PE) in [START_REF] Biler | Global existence versus blow up for some models of interacting particles[END_REF]) where the Fourier transform of solutions is analyzed: we prove that it tends to infinity in an appropriate sense in finite time, so that the solution itself necessarily becomes non-smooth. Thus, the technique used in the proof of Theorem 4.1 differs much from the usual approaches to nonlinear parabolic equations (see e.g. [START_REF] Quittner | Superlinear parabolic problems. Blow-up, global existence and steady states[END_REF]Ch. 17]). A more traditional approach involving moments of solutions is presented in Sec. 5 for the model (TM') considered in bounded domains.

Besides formal similarities of (TM), (TM') with (NLH) and (PE) when τ = 0, or with (PP) when τ > 0, there are rigorous results on singular perturbation limits τ 0 for small solutions of Keller-Segel systems for which blowup does not occur. Indeed, solutions with ϕ 0 = (-∆) -1 u 0 are recovered for τ = 0, in the sense of suitable convergences in [START_REF] Raczyński | Stability property of the two-dimensional Keller-Segel model[END_REF][START_REF] Biler | On the parabolic-elliptic limit of the doubly parabolic Keller-Segel system modelling chemotaxis[END_REF][START_REF] Lemarié-Rieusset | Small data in an optimal Banach space for the parabolic-parabolic and parabolic-elliptic Keller-Segel equations in the whole space[END_REF][START_REF] Kurokiba | Singular limit problem for the Keller-Segel system and drift-diffusion system in scaling critical spaces[END_REF] where various functional settings are proposed. Similar results hold also for the toy models but our main interest here is rather the behaviour for large τ than for τ close to 0.

Note however, that a natural conjecture that blowup phenomena should be continuous with respect to the parameter τ → 0 has not been, unfortunately, up to now rigorously proved, neither for (TM)→(NLH) nor (PP)→(PE).

In the global-in-time existence results below the initial data u 0 does not have to be positive. We make the assumption ϕ 0 = 0 which simplifies the presentation but it is not essential.

Our analysis illustrate that existence and well-posedness results for both toy models have nearly optimal character if one considers (both the regularity and) the size of the initial data. This in turn sheds new light on these issues for the classical parabolic-parabolic Keller-Segel system, as discussed in [START_REF] Biler | Large global solutions of the parabolic-parabolic Keller-Segel system in higher dimensions[END_REF].

Notation

In this paper we adopt the following notation and conventions. The expression A B, where A and B may depend on several parameters, means that there exists a constant c > 0, depending only on the space dimension, such that A ≤ c B. When both A B and B A we will write A ≈ B.

For a function f ∈ L 1 (R d ), the definition of the Fourier transform that we use is f (ξ) = f (x) exp(-iξ • x) dx. This definition is extended to S (R n ), the space of tempered distributions, in the usual way. The space of general distributions is denoted D (R d ).

In this paper we will deal with mild solutions. These are solutions of the integral formulation of (PP). The exact meaning of the integral must be understood in the specific functional setting.

Besov spaces

In this section, we study the local and the global solvability of (TM) in Besov-type spaces. The singularity of ∆ϕ is somehow too strong to directly perform the bilinear estimates arising from the the nonlinear term u∆ϕ in the usual way. For this reason, we will study (TM) as a perturbation of the original model (PP). Accordingly, the nonlinear term in will be decomposed as

u∆ϕ = ∇ • (u∇ϕ) -∇u • ∇ϕ. (2.1)
The first term in the right-hand is nothing but the nonlinearity of (PP). In our previous paper [START_REF] Biler | Large global solutions of the parabolic-parabolic Keller-Segel system in higher dimensions[END_REF], the solutions to (PP) were constructed in a ball of the space

E p := u ∈ L ∞ (0, ∞; L p (R d )), |||u||| p := ess sup t 1-d/(2p) u(t) p < ∞ . (2.2)
But the presence of the second term requires to have additional information on ∇u(t) p . For this reason, we will need to work in the subset of the space E p given by

F p := u ∈ L ∞ (0, ∞; W 1,p (R d )) : [[u]] p := |||u||| p + |||u||| 1,p < ∞ , (2.3)
where we define |||u||| 1,p := ess sup t>0 t 3/2-d/(2p) ∇u(t) p .

The purpose of this section is to establish the following:

Theorem 2.1. Let d ≥ 3. For 2d/3 < p < d and 2/d -1/p < 1/q < 1/d, there exist C q,d , C q,d > 0 such that if we denote γ := 1/2 -d/2(1/p -1/q) (2.4)

(when p → d, we have q → d, which forces γ → 1/2), and

u 0 Ḃ-(2-d/p) p,∞ < C q,d τ γ , (2.5)
then there exists a global mild solution u ∈ F p to (TM), such that

[[u]] p < C q,d u 0 Ḃ-(2-d/p) p,∞ . (2.6)
Such a solution is uniquely defined by condition (2.6). Moreover, u(t) -e t∆ u 0 ∈ BC(0, ∞; L d/2 ).

We will reduce the proof of this theorem to a few lemmas, established in the following subsections. For larger initial data the size condition (2.5) may not apply. In this case, local-in-time solutions can be proved to exist, under stronger regularity conditions on the initial data. Remark 2.2. As mentioned in the introduction, our main interest is for τ 1. Notice that the exponent in (2.4) is such that γ < 1/2. Thus, according to Theorem 2.1, the size of the admissible data for the global existence would be worse than O( √ τ ), as τ → +∞. This is not completely satisfactory, in view of next blowup result (Theorem 4.1 below): the latter will make evidence of a class of initial data of size O(τ ), whose solutions blow up in finite time. Thus, there is a substantial gap between our two results on global existence and on finite time blowup.

Such a gap is due to a technical limitation, related to the use of Besov spaces. Assuming that the initial data belong to slightly less rough spaces, we can almost completeley close this gap, up to a logarithmic factor in τ . This will motivate the analysis of (TM) in a different functional setting. See Section 3.

Review of known estimates in E p

In this short subsection we briefly recall without proof some results obtained in [START_REF] Biler | Large global solutions of the parabolic-parabolic Keller-Segel system in higher dimensions[END_REF] in order to study the original system (PP).

We recall the classical L p -L q inequalities for the heat semigroup (see e.g. [START_REF] Taylor | Partial Differential Equations III. Nonlinear Equations[END_REF]Ch. 15,(1.15)] or [START_REF] Giga | Navier-Stokes flow in R d with measures as initial vorticity and Morrey spaces[END_REF]), valid for 1 ≤ p ≤ q ≤ ∞ :

e t∆ f q ≤ C(d, p, q)t -d(1/p-1/q)/2 f p , ∇e t∆ f q ≤ C(d, p, q)t -1/2-d(1/p-1/q)/2 f p .
(2.7) Following [START_REF] Biler | Large global solutions of the parabolic-parabolic Keller-Segel system in higher dimensions[END_REF], let us introduce the (respectively, linear and bilinear) operators L and B given by:

Lz(t) := τ -1 t 0 ∇e τ -1 (t-s)∆ z(s) ds, (2.8) B(u, z)(t) := - t 0 ∇e (t-s)∆ • (u(s)Lz(s)) ds. (2.9)
Below, all constants are implicitly assumed to depend on d. They are also implicitly assumed to depend on p and q in the lemmas.

Lemma 2.3. If 0 ≤ 1/q ≤ 1/p < 1/q + 1/d and 0 < 1/p, then for t > 0

Lz(t) q ≤ Cτ -1/2+d/2(1/p-1/q) t -1/2+d/2q |||z||| p . Lemma 2.4. If 1/d < 1/p + 1/q ≤ 1, 0 ≤ 1/q < 1/d and the conditions of Lemma 2.3 hold, then we have |||B(u, z)||| p ≤ Cτ -1/2+d/2(1/p-1/q) |||u||| p |||z||| p .
In order to apply both lemmas, we should check the compatibility of the conditions of the exponents p and q. The two lemmas, combined, require that:

|1/p -1/d| < 1/q ≤ min{1/p, 1 -1/p}; 1/q < 1/d and 0 < 1/p, with d > 1.
This requires choosing p so that 1/d < 2/p < min{1 + 1/d, 4/d}. In other words, in order to find a suitable q to use the lemmas above, a necessary and sufficient assumption is:

d ≥ 2; max{d/2, 2d/(d + 1)} < p < 2d.
The last lemma in [START_REF] Biler | Large global solutions of the parabolic-parabolic Keller-Segel system in higher dimensions[END_REF] that we will need is the following: Lemma 2.5. If d ≥ 2 and 2d/3 < p ≤ d and u ∈ E p , we have:

B(u, u)(t) d/2 ≤ Cτ -3/2+d/p |||u||| 2 p .

New estimates in F p

Here we complete the results of the previous subsection with some additional estimates that we need to solve (TM). We will assume that d ≥ 3. Otherwise, we cannot find admissible parameters for the technical lemmas which we are using. In the same way as we did when studying the classical Keller-Segel model in [START_REF] Biler | Large global solutions of the parabolic-parabolic Keller-Segel system in higher dimensions[END_REF]Sec. 3], all constants are implicitly assumed to depend on d, as well as on p and q in the lemmas. The second term in the right-hand side of (2.1) motivates the introduction of a second bilinear operator. We define

B(u, z)(t) := t 0 e (t-s)∆ (∇u(s) • Lz(s)) ds.
(2.10)

We will proceed similarly to the proofs in [START_REF] Biler | Large global solutions of the parabolic-parabolic Keller-Segel system in higher dimensions[END_REF] of the results in the previous subsection, where we have already estimated |||B(u, z)||| p . Here we need to study the quantities

|||B(u, z)(t)||| 1,p + ||| B(u, z)(t)||| 1,p and ||| B(u, z)(t)||| p .
This will be achieved, respectively in Lemma 2.7 and Lemma 2.8 below. Both of these lemmas require a preliminary linear estimates for |||L∇v(t)||| q .

Lemma 2.6. Under the assumptions

0 ≤ 1/q ≤ 1/p < 1/q + 1/d; p < d, we have |||L∇v(t)||| q ≤ Cτ -1/2+d/2(1/p-1/q ) t -1+d/2q |||v||| 1,p .
Proof. Using (2.7) we get:

|||L∇v(t)||| q ≤ Cτ -1 t 0 [τ -1 (t -s)] -1/2-d/2(1/p-1/q ) ∇v(s) p ds ≤ Cτ -1/2+d/2(1/p-1/q ) |||v||| 1,p t 0 (t -s) -1/2-d/2(1/p-1/q ) s -3/2+d/2p ds ≤ Cτ -1/2+d/2(1/p-1/q ) t -1+d/2q |||v||| 1,p .
Lemma 2.7. If q satisfies the assumptions of Lemma 2.6 (and a fortiori those of Lemma 2.3 with q playing the role of q), and moreover

2/d < 1/p + 1/q ≤ 1; 0 ≤ 1/q < 1/d, then ∇B(u, z)(t) p + ∇ B(u, z)(t) p ≤ Cτ -1/2+d/2(1/p-1/q ) t -3/2+d/2p [[u]] p [[z]] p .
Proof. Since we have

∇B(u, z) = B(∇u, z) + B(u, ∇z) and ∇ B(u, z) = -B(∇u, z), it suffices to estimate B(∇u, z) p + B(u, ∇z) p .
First, we use (2.7) and Hölder's inequality and then, respectively, Lemma 2.3 and Lemma 2.6), ending by the change of variables s = t/s:

B(∇u, z) p ≤ C t 0 (t -s) -1/2-d/2q ∇u(s) p Lz(s) q ds ≤ Cτ -1/2+d/2(1/p-1/q ) t 0 (t -s) -1/2-d/2q s -3/2+d/2p s -1/2+d/2q ds • |||u||| 1,p |||z||| p ≤ Cτ -1/2+d/2(1/p-1/q ) t -3/2+d/2p |||u||| 1,p |||z||| p .
Then, similarly, we get:

B(u, ∇z) p ≤ C t 0 (t -s) -1/2-d/2q u(s) p L∇z(s) q ds ≤ Cτ -1/2+d/2(1/p-1/q ) t 0 (t -s) -1/2-d/2q s -1+d/2p s -1+d/2q ds • |||u||| p |||z||| 1,p ≤ Cτ -1/2+d/2(1/p-1/q ) t -3/2+d/2p |||u||| p |||z||| 1,p .
Lemma 2.8. With the same assumptions as in Lemma 2.3 for q playing the role of q, and if moreover

2/d < 1/p + 1/q ≤ 1; 0 ≤ 1/q < 2/d, we have B(u, z)(t) p ≤ Cτ -1/2+d/2(1/p-1/q ) t -1+d/2p [[u]] p [[z]] p .
Proof. Using (2.7), Hölder's inequality and then Lemma 2.3 we get

B(u, z)(t) p ≤ C t 0 (t -s) -d/2q ∇u(s) p Lz(s) q ds ≤ Cτ -1/2+d/2(1/p-1/q ) t 0 (t -s) -d/2q s -3/2+d/2p s -1/2+d/2q ds|||u||| 1,p |||z||| p ≤ Cτ -1/2+d/2(1/p-1/q ) t -1+d/2p |||u||| 1,p |||z||| p .
To find an admissible p < d, so that q from Lemma 2.4, q from Lemma 2.7 and q from Lemma 2.8 exist, a sufficient condition is

max{1/p -1/d, 2/d -1/p} < 1/q, 1/q , 1/q < min{1/d, 1 -1/p}.
(2.11)

(The conditions on q and q are in fact less restrictive than those needed for q ). This is possible for d ≥ 3 and d/2 < p < d.

Our last lemma will be useful to see that the fluctuation of the solution, i.e. the difference u -e •∆ u 0 , not only belongs to the space F p where u will be constructed, but also to

L ∞ (0, ∞; L d/2 ). Lemma 2.9. Let d ≥ 3 and d/2 < p < d. There exists a constant C τ > 0 such that B(u, u) d/2 ≤ C τ [[u]] 2 p . (2.12)
Proof.

Here and below we assume that d/2 < p < d, and r, q = q = q given by 1/r = 3/(2d) + 1/(2p); 1/p + 1/q = 1/r.

Therefore 2/d -1/p < 1/q = 3/(2d) -1/(2p) < 1/d ≤ 1 -1/p,
and moreover 1/q > 1/p -1/d, so the assumptions (2.11) are satisfied. In the four inequalities below, we use respectively (2.7), Hölder's inequality, Lemma 2.3 and the change of variables s = t/s:

B(u, u)(t) d/2 ≤ C t 0 (t -s) -d/2(1/r-2/d) ∇u(s)Lu(s) r ds ≤ C t 0 (t -s) 1/4-d/(4p) ∇u(s) p Lu(s) q ds ≤ C τ t 0 (t -s) 1/4-d/(4p) s -5/4+d/(4p) ds • [[u]] 2 p ≤ C τ [[u]] 2 p .
Proof of Theorem 2.1. To prove the first theorem, we observe that if

u 0 ∈ Ḃ-(2-d/p) p,∞ , then ∇u 0 ∈ Ḃ-2(3/2-d/2p) p,∞
, and the Besov norm of ∇u 0 is equivalent to that of u 0 . Applying twice [1, Theorem 2.34], first to u 0 , next to ∇u 0 , we deduce that e

[[B(u, v)]] p + [[ B(u, v)]] p ≤ Cτ -1/2+(d/2)(1/p-1/q) [[u]] p [[v]] p ,
that follows combining Lemmas 2.6, 2.7 and 2.8 with q = q = q . The solution u is thus constructed in a ball of F p . The fluctuation u -e •∆ u 0 also belongs to BC(0, ∞; L d/2 ). Indeed, for t > 0, we have the estimate,

u(t) -e t∆ u 0 d/2 ≤ B(u, u)(t) d/2 + B(u, u)(t) d/2 ≤ C τ [[u]] 2 p ,
that follows combining Lemma 2.5 and Lemma 2.9. The continuity with respect to t is standard and we skip it here (see [START_REF] Biler | Large global solutions of the parabolic-parabolic Keller-Segel system in higher dimensions[END_REF]).

Remark 2.10. For more regular, but possibly large, initial data, local-in-time existence of solutions to (TM) can be proved in several ways. For example, if d ≥ 3, u 0 ∈ L d/2 (R d ) and p is as in Theorem 2.1, then there exists T > 0 such that a solution u ∈ F p,T does exist. Here, F p,T is defined in the same way as F p , with the time interval (0, T ) replacing (0, ∞). The space F p,T is the subspace of F p,T of functions v satisfying the additional condition lim

t→0 + t 1-d/(2p) v(t) p + t 3/2-d/(2p) ∇v(t) p = 0.
A simple approximation argument (using the density of L p in L d/2 and the L d/2 -L p estimates for the heat hernel and its gradient) proves that if u 0 ∈ L d/2 (R d ), then e t∆ u 0 ∈ F p,T . Thus, the size condition on the data needed to apply the usual fixed point theorem can be ensured just taking T > 0 small enough (with a nontrivial dependence on u 0 and not just |u 0 | d/2 , due to the approximation procedure).

The solution that we obtain is easily proved to be in BC(0, T ; L d/2 ) and to be unique in F p,T .

Remark 2.11. An analogous analysis of (TM') is not possible since the bilinear term has a more singular stricture when using Besov spaces. This is one of the reasons why we need to introduce pseudomeasures.

Pseudomeasures

In this section we will study both toy models (TM) and (TM') in pseudomeasure spaces. We begin by recalling some definitions and statements of results from [START_REF] Biler | Large global solutions of the parabolic-parabolic Keller-Segel system in higher dimensions[END_REF]. The functional setting is more restrictive than in the previous section, but the size conditions will be less stringent when τ 1 and the proofs will be shorter. Indeed, when working in pseudomeasure spaces, it is no longer necessary to treat the nonlinearity of the toy models as a perturbation to that of (PP).

Let a ∈ R. We introduce the pseudomeasure space

PM a = {f ∈ S (R d ) : f PM a = ess sup ξ∈R d |ξ| a | f (ξ)| < ∞}, (3.1)
where f denotes the Fourier transform of the tempered distribution f . We will construct our solutions in the space

Y a = {u ∈ L ∞ loc (0, ∞; S (R d )) : u Ya = ess sup t>0, ξ∈R d t 1+(a-d)/2 |ξ| a | u(ξ, t)| < ∞}. (3.2)
When a = d -2, the space Y d-2 , agrees with the space

X = L ∞ (0, ∞; PM d-2 ),
already used in [6, Theorem 2.1] to establish a global existence result for the parabolic-elliptic Keller-Segel system (PE) for d ≥ 4 and small initial data in PM d-2 . When a = d -2, our space Y a is slightly larger than the space

Y a := Y a ∩ X considered in [6, Section 4] (d ≥ 3) or in [16] (d = 2).
We start recalling a result of [START_REF] Biler | Large global solutions of the parabolic-parabolic Keller-Segel system in higher dimensions[END_REF] for the classical parabolic-parabolic Keller-Segel system:

Theorem 3.1. Let d ≥ 2, τ > 0 and u 0 ∈ PM d-2 (R d ).
There exists κ d > 0, depending only on d, such that if

u 0 PM d-2 < κ d max 1, 3 3 τ (e ln τ ) 3 (3.3)
then (PP) possesses a global mild solution u ∈ X . More precisely:

-

If u 0 PM d-2 < κ d then u ∈ X ∩ Y d-4 3
, and this solution is unique in the ball {u :

u Y d-4 3 < 2κ d }.
-If τ ≥ e 3 , then under the weaker condition u 0 PM d-2 < 3 3 κ d τ /(e ln τ ) 3 we have u ∈ X ∩ Y d-4/ ln τ , and u is the unique solution in the ball {u : u Y d-4/ ln τ < 2κ d 3 3 τ /(e ln τ ) 3 }.

The analog of Theorem 3.1 for our two toy models is the following result:

Theorem 3.2. For d ≥ 3 the assertion of Theorem 3.1 holds for both systems (TM) and (TM').

The proof of Theorem 3.2 is conceptually the same as that of Theorem 3.1, but different technical restrictions appear on the parameters that are involved in the estimates. For reader's convenience we briefly outline the proof of Theorem 3.2. This will allow us to illustrate why the case d = 2 should be excluded for the toy model.

Lemma 3.3. Let 0 < α, β < d such that α + β > d. Then |x| -α * |x| -β = C(α, β, d)|x| -(α+β)+d , with C(α, β, d) = π d/2 Γ( d-α 2 )Γ( d-β 2 )Γ( α+β-d 2 ) Γ( α 2 )Γ( β 2 )Γ(d -α+β 2 ) . (3.4) Proof. See [6, Lemma 2.1].
Another Lemma used in [START_REF] Biler | Large global solutions of the parabolic-parabolic Keller-Segel system in higher dimensions[END_REF] (refining [16, Lemma 3.2]), is the following:

Lemma 3.4. Let s > 0, A > 0, δ > 0, 0 ≤ b ≤ 1 and δ * = min{δ, 1}. Then s 0 e -(s-σ)A σ -1+δ ds ≤ 4δ -1 * A -b s δ-b .

Pseudomeasures for the system (TM)

Taking the Fourier transform in (TM) we get

u(ξ, t) = e -t|ξ| 2 u 0 (ξ) - t 0 e -(t-s)|ξ| 2 u∆ϕ(ξ, s) ds = e -t|ξ| 2 u 0 (ξ) + (2π) -d t 0 R d e -(t-s)|ξ| 2 u(ξ -η, s)|η| 2 ϕ(η, s) dη ds. (3.5) Notice that ϕ(η, s) = τ -1 s 0 e -τ -1 (s-σ)|η| 2 u(η, σ) dσ. Therefore, u(ξ, t) = e -t|ξ| 2 u 0 (ξ) + (2π) -d t 0 s 0 R d |η| 2 τ e -(t-s)|ξ| 2 e -1 τ (s-σ)|η| 2 u(ξ -η, s) u(η, σ) dη dσ ds. (3.6)
By definition, by a solution to (TM) on (0, T ), with T > 0, we mean any function u on (0, ∞) with values in S (R d ) such that for a.e. 0 < t < T , u(t, •) is locally integrable in R d and for a.e. 0 < s < t < T and ξ ∈ R d the integrand in (3.6) is integrable in (0, t) × (0, s) × R d , and (3.6) holds for a.e. (ξ, t) ∈ R d × (0, T ).

Such solutions can be constructed via the standard fixed point algorithm as the limit u = lim k→∞ u k , where

u k+1 = e t∆ u 0 - t 0 e (t-s)∆ 1 τ u k (s)∆ s 0 e 1 τ (s-σ)∆ u k (σ) dσ ds,
for k = 1, 2, . . . , in an appropriate function space. We already gave an instance of such a construction for (PP), using pseudomeasure spaces in [5, Sec. 2].

Proof of Theorem 3.2. The case of (TM). We can write the problem (TM) in the integral form

u(t) = e t∆ u 0 + B (u, u)(t),
where B is the bilinear operator

B (u, v)(ξ, t) = (2π) -d t 0 s 0 R d |η| 2 τ e -(t-s)|ξ| 2 e -1 τ (s-σ)|η| 2 u(ξ -η, s) v(η, σ) dη dσ ds. (3.7)
Proceeding as in [START_REF] Biler | Large global solutions of the parabolic-parabolic Keller-Segel system in higher dimensions[END_REF], assuming for simplicity that u Ya = v Ya = 1 and applying Lemma 3.4, Lemma 3.

| B (u, v)|(ξ, t) ≤ (2π) -d τ -1 t 0 s 0 R d e -(t-s)|ξ| 2 e -1 τ (s-σ)|η| 2 s -1+(d-a)/2 σ -1+(d-a)/2 |ξ -η| -a |η| -a+2 dη dσ ds ≤ 8τ b-1 (2π) d (d -a) * t 0 R d e -(t-s)|ξ| 2 s -1+d-a-b |ξ -η| -a |η| -a+2-2b dη ds ≤ 8τ b-1 C(a, a -2 + 2b, d) (2π) d (d -a) * t 0 e -(t-s)|ξ| 2 s -1+d-a-b |ξ| -2a+2-2b+d ds ≤ K τ b-1 |ξ| -a t -1-(a-d)/2 ,
where,

K = K (a, b, d) = 32 C(a, a -2 + 2b, d) (2π) d (d -a) * (d -a -b) * . (3.8)
However, the conditions for the validity of this estimate are not the same as for the analogous estimate established for (PP) in [START_REF] Biler | Large global solutions of the parabolic-parabolic Keller-Segel system in higher dimensions[END_REF], because we needed to make a different choice for the parameters. Namely, we need here a < d and 0 ≤ b ≤ 1 in the first application of Lemma 3.4. We need also 0 < a < d, 0 < a + 2b -2 < d and 2a + 2b -2 > d for the application of Lemma 3.

B (u, v) Y d-4 3 b ≤ 1 4κ d b -3 τ b-1 u Y d-4 3 b v Y d-4 3 b (0 < b ≤ 1), (3.11)
for some constant κ d > 0 only depending on d ≥ 3. For the system (TM), the remainder of the proof is now carried just like in Theorem 3.1.

Remark 3.5. We do not know if the above result for (TM) holds for d = 2.

Pseudomeasures for the system (TM')

The goal of this subsection is to establish the proof, for the system (TM'), of Theorem 3.2.

Proof of Theorem 3.2. The case of (TM'). The bilinear operator associated with (TM') is

B (u, v)(t) = t 0 e (t-s)∆ (∆ϕ)(∆ψ)(s) ds with        ϕ(s) = 1 τ s 0 e (s-σ)τ -1 ∆ u(s) ds, ψ(s) = 1 τ s 0 e (s-λ)τ -1 ∆ v(λ) dλ.
(3.12)

It is convenient to rewrite B (u, v) in terms of its Fourier transform,

B (u, v)(ξ, t) = 1 (2π) d t 0 e -(t-s)|ξ| 2 (| • | 2 ϕ) * (| • | 2 ψ)(ξ, s) ds, or, more explicitly, B (u, v)(ξ, t) = τ -2 (2π) d t 0 s 0 s 0 R d e -(t-s)|ξ| 2 e -(s-σ)|ξ-η| 2 τ -1 e -(s-λ)|η| 2 τ -1 × |ξ -η| 2 |η| 2 u(ξ -η, σ) v(η,
| B (u, v)|(ξ, t) ≤ 64 τ -2+2b (2π) d (d -a) 2 * t 0 R d e -(t-s)|ξ| 2 |ξ -η| 2-a-2b |η| 2-a-2b s d-a-2b ds dη. (3.15)
Applying now Lemma 3.

3 with α = β = a + 2b -2 we get | B (u, v)|(ξ, t) ≤ 64 C(a + 2b -2, a + 2b -2, d) τ -2+2b (2π) d (d -a) 2 * |ξ| 4-2a-4b+d t 0 e -(t-s)|ξ| 2 s -1+(1+d-a-2b) ds.
The last step is the application of Lemma 3.4 with A = |ξ| 2 and δ = 1 + d -a -2b. We find, for any 0 ≤ γ ≤ 1,

| B (u, v)|(ξ, t) ≤ K (a, b, d) τ -2+2b |ξ| 4-2a-4b+d-2γ t 1+d-a-2b-γ , (3.16)
where

K (a, b, d) = 256 C(a + 2b -2, a + 2b -2, d) (2π) d (d -a) 2 * (1 + d -a -2b) * . ( 3.17)
is continuous. Indeed, for this, we just need to choose a different value for γ, namely

γ = 1 + d -a -2b.
This leads us to put the new condition 0 ≤ 1 + d -a -2b ≤ 1. But this new condition is already ensured by system (3.18). At this stage, one sees that the existence (and the uniqueness in a ball) of a solution to (TM') can be established in the same space as in Theorem 3.1 and under the same size condition (3.3).

4 Finite time blowup for both toy models: an iterative method involving the Fourier transform

In this section we establish the finite time blowup for solutions to both toy models, assuming that the initial data are large. For the blowup we will discuss to what extent the size conditions found in the previous section are sharp in the limit τ → ∞.

The first method leads to finite time blowup for (TM), for a class of initial data such that u 0 PM d-2 τ , for τ large enough. Therefore, the gap between the global existence and the blowup for such models is only logarithmic in τ .

The second method, a more "traditional" eigenfunction or moment one, will be presented in the next section for the system (TM') considered in bounded domains with the Dirichlet boundary conditions. It seems that other classical approaches to blowup questions in nonlinear parabolic equations (such as the energy method, convexity method, Fujita method), see [START_REF] Quittner | Superlinear parabolic problems. Blow-up, global existence and steady states[END_REF]Ch. 17] do not work for those systems because of the structure differences compared to the case of a single equation, and in particular, because of different diffusivities in the both equations.

The model (TM)

Let us recall that for, u 0 ∈ L d/2 , d ≥ 3, there exists a local-in-time solution to (TM), unique in the class precised in Remark 2.10. The following theorem gives a sufficient condition on u 0 , for such a solution, to blowup in finite time.

Our approach is closely related to that in [3, Theorem 3.1] which was inspired by the blowup result in [START_REF] Montgomery-Smith | Finite time blow up for a Navier-Stokes like equations[END_REF] for the so called "cheap" Navier-Stokes equations, producing lower bound estimates for the Fourier transform. The key idea is that the structure of the integral on the right hand side of (3.6) implies that, if u 0 (ξ) ≥ 0, then the positivity of the Fourier transform will be preserved by the sequence of approximate solutions u k , and so by the limit u of such sequence as k → ∞ when the latter does exist in an appropriate sense.

Consider

w 0 ∈ L 2 (R d ) defined by w 0 (ξ) = 1 I B 0 (ξ),
where 1 I E denotes the indicator function of a measurable set E, and B 0 is the ball with center (1, 0 . . . , 0) and radius 1 4 . Thus, the support of w 0 is contained in the annulus

E 0 = { 1 2 ≤ | • | ≤ 1}.
Notice that w 0 is a complex valued function, but there exist of course real valued functions u 0 ∈ L 2 (R d ), and even in S (R d ) such that u 0 ≥ w 0 : for example, one can consider a function proportional to a well-chosen Gaussian. Let t * be the maximal lifetime of the (unique) solution to (TM). There exists a constant κ d > 0 (only dependent on d) such that if A > κ d e 1/τ τ, (4.1)

then t * < 1.
Notice that the right-hand side in (4.1) behaves like τ as τ 1. Thus, the best possible size condition to be put on the initial data, in order to obtain the global existence for (TM), would be of the form u 0 τ, no matter the choice of the norm, and irrespectively of the functional setting where one constructs the solution.

Remark 4.2. In fact, our argument proves more than this. Indeed, we will actually prove the following assertion that is stronger than Theorem 4.1. Let τ > 0, A > 0, t * ≥ 1 and

u 0 (ξ) ≥ A w 0 (ξ). If u ∈ C((0, T ), L ∞ (R d )
) is a mild solution to (TM) on (0, T ) such that u ≥ 0, and if A is such that

(3t * -1 + e -4t * )A > κ d e t * /τ e t * τ, (4.2)
then we must have 0 < T ≤ t * . In particular, under condition (4.2) the L ∞ (R d )-norm of u(t) must blow up in a finite time.

Proof of Remark 4.2. By contradiction, assume T > t * and let u ∈ C((0, T ), L ∞ (R d )) be a solution to (TM) with u ≥ 0. Let us define, for k = 1, 2, . . . , the dyadic ball

B k = B k-1 + B k-1
and the dyadic annulus

E k = ξ ∈ R d : 2 k-1 ≤ |ξ| ≤ 2 k .
Let, for any integer k ≥ 1, w k = w 2 k 0 , in a such way that w k = w 2 k-1 . Then

w k = (2π) -d w k-1 * w k-1 ,
and therefore, supp

w k ⊂ B k ⊂ E k .
Lemma 4.3. Under the assumption of Theorem 4.1, for all k = 0, 1, 2, . . ., and t ∈ [0, T ], we have

u(ξ, t) ≥ β k e -2 k t 1 I {t k ≤t<t * } (t) w k (ξ), (4.3)
where (β k ) and (t k ) are two sequences defined below in (4.5) and (4.6).

Proof. For k = 0, the conclusion immediately follows from (3.7). Indeed, the first term gives us u(ξ, t) ≥ Ae -t|ξ| 2 w 0 (ξ), provided we choose β 0 = A and t 0 = 0. Let k ≥ 1. Assume that the inequality of the lemma holds with k -1 instead of k. Then, for all t k ≤ t < t * , considering only the bilinear part in (3.7), we get

u(ξ, t) ≥ (2π) -d t t k-1 s t k-1 R d |η| 2 τ e -(t-s)|ξ| 2 e -1 τ (s-σ)|η| 2 β 2 k-1 e -2 k-1 s e -2 k-1 σ × w k-1 (ξ -η) w k-1 (η) dη dσ ds ≥ (2π) -d t t k-1 R d (s -t k-1 ) 2 2k-4 τ e -(t-s)|ξ| 2 e -1 τ (t * -t k-1 )2 2k-2 β 2 k-1 e -2 k s × w k-1 (ξ -η) w k-1 (η) dη ds.
Thus, we can bound u(ξ, t) from below as follows

u(ξ, t) ≥ t t k-1 (s -t k-1 ) 2 2k-4 τ e -(t-s)2 2k e -1 τ (t * -t k-1 )2 2k-2 β 2 k-1 e -2 k s w k (ξ) ds ≥ t t k-1 (s -t k-1 )e -(t-s)2 2k ds 2 2k-4 τ e -1 τ (t * -t k-1 )2 2k-2 β 2 k-1 e -2 k t w k (ξ).
Integrating by parts then we get, always for

t k ≤ t < t * , u(ξ, t) ≥ t -t k-1 2 2k - 1 -e -(t-t k-1 )2 2k 2 4k 2 2k-4 τ e -1 τ (t * -t k-1 )2 2k-2 β 2 k-1 e -2 k t w k (ξ) ≥ (t k -t k-1 ) - 1 -e -(t * -t k-1 )2 2k 2 2k 2 -4 τ e -1 τ (t * -t k-1 )2 2k-2 β 2 k-1 e -2 k t w k (ξ).
We would like to establish the two bounds from below

(t k -t k-1 ) ≥ 3t * 2 -2k and e -1 τ (t * -t k-1 )2 2k-2 ≥ e -δ , for all k ≥ 1, (4.4)
for a suitable positive constant δ to be chosen later. This would imply that, for

t k ≤ t < t * , u(ξ, t) ≥ (3t * -1 + e -4δτ )2 -2k-4 τ -1 e -δ β 2 k-1 e -2 k t w k (ξ).
This last inequality is interesting only if the right hand side is positive. This requires 3t * -1 + e -4δτ > 0.

To ensure the validity of the first bound of (4.4) (and the condition t 0 = 0), a natural choice is

t k = t * 1 -1 4 k , (4.5)
so that t k t * . The second bound of (4.4) can be rewritten as (t * -t k )2 2k ≤ δτ ; we can ensure its validity for all k ∈ N as soon as δ is such that t * ≤ δτ.

The above conditions 3t * -1 + e -4δτ > 0 and δτ ≥ t * motivate us to restrict ourselves to t * ≥ 1.

We choose (β k ) in such a way that

β 0 = A, β k = (3t * -1 + e -4δτ )2 -2k-4 τ -1 e -δ β 2 k-1 , k = 1, 2, . . . .
This choice leads to inequality (4.3). The assertion of the lemma follows by induction.

In order to compute β k , we introduce M , such that

2 M = (3t * -1 + e -4δτ )e -δ 2 -4 τ -1 .
Notice that we have

β k = 2 M -2k β 2 k-1 for k ≥ 1.
Recalling that β 0 = A, we deduce that

β k = A 2 M -4 2 k 2 4-M +2k , k = 0, 1, . . . . (4.6)
The Fourier inversion formula and the positivity of u imply that

u(t) ∞ ≥ |u(0, t)| = (2π) -d u(t) 1 .
Moreover, by Lemma 4.3, for all k ∈ N,

u(t k ) 1 ≥ β k e -2 k t k w k 1 ≥ β k e -2 k t * w k 1 .
By the fact that w k ≥ 0 for all k ≥ 0, Fubini's theorem and the formula for the volume of the unit ball ω d = π d/2 Γ(1+d/2) , we have

w k 1 = (2π) -d w k-1 2 1 = . . . = ((2π) -d ) 2 k -1 w 0 2 k 1 = (2π) d (K d ) 2 k , with K d = 1 8 d π d/2 Γ(1 + d/2) .
We conclude that when

A2 M -4 e -t * K d > 1, then we have u(t k ) L ∞ → ∞ for k → ∞. But t k → t * and so lim sup t→t * u(t) ∞ = ∞.
This is in contradiction with the assumption that u ∈ C((0, T ), L ∞ (R d )) with T > t * . The size condition on A above can be rewritten in an equivalent form as

(3t * -1 + e -4δτ )A > κ d e δ τ e t *
where κ d > 0 is an explicit constant depending only on d. But δ ≤ t * /τ and so a sufficient blowup condition is (3t * -1 + e -4t * )A > κ d e t * (1+1/τ ) τ.

The model (TM')

Similarly to the subsection above, we want to prove explosion in finite time for large initial data. This time, we obtain a lower bound of the type τ 2 so there is a much larger discrepancy with respect to the upper estimates of the form τ (log τ ) -3 sufficient for the global existence in Sec. 3. Similar lower bounds appear in Proposition 5.1 for the problem in bounded domains.

To simplify notation, we assume from the beginning that τ ≥ 1/2, and we look for initial data which would ensure explosion at time t * = 1. As previously, our goal is to prove by induction that for some A(τ ) > 0, there is a sequence β k → ∞ such that for t k = 1 -4 -k , provided u 0 (ξ) ≥ A w 0 (ξ), we have u(ξ, t) ≥ β k e -2 k t 1 I {t k ≤t<t * } (t) w k (ξ), (4.7)

with the w k defined as above, which implies as previously u(t k ) ∞ → ∞. We begin by using the Fourier representation (3. e -(t-s)|ξ| 2 e -2 2k-2 (s-σ)/τ 2 2k-4 β k-1 w k-1 (ξ -η)e -2 k-1 t × e -2 2k-2 (s-λ)/τ 2 2k-4 β k-1 w k-1 (η)e -2 k-1 t dη dσ dλ ds ≥2 4k-8 τ -2 e -2 2k-1 (1-t k-1 )/τ t t k-1

(s -t k-1 ) 2 e (s-t)2 2k ds × β 2 k-1 e -2 k t w k (ξ) ≥2 4k-8 τ -2 e -2 2k (1+1/2τ )(1-t k-1 )

t t k-1 (s -t k-1 ) 2 ds × β 2 k-1 e -2 k t w k (ξ) ≥ 2 4k-8 3 × τ -2 e -2 2k (1+1/2τ )(1-t k-1 ) (t -t k-1 ) 3 β 2 k-1 e -2 k t w k (ξ) ≥ 2 4k-8 3 × τ -2 e -2 2k+1
(1-t k-1 ) (t k -t k-1 ) 3 β 2 k-1 e -2 k t w k (ξ).

In the last inequality, we used that τ ≥ 2. Now, with the same t k as in (4.5) (we recall that here t * = 1), we get for some K > 0 and for t k ≤ t < t * : u(ξ, t) ≥ (Kτ 2 ) -1 2 -2k e -2 k t β 2 k-1 w k (ξ).

So if we define β 0 := A; β k := (Kτ 2 ) -1 2 -2k β 2 k-1 , we have proved that u(ξ, t) ≥ β k e -2 k t 1 I {t k ≤t<t * } (t) w k (ξ).

As previously, we obtain that

β k = A/16Kτ 2 2 k 16Kτ 2 2 2k ,
which implies as above that we have blow up provided

A ≥ 16C d Kτ 2 ,
with C d > 0 only depending on the dimension.

Blowup for the toy model (TM') in bounded domains

The toy model (TM') features cross-diffusion terms which makes its existence and regularity of solutions delicate as we have already seen. Compared to the standard example of the square nonlinear

  3 and then again Lemma 3.4 with γ = -b + 1 + (d -a)/2 playing the role of b we get the estimate:

3 .For d = 3 and 3 2

 32 And we finally need d -a -b > 0 and 0 ≤ γ ≤ 1, with γ = -b + 1 + (d -a)/2, for the second application of Lemma 3.4. All these conditions are not compatible when d = 2. When d ≥ 3 they reduce to: < a < 3, or for d ≥ 4, d -2 ≤ a < d and a = 2, one can always find b satisfying such conditions. Hence the bilinear operator B : Y a × Y a → Y a is continuous when a is in such ranges. In this range, B is continuous also as an operator B : Y a × Y a → X . Indeed, to see this we just need to change the choice of the parameter γ, and to take γ = d -a -b: we should now replace the previous condition 0 ≤ -b + 1 + (d -a)/2 ≤ 1 with the new condition 0 ≤ d -a -b ≤ 1. But this new condition is satisfied given (3.10). Now, for any 0 < b ≤ 1, we can always choose, for example, a = d -4 3 b in a such way that the required conditions on a hold. Moreover, K (d -4 3 b, b, d) ≈ b -3 as b 0 and K (d -4 3 b, b, d) remains bounded as 0 < b ≤ 1 with b bounded away from zero. So, the following esimate holds:

2 (

 2 λ) ds dλ dσ dη.(3.13) When u and v both belong to Y a (with u Ya = v Ya = 1 for simplicity) we can estimate| B (u, v)|(ξ, t) ≤ τ -t-s)|ξ| 2 e -(s-σ)|ξ-η| 2 τ -1 e -(s-λ)|η| 2 τ -1× |ξ -η| 2-a |η| 2-a σ -1+(d-a)/2 λ -1+(d-a)/2 ds dλ dσ dη. (3.14) Applying twice Lemma 3.4 with A = |ξ -η| 2 τ -1 or A = |η| 2 τ -1 , and in both cases δ = (d -a)/2, we obtain, for any 0 ≤ b ≤ 1,

Theorem 4 . 1 .

 41 Let d ≥ 3, τ > 0, A > 0, and u 0 ∈ S (R d ), such that u 0 (ξ) ≥ A w 0 (ξ).

2 ( 2 t t k- 1 s t k- 1 s t k- 1 2 t

 221112 [START_REF] Lemarié-Rieusset | Small data in an optimal Banach space for the parabolic-parabolic and parabolic-elliptic Keller-Segel equations in the whole space[END_REF].u(t) = e t∆ u 0 + B (u, u)(t) with B (u, u)(ξ, t) = τ -e -(t-s)|ξ| 2 e -(s-σ)|ξ-η| 2 τ -1 e -(s-λ)|η| 2 τ -1 × |ξ -η| 2 |η| 2 u(ξ -η, σ) u(η, λ) ds dσ dλ dη.Again, we proceed by induction and assuming that (4.7) holds with k -1 instead of k, we deducethat for t k ≤ t < t * , u(ξ, t) =(2π) -d τ -2 t-s)|ξ| 2 e -(s-σ)|ξ-η| 2 τ -1 e -(s-λ)|η| 2 τ -1 × |ξ -η| 2 |η| 2 u(ξ -η, σ) u(η, λ) ds dσ dλ dη ≥(2π) -d τ -R d e -(t-s)|ξ| 2 e -(s-σ)|ξ-η| 2 τ -1 |ξ -η| 2 u(ξ -η, σ)× e -(s-λ)|η| 2 τ -1 |η| 2 u(η, λ) ds dσ dλ dη ≥(2π) -d τ -

  •∆ u 0 ∈ F p and

	e •∆ u 0 Fp	p,∞ u 0 Ḃ-(2-d/p)	.
	The usual fixed point lemma (see, e.g., [12, Theorem 13.2] or [2, Lemma 1.1.1]) applies in a ball of
	F p , thanks to the bilinear estimate

Acknowledgments

The authors thank Lucilla Corrias, Pierre-Gilles Lemarié, Miko laj Sierżȩga, Philippe Souplet for interesting conversations on the topic of that work.

The first named author would like to thank Institut Camille Jordan, Université Claude Bernard-Lyon 1 for hospitality during his sabbatical stay (Sep 2021-Jan 2022) as a fellow of Institut des Études Avancées -Collegium de Lyon, partially supported by the Polish NCN grant 2016/23/B/ST1/00434.

Let us collect all the conditions that we need on a and b for the applicability of the previous lemmas: we needed 0 ≤ b ≤ 1, d-a > 0 for the first application of Lemma 3.4, 0 < a+2b-2 < d, 2a+4b-4 > d for the application of Lemma 3.3, 1 + d -a -2b > 0 for the second application of Lemma 3.4 and finally 0 ≤ 4 + d -a -4b ≤ 2 (the latter condition correspond to the restriction 0 ≤ γ ≤ 1). Thus, we can collect all these conditions in the following system

Conditions (3.18) are not compatible for d = 1, 2. For d = 3 and 3 2 < a < 3, or otherwise for d ≥ 4 and d -2 ≤ a < d, a = 2, we can always find b satisfying the system. Therefore, we deduce that

For each 1 2 < b ≤ 1, a suitable choice for a which is compatible with system (3.18) is, e.g.,

(for d ≥ 4 we could also choose, e.g., a = d + 3 2 -3b). We have

Therefore we get

and that K (d + 4 3 -8 3 b, b, d) remains bounded for b in any interval of the form [ 1 2 + , 1], with 0 < ≤ 1 2 . Let us set β = 2b -1, so that 0 < β ≤ 1. Then, for d ≥ 3, we can find a constant κ d , depending only on d such that

for any 0 < β ≤ 1. We can also prove that, still for 0 < β ≤ 1, the operator

heat equation (NLH), the nonlinearity in (TM') is weaker that that in (NLH). Indeed, one may represent

Now it is clear that -∆ϕ(t) is somewhat smoother and smaller than 1 τ u(t) itself since (5.1) involves two smoothing (or rather "averaging") operators with respect to the both variables t and x, namely t 0 . . . ds and ∆e t∆ . Thus, one expects sufficient conditions for blowup of (TM') to be stronger with τ 0 than those for τ = 0, i.e. for equation (NLH). Note that blowup results for parabolic systems, in particular for various versions of doubly parabolic systems of chemotaxis, are scarce since there are no general comparison principles available for parabolic systems, even for linear ones. For some examples for the Keller-Segel system see, e.g., [START_REF] Winkler | Finite-time blow-up in the higher-dimensional parabolic-parabolic Keller-Segel system[END_REF][START_REF] Winkler | Blow-up profiles and life beyond blow-up in the fully parabolic Keller-Segel system[END_REF][START_REF] Winkler | Single-point blow-up in the Cauchy problem for the higher-dimensional Keller-Segel system[END_REF], where supplementary information on solutions is extracted from entropy functionals and other specific properties of those drift-diffusion systems.

Here, for simplicity of presentation, let us consider (TM') in a smooth bounded domain Ω ⊂ R d , supplemented with the homogeneous Dirichlet conditions for both u and ϕ u(x, t) = ϕ(x, t) = 0 for each x ∈ ∂Ω, t ≥ 0.

(5.2)

Proposition 5.1. There are no global-in-time classical solutions of system (TM') with τ ≥ 2 and suitably large u 0 ≥ 0, ϕ 0 ≥ 0 (of order τ 2 , τ , respectively). More precisely, if

where ψ ≥ 0 is the normalized eigenfunction of ∆ with the first eigenvalue λ, then classical solutions of system (TM') cannot be defined globally in time.

The proof of that result involves a generalization of the classical construction for nonlinear heat equations, see e.g. [15, Th. 17.1], namely the moment method invented about 50 years ago by S. Kaplan. This can be applied to sufficiently regular weak solutions. Indeed, a smooth function vanishing at the boundary (ψ below) can be used as a test function, and time-continuity in the weak sense is needed to perform the computations below.

Proof. Let ψ be the first eigenfunction of the Dirichlet Laplacian on the domain Ω ∆ψ + λψ = 0, ψ(x) = 0 for each x ∈ ∂Ω, which may be chosen so that ψ(x) > 0 for each x ∈ Ω and Ω ψ(x) dx = 1, here λ > 0 is the first eigenvalue of the Dirichlet Laplacian. Define the moments

(5.4) Solutions of the system (TM') satisfy a single, second order in time equation for ϕ τ ϕ tt = (τ + 1)∆ϕ t -∆ 2 ϕ + (∆ϕ) 2 .

(5.5)

Suppose that we have a global-in-time solution ϕ of (5.5). It is easy to see that after some integrations by parts we obtain

(5.6) Indeed, one has by the Jensen inequality

Now, the differential inequality (5.6) reads

(5.7)

Taking X(t) = e αt J(t), i.e. "zooming in" J(t), with α = λ 2 1 + 1 τ , we arrive at

Clearly, X(t) ≥ 0 but also Ẋ(t) ≥ 0 since sgn Ẋ = sgn J + αJ . Indeed, we have J =

for all τ ≥ 1. In particular, X(t) > 0 holds if 0 ≤ ϕ ≡ 0.

Therefore we obtain 2 Ẍ(t) Ẋ ≥ 2 λ 2 τ e -αt X 2 Ẋ ≥ 2 3 λ 2 τ d(e -αt X(t) 3 )/dt, so we get that

Suppose that we can guarantee that

(5.10)

Under this condition we continue with Ẋ(t)

and after an integration

Rearranging the terms we get 1

as long as the solution exists. Now, if τ ≥ 1 and

then we arrive at a contradiction for large times. From (5.9) and the inequality between the arithmetic and geometric means we obtain

Thus, to get (5.10), it suffices for τ ≥ 2 to have I(0) ≥ 2 3 λJ(0) 2 since J(0) = X(0). Now, conditions (5.10) and (5.12), sufficient for a finite time blowup, for τ ≥ 2 follow from

Thus, the condition (5.3) in Proposition 5.1 can be satisfied for suitably chosen J(0) and I(0) of order τ and τ 2 , respectively. Summarizing, with sufficiently large ϕ 0 (ϕ 0 τ ) and suitably large u 0 (u 0 τ 2 ) blowup of solutions occurs at a finite time.

Remark 5.2. An immediate generalization of the proof of Proposition 5.1 to the case of the whole space R d seems not possible since the relation -∆ψ ≥ λψ for some λ > 0 cannot be valid for a function 0 = ψ ≥ 0 where, for instance, ψ has either compact support or decays fast enough as x → ∞, and is sufficiently smooth to integrate ∆ 2 by parts four times. The above inequality fails near the boundary of the support of ψ.


Jeremy Donon

Jean-Xavier Bardaud

Shun-Ichi Valérie Brenner

Masaaki Ishiuchi

Eric Fujii
email: mfujii@res.titech.ac.jp

Gloaguen

Shun-Ichi Ishiuchi
email: ishiuchi.s.aa@m.titech.ac.jp

Masaaki Fujii

Eric Gloaguen
email: eric.gloaguen@cea.fr

Stepwise dissociation of ion pairs by water molecules: cationdependent separation mechanisms between carboxylate and alkali-earth metal ions †

come   L'archive ouverte pluridisciplinaire

Ion dissociation and pairing are ubiquitous phenomena occurring in electrolytes, 1 including biological media. 2 Despite numerous experimental studies in solution, identifying the supramolecular structure of paired ions together with their surrounding solvent molecules is often a challenging task. 3,4 Pairing between acetate and alkali cations in aqueous solution is a prototypical example, 5 where structure elucidation is often limited to the characterization of solvent-shared vs. contact ion pairs. [6][7][8] A few rare studies managed to push further the interpretation of experimental data at the cost of an advanced theoretical work with explicit solvation. 9,10 In addition to these studies in solution, the investigation of ion pairs in gas-phase experiments provide a unique opportunity to decompose, and thus better understand, the processes at play in ion pairing. [11][12][13][14][15][16] In particular, microsolvation experiments, 13,14,[16][17][18][19][20] where a controlled number of water molecules are attached to an ion pair, provide valuable information on the role played by the first surrounding water molecules on the supramolecular structure, and have been applied to investigate the pairing properties of the carboxylate group. [18][19][20] Recently, a series of IR spectra of deuterated propionate attached to both a doubly-charged calcium atom and an increasing number of water molecules (Ca 2+ , d-PrO -)(H2O)n=0-12 enabled Denton et al. to elegantly reveal the microsolvation scheme of this ion pair up to a relatively large number of water molecules. 20 However, these experiments are quite challenging and very few microsolvation schemes of ion pairs are currently known. Taking advantage of the potential of such an approach, the present work aims at revealing the role of water molecules in the ion dissociation involving cations of different nature, i.e. Ca 2+ and Ba 2+ . Their different ionic radius (100 vs. 135 pm respectively) 21 is expected to lead to different ion pairing properties as rationalized by Collins et al. who compared the surface charge density of several cations: Ca 2+ has a strong propensity to form contact ion pairs with carboxylate groups whereas Ba 2+ , with a smaller surface charge density similar to that of Na + , has a relatively smaller pairing inclination. 22 IR photodissociation spectra of H2-tagged calcium and barium acetate microsolvated clusters, i.e. (Ca 2+ , AcO -)(H2O)n=0-8 (Fig. 1) and (Ba 2+ , AcO -)(H2O)n=0-5 (Fig. 2) have been recorded in the carboxylate symmetric and antisymmetric stretches (s(CO2 -) and a(CO2 -)) as well as water bending δ(H2O) spectral ranges in order to compare their respective microsolvation schemes. While transitions above 1575 cm -1 can be readily assigned to δ(H2O), couplings of (CO2 -) with methyl bendings δ(CH3) make any assignment hazardous without the assistance of calculations. Indeed, (Ca 2+ , AcO -)(H2O)n=0-3 and (Ba 2+ , AcO -)(H2O)n=0-5 clusters present complex patterns in the 1350-1575 cm -1 spectral range. However, the isolated intense transition lying around 1520-1550 cm -1 for (Ca 2+ , AcO -)(H2O)n=4-8 clusters can be safely assigned to a(CO2 -) based on previous works. 19,20 By extrapolation, the blue-most intense transitions in this 1350-1575 cm -1 spectral range for the other clusters are tentatively assigned to a(CO2 -), which will be further confirmed by calculations. Remarkably, a blue-shift of these a(CO2 -) transitions is observed along both series (Fig. 3), but with significantly different behaviors: a rather constant blue-shift for the Ca 2+ series with a notable exception between n=4 and n=5, vs. a rapid increase followed by a plateau from n=2 for Ba 2+ . The number and relative intensities of δ(H2O) transitions appear also to be quite different for n=3-5 clusters. These observations evidence two significantly different microsolvation schemes, which need to be identified by a dedicated theoretical study. Experimental IR spectra of H2-tagged (Ca 2+ , AcO -)(H2O)n=0-8 clusters. For each spectrum, the total integrated intensity has been set to match the sum of the calculated intensities on the 1200-1700 cm -1 range (Fig. S2). Scaled harmonic frequencies (stick spectra) are also shown for comparison purposes. The vibrational mode assignment resulting from a normal mode analysis (Section S2) is shown by a color code: a(CO2 -) (red); s(CO2 -) (blue); δ(H2O) (cyan); other (black). (*) marks a doublet which may result from anharmonic couplings 23 or tunneling splitting. 24,25 In order to find theoretical structures capable to match the observations, we optimized the geometry and calculated the harmonic vibrational frequencies of 106 (resp. 82) structures for Ca 2+ -containing (resp. Ba 2+ ) clusters at the RI-B97-D3/dhf-TZVPP level, including H2 tags in our calculations (Section S2). Then, we assigned the experimental spectra to these structures according to two reliable criteria: (i) the frequency of a(CO2 -) and δ(H2O) transitions ; (ii) the number and intensity of the δ(H2O) transitions up to n=5, the decrease of the signal-to-noise ratio making this criterion less reliable from n=6. Interpretation of the coupled s(CO2 -) and δ(CH3) modes would demand, for their part, an appropriate theoretical treatment of anharmonic effects, [26][27][28] which is out of the scope of this study. Consequently, the corresponding spectral range displays obvious differences between experiments and simulations, and is thus not exploited for assignment purposes. The assigned structures are shown in Fig. 4 and5, and their calculated frequencies are reported in Fig. 1, 2 and 3 together with the vibrational mode assignment resulting from the vibrational analysis (Section S2). The agreement between both experimental and theoretical datasets is not only quite satisfactory for each structure according to the abovementioned criteria (Fig. 1 and2), but is also very good when considering the evolution along the series (Fig. 3). Focusing first on (Ca 2+ , AcO -)(H2O)n clusters (Fig. 4), the microsolvation scheme can be decomposed as a function of n. For n=1-4, one observes that water molecules fill the first solvation shell of calcium, and stay too far from the carboxylate group to make strong H-bonds (Fig. 6), forming at best an elongated H-bond for n=4 (229 pm). Then, the fifth water Please do not adjust margins Please do not adjust margins molecule starts to fill the second solvation shell of calcium. Finally, for n ≥ 6, water molecules fill both the second solvation shell of calcium and much better the first shell of carboxylate than the previous molecules (H-bonds shorter than 200 pm), thus making a H-bond network around the ion pair. As already observed for (Ca 2+ , d-PrO -)(H2O)n, 20 several spectroscopic trends accompany this microsolvation scheme, e.g. a red-shift of δ(H2O) and a blue-shift of a(CO2 -) transitions. Remarkably, the blue-shift induced by hydration has a spectacular effect on the coupling of a(CO2 -) with δ(CH3): while the intensity of a(CO2 -) is spread over several transitions for n=0-3 due to a strong coupling, this coupling suddenly disappears due to the on-going blue-shift of the a(CO2 -) transition, which becomes intense and sharp for n ≥ 4. a(CO2 -) has also proven to be an indicator of the ionic bond strength in bidentate species: (i) an increase of the splitting a-s(CO2 -) together with a blue-shift of a(CO2 -) have been observed for carboxylate groups paired with alkali cations of increasing size, revealing the relaxation of the carboxylate group resulting from the weakening of the ionic bond; 9 (ii) a blue-shift of a(CO2 -) is also observed when the cation is simultaneously involved in a competing cation- interaction, weakening the ionic bond. 15 Similarly, a(CO2 -) can be used in this microsolvation study to trace the effect of the carboxylate binding partners, i.e. both the cation and water molecules. The blue-shift observed for n=0-4 and n=5-8 (Fig. 3) reflects the gradual weakening of the ionic bond induced by hydration, which is also monitored by the increase of the cation-anion distance (Fig. 6). The fact that water molecules make efficient H-bonds (n=5-8) or not (n=0-4) has no real impact on the relaxation of the carboxylate group, which is dominated by the motion away of the cation. The passage from n=4 to n=5 is an exception, as the fifth water molecule goes to the second shell of the cation, which induces neither a blue-shift (Fig. 3) nor a significant increase of the cation-anion distance (Fig. 6). Interestingly, this characteristic "stop" of the blue-shift at n=5 has also been observed in (Ca 2+ , d-PrO -)(H2O)n experimental spectra, 20 but has not been assigned to the addition of the fifth water molecule to the second solvation shell of calcium. However, IR spectra of (Mg 2+ , AcO - )(H2O)n=0-5 Sup. Info. of ref. 19 in the OH stretch region show that strong H-bonds are indeed created at n=5, supporting a microsolvation scheme for this system similar to that proposed here for (Ca 2+ , AcO -)(H2O)n clusters. Regarding (Ba 2+ , AcO -)(H2O)n clusters (Fig. 5), the microsolvation scheme up to n=5 consists of adding molecules in the first solvation shell of barium. The first four water molecules make simultaneously elongated H-bonds with carboxylate (214 < This journal is © The Royal Society of Chemistry 20xx

Please do not adjust margins Please do not adjust margins d(CH3CO2 -•••HOH) < 237 pm, Fig. 6), and the fifth one lies opposite to acetate. Taking again a(CO2 -) as an indicator of the interactions of carboxylate with its binding partners, Fig. 3 suggests that the strong blue-shift induced by the first two water molecules correspond to a significant weakening of the ionic bond. However, additional water molecules (n=3-5) induce no additional blue-shift whereas the cation continues to move away from the carboxylate group (Fig. 6). As a consequence, the a(CO2 -) transition remains coupled to δ(CH3) up to n=5 for (Ba 2+ , AcO -)(H2O)n clusters, displaying a marked difference with (Ca 2+ , AcO -)(H2O)n clusters (Fig. 1 and2). The weakening of the ionic bond expected from the increase of the cation-anion distance (Fig. 6) thus appears in contradiction with the plateau observed for the a(CO2 -) transition. This fact can be explained if one considers that the water molecules binding to carboxylate compensate the motion away of the cation, resulting in a globally unrelaxed carboxylate group between n=2 and n=5. These microsolvation schemes provide clues about the respective role of water molecules in solvating these ion pairs. The primary role of the first water molecules is to solvate the divalent cations in both cases, the evolution of the H2O-cation distances being qualitatively the same as that observed for the hydrated cations. 29 However, water molecules interact with acetate stronger and earlier in the microsolvation scheme when paired with Ba 2+ compared to Ca 2+ (H-bond lengths to acetate, Fig. 6), which can be seen as a direct consequence of the larger ionic radius of the former cation: the relatively weaker electrostatic interaction with Ba 2+ release water molecules, which can bind more efficiently to acetate. In parallel, the cation-anion distance increases faster for Ba 2+ than for Ca 2+ (Fig. 6), showing that water molecules, by interacting simultaneously with both the cation and the anion, have a greater dissociation efficiency in the former case. Two mechanisms involved in the ion dissociation emerge from these microsolvation schemes. In the first one, water molecules mainly solvate the cation, which, in turn, gradually become a relatively worse partner for an ionic bond with acetate like for other water molecules 29 or any other ligand. Thus, water molecules play here an indirect role in the weakening of the ionic bond. This mechanism, not specific to ion pairs, is at play for both cations investigated, and appears as the first driving force of the first stages of ion dissociation. In the case of Ba 2+ , however, water molecules also bind simultaneously to the cation and the acetate group (Fig. 5), which may be seen as the beginning of the insertion of water in between the ions. This specific mechanism reveals a more direct role of water molecules in the ion dissociation, and has already been documented for other types of pairs such as alkali halide 13,17 or alkali acetate. 18 Interestingly, this mechanism is barely observed for Ca 2+ (Fig. 4), even when water strongly binds to acetate (n=6-8): none of these H-bond donor water molecules belong to the first solvation shell of Ca 2+ , highlighting the stability of this contact ion pair upon microsolvation. Furthermore, one should mention the non-intuitive reading of IR spectra. While the first mechanism is associated with a constant blue-shift of the a(CO2 -) transition, the second mechanism tends to cancel this effect, ending in a steady position of the a(CO2 -) transition. The IR spectra of e.g. n=5 clusters (Fig. 123) thus show a carboxylate group more strongly bound to its partners for Ba 2+ than for Ca 2+ , whereas the ions are more dissociated in the former case as illustrated by a steeper increase of the cation-anion distance between n=0 and n=5 (Fig. 6). This example illustrates once more the difficulty to disentangle the respective effects of the counter-ions and the surrounding water molecules on the spectroscopy of carboxylate. Finally, the frequency of the a(CO2 -) transition of free acetate ions, 1551.5 cm -1 , 10 is almost reached for (Ca 2+ , AcO -)(H2O)8 clusters (1547 cm -1 ), despite the fact that ions are far from being fully dissociated. This result echoes that of Denton et al. for (Ca 2+ , d-PrO -)(H2O)12 clusters, illustrating the difficulty to follow the complete dissociation process by IR spectroscopy, as already reported for alkali cations. 9,10 Consequently, the investigation of larger clusters by this approach should be carried out at an advanced theoretical level capable to reproduce experimental frequencies with an accuracy close to the wavenumber scale.

Figure 1 .

 1 Figure 1.Experimental IR spectra of H2-tagged (Ca 2+ , AcO -)(H2O)n=0-8 clusters. For each spectrum, the total integrated intensity has been set to match the sum of the calculated intensities on the 1200-1700 cm -1 range (Fig.S2). Scaled harmonic frequencies (stick spectra) are also shown for comparison purposes. The vibrational mode assignment resulting from a normal mode analysis (Section S2) is shown by a color code: a(CO2 -) (red); s(CO2 -) (blue); δ(H2O) (cyan); other (black). (*) marks a doublet which may result from anharmonic couplings23 or tunneling splitting.24,25

Figure 2 .

 2 Figure 2. Experimental IR spectra of H2-tagged (Ba 2+ , AcO -)(H2O)n=0-5 clusters. Same legend as Figure 1.

Figure 3 .

 3 Figure 3. Experimental (solid symbols) and theoretical (open symbols) wavenumbers  of a(CO2 -) for (Ca 2+ , AcO -)(H2O)n=0-8 (black squares) and (Ba 2+ , AcO -)(H2O)n=0-5 (red circles) clusters.

Figure 4 .

 4 Figure 4. RI-B97-D3/dhf-TZVPP theoretical structures of (Ca 2+ , AcO -)(H2O)n=0-8 clusters (H2 tags not shown). H-bonds (< 240 pm) are shown by black dotted lines, while the cation is linked to its first solvation shell by orange dotted lines.

Figure 5 .

 5 Figure 5. RI-B97-D3/dhf-TZVPP theoretical structures of (Ba 2+ , AcO -)(H2O)n=0-5 clusters (H2 tags not shown). Same legend as Figure 4.

Figure 6 .

 6 Figure 6. Evolution of intermolecular distances in (Ca 2+ , AcO -)(H2O)n=0-8 (black) and (Ba 2+ , AcO -)(H2O)n=0-5 (red) clusters. Left axis: metal ion M 2+ -CO2 -interatomic distance (cation-anion distance) ; right axis: H-bonds to acetate averaged over all the H-bonds).

Acknowledgments

JD, SII, MF and EG acknowledge the JSPS Summer program 2018. JD, JXB, VB and EG acknowledge ANR (grants IONPAIRS ANR-16-CE29-0017 and VAPOBIO ANR-20-CE29-0016). SII and MF acknowledge KAKENHI grants JP18H01938, JP19H05527, JP19K23624, JP20H00372 and JP21H04674, World Research Hub Initiative in Tokyo Institute of Technology, Cooperative Research Program of the ''Network Joint Research Center for Materials and Devices'' from the Ministry of Education, Culture, Sports, Science and Technology (MEXT).

Conflicts of interest

There are no conflicts to declare.


Samuel Kubler
email: samuel.kubler@pasteur.fr

Jean-Christophe Olivo-Marin

Thibault Lagache

STATISTICAL COUPLING BETWEEN TIME POINT-PROCESSES

Keywords:

The observation of physical phenomena often goes through the recording of discrete time series of events, that can be represented with marked point processes. The robust estimation of the correlation between two point processes can, therefore, unveil physical mechanisms underlying the observed phenomena. However, the robust estimation of correlation between two, or more, point-processes is hindered by the signal noise (leading to false and missing point detections), the important density of points, and possible time-shift between coupled points. We propose a statistical framework that uses hypothesis testing to estimate coupling between time pointprocesses. Using simulations, we show that our framework accurately estimates the coupling between two time point-processes even for noisy signal (with false point detections), for high density of points and in the presence of a time shift between coupled points. By applying our statistical framework to the recordings of neuron population activity in mouse visual cortex, we measure the functional coupling between individual neurons, and cluster neurons into functional ensembles.

INTRODUCTION

In many fields, the observation of physical events can only be done through discrete time series of events. This is the case, for example, of volcanic eruptions [START_REF] Dzierma | Eruption time series statistically examined: Probabilities of future eruptions at villarrica and llaima volcanoes, southern volcanic zone, chile[END_REF] and earthquakes [START_REF] Ogata | Seismicity analysis through point-process modeling: A review[END_REF] in geology, molecules arrival and departures form specific processes sites in cellular biology (e.g. endocytosis and pathogen entry [START_REF] Basquin | Membrane protrusion powers clathrin-independent endocytosis of interleukin-2 receptor[END_REF], and neuronal activity through the firing of action potentials from individual neurons [START_REF] Donald H Perkel | Neuronal spike trains and stochastic point processes: Ii. simultaneous spike trains[END_REF] which the study case of the experimental part of our article. The statistical characterization of relations (coupling) between two (or more) time series of events can unveil important mechanisms that underly the observed processes. For example, the observation of the sequential arrival of molecules at endocytic sites with fluorescence microscopy helped to unravel the mechanisms of cell entry [START_REF] Basquin | Membrane protrusion powers clathrin-independent endocytosis of interleukin-2 receptor[END_REF]. Another example is given by the observation of the firing of individual neurons within a population that provides information on neuronal communication and coding [START_REF] Carrillo-Reid | Controlling visually guided behavior by holographic recalling of cortical ensembles[END_REF]. Observed time series of events can often be modeled as marked point processes [START_REF] Descombes | Marked point process in image analysis[END_REF], with the point being the time location of an event and the mark its attributes (e.g. intensity, color, duration...). Therefore, the characterization of the correlation between different time series of events reduces to the estimation of the coupling between time point-processes. In the case of neuronal activity studies, the two main methods used to estimate the coupling between point-processes are either based on the estimation of underlying firing rates (i.e. the intensity of associated point-processes) [START_REF] Coeurjolly | Median-based estimation of the intensity of a spatial point process[END_REF] and the analysis of the correlation between estimated intensities, or coupling estimation is directly performed with colocalization analysis between discrete time point-processes [START_REF] Donald H Perkel | Neuronal spike trains and stochastic point processes: Ii. simultaneous spike trains[END_REF]. While the methods using the estimation of firing rates are more robust to the missing and false detections of single point events, they depend on the robust estimation of firing rates and are not well-suited for detecting the synchronization of single point events. On the other hand, the colocalization methods are sensitive to false and missed point detections. Moreover, high point density can lead to fortuitous point colocalization and overestimation of processes' coupling, whereas time shifts between coupled points can lead to coupling underestimation. To tackle these technical issues in colocalization analysis, we developed a statistical method to robustly evaluate the coupling between time point-processes, even in the presence of a time-shift between individual coupled points. Our method uses the multi-distance Ripley's K function [START_REF] Ripley | The second-order analysis of stationary point processes[END_REF] to measure the time-shifted accumulation of points from one time point-process relatively to the other. It corresponds to an adaptation of state-of-the-art statistics of pointprocesses to 1-dimensional temporal case of spiking events. This method has been developed to account for false inferred spike detections and potential time-shifts that results from spike deconvolution methods in calcium imaging [START_REF] Kubler | A robust and versatile framework to compare spike detection methods in calcium imaging of neuronal activity[END_REF]. As our method is based on statistical characterization of the Ripley's K function and hypothesis testing, it is robust to noise (false point detections) and remains accurate even for high point density. We assess the robustness of our method with synthetic simulations, and show that it outperforms state-of-the-art colocalization metrics. Finally, we use our method to compute the functional relations within a population of neurons in the visual cortex of mice, from their individual spiking activity.

METHOD

Measuring the coupling between time point-processes with the Ripley's K function

Ripley's K function introduced by Brian Ripley in the 70's [START_REF] Ripley | The second-order analysis of stationary point processes[END_REF] remains the gold-standard to measure the coupling between spatial point processes. The idea of our temporal adaptation is to measure the coupling between two time point-processes s1 = [t1, . . . , tn 1 ] and s2 = [t ′ 1 , ..., t ′ n 2 ] by creating a regular mesh grid whose topology is dependent on s1 point locations in time and by comparing the effective distribution of s2 points falling in each mesh to the distribution expected under a random distribution assumption through a statistical hypothesis rejection test. Thus, the Ripley's vector GN = [G0, . . . , Gi, . . . GN-1] embeds a 1-dimensional mathematical mesh grid implementation composed by an ensemble of N fixed size rings of radius [ri, ri+1] centered around s1 points. Thus, Ripley's function just corresponds to a statistical effective counter of s2 points that fall into rings around s1 objects with a boundary corrective term w which corrects for the potential underestimation of neighbors to points that are close to the starting-and endingpoints of the time study period. This correction is inspired by the 2-dimensional corrective term used by Ripley in [START_REF] Ripley | The second-order analysis of stationary point processes[END_REF].

Gi = |Ω| n1n2 t k ∈s 1 t ′ l ∈s 2 1 {r i ≤|t k -t ′ l |≤r i+1 } w(t k , t ′ l ) (1)
with

w(t k , t ′ l ) = 1 + 1 {|t k -t ′ l |>|t k -δΩ(t k )|} ( 2
)
where |Ω| is the length of the time period over which the two time point-processes are observed, 1 is an indicator function such that

1 {A} = 1 if A is True, 0 otherwise. δΩ is the coordinate of the clos- est boundary of the study domain ie δΩ(t k ) = 0 if t k ≤ |Ω|/2, Ω otherwise.
To detect a significant coupling between time series (point-processes) s1 and s2, we design an hypothesis testing approach. We compare how far is the effective Ripley's function counter in the rings regarding the number of points expected under a null hypothesis H0 of Complete Spatial Randmness where points are located according to an homogeneous Poisson distribution for s2 point-process. Analytical mean and standard deviation parameters are derived by calculating intersection of 1-dimensional volumes corresponding to the overlapping of rings. Under H0, Ripley's K function tends to a normal distribution [START_REF] Thibault Lagache | A statistical analysis of spatial colocalization using ripley's k function[END_REF] in accordance with the central limit theorem. Thus, the distribution of G = [Gi]i=0..N-1 is fully characterized by its mean MN = [E{Gi} = µi]i=0..N-1 and its standard deviation

ΣN = [E{G 2 i } -µ 2 i = σ 2 i ]i=0..N-1.
Using the CSR hypothesis for s2 time points and the boundary correction (Eq. 2), we compute that µi = 1 n1

t k ∈s 1 y∈Ω

1 {r i ≤|t k -y|≤r i+1 } 1 + 1 {|t k -y|>|t k -δΩ|} dy (3)
and

σ 2 i = |Ω| n 2 1 n2 t k ∈s 1 I21(t k , ri, ri+1) + t j ∈s 1 ,t j ̸ =t k I22(t k , tj, ri, ri+1) (4)
with

I21(t k , ri, ri+1) = -µ 2 i |Ω| + y∈Ω 1 {r i ≤|t k -y|≤r i+1 } 1 + 3 × 1 {|t k -y|>|t k -δΩ|} dy I22(t k , t l , ri, ri+1) ≈ - µi |Ω| + 1 {|t k -δΩ|>r i+1 } 1 {|t l -δΩ|>r i+1 } × y∈Ω 1 {r i ≤|t k -y|≤r i+1 } 1 {r i ≤|t l -y|≤r i+1 } dy
The size and number of the rings provide a maximum duration beyond which colocalization can no longer be detected and a temporal resolution to distinguish two close interactions.

Statistical test of time point-processes' coupling

To build a statistical test of time point-processes' coupling, we introduce the reduced statistics with A a correction matrix for ring's overlapping [START_REF] Thibault Lagache | Mapping molecular assemblies with fluorescence microscopy and object-based spatial statistics[END_REF]. Under the null hypothesis of s2 randomness, G is a standard normal vector (i.e. each of its component Gi ∼ N (0, 1). Therefore, a significantly high value of a vector component Gi would indicate an accumulation of coupled points around reference points corresponding to a positive coupling (while a low value would indicate a depletion of s2 points corresponding to a negative coupling).

G = A -1 G -MN ΣN (5)
Similarly to the statistical test introduced for spatial point-processes [START_REF] Thibault Lagache | Mapping molecular assemblies with fluorescence microscopy and object-based spatial statistics[END_REF], we use the maximum component of reduced Ripley vector Gmax = sup 0≤i≤N -1 Gi to test if there is at least one ring [ri; ri+1] where s2 time points accumulate significantly. To compute the pvalue associated with the observed maximum component Gmax, we compute that, ∀x > 0,

Pr{ Gmax ≥ x} = 1 -Pr{∀i ∈ [1..N -1], Gi < x} = 1 -(Pr{N (0, 1) < x}) N = 1 -cdf N (x),
where cdf(x) is the cumulative density function of the standard normal law. Finally we obtain the p -value

p-value = 1 -cdf N ( Gmax) (6)

Quantitative characterization of time point-processes' coupling

To further characterize the putative coupling between two time pointprocesses, we determine the components of Ripley's reduced vector G that are significantly high by using the universal threshold T (N ) = 2 log(N ) [START_REF] David | Ideal spatial adaptation by wavelet shrinkage[END_REF], which is widely used in image processing to determine the significant component of a signal corrupted with standard white noise. Thus, Gi > T (N ) indicates that there is a significant accumulation of s2 points at a time shift comprised between ri and ri+1 from s1 points. This allows the detection of coupling between two time point-processes at different distances and not only the detection of points co-occurence at the same time. Hence, our framework can handle coupling estimation with time shifts and varying delays. Finally, we convert the reduced Ripley vector components into a coupling probability between all the points (t k , t l ) of time point-processes s1 and s2

P (t k , t l ) = N -1 i=0
1 {r i ≤|t k -t l |<r i+1 } σi Gi1 { Gi >T (N )} Gi (7)
and extract a global coupling metric on the entire time series

GC(s1, s2) = 1 n1 t k ∈s 1 1 n l t l ∈s 2 P (t k ,t l )̸ =0 P (t k , t l ) (8)
where n l is the number of t l coupled points in s2, i.e. such that their coupling probability with t k in s1 is strictly positive. 3. RESULTS

Synthetic simulations

To validate our proposed statistical framework we use simulations where the coupling characteristics between time point-processes are known, and we compare the results of our method to standard measures of signals' correlation.

Robustness to variations of coupling level and time shift

Using simulations of time point-processes with varying coupling level and time shifts, we compare the accuracy of our statistical method with standard correlation metrics (see Appendix for details): 1) the Pearson correlation coefficient, 2) F1 score and 3) the Cluster Core Index (C3I) [START_REF] Mukherjee | Evaluating the stability of spatial keypoints via cluster core correspondence index[END_REF]. To simulate coupled time point-processes with effective coupling level pc, we first generate a reference homogeneous Poisson point-process s1, with n ∼ 30 points over Ω, with length |Ω| = 3700. Then, a proportion pcn of point process s2 are coupled to random s1 points with a time shift t shif t ∼ N (µs, σs), the other n(1 -pc) points of s2 being randomly distributed (homogeneous Poisson process over Ω). Finally, to simulate a video acquisition in biological imaging, we discretize the time period Ω = {Ωt} 1≤t≤T , where each Ωt is a time step with length ∆t = 1.

For increasing simulated levels of coupling and several time shifts ((µs, σs) = (0, 0), (1, 0), (0, 0.3), and (1, 0.3), we compare our statistical method with other standard methods (Figure 2). For µs = 0, i.e. when most coupled time points are co-occurring in same time steps Ωt, all methods provide a correct estimation of the effective coupling between time point-processes. We highlight that the non-linear, convex shape of C3I curve is due to the used index of correlation [START_REF] Mukherjee | Evaluating the stability of spatial keypoints via cluster core correspondence index[END_REF]. However, when the mean time shift increases to µs = 1, our statistical framework and C3I outperform classical correlation coefficients that are not well-equipped for time-shifted coupled point-processes. We conclude that our framework can handle coupling with or without time-shifts between coupled points and correctly estimates the simulated level of coupling.

Robustness to false and missed detections

A common issue when dealing with point-processes is the presence of false and missed detections in time point series. These artefacts are due to the presence of noise in the acquired biological signal. For example, when analyzing the spiking activity of individual neurons with calcium imaging, noise in calcium intensity traces and imperfections in deconvolution algorithms used for extracting neuron discrete spikes (time points) lead to false and missed detections (see section 3.2). In Figure 3, we measure the robustness of our method to increased levels of false and missed detections. We plot the F1 score of the proportion of couples estimated with our method.

For any simulated level of false detections (0%, 10%, 30% and 50 %), our method estimates an accurate number of points' couples (> 80% of couples are recovered). Conversely, our method is much more sensitive to missed detections. Indeed, an increased number of missed detections reduces the number of observable couples. Typically, for a percentage 0 ≤ p ≤ 1 of missed points in processes s1 and s2, the expected proportion of missed couples will be equal to 1 -p 2 . Finally, we also test the sensitivity of our method to the intensity λ0 of underlying point-processes. Indeed, an increased intensity (i.e. overall number of points over the observed period ∼ λ0|Ω|), can lead to an increased number of points' couples due to chance. We observe that, actually, our statistical method is very robust to an increased intensity of underlying point processes from λ0 = 0.1 Hz to 1 Hz.

From previous simulations, we conclude that our method outperforms other standard correlation methods when there is a time-shift between the coupled points of the two time point-processes, i.e. when coupled points are not necessarily co-occurring in the same time steps. Moreover, we assessed the robustness of our method when false detections (i.e. random points) are added to the original coupled time point-processes, as well as when the intensity of processes is overall increased. The only sensitive parameter is the level of missing points that leads to an expected decreased of detected points' couples.

Functional coupling between individual neurons

We apply our method to measure the functional coupling between individual neurons from their monitored spiking activity. We use the online dataset from [START_REF] Pérez-Ortega | Long-term stability of neuronal ensembles in mouse visual cortex[END_REF], corresponding to two-photon imaging of neuron activity in mouse visual cortex (File M1d1AS in the dataset). From calcium fluorescence traces, the exact spiking times can be obtained using spike inference techniques with variable accuracy and robustness [START_REF] Kubler | A robust and versatile framework to compare spike detection methods in calcium imaging of neuronal activity[END_REF]. A representative rasterplot of neuron spiking activity obtained with the constrained FOOPSI deconvolution algorithm [START_REF] Pnevmatikakis | Simultaneous Denoising, Deconvolution, and Demixing of Calcium Imaging Data[END_REF] is shown in Figure 1-A (|Ω| = 5 minutes, image acquisition rate = 12.3 Hz). We measure the coupling between individual neurons using the coupling index (eq. 8) (Figure 4-A). Size of used Ripley's vector (eq. 2) is equal to N = 4, with identical time-shifts ri+1 -ri = 1 frame (∼ 80 ms) for i = 0.. We identify neuronal ensembles with a Louvain clustering algorithm [START_REF] Vincent D Blondel | Fast unfolding of communities in large networks[END_REF] which maximizes graph modularity, and obtain 3 neuronal ensembles (with n ≥ 10 neurons), in agreement with [START_REF] Pérez-Ortega | Long-term stability of neuronal ensembles in mouse visual cortex[END_REF].

We conclude that our statistical framework allows the robust estimation of functional coupling between individual neurons from the calcium imaging of their spiking activity. Contrary to standard correlation metrics, it allows the estimation of points' coupling even in the presence of time-shifts, and does not require any thresholding of correlation coefficients for the network representation and clustering of neuronal connectivity.

CONCLUSION

We have proposed a statistical method for estimating the coupling between time point-processes that use the multi-distance Ripley's K function and hypothesis-testing framework. Our method is able to accurately estimate the coupling between time-shifted correlated point processes, and is robust to high intensity of point processes and false detections. The unique ability of our framework to compute the coupling between time-shifted point-processes is used to quantify the functional coupling between individual neurons imaged with fluorescent calcium indicators in mouse visual cortex [START_REF] Pérez-Ortega | Long-term stability of neuronal ensembles in mouse visual cortex[END_REF].

APPENDIX

Pearson correlation coefficient For a discretized time period Ω = {Ωt} 1≤t≤T , we introduce the indicator functions δt(s k ) = 1 {∃t l ∈ s k |t l ∈ Ωt} for each time-process {s k } k=1,2 that determines whether at least one point of each time point-process falls into the time step Ωt. The Pearson correlation coefficient is then given by r(s1, s2) = 1≤t≤T δt(s1)δt(s2) -T s1 s2

1≤t≤T δt(s1) -T s1 2 1≤t≤T δt(s2) -T s2 2

, with sk = T -1 1≤t≤T δt(s k ).

F1 score

To compute the F1 score between time point-processes s1 and s2, we set a tolerance for point-matching (tol = 2 frames) and define true positive (TP) as T P = 1≤t≤T δt(s1)δ t±tol (s2) with δ t±tol (s2) = 1 ∃t l ∈ s2|t l ∈ t-tol≤j≤t+tol {Ωj} . False positive (FP) and false negative (FN) are respectively given by F P = 1≤t≤T δt(s2)(1 -δ t±tol (s1)) and F N = 1≤t≤T δt(s1)(1δ t±tol (s2)). Finally F1 score is equal to F 1 = 2 × T P 2 × T P + (F P + F N ) .

COMPLIANCE WITH ETHICAL STANDARDS

This is a numerical simulation study for which no ethical approval was required.

Fig. 1 :

 1 Fig. 1: Statistical analysis of time point-processes A-Stochastic timepoint processes (e.g. rasterplot of individual neurons' spiking). B-Multitime-shift analysis of the coupling between two time point-processes. Depending on the accumulation of points from one point-process around the other for different time-shifts, a coupling probability is assigned to each pair of time points.

AF1Fig. 3 :

 3 Fig. 3: Robustness to false (A) and missed (B) detections, and to an increased intensity of the time point-process (C).

Fig. 4 :

 4 Fig. 4: Analyzing functional coupling between individual neurons in mouse visual cortex. A-Coupling index (eq. 8) between individual neurons computed from their extracted spiking activity (rasterplot). B-Network representation and Louvain clustering of neuronal ensembles.spike times are co-occuring in the same time step Ωt, while 1/3 are time-shifted by more than one frame. These latter, time-shifted coupled spikes are completely missed by standard correlation techniques. Using the coupling information between individual neurons, we represent neurons' couples with an undirected network graph, with edges corresponding to strictly positive coupling indexes (Figure4-B). We identify neuronal ensembles with a Louvain clustering algorithm[START_REF] Vincent D Blondel | Fast unfolding of communities in large networks[END_REF] which maximizes graph modularity, and obtain 3 neuronal ensembles (with n ≥ 10 neurons), in agreement with[START_REF] Pérez-Ortega | Long-term stability of neuronal ensembles in mouse visual cortex[END_REF]. We conclude that our statistical framework allows the robust estimation of functional coupling between individual neurons from the calcium imaging of their spiking activity. Contrary to standard correlation metrics, it allows the estimation of points' coupling even in the presence of time-shifts, and does not require any thresholding of correlation coefficients for the network representation and clustering of neuronal connectivity.

ACKNOWLEDGMENTS

No funding was received for conducting this study. The authors have no relevant financial or non-financial interests to disclose.