Elena Agliari

Francesco Alemanno

Adriano Barra

Michele Castellana

Daniele Lotito

Matthieu Piel

Inverse modeling of time-delayed interactions via the dynamic-entropy formalism

come

Introduction The ability to move is a fundamental property of many living systems, which allows them to perform key biological functions, such as searching for food and mates, or fulfilling the instructions encoded in their DNA , see e.g., [START_REF] Piel | Cell migration guided by long-lived spatial memory[END_REF][START_REF] Mishra | Cell interactions in collective cell migration[END_REF]. This occurs broadly, ranging on different scales, from animals to cells. As for the latter, advances in live imaging allowed scientists to collect big data for which statistical analysis has now become a robust and reliable tool [START_REF] Meijering | Methods for cell and particle tracking[END_REF][START_REF] Loeffler | Understanding cell fate control by continuous single-cell quantification[END_REF]. These investigations suggest that cells often migrate in groups and communicate as they move [START_REF] Suárez Korsnes | Single-cell tracking data aimed for big data analyses[END_REF]. Establishing whether such interactions are present is, in many cases, of utmost importance in order to possibly control and anticipate the evolution of a system. However, answering this question is notoriously difficult, especially in cases where unit-to-unit communication is not supported by the existence of physical bonds, e.g., in juxtacrine interactions, but it is rather induced by some signaling pathways, as in paracrine interactions. In order to inspect these possible interactions, several models and methods have been introduced and developed in the pas few decades. Among these, inferential techniques have been designed to take as input some tracks of moving units, such as birds in flocks or migrating cells, and output the parameters which provide information about the motion and the existence of interactions among the units, and/or between the units and an external source, such a predator for a bird flock [START_REF] Procaccini | Understanding cell fate control by continuous single-cell quantification[END_REF], or a cancer cell producing proteins which attract immune cells [START_REF] Alemanno | Quantifying heterogeneity to drug response in cancer-stroma kinetics[END_REF].

In particular, approaches based on the maximum entropy (ME) principle [START_REF] Jaynes | Information theory and statistical mechanics[END_REF] have been proven successful in a broad variety of contexts, and allowed one to unveil, for example, interaction patterns resulting from amino-acid sequences in protein families [START_REF] Seno | Maximum entropy approach for deducing amino acid interactions in proteins[END_REF][START_REF] Weigt | Identification of direct residue contacts in protein-protein interaction by message passing[END_REF], interaction structures of genetic networks [START_REF] Lezon | Using the principle of entropy maximization to infer genetic interaction networks from gene expression patterns[END_REF] and effective interactions in networks of neurons [START_REF] Bialek | Statistical mechanics for natural flocks of birds[END_REF] or in assemblies of immune cells [START_REF] Agliari | A statistical inference approach to reconstruct intercellular interactions in cell migration experiments[END_REF]. However, the state of the art of this inverse modeling still exhibits some limitations which make it unsatisfactory in some cases. Specifically, these models are usually based on instantaneous interactions among units: This assumption may not hold when the inter-action is mediated by some chemicals which take some time to move from the emitter to the receiver. For example, in cell-migration experiments, cell i may locally release a chemical compound along its migratory path, and cell j may cross the former path of i at a later time and thus feel a delayed interaction with i mediated by this compound. In this work, we aim at introducing, developing and testing a time-lagged ME framework that allows us to unveil not only time-lagged interactions, but also and more generally, interactions that arise indirectly as cell i alters the environment of cell j. Remarkably, this approach allows us to estimate the time lag of the interactions, and therefore speculate about the kind of messengers involved or the effectiveness of different embeddings.

Settings Let us consider the motion of N particles (i.e. units) in a D-dimensional space, and denote the set of trajectories by

x t i , 1 ≤ i ≤ N, 1 ≤ t ≤ N T + 1, (1)
where x t i is a D-dimensional position vector of particle i at time t, and the superscript index indicates instants of time separated by δt, for a total of N T + 1 temporal samples. The velocity and direction of motion of each particle i at time t then read, respectively,

v t i ≡ x t+1 i -x t i δt , s t i ≡ v t i |v t i | . (2)
The instantaneous mean alignment of the moving units at time t is given by

m(s t ) ≡ 1 N N i=1 s t i , (3)
which shall also be referred to as the 'magnetization', where s t ≡ {s t 1 , ..., s t N } are interpreted as a configuration of N Heisenberg spins. For practical purposes, we assume that the motion of these units may be temporally correlated up to a maximum temporal window whose length is τ M , with 0 ≤ τ M < N T . Hereafter we assume that the correct value for τ M is known a priori; in the following, once developed the necessary tools, we will provide a means for determining τ M from the data in a self-consistent manner. A natural observable to quantify this temporal correlation is the two-point correlation function R with delay τ R (τ, s) ≡

1 N T -τ M NT t=τM+1 m(s t )•m(s t-τ ), 1 ≤ τ ≤ τ M .
(4) It is also convenient to introduce the temporally averaged magnetization M

M [s] ≡ 1 N T -τ M NT t=τM+1 m(s t ), (5)
which quantifies the direction in which the units move on average.

Inference method Following the ME principle [START_REF] Jaynes | Information theory and statistical mechanics[END_REF][START_REF] Agliari | A statistical inference approach to reconstruct intercellular interactions in cell migration experiments[END_REF][START_REF] Cavagna | Dynamical maximum entropy approach to flocking[END_REF], we define the path entropy S as

S [P ] = -DsP [s] log P [s] (6)
where Ds ≡ NT t=1 N i=1 ds t i , ds t i is the surface element of the (D -1)-dimensional sphere. We seek P as the minimal probability measure whose average magnetization [START_REF] Suárez Korsnes | Single-cell tracking data aimed for big data analyses[END_REF] and two-point correlation function (4) match their respective empirical values [START_REF] Jaynes | Information theory and statistical mechanics[END_REF], referred to as M E and R E , respectively:

max P S[P ] (7)
subject to

DsP [s]M (s) = M E , (8)
DsP [s]R (τ, s) = R E (τ ), 1 ≤ τ ≤ τ M , (9)
DsP [s] = 1. (10)
To solve Eqs. [START_REF] Alemanno | Quantifying heterogeneity to drug response in cancer-stroma kinetics[END_REF] to [START_REF] Weigt | Identification of direct residue contacts in protein-protein interaction by message passing[END_REF], we introduce the Lagrangian multipliers J ≡ {J τ } τM τ =1 and H ≡ {H ℓ } D ℓ=1 , and obtain

P [s] = 1 Z exp N T t=τ M +1 N i=1 s t i • 1 N τ M τ =1 Jτ N j=1 s t-τ j + H , (11)
where Z = DsP [s] is a normalization constant, see Supplementary Material (SM) for details.

We notice that P [s] can be looked at as a Boltzmann-Gibbs distribution, and the argument of the exponential in it can be interpreted as a Hamiltonian: The direction of the ith unit at time t tends align with the direction of H-which is analogous to an external field in the direct problem-and, if J τ > 0 (J τ < 0), the direction of the ith particle at time t tends to aligned (misaligned) with the average direction at time t -τ . Now that we obtained an analytical expression for P [s], we infer the parameters J and H from a set of data. Specifically, let us suppose that we experimentally observe a system of N particles at regular time intervals for a total of N T + 1 timepoints. Then, according to (2), we know the velocity and the direction for each particle and each time:

s E = {s E t i |1 ≤ i ≤ N, 1 ≤ t ≤ N T }, (12)
where the superscript E stands for 'empirical', and we obtain

m E t ≡ m E (s E t ) = 1 N N i=1 s E t i , (13)
M E [s E ] = 1 N T -τ M NT t=τM+1 m E t , (14)
R E (τ, s E ) = 1 N T -τ M NT t=τM+1 m E t • m E t-τ . ( 15
)
By requiring that the distribution [START_REF] Lezon | Using the principle of entropy maximization to infer genetic interaction networks from gene expression patterns[END_REF] satisfies Eqs. ( 8) to [START_REF] Weigt | Identification of direct residue contacts in protein-protein interaction by message passing[END_REF], we may estimate the parameters J and H, and thus build a model for the dataset s E . However, this route (or analogous ones, e.g. those based on likelihood methods, see SM) entails several numerical issues as for the stability of the solution and the computational time, which would make the procedure cumbersome. Here, looking for an approach that can be routinely applied, we resort to an approximation that greatly simplifies and speeds up the calculation, whose consistency can be checked a posteriori. First, proceeding along the lines of the mean-field approximation in statistical mechanics [START_REF] Huang | Statistical Mechanics, 2nd Edition[END_REF], in Eq. ( 11) we replace the term in parentheses with its empirical counterpart

1 N τM τ =1 J τ N j=1 s t-τ j → τM τ =1 J τ m E t-τ . (16)
We denote the resulting distribution by P A to highlight that it is an approximation of P . Denoting by ⟨⟩ A the average performed with respect to P A , we get

⟨m t ⟩ A = NT t=τM+1 N i=1 ds t i P A [s] 1 N N i=1 s t i = = M τM τ =1 J τ m E t-τ + H if τ M ≤ t ≤ N T , m E t if 1 ≤ t ≤ τ M , (17)
where

M(x) ≡ x |x| I 1 (|x|) I 0 (|x|) , (18)
and I 0 , I 1 are the hyperbolic Bessel functions of order 0 and 1, respectively [START_REF] Abramowitz | Handbook of mathematical functions: with formulas, graphs and mathematical tables[END_REF].

Next, we notice that P A is factorized with respect to every variable s t i , thus the average of a product of magnetizations at different times equals the product of averages:

⟨m t m t ′ ⟩ A = ⟨m t ⟩ A ⟨m t ′ ⟩ A .
(19) This property allows us to recast Eqs. ( 8) and [START_REF] Seno | Maximum entropy approach for deducing amino acid interactions in proteins[END_REF] as

M E = 1 N T -τ M N T t=τ M +1 ⟨mt⟩ A , (20)
R E (τ ) = 1 N T -τ M N T t=τ M +1 ⟨mt⟩ A • ⟨m t-τ ⟩ A , 1 ≤ τ ≤ τ M . ( 21
)
These equations are simultaneously fulfilled by requiring that

⟨m t ⟩ A = m E t , τ M + 1 ≤ t ≤ N T , (22)
which, exploiting Eq. ( 17), can be written as

M τM τ =1 J τ m E t-τ + H = m E t , τ M + 1 ≤ t ≤ N T .
(23) Now, by introducing the inverse function

G E t ≡ M -1 m E t , we obtain the linear system τM τ =1 J τ m E t-τ + H = G E t , τ M + 1 ≤ t ≤ N T . ( 24
)
Before proceeding further, we stress that Eq. ( 22) constitutes a stricter condition than Eqs. ( 20) and ( 21): in order to avoid that ( 24) is overconstrained, we included an additional and controllable term in it, i.e.,

τM τ =1 J τ m E t-τ + H + σϵ t = G E t
, where ϵ t ∼ N (0 D , 1 D ) and σ ≥ 0 is an additional parameter to be determined. In the procedure above, we have modified the model in such a way that the relation [START_REF] Sáez | Atp promotes the fast migration of dendritic cells through the activity of pannexin 1 channels and p2x7 receptors[END_REF] does not hold deterministically. Thus, the related log-likelihood l A reads

l A (J , H, σ|s E ) ≡ - 1 2σ 2 NT t=τM+1 τM τ =1 J τ m E t-τ + H -G E t 2 -(N T -τ M ) D 2 log 2πσ 2 . ( 25
)
We now obtain the optimal estimates for J , H and σ.

Setting ∇ H l A = 0, we get

H = ⟨G E t ⟩ t - τM τ =1 J τ m E t-τ t (26)
and, proceeding analogously for ∇ Jτ l A = 0, we obtain

J = A -1 • B, (27)
with

A λτ ≡ m E t-τ • m E t-λ t -m E t-τ t • m E t-λ t ,(28) B λ ≡ ⟨G E t • m E t-λ ⟩ t -⟨G E t ⟩ t • ⟨m E t-λ ⟩ t . (29)
Finally, setting ∇ σ l A = 0, we obtain

σ 2 = 1 D H + τM τ =1 J τ m E t-τ -G E t 2 t . ( 30
)
where

⟨•⟩ t ≡ 1 N T -τ M NT t=τM+1 • (31)
denotes the time average.

As detailed in the SM, an estimate of the errors for the inferred parameters can be obtained by using the Fisher information for the log-likelihood [START_REF] Frieden | Science from Fisher Information: A Unification[END_REF]. In fact, one can prove that Var (J , H, σ) ≈ diag I -1 .

(32)

where the approximation holds as long as the sample size is large enough and I is the Fisher information matrix, which reads in the present setting

I = N T -τ M σ 2   m E t-τ • m E t-λ t m E t-τ t 0 m E t-λ t δ 0 0 0 2D   .
Inference of the interaction delay So far in our analysis, we assumed that τ M is a known constant, yet in most practical inference scenarios this quantity is not accessible and, actually, its estimate would provide valuable knowledge about the system under study. To estimate τ M , we rely on the Akaike Information Criterion (AIC) [START_REF] Banks | AIC under the framework of least squares estimation[END_REF], which assesses the quality of different models, each characterised by a different value of τ M , providing a tool for model selection. Given a dataset and a model with a mean squared error ϵ 2 , n p parameters and n o observations, the quantity

AIC ≡ 2n p n o + log ϵ 2 (33)
provides a measure of the quality of such model for the data. Notice that, the AIC favors models with a small error and a small number of parameters. In the analysis above, the parameters in Eq. (33) are

n o = D (N T -τ M ) , n p = τ M +D+1, ϵ 2 = N T N T -τ M σ 2 , thus AIC = 2 τ M + D + 1 D (N T -τ M ) + log N T σ 2 N T -τ M . (34)
In the following, after performing inference for multiple values of τ M and determining σ 2 from Eq. (30), we will evaluate Eq. (34) for each of the inferred models, and we will choose the model with the smallest value of AIC as the best model.

Testing the model We test our inference framework on both synthetic and real datasets. The statistical models used to generate the synthetic data is the Heisemberg-Kuramoto model (that is the direct model whose inverse one is exactly the ME) and the standard Vicsek model (VM) [START_REF] Vicsek | Novel type of phase transition in a system of self-driven particles[END_REF]. Each model under consideration has its own set of parameters, which we denote by the planted parameters. In the first setting inferred and planted parameters must coincide, indeed their agreement constitutes a positive test of the inference method. On the other hand, if the planted model is the VM, there is no one-to-one correspondence between inferred and planted parameters, yet we can still check that the inferred model generates a dataset statistically compatible with the planted one.

The validation of the ME model against itself is shown in Fig. 1. In particular Figure 1a shows that planted and inferred parameters match along the diagonal x = y, displaying an excellent agreement with each other. In addition, Fig. 1d shows that the inference framework is able to detect the exact number of interactions most of the times: in some cases, the number of interactions is overestimated, but we verified that the magnitude of the extra interactions is statistically compatible with zero within the error on the planted data. This can be seen from Fig. 1b: as J planted τ >τ planted ≡ 0, if J inferred τ >τ planted ̸ = 0, its magnitude corresponds to the absolute error and in these cases we reveal only small errors (∼ O(10 -2 )), while in Fig. 1c we report the Pearson correlation matrix between planted and inferred parameters.

Next, we consider the VM, which is defined by the

|M M E -M V | |M V | + |R V (1)-R M E (1)| R V (1)
as functions of H x and T , d) Logarithm of the total number of inferred interactions τ M for the ME model as a function of T and H x .

following dynamical equations [START_REF] Vicsek | Novel type of phase transition in a system of self-driven particles[END_REF][START_REF] Ginelli | Relevance of metric-free interactions in flocking phenomena[END_REF]

v t+1 k = v 0 Γ 1 |n t k | j∈n t k v t j +H + √ 2T η t k , η t k ∼ N (0, 1) (35)
x t+1 k = x t k + ∆t v t+1 k , k = 1, • • • , N, (36)
where v t i , x t i are the velocity and position of the ith particle at time t, respectively. Also, n c and n t 1 , • • • , n t nc are the total number of topological neighbors and the set of topological neighbors of the ith particle, including the i-th particle itself, at time t, respectively, further

Γ(v) ≡ v |v| , (37)
H is the external field applied to each particle, and T , which is analogous to a temperature, sets the noise strength. In this model, the particles are self-propelled, i.e., v t j = v 0 . Note that, as we discussed above, we denoted the planted parameter H with the same symbol as the magnetic field in the ME distribution [START_REF] Lezon | Using the principle of entropy maximization to infer genetic interaction networks from gene expression patterns[END_REF], but they need not be the same quantity. Finally, the magnetiza-tion of the VM is given by the expression

m t = 1 N v 0 N k=1 v t k , (38)
which is analogous to Eq. ( 3). The validation procedure for the ME inference in the case of the VM is different from the previous one. In what follows, we set H y = 0, i.e., the parameters of the VM are T and H x . As these parameters are different from those of the ME model, a comparison between planted and inferred parameters would be pointless.

Instead, a method that allows for validation in this case is to compare the lowest-order moments of some observables for the ME model and the VM. In particular, we measure the magnetization M and the correlation function R(1) = R(τ = 1, s) in both models, and check if their values are statistically compatible: if they do, the inference procedure will prove successful. Otherwise our procedure does not yield a valid ME model for the given configuration of the parameters T and H. Figure 2a, b andc provide a positive test of the inference method, indicating that the values of the correlations R( 1) and magnetization M from the inferred model match those of the planted model for most realisations of the planted parameters. Given that the coordinates at any time step in the VM depend only on the ones at the preceding step, we expect that, among the inferred couplings J(τ ), only J(1) is significantly different from zero: this prediction is mostly confirmed in Fig. 2d, where, apart from a small region in the parameter space (T, H x ), τ M is always equal to one.

Finally, to highlight the inconsistency that would result from neglecting interaction delays, we consider real datasets obtained in two-dimensional migration experiments where leukocytes move in the presence of a chemokine gradient, that acts as a chemoattractant [START_REF] Agliari | A statistical inference approach to reconstruct intercellular interactions in cell migration experiments[END_REF]: note that all the aspects of biological relevance will be discussed elsewhere, here we just focus on the methodology. The positions of cells were recorded at regular time intervals (2 minutes) and we studied two regions of the experiment, the former (region A) where the gradient is low and the latter (region B) where the gradient is high, providing two datasets suitable to be addressed by our approach: we refer to the last section of SM for technical details. We estimate

δAIC(τ M ) ≡ AIC(τ M ) -min τ ′ M AIC(τ ′ M ),
i.e., the difference between AIC for a given τ M and AIC for the optimal value of τ M , as well as the inferred parameters J(τ ), H, σ and the correlation function

C(τ ) ≡ m(s t ) • m(s t-τ ) t -m(s t ) t • m(s t-τ ) t (39)
for the optimal value of τ M . These results are reported in Figure 3 for both datasets A and B, highlighting marked differences between the prediction by the standard ME and our delayed inference technique. Indeed, as shown in Fig. 3a, the optimal τ is 13 for region A and 14 for region B, and, as shown in Fig. 3b, the corresponding couplings J(τ = 13) and J(τ = 14) are both positive, while in the τ → 0 region of these plots the interactions are actually repulsive: accordingly the connected correlation functions provided in Fig. 3d are small but systematically different from zero for τ < τ M , then they vanish.

As expected, the method correctly returns a higher field in region B, w.r.t. A, as shown in Fig. 3c.

To summarize, in this letter we revised the standard Maximum Entropy (ME) inferential protocol to extend its usage to scenarios where communication among units are not instantaneous. While from a gnoseological perspective such a general setting should always be preferable to versions where equal-time quantities are constrained, when the communication among units is reasonably faster than the fastest timescale of their dynamics, the standard ME and our method return the same outcome. However, when nothing is known on the speed of the interaction or it must be assumed as slow, here we have demonstrated how it should be mandatory the usage of an inference setting which incorporates temporal delays.

We acknowledge financial support from PNRR MUR (PE0000013-FAIR), from MAECI (BulBul, F85F21006230001), from PRIN MUR (20229T9EAT), from INFN (FieldTurb) and from Sapienza University of Rome (RM120172B8066CB0, RM12117A8590B3FA). The Authors are grateful to Raphael Voituriez and Theresa Jakuszeit for useful discussions. In this Supplementary Material we aim to provide a detailed mathematical formulation of the problem of inferring delayed interactions among dynamical units (i.e., an ensemble of migrating cells) by providing to the inferential protocol solely their trajectories: we want to solve this problem via the time-delayed dynamical maximum entropy approach and to present results we achieved with this tool on extensive synthetic and real tests. We naturally split this document in two main sections: in the former, Section S1, we provide a detailed derivation of the time-delayed Dynamic Entropy Formalism, in the latter, Section S2, we vaildate the model on three different datasets. The exposition of the method is split in Sec. S1.1, where we formulate the problem we aim to address as a maximum entropy extremization and in Sec. S1.2 where we expose our strategy of solution, based on log-likelihood estimation. In particular in Sec. S1.2.1 we prove an analytical estimate for the parameters inferred by solving the related gradient equations for the score, in Sec. S1.2.2 we rely on the Fisher Information (i.e., we inspect the Hessian of the log-likelihood) to quantify errors in the inference procedure and finally in Sec. S1.2.3 we rely on Akaike Information Criterion to select a cufoff on the past evolution to include in the inference to make it optimal (w.r.t. the Akaike criterion), providing an effective estimate for the elapsed time involved to let the interaction work. As for the datasets, the first -analyzed in Sec. S2-is a synthetic dataset storing trajectories of N particles generated by the Heisember-Kuramoto model [START_REF] Agliari | A statistical inference approach to reconstruct intercellular interactions in cell migration experiments[END_REF] with delay (namely exactly the direct problem whose inverse modelling we tackled in the previous Section), the second -analyzed in Sec. S2.1-is a synthetid dataset whose trajectories are generated accordingly to the Vicsek model [START_REF] Vicsek | Novel type of phase transition in a system of self-driven particles[END_REF] (a Markovian self-propelled model widely used to mimic cell's behavior) and the third -analyzed in Sec. S2.3-is a real dataset containing the trajectories of dendritic cells globally migrating toward a chemoattractive source [START_REF] Agliari | A statistical inference approach to reconstruct intercellular interactions in cell migration experiments[END_REF].

Dataset [A] Dataset [B]

S1 Dynamic maximum-entropy formalism

S1.1 Problem formulation

Let us consider the motion of N particles in a D-dimensional space, and denote by

x t i , 1 ≤ i ≤ N, 1 ≤ t ≤ N T + 1 (S1)
the set of trajectories, where x t i is the position of the particle i at time t and the superscript index indicates instants of time separated by δt, for a total of N T + 1 temporal samples. By using these collected positions x t i , we can evaluate the velocity v and direction of motion s of each particle i at time t:

v t i ≡ x t+1 i -x t i δt , s t i ≡ v t i |v t i | , 1 ≤ i ≤ N, 1 ≤ t ≤ N T . (S2)
Focusing on the collection of variables s ≡ {s t i }, 1 ≤ i ≤ N, 1 ≤ t ≤ N T , a natural framework to determine a probability measure for these units is the path-integral formalism [START_REF] Agliari | A statistical inference approach to reconstruct intercellular interactions in cell migration experiments[END_REF][START_REF] Cavagna | Dynamical maximum entropy approach to flocking[END_REF]: We define the path entropy S as

S [P ] = -DsP [s] log P [s] (S3)
where Ds ≡ NT t=1 N i=1 ds t i is the integration measure and each ds t i is a surface element of the (D -1)-dimensional sphere.

The average alignment of the moving units at time t is given by

m(s t ) ≡ 1 N N i=1 s t i , (S4)
where

s t ≡ {s t i }, 1 ≤ i ≤ N. (S5)
The latter can be interpreted as a configuration of soft spins that evolve in time, thus, we will refer to m(s t ) as the magnetization at time t. Furthermore, to quantify the direction towards which the units move on average, we introduce the temporally averaged magnetization M :

M [s] ≡ 1 N T -τ M NT t=τM+1 m(s t ). ( S6
)
In our analysis, we imagine that the motion of the units above can be temporally correlated up to a maximum temporal window whose length is τ M with 0 ≤ τ M < N T . A natural observable to quantify this persistence is the two-point correlation function R with delay τ

R (τ, s) ≡ 1 N T -τ M NT t=τM+1 m s t • m s t-τ , 1 ≤ τ ≤ τ M . (S7)
For the moment we assume that the correct value for τ M is known a priori; in the following sections, after having developed the necessary tools, we will provide a means for selecting τ M in a self-consistent manner for any dataset. Notice that the expression in (S7) implements a "mean-field" correlation as the direction of each unit i is related to the direction of any other unit at the previous τ M time steps, regardless of their spatial distance.

Following the maximum-entropy (ME) principle [START_REF] Jaynes | Information theory and statistical mechanics[END_REF], we seek P as the minimal probability measure whose average magnetization (S6) and two-point correlation function (S7) match their respective empirical values:

max P S[P ] (S8) subject to DsP [s]M (s) = M E , (S9) DsP [s]R (τ, s) = R E (τ ), 1 ≤ τ ≤ τ M , (S10)
DsP [s] = 1, (S11)
where M E and R E are the empirical values of the average magnetization and two-point correlation, respectively.

To solve Eq. (S8), we use the method of Lagrangian multipliers. The Lagrangian S ⋆ reads

S ⋆ [P, J, H] = - 1 N (N T -τ M ) DsP [s] log P [s]+ + τM τ =1 J τ DsP [s]R (τ, s) -R E (τ ) + + H • DsP [s]M (s) -M E + + ζ N (N T -τ M ) DsP [s] -1 . (S12)
where

J ≡ {J τ } τM τ =1 , H ≡ {H ℓ } D ℓ=1
, ζ play as Lagrangian multipliers and the normalization factor

1 N (N T -τ M )
ensures that S ⋆ is an intensive quantity with respect to N and N T , and that the final expression for P [s] is well defined (vide infra). The Lagrangian can be easily extremised with respect to P :

δS ⋆ δP = - 1 N (N T -τ M ) (log P [s] + 1 -ζ) + τM τ =1 J τ R (τ, s) + H • M (s) = 0 (S13) yielding P [s] = 1 Z exp N (N T -τ M ) τM τ =1 J τ R (τ, s) + H • M (s) (S14)
where the Lagrangian multipliers J , H are determined implicitly by Eqs. (S9) to (S11), and

Z ≡ exp{(1 - ζ)/[N (N T -τ M )]}.
By substituting Eqs. (S6) and (S7) into Eq. (S14), we get

P [s] = 1 Z exp N NT t=τM+1 m s t • τM τ =1 J τ m s t-τ + H . ( S15
)
Finally by substituting Eq. (S4) into Eq. (S15), we obtain

P [s] = 1 Z exp   NT t=τM+1 N i=1 s t i •   1 N τM τ =1 J τ N j=1 s t-τ j + H     , (S16)
where the normalisation constant Z reads

Z = Ds exp   NT t=τM+1 N i=1 s t i •   1 N τM τ =1 J τ N j=1 s t-τ j + H     . (S17)
We emphasize that the resulting P [s] can be looked at as a Boltzmann-Gibbs distribution in such a way that the argument of the exponential plays as an energy function and this, consistently with thermodynamic principles, scales linearly with the effective number of degrees of freedom N and N T . More generally, the structure of P [s] suggests that the direction of the i-th particle at time t tends to be aligned with the direction given by H -which plays as an external field -and that, if J τ > 0 (resp. J τ < 0), the direction of the i-th particle at time t tends to be aligned (resp. misaligned) with the average direction at time (t -τ ) -which plays as an internal field.

S1.2 Solution method

In this section we present our strategy to obtain an estimate for the parameters J and H starting from a sample of experimental data. As we will detail, we first recast the problem into a maximum likelihood setting, which allows us to implement a simplification in the expression for P [s], yielding to an approximation expression P A [s]. The latter can then be handled analytically in such a way we can finally obtain a set of algebraic equations for J and H which can be solved numerically and straightforwardly.

Suppose that we experimentally observe a system of N particles at regular time intervals for a total of N T + 1 timepoints, then according to Eq. (S2) we can evaluate the velocity v and the direction s for each particle and time point, overall collecting the dataset:

s E = {s E,t i |1 ≤ i ≤ N, 1 ≤ t ≤ N T }, (S18)
whence we can obtain

m E (s E,t ) = 1 N N i=1 s E,t i , ( S19
)

M E [s E ] = 1 N T -τ M NT t=τM+1 m E (s E,t ) (S20) R E (τ, s E ) = 1 N T -τ M NT t=τM+1 m E (s E,t ) • m E (s E,t-τ ) (S21)
In order to estimate the parameters J and H appearing in Eq. (S16), and thus build a model for the dataset s E , we can use the method of maximum likelihood [START_REF] Cox | Principles of Statistical Inference[END_REF]. The log-likelihood l of the model (S16) given the dataset s E is

l(J , H|s E ) ≡ log P [s E |J , H] = NT t=τM+1 N i=1 s E,t i •   1 N τM τ =1 J τ N j=1 s E,t-τ j + H   (S22) -log Ds exp NT t=τM+1 N i=1 s t i • 1 N τM τ =1 J τ N j=1 s t-τ j + H .
The value of J and H that maximize l(J , H|s E ) represent the maximum likelihood estimate for the parameters. We therefore derive l(J , H|s E ) in Eq. (S22) and solve for

∂l ∂J = 0, ∂l ∂H = 0. (S23)
As well known, with this operation we retrieve the constraints (S9) and (S10). In other words, the estimate of J and H obtained by this maximum likelihood framework, namely by solving Eq. (S23), coincides with the one we would obtain by the maximum entropy framework, namely requiring the fulfilment of Eqs.(S9)-(S11) with the distribution (S16). However, for both these routes the solution requires the evaluation of high-dimensional multivariable integrals, like the one appearing in the last term of Eq. (S22), which entail several numerical issues as for the stability of the solution and the computational time (e.g., in our datasets the number of single-variable integrals to perform is O(10 4 ), being N T = 150 and N = 300). Thus, we resort to an approximation: First, in Eq. (S16) we replace the term in parentheses with its empirical counterpart

1 N τM τ =1 J τ N j=1 s t-τ j → τM τ =1 J τ m E t-τ , (S24)
along the lines of the mean-field approximation in statistical mechanics [START_REF] Huang | Statistical Mechanics, 2nd Edition[END_REF].

With the substitution (S24), the dependence of P [s] on the variables s t i , t = 1, • • • , N T -1 through the term in parentheses is dropped, thus, although P [s] remains defined over the whole set of variables, it no longer depends on s t i , t = 1, • • • , τ M . In order to restore this dependence, at least formally, we introduce some constraints that set the missing variables to their experimental values:

P [s] → P A [s] ≡ 1 Z A N i=1 NT t=τM+1 exp s t i • τM τ =1 J τ m E t-τ + H τM t=1 δ s t i -s E,t i , (S25)
with

Z A ≡ NT t=τM+1 N i=1 ds t i exp s t i • τM τ =1 J τ m E t-τ + H , ( S26
)
where the subscript A stands for approximated. In this way, we are still able to evaluate all the momenta of the variables s, including those involving s t i , t = 1, • • • , τ M . Denoting by ⟨⟩ A the average performed with respect to P A , i.e., for the generic observable F (s),

⟨F (s)⟩ A ≡ DsP A [s]F (s), ( S27
)
we notice that

⟨m t ⟩ A = NT t=τM+1 N i=1 ds t i P A [s] 1 N N i=1 s t i = = M τM τ =1 J τ m E t-τ + H τ M ≤ t ≤ N T m E t 1 ≤ t ≤ τ M (S28)
where the vector function M is defined as

M(x) ≡ x |x| I 1 (|x|) I 0 (|x|) (S29)
and I 0 , I 1 are hyperbolic Bessel functions (also known as modified Bessel function of the first kind) of order 0 and 1, respectively [START_REF] Abramowitz | Handbook of mathematical functions: with formulas, graphs and mathematical tables[END_REF].

In addition, we notice that P A is factored with respect to every variable s t i , thus any average of products of magnetizations at different times t and t ′ equals a product of averages:

⟨m t m t ′ ⟩ A = ⟨m t ⟩ A ⟨m t ′ ⟩ A . (S30)
The property (S30) allows us to recast Eqs. (S9) and (S10) in the following form:

1 N T -τ M NT t=τM+1 ⟨m t ⟩ A = M E = 1 N T -τ M NT t=τM+1 m E t , (S31) 1 N T -τ M NT t=τM+1 ⟨m t ⟩ A • ⟨m t-τ ⟩ A = R E (τ ) = 1 N T -τ M NT t=τM+1 m E t • m E t-τ , 1 ≤ τ ≤ τ M . ( S32
)
These equations are simultaneously fulfilled by requiring that

⟨m t ⟩ A = m E t , τ M + 1 ≤ t ≤ N T . (S33)
In fact, given that Eq. (S33) is automatically satisfied for 1 ≤ t ≤ τ M , see Eqs. (S25) and (S26), if Eq. (S33) holds, then Eqs. (S31) and (S32) are satisfied.

We stress that Eq. ( S33) is only a sufficient condition for Eqs. (S31) and (S32) to hold and this condition introduces stricter constraints than those stemming from Eqs. (S31) and (S32). However, by moving from Eqs. (S31) and (S32) to Eq. (S33) we retain the very same parameters, that are H and J . Specifically, in this problem, the amount of data available N T × N is typically much larger than the number of parameters D + τ M . Thus, this stricter condition does not imply any overfitting (we are introducing constraints and not parameters) which, on the other hand, could be controlled by tuning τ M . By exploiting Eq. (S28), Equation (S33) can be rewritten as that

M τM τ =1 J τ m E t-τ + H = m E t , τ M + 1 ≤ t ≤ N T . (S34)
By noticing that M is a monotone function, its inverse function is well defined and we can rewrite Eq. (S34) as

G E t ≡ M -1 m E t , (S35) τM τ =1 J τ m E t-τ + H = G E t , τ M + 1 ≤ t ≤ N T (S36)
thus obtaining that the operator G E t is linear in the variables J , H: this last step will turn out to be crucial for the following developments. Let us underline that Eq. (S36) constitutes a linear system, whose solution provides us with the parameters J , H that we are looking for, but this system exhibits a number of constraints that is larger than the number of parameters. More precisely, the number of parameters is still D + τ M , and the number of constraints is (N T -τ M ) × D. This unbalance results from (S33) We can therefore enrich the model by stating that the relation (S36) does not hold deterministically and by introducing a source of noise represented by a standard Gaussian variable tuned by the additional parameter σ. This stochastic term turns the linear system Eq. (S36) into an auto-regressive model

1 of order τ M τM τ =1 J τ m E t-τ + H + σϵ t = G E t , τ M + 1 ≤ t ≤ N T , (S37)
ϵ t ∼ N (0 D , 1 D ) , σ ≥ 0. (S38)
In this setting we are saying that, given the parameters J , H and σ, the magnetization m t of model at time t depends on the magnetization exhibited in the previous τ M times according to the evolution rule:

m t = M H + τM τ =1 J τ m t-τ + σϵ t . (S39)
Now, the log-likelihood l A related to the auto-regressive model is

l A (J , H, σ|s E ) ≡ - 1 2σ 2 NT t=τM+1 τM τ =1 J τ m E t-τ + H -G E t 2 -(N T -τ M ) D 2 log 2πσ 2 . ( S40
)
and it is built by assuming that the ϵ t are standard Gaussians. In the next section we will handle l A (J , H, σ|s E ) to obtain an estimate for J , H, and σ.

S1.2.1 Analytical solution for the estimated parameters

As in maximum likelihood estimations, we evaluate the gradient of Eq. (S40) and set it to zero in order to recover the optimal estimates for J , H and σ.

Computing the derivatives of Eq. (S40), we obtain

∇ H l A = - N T -τ M σ 2 H + τM τ =1 J τ m E t-τ -G E t t , (S41)
∇ Jτ l A = - N T -τ M σ 2 H + τM λ=1 J λ m E t-λ -G E t • m E t-τ t , (S42)
∇ σ l A = N T -τ M σ 3 τM τ =1 J τ m E t-τ + H -G E t 2 t -(N T -τ M ) D σ , (S43)
where

⟨•⟩ t ≡ 1 N T -τ M NT t=τM+1 • (S44)
denotes the time average.

Setting Eq. (S41) to zero and solving for H, we get

H = ⟨G E t ⟩ t - τM τ =1 J τ m E t-τ t . (S45)
Proceeding along the same lines, from Eq. (S42), for J we obtain

τM τ =1 A λτ J τ = B τ ↔ J ≡ A -1 • B, (S46)
where the τ M × τ M matrix A λτ and the τ M -dimensional vector B λ are

A λτ ≡ m E t-τ • m E t-λ t -m E t-τ t • m E t-λ t , (S47)
B λ ≡ ⟨G E t • m E t-λ ⟩ t -⟨G E t ⟩ t • ⟨m E t-λ ⟩ t . (S48)
Finally, setting ∇ σ l A = 0, we obtain the variance σ 2 of the auto-regressive process (S39):

σ 2 = 1 D H + τM τ =1 J τ m E t-τ -G E t 2 t . (S49)

S1.2.2 Errors to inferred measurements

In order to estimate the errors for the inferred parameters in maximum likelihood estimation, we exploit the Fisher information for the log-likelihood (see e.g., [START_REF] Frieden | Science from Fisher Information: A Unification[END_REF]). We recall that, given some dataset X where each observation x i is assumed to be identically and independently distributed according to a true underlying distribution, and being f θ (x) a model probability density function parametrized by θ, we can write the log-likelihood function as

l(θ|x) = n i=1 log f θ (x i ),
where n is the sample size. Then, the (empirical) Fisher information matrix I has elements given by

I θa,θ b ≡ -exp[l(θ|x)] ∂ 2 l(θ|x) ∂θ a ∂θ b dx, (S50)
namely its elements correspond to the expectation of the elements of the Hessian matrix of the log-likelihood. Let us suppose that true parameter is θ 0 , and that the maximum-likelihood estimate of θ 0 is θ * = argmax θ l(θ|x). Then, one can prove that θ * ∼ N (θ 0 , I -1 (θ 0 )). Since when the sample size approaches infinity, the maximum-likelihood estimate approaches the true parameter (this is also known as the consistency property), we can write that the covariance matrix for the estimated parameters is just

Σ θ * ≈ n≫1 I -1 (θ * ). (S51)
We will use this result to provide the error estimates for our model. By Eq. (S40) the matrix elements of the Fisher information read

I Ha,H b = N T -τ M σ 2 δ ab , (S52)
I Jτ ,H = N T -τ M σ 2 m E t-τ t , (S53)
I Jτ ,J λ = N T -τ M σ 2 m E t-τ • m E t-λ t , (S54)
I σ,σ = N T -τ M σ 2 2D, (S55)
I σ,Jτ = I σ,H = 0, (S56)
and Eq. (S51) becomes

I = N T -τ M σ 2   m E t-τ • m E t-λ t m E t-τ t 0 m E t-λ t δ 0 0 0 2D   . (S57)
Finally, by substituting Eq. (S57) into Eq. (S51), we obtain the error estimates for J , H and σ:

Var (J , H, σ) = diag (Σ θ * ) . (S58)

S1.2.3 Effective number of delayed interactions

So far in our analysis, we assumed that the parameter τ M is a known constant. In most practical inference scenarios, where one needs to find from scratch the best model for a given dataset, the number of the effective delayed interactions τ M that must be considered is unknown: in what follows, we will provide a way to estimate it.

In order to tackle this problem, we rely on the Akaike Information Criterion (AIC) [START_REF] Banks | AIC under the framework of least squares estimation[END_REF], which is used in the literature to probe the flexibility of an auto-regressive model. Given a set of models for the data, the AIC estimates the quality of each model relative to all others, thus providing a tool for model selection. Namely, given a dataset and a model with a mean squared error ϵ 2 , n p parameters and n o observations, the quantity

AIC ≡ 2n p n o + log ϵ 2 (S59)
estimates the prediction error of the model and, consequently, the quality of such model for the data. Notice that, by definition, the AIC favours those models yielding a small error and a small number of parameters.

In the analysis above, the parameters in Eq. (S59) are

n o = D (N T -τ M ) , n p = τ M + D + 1, ϵ 2 = N T N T -τ M σ 2 , (S60) thus AIC = 2 τ M + D + 1 D (N T -τ M ) + log N T σ 2 N T -τ M . (S61)
After performing inference of Sections S1.2.1 and S1.2.2 for multiple values of τ M and determining σ 2 from Eq. (S49), we will evaluate Eq. (S61) for each of the inferred models, and we will choose the model with the smallest value of AIC as the best model.

S2 Model validation

In what follows, we test the inference framework of Section S1. We consider two different physical models, the Heisember-Kuramoto model (addressed in Sec. Eq. (S25)) and the Vicselk model (addressed in Sec. Sections S2.1 and S2.2), that we use to generate controlled synthetic datasets and a biological real-dataset collected via time-lapse microscopy on dendritic cell migration in a chemoattractant field we run the inference framework on such datasets and obtain the parameters J and H, which we will call inferred parameters and thus evaluate the goodness of the inference process. For each model validation, the model under consideration has its own set of parameters, which we will denote by the planted parameters. If the planted model is the first, that model is directly the time-delayed dynamical maximum entropy model itself (i.e. it is the direct formulation whose inverse problem is the algorithm developed in this paper), then inferred and planted parameters match hence they must numerically coincide (and their agreement indicates that the inference framework works, see Figure S1). If the model is instead the VM, there is no one-to-one correspondence between inferred and planted parameters, however, we can still check that the inferred model generates a new dataset statistically compatible with the planted dataset; futher, as trajectories described by the Vicsek model are essentially Markov chains (vide infra, see (S62), (S63)) we expect that the timewindow τ M that the protocol decides to explore to infer the couplings and fields is only one temporal step back, as confirmed in Figure S1). Finally for the biological datasets we run the delayed-inferenece protocol to obtain a remarkable result: the trajectory of a given cell at time t is influenced by those of the other cells for a past temporal window spanning back to τ M ∼ O(10) time-steps, thus resulting in highly non-Markovian dynamics (Figure S4). To inspect the validity of this result, once obtained the temporally ordered magnetizations, we reshuffled them to produce a new dataset of fake (Bernoullian) trajectories that we gave back to the protocol for obtaining in this case that solely the detection of the external field H (secreted by the chemoattractant) persisted, while both τ M = 0 as well as the coupling J = 0 (Figure S5) providing an indirect positive test on this new inference method.

S2.1 Synthetic dataset: the Heisemberg-Kuramoto model

To perform the validation of the ME model against itself, we will perform the following steps:

1. We set D = 2, τ M = 4, σ = 0.1, and draw the planted parameters J 1 , • • • , J τM , H x , H y independently and identically from U (-1, 1), where U (a, b) denotes the uniform probability distribution between a and b.

2. With the parameters just set, we simulate the ME model using Eq. (S39) from t = 0 to t = 5 × 10 3 , and store the magnetizations from t = 1 × 10 3 to t = 5 × 10 3 in order to consider equilibrated snapshots.

3. With the newly stored magnetizations we carry out the inference process according to the prescription explained in Section S1, obtaining the estimates: σ, J 1 , • • • , J τM , H x , H y , and τ M .

4. We store the corresponding pairs of inferred and planted parameters.

This process is repeated 150 times: with the collected data, the comparisons in Fig. S1 were made. In particular, Fig. S1a shows that the planted and inferred parameters are proportional along the diagonal x = y, implying an excellent agreement between them. In addition, Fig. S1d shows that the inference framework is able to detect the exact number of interactions most of the times. In some cases we notice that the number of interactions is overestimated, but we verified that the magnitude of the extra interactions is statistically compatible with zero within the error on the planted data. This can be seen from Fig. S1b: as J planted τ >τ planted ≡ 0, if J inferred τ >τ planted ̸ = 0, its magnitude corresponds to the absolute error and in these cases we reveal only small errors (∼ O(10 -2 )).

S2.2 Synthetic dataset: the Vicsek model

In this Section, we consider a topological 2 variant of the VM [START_REF] Vicsek | Novel type of phase transition in a system of self-driven particles[END_REF][START_REF] Ginelli | Relevance of metric-free interactions in flocking phenomena[END_REF]. Such model is defined by the following dynamical equations

v t+1 k = v 0 Γ 1 |n t k | j∈n t k v t j + H + √ 2T η t k , η t k ∼ N (0, 1) (S62) x t+1 k = x t k + ∆tv t+1 k , k = 1, • • • , N, (S63)
where v t i , x t i are the velocity and position of the ith particle at time t, respectively. Also, n c and n t 1 , • • • , n t nc are the total number of topological neighbors and the set of topological neighbors of the ith particle, including the i-th particle, at time t, respectively,

Γ(v) ≡ v |v| , (S64)
H is the external field applied to each particle, and T sets the strength of the dynamical noise. In this model, the particles are self-propelled, i.e., v t j = v 0 . Note that, as we discussed above, we denoted the planted parameter H with the same symbol as the magnetic field in the ME distribution (S15), but they need not be the same quantity. Finally, the magnetization of the VM is given by the expression

m t = 1 N v 0 N k=1 v t k , (S65)
which is analogous to Eq. (S4). The validation procedure for ME inference in the case of the VM is different from that of Section S2.1. In the case of the VM, we have T, H x as free parameters: These parameters are different from those of the ME model, therefore any comparison between planted and inferred parameters would be useless. Instead, a method that allows for validation in this case is to compare the lowest-order moments of some observables across the ME model and the VM. In particular, we will measure the magnetization M and the correlation function R(1) in both models, and check if their values are statistically compatible: If they do, the inference procedure will prove successful. Otherwise, with the particular configuration of the parameters (T, H x ), the procedure does not yeld infer a valid ME model. Concretely, the validation procedure was carried out by means of the following algorithm:

J 1 J 2 J 3 J 4 H x H y planted J 1 J 2 J 3 J 4 H x H y inferred c) -1 0
1. Draw H x according to U (0, 3) and T according to U (0, 0.5).

2. Run the dynamics of Eqs. (S62) and (S63) from t = 0 to t = 2 × 10 3 , storing the magnetization at each step (excluding the first 10 3 that are used for equilibration).

3. Perform the ME inference process, using the magnetizations of Item 2 as input.

4. From the parameters inferred in Item 3, simulate the ME model from t = 0 to t = 2 × 10 3 , collecting the magnetization at each step (excluding the first 10 3 that are used for equilibration).

5. For both models, calculate the moments M and R(1), their relative difference and the total number of interactions deduced in [START_REF] Loeffler | Understanding cell fate control by continuous single-cell quantification[END_REF].

The results of this analysis are shown in Fig. S2. Figure S2a, b and c yield a positive test of the inference method, indicating that the values of the correlations R( 1) and magnetization M from the inferred model match those of the planted model for most realisations of the planted parameters. Given that the coordinates at any time step in the VM depend only on the ones at the preceding step, we expect that only J(1) is significantly different from zero among the inferred couplings J(τ ) of the ME model: this prediction is mostly confirmed in Fig. S2d, where, apart from a small region in the parameter space (T, H x ), τ M is always equal to 1.

S2.3 Real dataset: dendritic cells in a chemokine gradient

Once validated the ME inference method on synthetic datasets, we now apply it as a proof of concept to an experimental dataset.

In this Section, we consider a two-dimensional tracking experiment where cells move in the presence of chemokines, which guide cell migration by acting as a chemoattractant, see Figure S3. We use datasets derived from [START_REF] Agliari | A statistical inference approach to reconstruct intercellular interactions in cell migration experiments[END_REF], where bone marrow-derived dendritic cells (BMDCs) were prepared as previously described in [START_REF] Sáez | Atp promotes the fast migration of dendritic cells through the activity of pannexin 1 channels and p2x7 receptors[END_REF]. Briefly, mouse bone marrow precursors were obtained from wild-type C57/B6 mice and were differentiated in vitro for 10 days with granulocyte-macrophage colony-stimulating factor-containing culture medium, obtained from transfected J558 cells. At day 10, cells were stimulated with a pulse bacterial lipopolysaccharide (LPS; 100 ng/ml), as previously described in [START_REF] Bretou | Lysosome signaling controls the migration of dendritic cells[END_REF]. The migration experiments in collagen gels were performed as previously described in [START_REF] Sáez | Leukocyte migration and deformation in collagen gels and microfabricated constrictions[END_REF]. LPSstimulated BMDCs were mixed with bovine collagen type I at 3 mg/ml and loaded in a custom-made chamber of polydimethylsiloxane. After 30 min of incubation at 37 • , gel polymerization was reached and samples were bathed with a medium containing CCL19 (200 ng/ml). CCL19 chemokines are not uniform distributed in the medium of the experiment, rather they form a gradient and the system under investigation has been observed across multiple instants of time, and the positions of all cells at each time were recorded. Data was collected for two regions, A (low values of chemokine gradient) and B (high values of chemokine gradient), thus yielding two distinct datasets.

We apply the inferential ME framework to datasets A and B. For the inferred ME model, we estimate δAIC(τ M ) ≡ AIC(τ M ) -min

τ ′ M AIC(τ ′ M ),
i.e., the difference between AIC for a given τ M and AIC for the optimal value of τ M , as well as the inferred parameters J(τ ), H x , H y and σ and the correlation function

C(τ ) ≡ m(s t ) • m(s t-τ ) t -m(s t ) t • m(s t-τ ) t , 1 ≤ τ ≤ τ M , (

Figure 1 :

 1 Figure 1: Test of the statistical-inference method, where the Heisemberg-Kuramoto model is chosen as planted model. Here D = 2, τ M = 4, σ = 0.1, and we drew 150 planted parameters J 1 , • • • , J τM , H x , H y independently and identically from a uniform distribution over the range (-1, 1). a) Scatter plot of planted vs inferred parameters: b) Histogram of the absolute error between J , H planted and J , H inferred: for each realisation of the process we calculate |J inferred -J planted | and |H inferred -H planted | and the related 150 × (D + τ M ) components make up the sample; the black dot represents the average and the black horizontal line represents a confidence interval of 68%; the average error is around ∼ 10 -2 , which is the expected amount of error given the size of the planted dataset. c) Pearson correlation matrix between planted and inferred parameters: as expected distinct kinds of parameters, e.g., H and J , are uncorrelated. d) Histogram of τ inferred M

Figure 2 :

 2 Figure 2: Validation of our inference method with the Vicsek model (VM) simulated in a periodic, twodimensional square lattice with size L = 150, with N = 100 particles and ∆t = 1, v 0 = 1, n c = 4, H = (H x , 0), where H x and T are drawn independently and identically from a uniform distribution over the range (0, 3) and (0, 0.5), respectively. The superscripts V and M E correspond to quantities calculated for the VM and the ME model, respectively. a) Observables M and R(1) as functions of H x , for the both models, with T = 0.25. b) Observables M and R(1) as functions of T , for both models, with H x = 1.375. c) Sum of the relative errors

Figure 3 :

 3 Figure 3: Inference for the chemokine-gradient dataset, for both regions A and B. a) Difference between the Akaike Information Criterion (AIC) and its minimal value, as a function of the total number of interactions τ M b) Inferred values of the delayed interaction J(τ ) and corresponding errors. c) Inferred values for the two components of the field, H x and H y , and inferred value of σ. d) Correlation function C(τ ), see Eq. (39), from the experimental data and from the ME model.

Figure S1 :

 S1 Figure S1: Test of the statistical-inference method: comparison between the (inferred) parameters returned by time-delayed dynamical maximum entropy protocol and the (planted) parameters of the Heisemberg-Kuramoto models used to generate the trajectories provided to the inferential machinery. a) Scatter plot of planted vs inferred parameters: b) Histogram of the absolute error between J , H planted and J , H inferred : for each realisation of the process we calculate |J inferred -J planted | and |H inferred -H planted | and the related 150 × (D + τ M ) components make up the sample; the black dot represents the average and the black horizontal line represents a confidence interval of 68%; the average error is around ∼ 10 -2 , which is the expected amount of error given the size of the planted dataset. c) Pearson correlation matrix between planted and inferred parameters: as expected distinct kinds of parameters (e.g. H and J ) are uncorrelated. d) Histogram of τ inferred M

Figure S2 :

 S2 Figure S2: Validation of the ME inference method with the Vicsek model (VM). The superscripts V and M E correspond to quantities calculated for the VM and the ME model, respectively. a) Observables M and R(1) as functions of H x , for the both models, with T = 0.25. b) Observables M and R(1) as functions of T , for both models, with H x = 1.375. c) Sum of the relative errors |M M E -M V | |M V |

  Supplementary Material for 'Inverse modeling of time-delayed interactions via the dynamic-entropy formalism' Elena Agliari a , Francesco Alemanno b , Adriano Barra c * , Michele Castellana d,e , Daniele Lotito c,f , Matthieu Piel g a Dipartimento di Matematica "Guido Castelnuovo", Sapienza Università di Roma, Roma, Italy b Dipartimento di Matematica, Università di Bologna, Bologna, Italy c Dipartimento di Matematica e Fisica, Università del Salento, Lecce, Italy d Institut Curie, PSL Research University, CNRS UMR168, 75005 Paris, France e Sorbonne Universités, Université Pierre et Marie Curie Université Paris 06, 75006 Paris, France

f Dipartimento di Informatica, Università di Pisa, Lungarno Antonio Pacinotti, 43, 56126 Pisa, Italy g Laboratoire de Biologie Cellulaire et Cancer, Institut Curie, CNRS UMR 168, 75005 Paris, France Contents S1 Dynamic maximum-entropy formalism 9 S1.1 Problem formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 S1.2 Solution method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 S1.2.1 Analytical solution for the estimated parameters . . . . . . . . . . . . . . . . . . . . . . . . . 13 S1.2.2 Errors to inferred measurements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 S1.2.3 Effective number of delayed interactions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 S2 Model validation 15 S2.1 Synthetic dataset: the Heisemberg-Kuramoto model . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 S2.2 Synthetic dataset: the Vicsek model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 S2.3 Real dataset: dendritic cells in a chemokine gradient . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

In an autoregression model, we forecast the variable of interest using a linear combination of past values of the variable; the term autoregression indicates that it is a regression of the variable against itself, see e.g.,[START_REF] Theodoridis | Machine Learning, A Bayesian and Optimization Perspective[END_REF].

The topological variant ensures that, at any time, each unit is always interacting with a certain number of neighbours, in contrast, in the metric variant, the number of neighbours is a random variable which strongly depend on the initial configuration and on the density of the system. The former is therefore preferred as, being more controllable, it is more suitable as a reference model for validation.

for the optimal value of τ M . These results are shown in Figure S4 for both datasets A and B. The outcomes obtained for the two datasets are consistent and the unique picture that emerges can be streamlined as follows: the dynamics of a given leukocyte at a certain time is influenced by the past actions of its peers for a long time window (in these cases τ M ∼ 13, 14, see panels a of Figure S4, one per dataset), hence the trajectories that these cells paint while migrating seem far from the Markovian limit. Along the dynamics, the history of the coupling is hetergeneous, specifically, it is (mainly) negative in the short past (τ = 1, 2) and (mainly) positive in the long past (τ ≥ 3) while it gets irrelevant -in the sense of the AIC criterion -for τ ≥ 13, 14 as shown in panels b of Figure S4. Also, as expected, the inferential procedure correctly returns a (chemoattractant) field larger in the region B than in the region A (i.e., H x ∼ 0.27 in region A and H x ∼ 0.46 in region B) as shown in panels c of Figure S4, while in panels d we report the connected correlation function that, up to τ M ∼ 13, 14 is small but systematically different from zero, while it gets zero for larger values of τ . To test a posteriori the soundness of the output of our inference protocol on these real datasets, we abandon the collected temporally ordered magnetizations of the various frames, and we reshuffle them to produce two new randomized datasets (one per region) that we used as input for our protocol: results, reported in Figure S5, correctly return the same external field, panels c, as this is a one-point correlation information that is not destroyed by the permutation of frames, while the optimal temporal delay to infer the coupling is now τ M = 0, panel a, that is again correct as it captures the Bernoulli coin we generate by reshuffling and, consistently, a null value for the coupling, panel b, associated to a null connected correlation function, panel d. Finally, we would be overall tempted to speculate that a cell, during migration, integrates multiple signals in time to decide its movement. Yet, in the present setting, with τ M ∼ O(15) and an acquisition time of two minutes per frame, this implies that cells would sense for a time window ∼ 30 minutes, that is a rather large time lapse: several technical experimental reasons may lie behind this result, ranging from chemokine gradient temporal instability to collagen gel local heterogeneity or the way the chemokine affects the cell or binds on the collagen. Clearly, a systematic study has to deepen these biological aspects, but this will be faced in a forthcoming paper.

Dataset [A] Dataset [B]

Dataset [A] Dataset [B]

Figure S5: Test on the validity of the delayed inference procedure: we reshuffled trajectories in the chemokinegradient dataset, for both regions A and B, and gave this random dataset back to the inferential protocol to inspect its outcomes. a) Difference between the Akaike Information Criterion (AIC) and its minimal value, as a function of the total number of interactions τ M : note that this time the optimal τ M is zero, in agreement with a Bernoullian process. b) Inferred values of the delayed interaction J(τ ) and corresponding errors: note that, as expected, in this case there is no inferred coupling. c) Inferred values for the two components of the field, H x and H y , and inferred value of σ. d) Correlation function C(τ ), see Eq. (S66), from the experimental data and from the ME model. As expected, while permutation invariance over the frames destroyed the temporal correlation of the motion (induced by the couplings) it did not erase the one-point gradient of the chemoattractant, that is still well inferred.


Adrien Blanchet

Philippe Laurenc ¸ot

FINITE MASS SELF-SIMILAR BLOWING-UP SOLUTIONS OF A CHEMOTAXIS SYSTEM WITH NON-LINEAR DIFFUSION

Keywords: . 2000 Mathematics Subject Classification. Primary: 35K65, 34C10, Secondary: 92B99 Backward self-similar solutions, blowup, chemotaxis, Patlak-Keller-Segel model, degenerate diffusion

For a specific choice of the diffusion, the parabolic-elliptic Patlak-Keller-Segel system with non-linear diffusion (also referred to as the quasi-linear Smoluchowski-Poisson equation) exhibits an interesting threshold phenomenon: there is a critical mass M c > 0 such that all the solutions with initial data of mass smaller or equal to M c exist globally while the solution blows up in finite time for a large class of initial data with mass greater than M c . Unlike in space dimension 2, finite mass self-similar blowing-up solutions are shown to exist in space dimension d ≥ 3.

Introduction

In space dimension d = 2, the parabolic-elliptic Patlak-Keller-Segel (PKS) system is a simplified model which describes the collective motion of cells in the following situation: cells diffuse in space and emit a chemical signal, the chemo-attractant, which results in the cells attracting each other. If ρ denotes the density of cells and c the concentration of the chemo-attractant, the PKS system reads [START_REF] Keller | Initiation of slide mold aggregation viewed as an instability[END_REF][START_REF] Patlak | Random walk with persistence and external bias[END_REF] (1)

   ∂ t ρ(t, x) = div [∇ρ(t, x) -ρ(t, x)∇c(t, x)] , c(t, x) = (E 2 ⋆ ρ)(t, x) , E 2 (x) = - 1 2π ln |x| , (t, x) ∈ [0, ∞) × R 2 .
This model may be seen as an elementary brick to understand the aggregation of cells in mathematical biology as it exhibits the following interesting and biologically relevant feature: there is a critical mass above which the density of cells is expected to concentrate near isolated points after a finite time, a property which is related to the formation of fruiting bodies in the slime mold Dictyostelium discoideum. Such a phenomenon does not take place if the density of cells is too low. More precisely, given a non-negative integrable initial condition ρ 0 with finite second moment, the system (1) has a unique maximal classical solution (ρ, c) defined on some maximal time interval [0, T ), T ∈ (0, ∞]. Its first component ρ is non-negative and the mass of ρ (that is, its L 1 -norm) remains constant through time evolution

ρ(t) 1 = M := ρ 0 1 , t ∈ [0, T ) .
It is well-known that, if M < 8 π, the solution to (1) exists globally in time while it blows up in finite time if M > 8 π, see [START_REF] Blanchet | Two-dimensional Keller-Segel model: optimal critical mass and qualitative properties of the solutions[END_REF][START_REF] Dolbeault | Optimal critical mass in the two-dimensional Keller-Segel model in R 2[END_REF][START_REF] Horstmann | From 1970 until present: the Keller-Segel model in chemotaxis and its consequences. I[END_REF][START_REF] Jäger | On explosions of solutions to a system of partial differential equations modelling chemotaxis[END_REF] and the references therein.

More recently, it was shown that there is global existence as well for the critical mass M = 8 π, the blowup occurring in infinite time with a profile being a Dirac mass of mass 8π [START_REF] Blanchet | Infinite time aggregation for the critical two-dimensional Patlak-Keller-Segel model[END_REF]. When the mass M is above 8π, the shape of the finite time blowup is not self-similar according to asymptotic expansions computed in [START_REF] Chavanis | Virial theorem and dynamical evolution of self-gravitating Brownian particles in an unbounded domain. I. Overdamped models[END_REF][START_REF] Lushnikov | Critical chemotactic collapse[END_REF] (see also [START_REF] Herrero | Singularity patterns in a chemotaxis model[END_REF] for a related problem in a bounded domain). In addition, there is no integrable and radially symmetric blowing-up self-similar solution to (1) [START_REF] Naito | Self-similarity in chemotaxis systems[END_REF]Theorem 8].

In space dimension d ≥ 3, the system (1) seems to be less relevant from the biological point of view as blowup may occur whatever the value of M [START_REF] Herrero | Finite-time aggregation into a single point in a reaction-diffusion system[END_REF][START_REF] Nagai | Behavior of solutions to a parabolic-elliptic system modelling chemotaxis[END_REF]. This means that the diffusion is too weak to balance the aggregation resulting from the chemotactic term. It is however well-known that one can enhance the effect of diffusion to prevent crowding by considering a diffusion of porous medium type which increases the diffusion of the cells when their density ρ is large. This is the generalised version of the Patlak-Keller-Segel model considered in, e.g., [START_REF] Blanchet | Critical mass for a Patlak-Keller-Segel model with degenerate diffusion in higher dimensions[END_REF][START_REF] Chavanis | Anomalous diffusion and collapse of self-gravitating Langevin particles in D dimensions[END_REF][START_REF] Sugiyama | Global existence in sub-critical cases and finite time blow-up in super-critical cases to degenerate Keller-Segel systems[END_REF][START_REF] Sugiyama | Application of the best constant of the Sobolev inequality to degenerate Keller-Segel models[END_REF][START_REF] Suzuki | Degenerate parabolic equation with critical exponent derived from the kinetic theory, I, generation of the weak solution[END_REF]:

(2) ∂ t ρ(t, x) = div (∇ [ρ m (t, x)] -ρ(t, x)∇c(t, x)) , c(t, x) = (E d ⋆ ρ)(t, x) , E d (x) = c d |x| 2-d , (t, x) ∈ [0, ∞) × R d , where m > 1, c d := 1/((d -2) σ d ), and
σ d := 2 π d/2 /Γ(d/2
) denotes the surface area of the sphere S d-1 of R d . The system (2) also arises in astrophysics [START_REF] Chavanis | Anomalous diffusion and collapse of self-gravitating Langevin particles in D dimensions[END_REF] (being then referred to as the generalised Smoluchowski-Poisson equation), and ρ and c denote the density of particles and the gravitational potential, respectively. For (2), it turns out that there is only one critical exponent of the non-linear diffusion, namely m d := 2(d -1)/d, such that the mass plays a similar role to that in [START_REF] Blanchet | Infinite time aggregation for the critical two-dimensional Patlak-Keller-Segel model[END_REF]. Indeed, if m > m d the diffusion enhancement is too strong and the solutions always exist globally in time whereas if m < m d the diffusion is not strong enough to compensate the aggregation term and there are solutions blowing up in finite time whatever the value of the mass [START_REF] Sugiyama | Global existence in sub-critical cases and finite time blow-up in super-critical cases to degenerate Keller-Segel systems[END_REF][START_REF] Sugiyama | Application of the best constant of the Sobolev inequality to degenerate Keller-Segel models[END_REF]. The relevant diffusion is thus achieved in the case when m = m d . In this case, it was proved in [START_REF] Blanchet | Critical mass for a Patlak-Keller-Segel model with degenerate diffusion in higher dimensions[END_REF] that there is a unique threshold mass M c > 0 with the following properties: if the mass M = ρ 0 1 of the initial condition ρ 0 is less or equal to M c , then the corresponding solution to (2) exists globally in time, whereas given any M > M c there are initial data ρ 0 with mass M such that the corresponding solution blows up in finite time. Thus, for the peculiar choice m = m d and d ≥ 3, the system (2) exhibits the same qualitative behaviour as the PKS system (1) in space dimension 2. Still, there is a fundamental difference as the latter has no fast-decaying stationary solution with mass 8π while the former has a two-parameter family of non-negative, integrable, and compactly supported stationary solutions with mass M c for each d ≥ 3 [2, Section 3].

It is then tempting to figure out whether this striking difference extends above the critical mass M c and this leads us to investigate the existence of blowing-up (or backward) self-similar solutions with finite mass. More precisely, since mass remains unchanged throughout time evolution, we look for solutions (ρ, c) to (2) with m = m d and d ≥ 3 of the form

(3) ρ(t, x) = 1 s(t) d Φ x s(t) and c(t, x) = 1 s(t) d-2 Ψ x s(t) for (t, x) ∈ [0, T ) × R d with s(t) := [d(T -t)] 1/d ,
the time T being an arbitrary positive real number. Note that s(t) converges to zero as t increases to the blowup time T .

Our main result is then the following:

Theorem 1 (Existence of finite mass self-similar blowing-up solutions). There exists M 2 ∈ (M c , ∞) such that, for any M in (M c , M 2 ], there exists at least a nonnegative self-similar blowing-up solution (ρ M , c M ) to (2) of the form (3) with a radially symmetric, compactly supported, and non-increasing profile

Φ M satisfying ρ M (t) 1 = Φ M 1 = M for t ∈ [0, T ) and ρ M (t) ∞ → ∞ as t → T .
As a consequence of Theorem 1, we realize that non-negative, integrable, and radially symmetric self-similar blowing-up solutions to (2) with a non-increasing profile only exist below a threshold mass. Another by-product of our analysis is the existence of non-negative and non-integrable self-similar blowing-up solutions to [START_REF] Blanchet | Critical mass for a Patlak-Keller-Segel model with degenerate diffusion in higher dimensions[END_REF], see Proposition 8 below.

Blowing-up self-similar profiles

From now on,

d ≥ 3 and m = m d = 2(d -1) d ,
and we look for a solution (ρ, c) to (2) of the form

(4) ρ(t, x) = 1 s(t) d Φ x s(t) and c(t, x) = 1 s(t) d-2 Ψ x s(t)
with s(t) = [d(Tt)] 1/d and (t, x) ∈ [0, T ) × R d for some given T > 0. We further assume that Φ enjoys the following properties:

   Φ ∈ C(R d ) ∩ L 1 (R d ) is radially symmetric and non-negative, Φ m-1 ∈ W 1,∞ (R d ) .
Inserting the ansatz (4) in [START_REF] Blanchet | Critical mass for a Patlak-Keller-Segel model with degenerate diffusion in higher dimensions[END_REF] gives that (Φ, Ψ) solves div (∇ [Φ m (y)] -Φ(y)∇Ψ(y) -Φ(y) y) = 0 ,

Ψ(y) = (E d ⋆ Φ)(y) , for y ∈ R d . Since Ψ = E d ⋆ Φ
, the radial symmetry of Φ ensures that of Ψ and, introducing the profiles (ϕ, ψ) of (Φ, Ψ)

(5) Φ(y) = ϕ(|y|) , Ψ(y) = ψ(|y|) , y ∈ R d .
By [14, Theorem 9.7, Formula ( 5)], we have ( 6)

ψ(r) = 1 (d -2)r d-2 r 0 ϕ(s) s d-1 ds + 1 d -2 ∞ r ϕ(s) s ds
for r ≥ 0. We can also write the equation for ϕ as

(7) ∂ r r d-1 ϕ(r) ∂ r J(r) = 0 with J(r) := 2(d -1) d -2 ϕ (d-2)/d (r) -ψ(r) - r 2 2
for r ∈ (0, ∞). Since we are looking for an integrable profile, we formally conclude that

(8) ∂ r J(r) = 0 for r ∈ P ϕ := {s ∈ (0, ∞) : ϕ(s) > 0} .
In particular, J is constant on any connected component of P ϕ . But, if C is a connected component of P ϕ , we have either

(9) C = (0, R s ) for some R s ∈ (0, ∞] , or C = (R i , R s ) for some R i ∈ (0, ∞) and R s ∈ (0, ∞] .
Remark 2. If we additionally assume that the profile ϕ is non-increasing then P ϕ has only one connected component which is necessarily of the form (9). Now, take a connected component C of P ϕ . It follows from ( 8) that there is µ ∈ R such that ( 10)

J(r) = 2(d -1) d -2 ϕ (d-2)/d (r) -ψ(r) - r 2 2 = -µ for r ∈ C .
Owing to the assumed integrability of Φ, the function r → r d-1 ϕ(r) belongs to L 1 (0, ∞) and it follows from ( 6) that the function r → r d-2 ψ(r) is bounded in C. Therefore [START_REF] Herrero | Singularity patterns in a chemotaxis model[END_REF] only complies with the integrability of Φ if R s < ∞ which implies the boundedness of C. Introducing Ξ := ϕ (d-2)/d and taking the Laplacian of both sides of (10) yield that Ξ is a positive solution to [START_REF] Horstmann | From 1970 until present: the Keller-Segel model in chemotaxis and its consequences. I[END_REF] -

d 2 Ξ dr 2 (r) - d -1 r dΞ dr (r) = d -2 2(d -1) Ξ(r) d/(d-2) -d in C , with either (12) ∂ r Ξ(0) = Ξ(R s ) = 0 if C = (0, R s ) or (13) Ξ(R i ) = Ξ(R s ) = 0 if C = (R i , R s ) .
A final change of scale, namely

η(r) := 1 λ d Ξ r µ d , λ d := d (d-2)/d , µ d := d 1/d d -2 2(d -1) 1/2
, leads us to the following boundary-value problem for η: either ( 14)

       d 2 η dr 2 (r) + d -1 r dη dr (r) + η(r) d/(d-2) -1 = 0 , r ∈ (0, µ d R s ) , dη dr (0) = 0 , η(µ d R s ) = 0 , or (15)
     d 2 η dr 2 (r) + d -1 r dη dr (r) + η(r) d/(d-2) -1 = 0 , r ∈ (µ d R i , µ d R s ) , η(µ d R i ) = 0 , η(µ d R s ) = 0 .
We have thus reduced our study to one or several boundary-value problems (depending on the number of connected components of P ϕ ) for a nonlinear second order differential equation. The purpose of the next section is then a precise study of this ordinary differential equation. However, before going on, let us point out that [START_REF] Horstmann | From 1970 until present: the Keller-Segel model in chemotaxis and its consequences. I[END_REF] is not equivalent to [START_REF] Herrero | Singularity patterns in a chemotaxis model[END_REF]. Indeed, since [START_REF] Felmer | Monotonicity properties for ground states of the scalar field equation[END_REF], the fact that Ξ is a solution to [START_REF] Horstmann | From 1970 until present: the Keller-Segel model in chemotaxis and its consequences. I[END_REF] only guarantees that ∂ r (r d-1 ∂ r J(r)) = 0 for r ∈ C. Consequently, there are constants C 1 and C 2 such that

∂ r J(r) = 2(d -1) d -2 ∂ r Ξ(r) + 1 r d-1 r 0 Ξ(s) d/(d-2) s d-1 ds -r , r ∈ C , by
∂ r J(r) = - (d -2) C 1 r d-1 , J(r) = C 1 r d-2 + C 2 , r ∈ C , from which (10) follows only if C 1 = 0. On the one hand, if C = (R i , R s ) with 0 < R i < R s ,
it is yet unclear whether the boundary conditions [START_REF] Keller | Initiation of slide mold aggregation viewed as an instability[END_REF] might imply this property. On the other hand, if C = (0, R s ), the boundary conditions [START_REF] Jäger | On explosions of solutions to a system of partial differential equations modelling chemotaxis[END_REF] ensure that ∂ r J(0) = 0 and thus C 1 = 0. We shall only deal with this case in the remaining of this paper and thus focus on the non-increasing profiles ϕ.

An auxiliary ordinary differential equation

For a ∈ R, let u(., a) ∈ C 1 ([0, r max (a))) denote the maximal solution to the Cauchy problem ( 16)

     u ′′ (r, a) + d -1 r u ′ (r, a) + |u(r, a)| p-1 u(r, a) -1 = 0 , r ∈ [0, r max (a)) , u(0, a) = a , u ′ (0, a) = 0 , with r max (a) ∈ (0, ∞] and p = d/(d -2).
Clearly, if a = 1 then u(., 1) ≡ 1 is a stationary solution and r max (1) = ∞. We first show that u(., a) is global for all a ∈ R and oscillates around the value 1 if a = 1. Lemma 3. For each a ∈ R \ {1}, r max (a) = ∞, and the solution u(., a) to ( 16) is an oscillatory function in (0, ∞). More precisely,

• if a > 1, there is an increasing sequence (r i (a)) i≥0 of real numbers such that r 0 (a) = 0,

   u ′ (r i (a), a) = 0 , (-1) i u ′ (r, a) < 0 for r ∈ (r i (a), r i+1 (a)) , u(r 2i (a), a) > u(r 2i+2 (a), a) > 1 > u(r 2i+3 (a), a) > u(r 2i+1 (a), a) for i ≥ 0, • if a < 1, there is an increasing sequence (r i (a)) i≥1 of real numbers such that r 1 (a) = 0    u ′ (r i (a), a) = 0 , (-1) i u ′ (r, a) < 0 for r ∈ (r i (a), r i+1 (a)) , u(r 2i (a), a) > u(r 2i+2 (a), a) > 1 > u(r 2i+1 (a), a) > u(r 2i-1 (a), a) for i ≥ 1.
These properties are illustrated in Figure 1. Notice that, for a = 7, u(., 7) vanishes at a finite r and thus provides a solution to [START_REF] Lieb | Analysis[END_REF]. By [START_REF] Kwong | Uniqueness of positive solutions of ∆uu + u p = 0 in R N[END_REF], for all r ∈ [0, r max (a))

(18) dE dr (r, a) = - d -1 r |u ′ (r, a)| 2 ≤ 0 ,
Obviously E(r, a) ≥ -p/(p + 1). Owing to [START_REF] Naito | Self-similarity in chemotaxis systems[END_REF], E(r, a) ∈ [-p/(p + 1), E(0, a)] for r ∈ [0, r max (a)) which prevents u(., a) of becoming unbounded at a finite value of r, thereby implying that r max (a) = ∞. We next argue using Sturm's oscillations theorem as in [16, Lemma 9], to establish the oscillatory behaviour of u(., a) for a = 1.

According to [START_REF] Lieb | Analysis[END_REF], we are interested in finding solutions to the initial value problem [START_REF] Kwong | Uniqueness of positive solutions of ∆uu + u p = 0 in R N[END_REF] which are positive and vanish at a finite value of r. We thus focus on the case a > 0 and investigate the positivity properties of u(., a). Lemma 4. There is a constant a c > 1 such that

• if a ∈ (0, a c ), then u(r, a) > 0 for all r ≥ 0,

• if a = a c , then there is R(a c ) > 0 such that

     u(R(a c ), a c ) = 0 u ′ (R(a c ), a c ) = 0 u(r, a c ) > 0 for r ∈ [0, R(a c )), • if a ∈ (a c , ∞), then there is R(a) > 0 such that      u(R(a), a) = 0 u ′ (R(a), a) < 0 u(r, a) > 0 for r ∈ [0, R(a)).
These three possibilities are drawn in Figure 2. Clearly, P ∪ N ∪ N 0 = (0, ∞) and 1 ∈ P. Actually, if a ∈ (0, (p + 1) 1/p ), then E(0, a) < 0 and the monotonicity ( 18) of E entails that E(r, a) < 0 for all r ≥ 0. But, if R(a) < ∞, it readily follows from the definition (17) of the functional E that E(R(a), a) ≥ 0 whence a contradiction. Therefore, R(a) = ∞ for any a ∈ (0, (p + 1) 1/p ) so that (19) 0, (p + 1) 1/p ⊂ P .

Consider now a ∈ N 0 . Then U(x) := u(|x|, a) is a radial positive solution to the homogeneous Dirichlet-Neumann free boundary problem ∆U + U p -1 = 0 in B(0, R(a)) with U = ∂ ν U = 0 on ∂B(0, R(a)). According to [20, Theorem 3 (iii)], there is only one value of a for which this solution has a positive radial solution and it is unique. Consequently, there is a unique a c > 0 such that N 0 = {a c }.

Consider next a ∈ N ∪N 0 and recall that a > 1 by [START_REF] Patlak | Random walk with persistence and external bias[END_REF]. Following [START_REF] Kwong | Uniqueness of positive solutions of ∆uu + u p = 0 in R N[END_REF]Lemma 11], let us assume for contradiction that there is ̺ ∈ (0, R(a)) such that u ′ (̺, a) = 0. Either ̺ ≤ 1 and we infer from the definition, the monotonicity of E, see ( 17)- [START_REF] Naito | Self-similarity in chemotaxis systems[END_REF], and the definition of R(a) that 0 > E(̺, a) ≥ E(R(a), a) ≥ 0 which is a contradiction. Or ρ > 1 and the oscillating behaviour of the solutions implies, using the notation of Lemma 3, that ̺ ≥ r 2 (a). This implies that r 1 (a) < R(a). Then u(r 1 (a), a) ∈ (0, 1) and using again ( 17), [START_REF] Naito | Self-similarity in chemotaxis systems[END_REF], and the definition of R(a), we conclude that 0 > E(r 1 (a), a) ≥ E(R(a), a) ≥ 0, hence a contradiction. Therefore, [START_REF] Serrin | Uniqueness of ground states for quasilinear elliptic equations[END_REF] u ′ (r, a) < 0 for r ∈ (0, R(a)) if a ∈ N ∪ N 0 .

Let us now prove that P and N are open subsets of (0, ∞). We first consider a ∈ N : by [START_REF] Serrin | Uniqueness of ground states for quasilinear elliptic equations[END_REF] there are ̺ > R(a) and ε > 0 such that u(̺, a) < 0 and u ′ (r, a) < -2ε for r ∈ (0, ̺). By continuous dependence, there is δ ∈ (0, a) such that u(̺, b) < 0 and u ′ (r, b) < -ε for r ∈ (0, ̺) and b ∈ (aδ, a + δ). Since u(0, b) = b > 0, we readily deduce that, for each b ∈ (aδ, a + δ), we have R(b) ∈ (0, ̺) with u ′ (R(b), b) < -ε < 0. Consequently, (aδ, a + δ) and N is open in (0, ∞). Consider next a ∈ P, a > 1. By Lemma 3 and (17), we have u(r, a) ≥ u(r 1 (a), a) ∈ (0, 1) for r ∈ [0, r 1 (a)] and E(r 1 (a), a) < 0. By continuous dependence, there is δ > a such that u(r, b) ≥ u(r 1 (a), a)/2 > 0 for r ∈ [0, r 1 (a)], u(r 1 (a), b) ∈ (0, 1), and E(r 1 (a), b) < 0 for b ∈ (aδ, a + δ). Assume now for contradiction that there is b ∈ (aδ, a + δ) such that R(b) < ∞. Owing to [START_REF] Nagai | Behavior of solutions to a parabolic-elliptic system modelling chemotaxis[END_REF], [START_REF] Naito | Self-similarity in chemotaxis systems[END_REF], and the definition of R(b), we obtain 0 > E(r 1 (a), b) > E(R(b), b) ≥ 0 and a contradiction. Consequently, (aδ, a + δ) ⊂ P and P is also open in (0, ∞).

We finally argue as in [START_REF] Kwong | Uniqueness of positive solutions of ∆uu + u p = 0 in R N[END_REF]Lemma 15] to show that there is A > 0 such that (A, ∞) ⊂ N .

Since P and N are open subsets of (0, ∞), N 0 = {a c }, (0, (p + 1) 1/p ) ⊂ P, and (A, ∞) ⊂ N , we readily conclude that P = (0, a c ) and N = (a c , ∞).

We next study the properties of the map a → R(a). An efficient tool for that purpose is the variation of u(., a) with respect to a defined by

ϑ(r, a) := ∂u ∂a (r, a) , (r, a) ∈ [0, ∞) × (0, ∞),
which solves the second order linear differential equation ( 21)

ϑ ′′ (r, a) + d -1 r ϑ ′ (r, a) + p u(r, a) p-1 ϑ(r, a) = 0 , r ∈ [0, ∞) , ϑ(0, a) = 1 , ϑ ′ (0, a) = 0 ,
We argue as in [START_REF] Felmer | Monotonicity properties for ground states of the scalar field equation[END_REF][START_REF] Tang | Uniqueness of positive radial solutions for ∆uu + u p = 0 on an annulus[END_REF] to prove the following lemma. In addition, u(z(a), a) > 1 and the ratio ϑ(., a)/u(., a) is a decreasing function of r on (0, R(a)).

Proof of Lemma 5. Since the proof follows rather closely that of [START_REF] Tang | Uniqueness of positive radial solutions for ∆uu + u p = 0 on an annulus[END_REF] and [7, Lemma 2.1], we sketch it briefly for the sake of completeness. Fix a > a c and set u = u(., a) and ϑ = ϑ(., a) to simplify notations. We first argue as in [START_REF] Kwong | Uniqueness of positive solutions of ∆uu + u p = 0 in R N[END_REF]Lemma 17] to show that ϑ vanishes at least once in the interval (0, z 1 (a)), where z 1 (a) denotes the unique zero in (0, R(a)) of u -1. Indeed, ( 16) also reads

(u(r) -1) ′′ + d -1 r (u(r) -1) ′ + u(r) p -1 u(r) -1 (u(r) -1) = 0 , r ∈ [0, ∞)
and (u(r) p -1)/(u(r) -1) ≤ p u(r) p-1 for r ∈ [0, z 1 (a)). It then follows from Sturm's comparison theorem that ϑ vanishes at least once in the interval (0, z 1 (a)). Let z ∈ (0, z 1 (a)) denote the first zero of ϑ.

We now aim at showing that ϑ cannot vanish once more in the interval (z, R(a)). To this end, we define

ξ(r) := r d-1 [u ′ (r) ϑ(r) -u(r) ϑ ′ (r)] = -r d-1 u(r) 2 ϑ u ′ (r) , r ∈ [0, R(a)) ,
which encodes the monotonicity of ϑ/u. It follows from ( 16) and ( 21) that

(22) ξ ′ (r) = r d-1 ((p -1) u p (r) + 1) ϑ(r) , r ∈ [0, R(a)) .
Clearly, ξ ′ (r) > 0 for r ∈ (0, z) and ξ(0) = 0, so that ξ(r) > 0 for r ∈ (0, z]. Assume now for contradiction that there is ̺ ∈ (z, R(a)) such that ξ(r) > 0 for r ∈ (0, ̺) and ξ(̺) = 0 .

Observing that ϑ ′ (z) < 0, we realize that, if ϑ(̺) ≥ 0, there is σ ∈ (z, ̺] such that ϑ(r) < 0 for r ∈ (z, σ) and ϑ(σ) = 0. In that case, ϑ ′ (σ) ≥ 0 and thus ξ(σ) = -σ d-1 u(σ) ϑ ′ (σ) ≤ 0, leading us to a contradiction. Consequently, [START_REF] Sugiyama | Application of the best constant of the Sobolev inequality to degenerate Keller-Segel models[END_REF] ϑ(̺) < 0 .

We next introduce the functions

T (r) := 2 (u(r) p -1) (p -1) u(r) p + 1 ξ(r) -ζ(r) , ζ(r) := r d [u ′ (r) ϑ ′ (r) + (u(r) p -1) ϑ(r)] + (d -2) r d-1 u ′ (r) ϑ(r) ,
for r ∈ [0, R(a)) and use ( 16), [START_REF] Slepčev | Selfsimilar blowup of unstable thin-film equations[END_REF], and ( 22) to obtain

ζ ′ (r) = 2 r d-1 (u(r) p -1) ϑ(r) , T ′ (r) = 2p 2 u(r) p-1 [(p -1) u(r) p + 1] 2 u ′ (r) ξ(r) , (24)
for r ∈ [0, R(a)). Integrating ( 24) over (0, ̺) and using the negativity of u ′ and the positivity of ξ on this interval give [START_REF] Tang | Uniqueness of positive radial solutions for ∆uu + u p = 0 on an annulus[END_REF] ζ(̺) = -T (̺) > 0 .

Since ξ(̺) = 0, we have u(̺) ϑ ′ (̺) = u ′ (̺) ϑ(̺) and we have

ζ(̺) = Q(̺) ϑ(̺) u(̺) ,
where

Q(r) := r d u ′ (r) 2 + u(r) p+1 -u(r) + (d -2) r d-1 u(r) u ′ (r) , r ∈ [0, R(a)) .
It then follows from ( 23), [START_REF] Tang | Uniqueness of positive radial solutions for ∆uu + u p = 0 on an annulus[END_REF], and the positivity of u that (26) Q(̺) < 0 .

Finally, define

P (r) := r d u ′ (r) 2 + 2 u(r) p+1 p + 1 -2 u(r) + (d -2) r d-1 u(r) u ′ (r)
for r ∈ [0, R(a)). On the one hand, we notice that (27)

P (r) = Q(r) -u(r) - p -1 p + 1 u(r) p+1 < Q(r) , r ∈ [0, R(a)) .
On the other hand, we deduce from ( 16) and ( 18) that

P ′ (r) = r d-1 u(r) d -2 d -1 u(r) p -(d + 2) , r ∈ [0, R(a)) .
At this point, we realize that we have necessarily a > (d + 2)(d -1)/(d -2) and that there is s ∈ (0, R(a)) such that P ′ (r) > 0 if r ∈ (0, s) and P ′ (r) < 0 if r ∈ (s, R(a)).

Since P (0) = 0 and P (R(a)) > 0, we conclude that P (̺) > 0 and then Q(̺) > 0 by (27). But this contradicts (26). We have thus established that ξ is positive in (0, R(a)) from which Lemma 5 follows.

We are now in a position to state and prove some properties of the map a → R(a).

Proposition 6. The map a → R(a) is a decreasing function on (a c , ∞) and there is

z 1 > 0 such that (28) lim aցac R(a) = R(a c ) and lim a→∞ a (p-1)/2 R(a) = z 1 .
The monotonicity of a → R(a) is shown in Figure 3. According to numerical simulations, the function a → a (p-1)/2 R(a) also seems to be a decreasing function of a ∈ [a c , ∞), see Figure 3.

     w ′′ (r) + d -1 r w ′ (r) + |w(r)| p-1 w(r) = 0 , r ∈ [0, ∞) , w(0) = 1 , w ′ (0) = 0 .
By [START_REF] Fowler | Further studies of Emden's and similar differential equations[END_REF], there is z 1 > 0 such that (33) w(r) > 0 and w ′ (r) < 0 for r ∈ [0, z 1 ) , w(z 1 ) = 0 , w ′ (z 1 ) < 0 .

Owing to (33), there is δ > 0 such that w(r) < 0 for r ∈ (z 1 , z 1 + δ). It then follows from (31) that, given r ∈ (z 1 , z 1 + δ), v(r, a) < 0 for a large enough (depending on r), whence a (p-1)/2 R(a) ≤ r for a large enough by (30). Letting r → z 1 guarantees that lim sup a→∞ a (p-1)/2 R(a) ≤ z 1 .

Next, if ̺ ∈ (0, γ), we have w(r) > w(̺) > 0 for r ∈ [0, ̺] and we infer from (31) that v(r, a) > w(̺)/2 > 0 for r ∈ [0, ̺] and a large enough. Consequently, ̺ < a (p-1)/2 R(a) for a large enough, from which we conclude that lim inf a→∞ a (p-1)/2 R(a) ≥ z 1 .

Combining the above two inequalities completes the proof of Proposition 6.

The above information allow us to estimate from above and from below a specific integral of u(., a). The upper bound (35) is then a straightforward consequence of (28) and the above inequality.

Next, recalling that v(., a) is defined by (29), we have

M(a) = d |B(0, 1)| a (p-1)/2 R(a) 0 v(r, a) p r d-1 dr ,
and we infer from ( 28) and ( 31) that (34) holds true.

Proof of Theorem 1

Thanks to the analysis done in the previous sections, we are now in a position to construct self-similar blowing-up solutions to (2) having either finite or infinite mass. Proposition 8. Given a > 0 and T > 0, define M 2 .

ϕ(r) := λ d/(d-2) d u(µ d r, a) d/(d-2) for r ∈ [0, ∞) if a ∈ (0, a c ) , and
ϕ(r) :=    λ d/(d-2) d u(µ d r, a) d/(d-2) for r ∈ [0, R(a)/µ d ] 0 for r ≥ R(a)/µ d , if a ∈ [a c , ∞) .

Discussion

We have proved the existence of non-negative, integrable, and radially symmetric self-similar blowing-up solutions for [START_REF] Blanchet | Critical mass for a Patlak-Keller-Segel model with degenerate diffusion in higher dimensions[END_REF]. The profile ϕ of these self-similar solutions is compactly supported and non-increasing, and the mass of the corresponding selfsimilar solution ranges in the bounded interval (M c , M 2 ], the threshold mass M c corresponding to the onset of blowup found in [START_REF] Blanchet | Critical mass for a Patlak-Keller-Segel model with degenerate diffusion in higher dimensions[END_REF]. Our analysis thus reveals the existence of a second threshold value M 2 > M c of the mass above which no radially symmetric and non-increasing self-similar blowing-up solution exist. The meaning of this second threshold value for the mass is yet unclear. It is worth mentioning at this point that a related situation was uncovered for the critical unstable thin-film equation

∂ t u = -∂ x u n ∂ 3 x u + u n+2 ∂ x u , (t, x) ∈ [0, ∞) × R ,
in [START_REF] Slepčev | Selfsimilar blowup of unstable thin-film equations[END_REF] for n ∈ (0, 3/2). It is likely that, given M ∈ (M c , M 2 ], there is only a unique radially symmetric and non-increasing self-similar blowing-up solution with mass M and Figure 4 provides some numerical evidence of this fact. Besides this uniqueness question, the question of stability of these blowing-up solutions is also of interest. Another challenging question is the existence (or non-existence) of integrable profiles ϕ with a non-connected positivity set as discussed in Section 2. Figure 5 provides numerical evidence that, if a > a c is large enough, u(., a) may have several zeroes and each positive "hump" actually corresponds to a solution of (15) for suitable values of R i and R s . Whether the additional constraint (10) may be satisfied does not seem to be clear.

Figure 1 .

 1 Figure 1. Various oscillating behaviours of u(., a) for a ∈ {0.2, 1, 3, 7}.

Figure 2 .

 2 Figure 2. Behaviour of u(., a) for a > a c , a = a c and a < a c .

Lemma 5 .

 5 If a > a c , there is a unique z(a) ∈ (0, R(a)) such that , a) > 0 for r ∈ [0, z(a)) , ϑ(z(a), a) = 0 ϑ(r, a) < 0 for r ∈ (z(a), R(a)] .

Figure 3 .v

 3 Figure 3. Monotonicity of the radius R and a → a (p-1)/2 R(a) (d = 3).

Proposition 7 . 1 0

 71 For a ∈ [a c , ∞), we define M(a) := d |B(0, 1)| R(a) 0 u(r, a) p r d-1 dr .Recalling that w is the solution to (32) and z 1 is its first positive zero, we havelim a→∞ M(a) = M c := d |B(0, 1)| z w(r) p r d-1 dr ,(34)M 2 := sup a∈[ac,∞) M(a) < ∞ . (35) Proof of Proposition 7. Let a ≥ a c . Since u(0, a) = a, it follows from the monotonicity of u(., a) that M(a) ≤ d |B(0, 1)| R(a) 0 a p r d-1 dr = |B(0, 1)| a (p-1)/2 R(a) d .

2 (d- 1 )/ 2 M 2 (

 2122 Define next ψ, Φ, and Ψ by (6) and (5), respectively. Then the functions (ρ, c) defined by (4) in (0, T ) × R d with s(t) = [d(Tt)] 1/d is a non-negative self-similar blowing-up solution to (2) with finite mass if a ≥ a c and infinite mass if a ∈ (0, a c ). The proof of Proposition 8 readily follows from the analysis performed in Sections 2 and 3. As for Theorem 1, it is a straightforward consequence of Proposition 8, the threshold values M c and M 2 being given by M c := d 1/d 2(d -1) d -c and M 2 := d 1/d 2(d -1) d -

Figure 4 .

 4 Figure 4. Monotonicity of the mass a → M(a).

Figure 5 .

 5 Figure 5. Positivity set of u(., a) with two (a = 50, left) and three (a =

GREMAQ, CNRS UMR 5604, INRA UMR 1291, Université de Toulouse,

Allée de Brienne, F-31000 Toulouse, France E-mail address: Adrien.Blanchet@univ-tlse1.fr2 Institut de Mathématiques de Toulouse, CNRS UMR 5219, Université de Toulouse, F-31062 Toulouse cedex 9, France E-mail address: laurenco@math.univ-tlse.fr

Acknowledgements

We thank José Antonio Carrillo, Jean Dolbeault, and Dejan Slepčev for stimulating discussions. Part of this work was done while the authors enjoyed the hospitality and support of the Centro de Ciencias Pedro Pascuale de Benasque. c 2008 by the authors. This paper is under the Creative Commons licence Attribution-NonCommercial-ShareAlike 2.5.

Partially supported by the project EVaMEF ANR-09-JCJC-0096-01.