Dan Pineau

Fran¸cois Orieux

Alain Abergel

Fast Multispectral and Hyperspectral Image Fusion via Hessian Inversion

Keywords: Inverse problems, data fusion, MM algorithm, pansharpening, multispectral, hyperspectral, variant blur, deconvolution, close form, infrared, JWST

Multispectral and hyperspectral data fusion allows data restoration with increased spatial and spectral resolutions. A common approach, in particular with the presence of variant blur, is to solve an ill-posed inverse problem by minimizing a mixed criterion. This minimization usually requires an iterative gradient-based method. Still, this paper demonstrates the existence of a reachable explicit solution without the need to solve a Sylvester equation, allowing for a three order of magnitude timesaving factor wrt. iterative algorithms. We exploit this explicit resolution to accelerate the minimization by a Majorize-Minimize algorithm of a semi-quadratic convex edge-preserving criterion. We conduct experiments on realistic synthetic measurements for the James Webb Space Telescope and show that our proposed solution outperforms the state-of-the-art in computation time and quality of reconstruction.

I. INTRODUCTION

M ANY fields, such as remote sensing, astrophysics, or Earth observation, use hyperspectral data for spectra observation. Several applications allow the use of integral field spectrographs to acquire spectra in a given field of view, but the images produced are generally poorly sampled. Parallel acquisition of high spatial resolution and well-sampled data, such as multispectral data, opens the field of data fusion to reconstruct highly spatially and spectrally resolved data [START_REF] Yokoya | Hyperspectral and multispectral data fusion: A comparative review of the recent literature[END_REF].

Usually, multispectral (MS) and hyperspectral (HS) datasets are represented as datacubes, with two spatial dimensions and one spectral dimension, that have been spatially or spectrally degraded by the instruments' effects like transmission of optics, filters, sampling of the detector, etc. Therefore, the MS/HS fusion problem consists of exploiting the information provided by these two degraded datasets to reconstruct a hyperspectral datacube corrected from instrumental effects.

Different approaches were explored in the literature to solve the MS/HS fusion problem like matrix factorization and tensor-based approaches [START_REF] Yokoya | Coupled nonnegative matrix factorization unmixing for hyperspectral and multispectral data fusion[END_REF]- [START_REF] Borsoi | Coupled tensor decomposition for hyperspectral and multispectral image fusion with inter-image variability[END_REF], deep neural networks approaches [START_REF] Palsson | Multispectral and hyperspectral image fusion using a 3-d-convolutional neural network[END_REF]- [START_REF] Hu | Fusformer: A transformer-based fusion network for hyperspectral image super-resolution[END_REF], or bayesian based approaches [START_REF] Zhang | Noise-resistant waveletbased bayesian fusion of multispectral and hyperspectral images[END_REF], [START_REF] Simoes | A convex formulation for hyperspectral image superresolution via subspace-based regularization[END_REF].

Nonetheless, some cases of MS/HS fusion problems can require using computationally burdening algorithms, especially when fusing high-dimensional datasets. The development of computationally efficient solutions to this fusion problem was attempted in the literature. In Wei et al. [START_REF] Wei | Fast fusion of multi-band images based on solving a sylvester equation[END_REF], the MS/HS data fusion problem is formulated through a Sylvester equation, and an explicit solution is obtained, allowing a significant reduction of the computational cost, unlike iterative approaches. Methods based on Tucker decompositions proposed by Prévost et. al. [START_REF] Prévost | Hyperspectral superresolution with coupled tucker approximation: Recoverability and svdbased algorithms[END_REF] and Borsoi et. al. [START_REF] Borsoi | Coupled tensor decomposition for hyperspectral and multispectral image fusion with inter-image variability[END_REF] also provide explicit solutions to the MS/HS fusion problem.

However, the problem tackled in this paper involves the fusion of datasets where instrumental effects include significant and spectrally varying spatial blurs, acquired on a relatively wide mid-infrared waveband for instance. The previously presented methods [START_REF] Wei | Fast fusion of multi-band images based on solving a sylvester equation[END_REF] or [START_REF] Borsoi | Coupled tensor decomposition for hyperspectral and multispectral image fusion with inter-image variability[END_REF] did not take into account such spectrally dependent effects, making them inadequate for the problem of this paper.

In contrast, Guilloteau et. al. [START_REF] Guilloteau | Hyperspectral and multispectral image fusion under spectrally varying spatial blurs-application to high dimensional infrared astronomical imaging[END_REF] proposes the solution of an inverse problem with a spectrally non-stationary spatial blur by minimizing a least squares criterion including data adequacy terms and a quadratic regularization. However, the minimization of this criterion required using an iterative algorithm instead of an explicit solution.

This paper approaches the MS/HS fusion problem with variational approaches, as done in [START_REF] Guilloteau | Hyperspectral and multispectral image fusion under spectrally varying spatial blurs-application to high dimensional infrared astronomical imaging[END_REF]. The first contribution is the demonstration that a close form and explicit solution of a regularized least square is reachable for the same inverse problem [START_REF] Guilloteau | Hyperspectral and multispectral image fusion under spectrally varying spatial blurs-application to high dimensional infrared astronomical imaging[END_REF], instead of the previously proposed iterative approach, and without involving the resolution of a Sylvester equation like in [START_REF] Wei | Fast fusion of multi-band images based on solving a sylvester equation[END_REF]. This proposed procedure allows for a significant time-saving factor of up to 7 000 compared to commonly used iterative approaches. A second contribution is the development of a fast convex edge-preserving solution for the MS/HS fusion problem. This method uses a semi-quadratic regularization [START_REF] Abi-Rizk | Super-resolution hyperspectral reconstruction with majorization-minimization algorithm and low-rank approximation[END_REF] and a Majorize-Minimize algorithm [START_REF] Chouzenoux | A majorize-minimize strategy for subspace optimization applied to image restoration[END_REF] to improve the quality of the reconstruction of objects featuring high spatial gradients in comparison to [START_REF] Guilloteau | Hyperspectral and multispectral image fusion under spectrally varying spatial blurs-application to high dimensional infrared astronomical imaging[END_REF] (NRMSE from 27 × 10 -3 to 22 × 10 -3 ). These advances exploit a new efficient way to calculate, store, and apply the inverse of a Hessian matrix by exploiting its structure of diagonal blocks, later described in section III.

Section II describes the model used to represent the object observed by the instruments and the data formation processes or forward models: one for an imaging instrument, and another for an integral field spectrograph. Section III formulates the problem and details solutions for two different regularizations: quadratic and convex differentiable. Finally, section IV presents an application of our method on simulated astronomical data acquired by instruments similar to MIRI, the mid-infrared instrument of the James Webb Space Telescope (JWST). The results are presented in terms of reconstruction quality and time savings against iterative approaches and deep neural network methods from the literature.

II. MODELS DESCRIPTION

A. Observed object model

The observed sky region is described as a discrete hyperspectral cube x with two spatial coordinates i ∈ [1, I] and j ∈ [1, J], and a spectral coordinate l ∈ [1, L]. As proposed in [START_REF] Guilloteau | Hyperspectral and multispectral image fusion under spectrally varying spatial blurs-application to high dimensional infrared astronomical imaging[END_REF], [START_REF] Abi-Rizk | Super-resolution hyperspectral reconstruction with majorization-minimization algorithm and low-rank approximation[END_REF], [START_REF] Hadj-Youcef | Spatio spectral reconstruction from low resolution multispectral data: application to the mid-infrared instrument of the james webb space telescope[END_REF], spectral correlations are introduced within the observed object x, such that the spectral content of any location (i, j) is a linear combination of T known spectra s t . The observed object x thus writes as

x[i, j, l] = T t=1 a t [i, j] s t [l] (1)
where a t [i, j] gives the abundance of the spectrum s t at (i, j). This assumption of a Linear Mixing Model offers several advantages. As the spectra s t can be defined for all common wavelengths covered by the instruments, they enable efficient joint processing of all available data. Therefore, the number of unknowns describing the object becomes P × T instead of P × L, where P = I × J is the number of pixels in an image. Thus, choosing T ≪ L significantly reduces the number of unknowns. Eq. 1 can be rewritten in a matrix form as

x = T a (2)
where x ∈ R LP ×P is the object with stacked monochromatic images, T ∈ R LP ×T P and a ∈ R T P ×P the stacked abundance images. The linear mixing matrix T is written as

T =     s 1 [1]I P . . . s T [1]I P . . . . . . s 1 [L]I P . . . s T [L]I P     (3)
with I P the identity matrix.

B. Imager model

Imagers contain several components, such as mirrors to focus light, filters to select the waveband to observe, and a detector to sample the light spatially. Due to the diffraction on the finite aperture of the instrument, the spectrally varying response blurs the images. Moreover, fields of view are considered relatively small, so the blurring effect can be modeled as a spatial convolution with spectrally varying impulse responses.

It is assumed that each wideband image used for the multispectral dataset can be produced from different imagers with different blurring responses. Let C be the total number of imagers. The blurred object x c hm writes

x c hm [i, j, l] = (x * i,j h c m )[i, j, l]
where * i,j

is the 2D spatial convolution operator on the axes of i and j with h c m the spectrally varying impulse response of

imager c ∈ [1, C].
Then, the spectral filter profiles and the detector's quantum efficiency of the imager must be taken into account. These effects write as a spectral weighting w c m with

x c w [i, j, l] = x c hm [i, j, l] w c m [l]. (4)
Finally, wideband images result from the integration of the incoming light on a given waveband

x c m [i, j] = l x c w [i, j, l].
Thus, transformations of the observed object x to obtain the collection of C wideband images x m can be written in a matrix form as

x m = W m C m x (5)
where C m ∈ R CLP ×LP writes as

C m =                     C 1 m,1 . . . C 1 m,L     . . .     C C m,1 . . . C C m,L                     (6)
and where C c m,l ∈ R P ×P is a convolution matrix with the impulse response h c m at the wavelength l. The matrix W m ∈ R CP ×CLP operates the spectral response and spectral integration for each of the C imagers, and is written as

W m = W 1 m . . . W C m (7)
where

W c m ∈ R CP ×LP writes W c m =           0 P . . . 0 P . . . . . . w c m [1]I P . . . w c m [L]I P . . . . . . 0 P . . . 0 P           (8)
and where the row of non-zero matrices is the row c. Using the Linear Mixing Model for Eq. 5 yields

x m = W m C m T a (9a) = M a (9b)
where M ∈ R CP ×T P is the matrix form of the imager forward model, with each block calculated as weighted sums of convolution matrices

M =     l w 1 m [l]s 1 [l]C 1 m,l . . . l w 1 m [l]s T [l]C 1 m,l . . . . . . l w C m [l]s 1 [l]C C m,l . . . l w C m [l]s T [l]C C m,l     . (10)

C. Spectrometer model

Common design choices exist between an integral field spectrograph and an imager, such as mirrors to focus light and detectors to sample the flux. Hypotheses do not differ from the imager's blur, as it varies spectrally but is spatially invariant. Therefore, the blurred observed object x hh is a convolution between the optical impulse response h h of the instrument and x as

x hh [i, j, l] = (x * i,j h h )[i, j, l].
Then, the diffraction grating does not transmit flux for all wavelengths equally. This effect, combined with the quantum efficiency of the detector, leads to the spectral weighting of the spectrometer w h , which must be taken into account as

x w [i, j, l] = x hh [i, j, l]w h [l].
The last effect is the detector spatial sampling of the incoming light with a spatial sampling step greater than the one used for the imager. It is computed as a spatial integration of every region of d i ×d j pixels of the object x w , d i , d j ∈ N * . Consequently, the resulting images are spatially decimated from P to P ′ = P/(d i d j ) pixels per image, and each pixel value is calculated as

x h [ ī, j, l] = ( ī+1)di i= īdi ( j+1)dj j= jdj x w [i, j, l] (11)
where ī and j are the new pixels coordinates in the spectrometer images.

Like the imager's model, the formation of images from the spectrometer model can be written in a matrix form as

x h = SΣW h C h x (12)
where the matrix C h ∈ R LP ×LP contains L convolution matrices C h,l ∈ R P ×P on its diagonal, applying the spectrally variant blur h h for each wavelength. The matrix W h ∈ R LP ×LP operates the spectral response of the spectrometer, written as

W h =     w h [1]I P . . . w h [L]I P     . ( 13
)
The spatial integration of the spectrometer's detector, described in Eq. 11, is modeled by combining two matrix operators: a summing operator Σ ∈ R LP ×LP and a subsampling operator S ∈ R LP ′ ×LP . The first operator replaces the value of each pixel of every image with the sum of the region of d i × d j pixels surrounding them and is written as

Σ = I L ⊗ Σ ( 14
)
where ⊗ stands for the Kronecker product, and Σ ∈ R P ×P is a convolution matrix with a 2D kernel of size d i × d j and full of ones.

The second operator then selects one pixel for every region of d i × d j pixels, hence a reduction of the image size from P to P ′ pixels, and writes as

S = I L ⊗ S ( 15
)
where each matrix S ∈ R P ′ ×P subsamples one image of the input cube. Using Eq. 3 of the Linear Mixing Model allows to write Eq. 12 as

x h = SΣW h C h T a (16a) = Ha (16b)
where H ∈ R LP ′ ×T P is the matrix form of the spectrometer forward model and is written as

H = SΣW h C h T =     w h [1]s 1 [1]SC Σ,1 . . . w h [1]s T [1]SC Σ,1 . . . . . . w h [L]s 1 [L]SC Σ,l . . . w h [L]s T [L]SC Σ,l     (17)
with C Σ,l = ΣC h,l the joint convolution operator for the wavelength l.

III. METHODOLOGY

This section formulates the inverse problem for the reconstruction of the observed object described with abundance maps a. Then, fast approaches are detailed to solve the problem for two regularization cases.

A. Formulation of the problem

The abundance maps describing the observed object are defined as the minimizer of a criterion a = arg min a J (a) [START_REF] Hunt | A matrix theory proof of the discrete convolution theorem[END_REF] composed of two data adequacy terms and a regularization term as

J (a) = µ m ∥y m -M a∥ 2 2 + µ h ∥y h -Ha∥ 2 2 + R(a) (19
) where µ m = 1/2σ 2 m , µ h = 1/2σ 2 h , and σ m and σ h are the standard deviations of additive white Gaussian noises. Because of the spatial convolution, a regularization R(a), detailed later, is used to stabilize the solution. An explicit minimizer of this criteria with a quadratic regularization identical to [START_REF] Guilloteau | Hyperspectral and multispectral image fusion under spectrally varying spatial blurs-application to high dimensional infrared astronomical imaging[END_REF] is given in section III-B, and a fast minimization procedure for a convex semi-quadratic regularization is described in section III-E.

B. Quadratic regularization

A formulation of Eq. 19 with a quadratic regularization (ℓ 2norm), as proposed in Guilloteau et al. [START_REF] Guilloteau | Hyperspectral and multispectral image fusion under spectrally varying spatial blurs-application to high dimensional infrared astronomical imaging[END_REF], writes as

J ℓ2 (a) = µ m ∥y m -M a∥ 2 2 + µ h ∥y h -Ha∥ 2 2 + µ r ∥D r a∥ 2 2 + ∥D c a∥ 2 2 ( 20
)
where

D r = I T ⊗ D r , (21a)
D c = I T ⊗ D c , (21b)
and D r , D c ∈ R T P ×T P are first order 2D difference matrices along the rows and columns respectively.

Minimizing J ℓ2 is equivalent to solving ∇J ℓ2 = 0, which yields

a = Q -1 ℓ2 q ℓ2 (22)
with

Q ℓ2 = µ m Q M + µ h Q H + µ r (Q Dr + Q Dc ), (23)
q ℓ2 = µ m M H y m + µ h H H y h , (24)
where • H is the conjugate transpose operator, Q ℓ2 ∈ R T P ×T P , q ℓ2 ∈ R T P ×P , and

Q M = M H M (25)
Q H = H H H (26)
Q Dr = D H r D r ( 27
)
Q Dc = D H c D c (28)
In [START_REF] Guilloteau | Hyperspectral and multispectral image fusion under spectrally varying spatial blurs-application to high dimensional infrared astronomical imaging[END_REF], authors propose a gradient based descent algorithm to compute â (Eq. 22). However, the minimizer a of J ℓ2 can be calculated explicitly if the inversion of Q ℓ2 is accessible. Section III-D demonstrates that the structure of Q ℓ2 , detailed in III-C, can be exploited for its fast inversion, that Q -1 ℓ2 is easily storable and its application to q ℓ2 is also fast.

C. Structure of Q ℓ2

Lemma 1: If convolution matrices C * and D * are circulant, then the matrices Q M , Q H , and Q D * are square with circulant blocks of equal size. Consequently, Q ℓ2 is also made of circulant blocks.

1) Structure of Q M : As shown in Eq. 10, the matrix M contains C ×T blocks M c,t that are sums of circulant matrices

M c,t = l w c m [l]s t [l]C m,l . (29)
As a result, the blocks M c,t are circulant themselves. Then, the matrix Q M is a square matrix of T × T blocks that writes

(Q M ) t,t ′ = c M H c,t M c,t ′ (30)
that are also circulant [START_REF] Hunt | A matrix theory proof of the discrete convolution theorem[END_REF].

2) Structure of Q H : From Eq. 17, the L × T blocks of H writes

H l,t = w h [l]s t [l]SC Σ,l (31)
Therefore, Q H is a square matrix of T × T blocks that writes as

(Q H ) t,t ′ = l w h [l] 2 s t [l]s t ′ [l]C H Σ,l S H SC Σ,l . (32)
However, the structure of these blocks in the Fourier space leads to exploitable matrix structures. Indeed, writing

C Σ,l = F H Λ l F gives C H Σ,l S H SC Σ,l = F H Λ H l F S H SF H Λ l F . ( 33
)
As S carries a subsampling while S H fills images with zeros, the resulting operation of S H S is a spatial weighing. More precisely, the matrix S H S is diagonal with 1 every d i d j coefficients and zeros elsewhere, i.e. the diagonal is a Dirac comb of step d i d j . Therefore, F S H SF H is a circulant convolution matrix with a Dirac comb as kernel, and as presented in Wei et al. [START_REF] Wei | Fast fusion of multi-band images based on solving a sylvester equation[END_REF], writes

F S H SF H = 1 d i d j J didj ⊗ I P ′ (34)
where J didj ∈ R didj ×didj is full of "1". Multiplying Eq. 34 by Λ l on the right gives

F S H SF H Λ l = 1 d i d j     I P ′ . . . I P ′ . . . . . . I P ′ . . . I P ′         Λ l,1 . . . Λ l,didj     = 1 d i d j     Λ l,1 . . . Λ l,didj . . . . . . Λ l,1 . . . Λ l,didj     (35)
where

Λ l,d ′ , d ′ ∈ [1, d i d j ]
, are subdivisions of the diagonal of the initial matrix Λ l . Multiplying Eq. 35 by Λ H l on the left gives

Λ H l F S H SF H Λ l = 1 d i d j     Λ H l,1 Λ l,1 . . . Λ H l,1 Λ l,didj . . . . . . Λ H l,didj Λ l,1 . . . Λ H l,didj Λ l,didj     . (36)
As each block of Λ H l F S H SF H Λ l is diagonal, the matrices C H Σ,l S H SC Σ,l defined at Eq. 33 are thus made circulant blocks, which turns also true for the blocks of Q H (Eq. 32).

3) Structures of Q D * : Using Eqs. 21a and 21b gives expressions for Q Dr and Q Dc as

Q Dr = I T ⊗ D H r D r and (37)
Q Dc = I T ⊗ D H c D c . (38)
Since D * are circulant, D H * D * are circulant and Q D * are block matrices with circulant blocks.

D. Inversion of

Q ℓ2 Proposition 1: If ker Q M ∩ ker Q H ∩ ker Q Dr ∩ ker Q Dc = ∅ then Q -1
ℓ2 exists and is a block matrix with circulant blocks. Moreover, thanks to Fourier transform, Q -1 ℓ2 is tractable, storable in Fourier space and applicable. The proof is described Appendix A.

E. Semi-quadratic regularization

The quadratic penalization of the previous regularization may turn inadequate when restoring high spatial gradients in images and can induce ringing artifacts. To solve this issue, non-quadratic regularizations can be convenient alternatives [START_REF] Idier | Convex half-quadratic criteria and interacting auxiliary variables for image restoration[END_REF], [START_REF] Geman | Nonlinear image recovery with half-quadratic regularization[END_REF]. A common choice is to use a strictly differentiable convex function φ within the regularization function R, such that φ is

[16] 1) φ is C 1 , even, coercive, 2) ∀u ∈ R + , φ( √ u) is concave, 3) ∀u ∈ R, 0 < φ(u)/u < +∞.
Therefore, the new convex criterion to minimize writes as

J ℓ2,1 (a) = µ m ∥y m -M a∥ 2 2 + µ h ∥y h -Ha∥ 2 2 + µ r t,i,j φ(D r a t )[i, j] + φ(D c a t )[i, j] (39)
As proposed by Geman and Yang (GY) [START_REF] Geman | Nonlinear image recovery with half-quadratic regularization[END_REF] and used in [START_REF] Hadj-Youcef | Spatio spectral reconstruction from low resolution multispectral data: application to the mid-infrared instrument of the james webb space telescope[END_REF], the half-quadratic strategy aims to minimize J ℓ2,1 (a) by defining an augmented criterion

K(a, b r , b c ) such that min br,bc K(a, b r , b c ) = J ℓ2,1 (a) (40)
where b r and b c are auxiliary variables. Therefore, under light hypothesis [START_REF] Idier | Convex half-quadratic criteria and interacting auxiliary variables for image restoration[END_REF] K is convex and a relaxed version of J with common minimizer, i.e.

â, br , bc = arg min a,br,bc

K and â = arg min a J ℓ2,1 . (41)
This convex [START_REF] Idier | Convex half-quadratic criteria and interacting auxiliary variables for image restoration[END_REF] surrogate criterion writes as

K(a, b r , b c ) = µ m ∥y m -M a∥ 2 2 + µ h ∥y h -Ha∥ 2 2 + µ r ∥D r a -b r ∥ 2 2 + ζ(b r ) + ∥D c a -b c ∥ 2 2 +ζ(b c ) (42)
where

ζ(u) = i ζ(u i )
, with ζ linked to φ through convex duality [START_REF] Idier | Convex half-quadratic criteria and interacting auxiliary variables for image restoration[END_REF]. The benefit of this strategy is to transform a highdimensional convex minimization with correlation problem into two separate minimization problems: a quadratic one with correlation (minimizing K with b r and b c fixed), and a separable one wrt. to the auxiliary variables (minimizing K with a fixed), as

       a (k) = arg min a K a, b (k) r , b (k) c , b (k+1) r , b (k+1) c = arg min br,bc K a (k) , b r , b c . (43) (44)
The alternative resolution of these minimization problems leads to finding the triplet ( a, b r , b c ) that minimizes K, and according to Eq. 41, a will minimize the original criterion J ℓ2,1 .

Solving the first problem (Eq. 43) yields

a = Q -1 ℓ2,1 q ℓ2,1 (45)
where

Q ℓ2,1 = µ m Q M + µ h Q H + µ r Q Dr + Q Dc (46) q ℓ2,1 = µ m M H y m + µ h H H y h + µ r D H r b r + D H c b c (47)
and where expressions of Q M , Q H , Q Dr and Q Dc are found from Eq. 25 to 28. Noticing that Q ℓ2,1 = Q ℓ2 , where the latter matrix was defined in section III-B for the resolution of the quadratic criterion J ℓ2 , leads to invert Q ℓ2,1 as it was done for Q ℓ2 in section III-D. As Q -1 ℓ2,1 does not depend on a, b r or b c , it can only be calculated once for a given problem, at the start of the alternative minimization procedure. As a consequence, the solution to Eq. 43 is explicit and fast.

The second minimization problem Eq. 44, used to update the auxiliary variables b r and b c , can be reformulated in two distinct problems with

         b r = arg min br 1 2 ∥D r a -b r ∥ 2 2 + ζ(b r ), b c = arg min bc 1 2 ∥D c a -b c ∥ 2 2 + ζ(b c ). ( 48
) (49)
Cancelling the gradients of both of these criteria gives formulations for the minimizers [START_REF] Hadj-Youcef | Spatio spectral reconstruction from low resolution multispectral data: application to the mid-infrared instrument of the james webb space telescope[END_REF] as

b r = D r a - 1 2 φ ′ D r a , ( 50
) b c = D c a - 1 2 φ ′ D c a (51)
where φ ′ independently applies φ ′ , the first derivative of φ, to every element of the input vector.

As illustrated, both minimization problems of Eqs. 43 and 44 are solved explicitly. Consequently, this procedure is expected to converge quickly in many test cases. This algorithm is described in Alg. 1, where the stopping criterion can be chosen to exploit the least square distance between two iterates, the value of the criterion J ℓ2,1 of Eq. 39, or the number of iterations.

Algorithm 1: Convex ℓ 2,1 HS-MS fusion

Input : µ m , y m , M , µ h , y h , H, µ r , D r , D c , φ ′ P ← Q -1 ℓ2,1 ; # See Eq. (46) and Appendix A ȳ ← µ m M H y m + µ h H H y h ; while Stopping criterion is not met do b r ← D r a -1 2 φ ′ D r a b c ← D c a -1 2 φ ′ D c a a ← P ȳ + µ r D H r b r + D H c b c end Output: a IV. EXPERIMENTAL RESULTS
The proposed method for the fusion of multispectral and hyperspectral data can be applied to a relatively wide variety of cases. Indeed, the impulse response can differ between the imager and the spectrometer model. Moreover, the wideband images obtained in the multispectral dataset can originate from different imagers (with different impulse responses and spectral responses) as long as their fields of view and spatial sampling steps are the same. Note that the fusion with imagery data with different sampling spatial steps or the fusion with more than one hyperspectral datacubes is possible by adding more data attachment terms in the originally formulated criterion of Eq. 19, but this was not explored in this paper.

In this section, the proposed algorithms' performances were evaluated by simulating the fusion of MS and HS data from the mid-infrared instrument MIRI of the James Webb Space Telescope [START_REF] Stsci | JWST documentation website[END_REF]. Instrument models described in sections II-B and II-C have been adapted for the imager MIRIM [START_REF] Bouchet | The mid-infrared instrument for the james webb space telescope, iii: Mirim, the miri imager[END_REF] and the Mid Resolution Spectrometer (MRS) [START_REF] Wells | The mid-infrared instrument for the james webb space telescope, vi: The medium resolution spectrometer[END_REF] belonging to MIRI.

The instrument models chosen for this study are described in section IV-A and the synthetic observation datasets in section IV-B. The performances of data fusion approaches are compared with reconstruction quality metrics, explained in section IV-C. Results obtained with the quadratic and semiquadratic regularizations, introduced in sections III-B and III-E, are presented in terms of time benefits and quality of reconstruction in sections IV-D and IV-E. Finally, inverse problem approaches are compared with a neural network method of the literature in section IV-F.

A. Study case instrument models

The PSF of the JWST, obtained with webbpsf [START_REF] Perrin | Simulating point spread functions for the james webb space telescope with webbpsf[END_REF], was used to simulate the impulse response used in the imager and spectrometer models of sections II-B and II-C. The FWHM of this PSF is proportionally linear to the wavelength, as illustrated in Fig. 1. The spectral response of the imager uses the real Photon Conversion Efficiency (PCE) curves of MIRIM [START_REF] Stsci | JWST documentation website[END_REF], where C = 9 different spectral channels are considered [START_REF] Bouchet | The mid-infrared instrument for the james webb space telescope, iii: Mirim, the miri imager[END_REF]. The spectral response induced by the MRS is composed of 4 channels with 3 sub-channels each [START_REF] Wells | The mid-infrared instrument for the james webb space telescope, vi: The medium resolution spectrometer[END_REF]. However, a simplified spectrometer model is considered with one unique channel and no sub-channel, with a spectral response taken to be the maximum of the concatenated PCE curves of the MRS [START_REF] Stsci | JWST documentation website[END_REF].

This model applies a spatial decimation with d i = d j = 4, i.e. integrates square region of 4 × 4 pixels into one, which induces a spatial subsampling for wavelength shorter than 20 µm [START_REF] Stsci | JWST documentation website[END_REF].

B. Synthetic data

Typical MS and HS data that MIRI/JWST could acquire were created from a set of T = 5 abundance maps a derived from real observations [START_REF] Guilloteau | Simulated jwst data sets for multispectral and hyperspectral image fusion[END_REF], [START_REF] Habart | High angular resolution near-ir view of the Orion bar revealed by keck/nirc2[END_REF], shown Fig. 3. These maps are associated with 5 spectra s [START_REF] Team | Pdrs4all: A jwst early release science program on radiative feedback from massive stars[END_REF] of length L = 300 covering the MIRI spectral range (from 4.6 to 28.6 µm), shown Fig. 2. The maps a and spectra s have been used as first inputs for the imager and spectrometer model to simulate multispectral and hyperspectral data.

C. Reconstruction quality metrics

The quality of the reconstructed cubes was evaluated using four different criteria : the Normalized Root Mean Squared Error (NRMSE), the Structural Dissimilarity (aDSSIM), the Spectral Angle Mapper (aSAM), and the Peak Signal-to-Noise Ratio (PSNR). For all criteria except the PSNR, the lower the value, the better.

The Normalized Root Mean Squared Error (NRMSE) [START_REF] Fienup | Invariant error metrics for image reconstruction[END_REF] measures the cumulative squared error between each pixel of the reconstructed cube x and the true cube x, and writes as

NRMSE(x, x) = ∥x -x∥ 2 ∥x∥ 2 (53)
with a value between 0 and 1, the lower the better the reconstruction.

The Structural SIMilarity (SSIM) [START_REF] Wang | Image quality assessment: from error visibility to structural similarity[END_REF] evaluates the quality of the perceived spatial quality of a reconstructed image w.r.t. the original one. To evaluate the spatial quality of the full reconstructed cube, the average SSIM between all images of both cubes x and x is calculated and rescaled between 0 to 1, where 0 means perfect spatial similarity and 1 perfect anti-correlation. This value is called the average Structural Dissimilarity (aDSSIM) and is written as

aDSSIM(x, x) = 1 2 1 - 1 L L l=1 SSIM(x l , x l ) . (54)
The average Spectral Angle Mapper (aSAM) [START_REF] Yuhas | Discrimination among semi-arid landscape endmembers using the spectral angle mapper (sam) algorithm[END_REF] measures the spectral distorsion between x and x. It calculates the average angular distance between the spectra of both cubes at all pixels. It is formulated as

aSAM(x, x) = 1 P P p=1 arccos ⟨x p , x p ⟩ ||x p || 2 || x p || 2 (55)
and returns a value between 0 and π. The lower this value is, the better the spectral content of the original cube x was recovered in x.

The Peak Signal-to-Noise Ratio [START_REF] Instruments | Peak signal-to-noise ratio as an image quality metric[END_REF] evaluates the quality of the reconstruction by calculating the maximum ratio of a pixel value to the noise value. It is calculated as PSNR(x, x) = 20 log 10 max(x)

MSE(x, x) (56)
and the higher it is, the higher the quality of reconstruction is.

D. Reconstruction with ℓ 2 approach

The explicit procedure for the minimization of the quadratically regularized criterion J ℓ2 proposed in section III-B, now called the ℓ 2 approach, was used to reconstruct abundance maps a ℓ2 from both MS and HS datasets.

Its performance was compared with a gradient-based iterative approach, as proposed in [START_REF] Guilloteau | Hyperspectral and multispectral image fusion under spectrally varying spatial blurs-application to high dimensional infrared astronomical imaging[END_REF], regarding the quality of reconstruction and computational load. Although the exact algorithm used in [START_REF] Guilloteau | Hyperspectral and multispectral image fusion under spectrally varying spatial blurs-application to high dimensional infrared astronomical imaging[END_REF] is not mentioned, the Linear Conjugate Gradient (LCG) [START_REF] Nocedal | Numerical Optimization[END_REF] optimization algorithm was specifically used for the comparison with the proposed ℓ 2 explicit approach of this paper. The value for the only hyperparameter µ r was chosen for both methods to maximize the quality of the reconstructed cube x ℓ2 wrt. the original cube x, according to the NRMSE, defined at Eq. 53.

Inverse problem approaches were also compared with more "naive" reconstruction methods, such as the coaddition, which consists of manually inverting the spectral response for both MS and HS datasets and reformating them to allow the calculation of the average cube. To do so, images of the MS dataset were duplicated for each wavelength on their corresponding spectral band, and images of the HS dataset were upsampled to match the dimensions of MS images.

Fig. 4 compares the reconstructions of the coaddition method, the iterative inverse problem approach (LCG) [START_REF] Guilloteau | Hyperspectral and multispectral image fusion under spectrally varying spatial blurs-application to high dimensional infrared astronomical imaging[END_REF], and the exact solution found with the proposed procedure at a relatively long wavelength (25.5 µm) where instruments data are blurry. Two images are shown for the iterative approach: one after 0.05 seconds of calculations (the computation time of the explicit approach, which corresponds to one iteration), and one after convergence of the algorithm (approximately 6 minutes, 7 000 iterations), with a convergence defined as reaching 0.1% of min(J ℓ2 )1 . Note that precalculations, such as the calculation of Q (Eq. 23), Q -1 (section III-D) and q ℓ2 (Eq. 24) in the Fourier space were not included in these time measurements. For information, in that case, precalculations required 15 seconds for LCG and 17 seconds for the proposed method (longer time due to the inversion of Q, not required for LCG). Tab. I compares the reconstruction quality measurements between the different reconstruction methods.

As expected, results for the coaddition method are relatively poor as the method does not deconvolve images. The intensities of pixels are slightly closer to the original than those given in the instrument's data, but much more precisely recovered with inverse problem approaches. Efficient deconvolution results of the latter approaches are explained by their ability to take into account instrument models (sections II-B and II-C) and exploit the high spatial resolution of short wavelengths to recover high-frequency features at long wavelengths through the spectral correlations introduced with the Linear Mixing Model (section II-A).

At equal computation time, the proposed method performs a significantly finer restoration of the signal than the iterative approach. When allowed to converge, the iterative approach performs visually equivalently to the proposed method, as expected. However, the reconstruction with the explicit approach is 7 000 times faster while demonstrating slightly superior reconstruction quality, as shown in Tab. I. Indeed, for the reconstruction of T = 5 abundance maps of size 124×248 from an MS dataset of size 124×248×9 and an HS dataset of size 31×62×300, the proposed explicit ℓ 2 approach required only 0.05 seconds for a ℓ2 = Q -1 ℓ2 q ℓ2 on a processor i9-10885H, 2.40 GHz (x 16), 64 Go RAM. In contrast, the iterative approach requires a computational time exponentially proportional to the signal-to-noise ratio SNR dB , as shown in Fig. 5. Therefore, for a low noise case (SNR dB = 100 dB), approximately 6 minutes were necessary with the LCG algorithm to reach 0.1% of the minimum of the criteria min(J ℓ2 ). Indeed, the less noisy the observation data, i.e. the higher the SNR dB , the more high frequencies need to be restored, thus requiring a more significant number of iterations. This issue is avoided with the proposed explicit ℓ 2 approach.

It is worth noting that the time-saving factor relative to the LCG algorithm could be higher by setting a convergence threshold stricter than 0.1% of min(J ℓ2 ). In any case, such a time-saving factor could save hours of computational time for the joint processing of heavier MS and HS datasets, e.g. real Required time (seconds)

LCG : 10% of min( ) LCG : 1% x min( ) LCG : 0.1% of min( ) Proposed explicit approach Fig. 5: Time required for the Linear Conjugate Gradient (LCG) algorithm [START_REF] Guilloteau | Hyperspectral and multispectral image fusion under spectrally varying spatial blurs-application to high dimensional infrared astronomical imaging[END_REF] to reach 10%, 1% and 0.1% of min(J ℓ2 ) (Eq. 20), compared to the time required for the proposed explicit solution as a function of the signal-to-noise ratio SNR dB of the MS and HS datasets. Required time for the proposed explicit approach: 0.05 seconds. Precalculations are not taken into account.

observation datasets. Nonetheless, our explicit ℓ 2 approach was tested for the fusion of datasets with dimensions equal to the ones used in [START_REF] Guilloteau | Hyperspectral and multispectral image fusion under spectrally varying spatial blurs-application to high dimensional infrared astronomical imaging[END_REF], i.e. 11 images of size 90 × 900 pixels for the MS dataset, and 4974 images of 30 × 300 pixels for the HS dataset, with only T = 4 spectra templates s used for the Linear Mixing Model. In that case, the fastest approach proposed by [START_REF] Guilloteau | Hyperspectral and multispectral image fusion under spectrally varying spatial blurs-application to high dimensional infrared astronomical imaging[END_REF], the "vectorized approach", allowed a reconstruction after approximately 36 minutes of preprocessing and 20 seconds of processing. The number of iterations or time per iteration was not mentioned in the paper. In comparison, our proposed explicit approach provided the exact minimizer of J ℓ2 with 4 minutes of preprocessing and again only 0.05 seconds for the calculation of a = Q -1 q ℓ2,1 . Indeed, the application of the Hessian Q -1 is not influenced by the number of wavelengths of the input data thanks to the subspace approximation.

E. Reconstruction with ℓ 2,1 approach

The second proposed algorithm aiming to minimize the semi-quadratically regularized criterion J ℓ2,1 of Eq.39, now called the ℓ 2,1 approach, has also been tested for the reconstruction of abundance maps and its performance compared with the ℓ 2 approach.

The chosen convex non quadratic function φ for the penalization term in J ℓ2,1 is the Huber function, defined as

φ(δ) = δ 2 if |δ| < θ, 2θ|δ| -θ 2 otherwise, (57)
where θ ∈ R + is the threshold at which the Huber function switches from quadratic to linear. Thus, higher gradients are less penalized than with a quadratic regularization (ℓ 2 approach), which allows for a finer reconstruction of edges. The ℓ 2,1 approach requires setting two hyperparameters : the regularization parameter µ r , and the Huber threshold θ. Again, hyperparameters were chosen to minimize NRMSE(x, x ℓ2,1 ), Comparison of reconstruction quality measurements between the coaddition, LCG at an equal time with the proposed method and when it converged, and the proposed explicit method for the minimization of the quadratic criterion J ℓ2 defined at Eq. 20.

Methods NRMSE ×10 -3 dSSIM ×10 -5 SAM ×10 -3 PSNR [dB] Time [s]
(a) Original abundance maps a (b) Exact minimum of J ℓ 2 [START_REF] Guilloteau | Hyperspectral and multispectral image fusion under spectrally varying spatial blurs-application to high dimensional infrared astronomical imaging[END_REF] (c) Minimum of J ℓ 2,1

Fig. 6: Comparison between original abundance maps a, on the left, and the ones reconstructed with the explicit minimization of J ℓ2 and with the fast minimization of J ℓ2,1 , where SNR dB = 30 dB for both MS and HS input datasets. Each of these abundance maps is associated in order with the spectra of Fig. 2. All maps share the same dynamic range.

where x ℓ2,1 is the hyperspectral cube reconstructed from the abundance maps a ℓ2,1 .

Ideally, the reconstructed abundance maps should perfectly separate ("unmix") features with different spectra originally superposed in the instrument data. To verify the unmixing for both proposed approaches, the reconstructed maps from both ℓ 2 and ℓ 2,1 approaches are shown in figure 6 for MS and HS datasets with a signal-to-noise ratio of 30 dB. Even though low-frequency features are well recovered with both approaches, the unmixing of high-frequency features (or edges) with the ℓ 2 approach is limited. In fact, a higher noise forces a higher quadratic regularization to stabilize the solution, with the side effect of strongly penalizing any highfrequency feature. Thus, the solution minimizing J ℓ2 in a high noise case are maps sharing these higher frequencies to avoid the presence of sharp edges. In contrast, the ℓ 2,1 approach performs a visibly superior unmixing, and the recovery of high-frequency features is also enhanced.

Images from the cubes reconstructed from the abundance maps of figure 6 are shown on top of Fig. 8. Even if the ℓ 2 approach shows clear deblurring and denoising capabilities despite its inaccurate spectral unmixing, the ℓ 2,1 approach performs a superior edge reconstruction, as expected. On top of an improved spatial reconstruction, the ℓ 2,1 approach performed a significantly finer reconstruction of spectra at longer wavelengths than the ℓ 2 approach, especially in high gradient regions as shown figure 7. Errors on the reconstructed spectra with the coaddition method also tend to be inaccurate edges and smooth gradients were created, and the resulting reconstructions are shown at the bottom of Fig. 8. The optimal reconstruction with the ℓ 2 approach is noisy, and the edges are poorly reconstructed. Conversely, the ℓ 2,1 approach still demonstrates high denoising and edge-preserving capabilities (see Tab. II for a comparison of reconstruction quality measurements).

Compared to the iterative gradient-based ℓ 2 approach from the literature [START_REF] Guilloteau | Hyperspectral and multispectral image fusion under spectrally varying spatial blurs-application to high dimensional infrared astronomical imaging[END_REF], the proposed ℓ 2,1 approach is superior in both quality of reconstruction and computational time. Indeed, for every case tested in this work, the output of the proposed approach at any given computational time had a lower NRMSE to the true object than the output obtained with the literature.

However, the ℓ 2,1 approach can present some limitations compared to the ℓ 2 approach. Indeed, it is iterative, thus necessarily slower, and requires setting convergence thresholds. Moreover, it involves setting 2 hyperparameters (µ r and θ) instead of 1 for the ℓ 2 approach (µ r ), which complicates the search for an optimal solution.

As demonstrated in this study, the ℓ 2 approach could be suited for the fusion of low-noise datasets or datasets mainly composed of low-frequency features. Other situations would require using the ℓ 2,1 approach for an improved reconstruction of high-frequency features, superior spectral unmixing, and resilience to noise. Fig. 9: Example of spectra reconstructed with the coaddition, the ℓ 2 [START_REF] Guilloteau | Hyperspectral and multispectral image fusion under spectrally varying spatial blurs-application to high dimensional infrared astronomical imaging[END_REF] and the ℓ 2,1 approaches, and Fusformer [START_REF] Hu | Fusformer: A transformer-based fusion network for hyperspectral image super-resolution[END_REF] for the same spatial position, taken in a high gradient region. SNR of input datasets: 40 dB. The study was extended with a comparison with neural network approaches as they have been increasingly used for MS/HS fusion. As demonstrated in [START_REF] Hu | Fusformer: A transformer-based fusion network for hyperspectral image super-resolution[END_REF], the self-attention feature of a transformer embedded in a greater neural network architecture allowed it to outperform other known methods of the literature in terms of quality of reconstruction. Indeed, the transformer enabled to inspect the global relationships between features, whereas convolution kernels, for example in CNNbased methods, could only account for local relationships of features [START_REF] Hu | Fusformer: A transformer-based fusion network for hyperspectral image super-resolution[END_REF]. This architecture, called Fusformer, was compared with the proposed inverse problem approaches for the fusion of astronomical data.

F. Comparison with neural network approaches

A trained model of Fusformer provided by its authors was used for this purpose. New astronomical datasets were created to satisfy the model requirements on the input dimensions : 3 × 48 × 48 pixels for the MS dataset, and 31 × 12 × 12 for the HS dataset. These reduced datasets were noised so their SNR is 40 dB, and were the ones used as inputs for all compared methods. Moreover, a normalization of these input datasets was required to set their values between 0 and 1. Then, the output cube of the network was denormalized such that its boundary values were the same as for the input datasets.

Images of the reconstructed cubes with the coaddition, ℓ 2 and ℓ 2,1 approaches and Fusformer are shown Fig. 10. It can be seen that Fusformer successfully recovered some highfrequency features, which is assumed related to the training of the network in the high-pass domain [START_REF] Hu | Fusformer: A transformer-based fusion network for hyperspectral image super-resolution[END_REF]. However, low frequencies were not well recovered, i.e. background values are too intense relative to high-frequency features compared to the original image. Furthermore, a substantial drawback of this network (and most deep neural networks) is that output values are intrinsically bounded between 0 and 1 before denormalization, yet such inverse or deconvolution problems usually lead to restoring unknown values of flux potentially falling outside of a predefined interval. Consequently, naive denormalization cannot restore correct values of flux, as seen in Fig. 10f. Therefore, even though spatial information can be well recovered with this method, it turns inadequate for spectra reconstruction, as shown with reconstructed spectra of Fig. 9.

V. CONCLUSION Typical imager and spectrometer models have been described considering a wavelength-dependent spatial blur, different spectral responses, spectral integrations for the imager, and spatial integrations and subsampling for the spectrometer. Efforts were made to extend the applicability of these models such that the impulse response can be different between the imager and spectrometer model, and such that the wideband images themselves originate from different imager models.

The data fusion problem was first posed through the formulation of a regularized least squares criterion, for which a calculable explicit solution has been demonstrated by exploiting an efficient way to calculate the Hessian matrix. Then, the problem was formulated with an edge-preserving semi-quadratic regularization. This latter problem involved an alternating minimization problem, which has been accelerated with the same Hessian inversion procedure.

These data fusion methods have been applied to simulated astronomical images from the James Webb Space Telescope. The benefits of the explicit solution for the exact minimization of the quadratic criterion against an iterative approach have been proven with a time-saving factor of 7 000 for the fusion of low noise data (SNR dB = 100 dB).

A comparison of the approaches with quadratic ℓ 2 and semi-quadratic ℓ 2,1 regularizations has also been carried out and highlighted the improved ability of the latter approach to unmix spectral information, to better reconstruct spatial information (especially in high gradient areas) and spectra, which also led to superior denoising capabilities.

Finally, a performance comparison with a deep neural network architecture with an embedded transformer trained for the fusion of MS and HS datasets has been conducted and led to the conclusion that the proposed inverse problem approaches outperformed such a method in this context of data fusion.

Further work could include estimations of optimal values of hyperparameters to develop a non-supervised reconstruction procedure. This would result in drastic time savings as iterative searches for optimal hyperparameter values would be avoided. Other works could be conducted to improve the quality of the reconstruction, for example by imposing a positivity constraint, or by using wavelet or deep neural network-based regularizations. where the blocks ∆ t,t ′ are diagonal, and the blocks ∆ ′ p , p ∈ [1, P ], are matrices of size T × T , built so that the coefficient of the row t and column t ′ is the p th coefficient of the diagonal of ∆ t,t ′ ∆ ′ p t,t ′ = (∆ t,t ′ ) p,p .

Such permutations allow for the creation of a new block diagonal matrix where all Fourier coefficients have been sorted by frequency (one block per frequency). As done in [START_REF] Hadj-Youcef | Fast joint multiband reconstruction from wideband images based on low-rank approximation[END_REF], this matrix is directly invertible. Using the orthogonality of permutations matrices, i.e. P -1 = P T , the inversion of Eq. 59 writes as where δ t,t ′ have the same size and structure than ∆ t,t ′ , i.e. diagonal matrices of size P × P .

The inversion of Q ℓ2 is then calculated as which proves that blocks of Q -1 ℓ2 indeed follow a diagonal structure in the Fourier space, as for Q ℓ2 . This property allows for fast computation and application of Q -1 ℓ2 on q ℓ2 in the Fourier space to explicitly obtain the minimizer of J ℓ2 with Eq. 22.

Q -1 ℓ2 =     F H     ∆ 1,
Therefore, the inversion of the initial matrix Q ℓ2 require T 2 fft to diagonalize the blocks, and P inversion of matrix of size T ×T . In the case of this paper, P is the number of pixels, and T is the number of spectral templates. Such configuration can be solved efficiently by an actual personal computer and linear algebra libraries.

Fig. 1 :

 1 Fig. 1: PSF of MIRI/JWST.

Fig. 2 :

 2 Fig. 2: Spectra for the Linear Mixing Model.

Fig. 3 : 2 2

 32 Fig. 3: Original abundance maps a

Fig. 4 :

 4 Fig. 4: Imagery data y m and spectrometry data y h with SNR dB = 100 dB on top, and reconstruction with the coaddition method Fig. 4c. The next figures show results of the minimization of J ℓ2 with the LCG algorithm in 0.05 second and 6 minutes (value of criterion at 0.1% of min(J ℓ2 )), and with the demonstrated explicit method (section III-B) in 0.05 seconds. The last figure is the original image. Central wavelength for all images is 25.5 µm. Unit is mJy arcsec -2 . The dynamical scales for the three last figures are identical.

Fig. 7 :Fig. 8 :

 78 Fig.7: Example of spectra reconstructed with the coaddition, the ℓ 2[START_REF] Guilloteau | Hyperspectral and multispectral image fusion under spectrally varying spatial blurs-application to high dimensional infrared astronomical imaging[END_REF] and the ℓ 2,1 approaches for the same spatial position, taken in a high gradient region. SNR of input dataset: 30 dB.

Fig. 10 :

 10 Fig. 10: Imagery data y m and spectrometry data y h with SNR dB = 40 dB on top, and below are reconstructions with the coaddition method, the ℓ 2 and ℓ 2,1 approaches, Fusformer, and the original image. All results were produced from the detailed reduced MS and HS datasets. The central wavelength for all images is 25.5 µm. Unit is mJy arcsec -2 . Same dynamical scale for figures (d), (e), and (g).

  APPENDIX A PROOF OF Q ℓ2 INVERSIONWe consider a square matrix made of diagonal blocks ∆ t,t ′

  1 . . . δ 1,T . . . . . . δ T,1 . . . δ T,T

1 .

 1

TABLE I :

 I

TABLE II :

 II Comparison of reconstruction quality measurements between the coaddition, the ℓ 2 and the ℓ 2,1 approaches for input MS and HS datasets with SNR dB = 30 dB.

Université Paris-Saclay, CNRS, CentraleSupélec, Laboratoire des Signaux et Systèmes, 3 rue Joliot Curie, 91190 Gif-sur-Yvette, France

Institut d'Astrophysique Spatiale (IAS), CNRS, UMR 8617, Université Paris-Saclay, 91405 Orsay, France

Innovative Common Laboratory For Space Spectroscopy (INCLASS), IAS, ACRI-ST, CNRS, Université Paris-Saclay, 91405 Orsay, France

The value of min(J ℓ

) comes from the proposed explicit approach.

ACKNOWLEDGMENTS

This work is funded by the Centre National d' Études Spatiales (CNES) and the Agence Nationale de la Recherche (ANR), and was supported by a collaboration with ACRI-ST as part of the INCLASS joint laboratory.


domain: math sdv

text:
Nathalie Eisenbaum
email: nathalie.eisenbaum@parisdescartes.fr

Nicolas Meunier
email: nicolas.meunier@univ-evry.fr

A stochastic lipid structured model for macrophage dynamics in atherosclerotic plaques

, is to include offloading of lipids through the dumping term. We find explicit necessary and sufficient conditions for macrophage trapping in the locally enriched region.

Introduction

Atherosclerosis is a major cause of death in industrialized societies since it is the primary cause of heart attack (acute myocardial infarction) and stroke (cerebrovascular accident). It is now accepted that atherosclerosis is a chronic inflammatory disease which starts within the intima, the innermost layer of an artery. It is driven by the accumulation of macrophage cells within the intima and promoted by modified low density lipoprotein (LDL) particles [START_REF] Moore | Macrophages in the pathogenesis of atherosclerosis[END_REF]. Inflammation occurs at sites within the arterial wall where modified low-density lipoproteins (LDL) accumulate after penetrating the wall from the bloodstream [START_REF] Williams | The response-to-retention hypothesis of early atherogenesis[END_REF]. The immune response attracts circulating monocytes to these sites [START_REF] Van Furth | The origin and kinetics of mononuclear phagocytes[END_REF]. After crossing the endothelium (the layer of cells lining the lumen of the blood vessel), monocytes enter the blood vessel wall and usually differentiate into macrophages that remove LDL by phagocytosis [START_REF] Moore | Macrophages in atherosclerosis: a dynamic balance[END_REF]. When the macrophages accumulate internalised lipid, they produce inflammatory signals that promote the recruitment of additional monocytes [START_REF] Tall | Cholesterol, inflammation and innate immunity[END_REF][START_REF] Meunier | Mathematical Study of an Inflammatory Model for Atherosclerosis: A Nonlinear Renewal Equation[END_REF]. In this way, the persistent influx of LDL and the influx of monocytes into the arterial wall produce large numbers of lipid-laden macrophages. This positive feedback loop creates an unresolved inflammatory response [START_REF] Nathan | Nonresolving inflammation[END_REF] that thickens and enlarges the arterial wall, to produce an atherosclerotic plaque. Plaques can continue to grow throughout life and accumulate lipid-rich necrotic debris (designated as the necrotic core) through necrotic cell death of live and apoptotic macrophages [START_REF] Moore | Macrophages in the pathogenesis of atherosclerosis[END_REF]. Plaques may rupture due to the deterioration of the arterial wall, resulting in the release of thrombotic (blood clot-forming) material. If these blood clots block arteries, they can cause heart attacks or strokes. The intensity of inflammation, in general, depends on the number of LDL particles and macrophages located within the arterial wall. The number of monocyte-derived macrophages inside the plaque is determined by the balance between rapid monocyte recruitment [START_REF] Van Furth | The origin and kinetics of mononuclear phagocytes[END_REF][START_REF] Moore | Macrophages in atherosclerosis: a dynamic balance[END_REF], macrophage proliferation [START_REF] Jenkins | Local macrophage proliferation, rather than recruitment from the blood, is a signature of th2 inflammation[END_REF], programmed cell death (apoptosis) [START_REF] Moore | Macrophages in the pathogenesis of atherosclerosis[END_REF], or macrophage emigration from the arterial wall [START_REF] Llodrá | Emigration of monocyte-derived cells from atherosclerotic lesions characterizes regressive, but not progress ive, plaques[END_REF]. Macrophages accumulate lipids by phagocytosis of LDL and apoptotic cells. The recent demonstration of the dynamic nature of inflammatory macrophage populations (see [START_REF] Moore | Macrophages in atherosclerosis: a dynamic balance[END_REF]) has led to the development of a new line of mathematical modeling of atherosclerotic plaque immunology (eg [START_REF] Ford | A lipid-structured model for macrophage populations in atherosclerotic plaques[END_REF]). Many mathematical models divide the macrophage population into "macrophages" and "foam cells". This is true in ODE models [START_REF] Cohen | Athero-protective effects of high density lipoproteins (hdl): an ode model of the early stages of atherosclerosis[END_REF], spatially resolved PDE models [START_REF] Chalmers | Bifurcation and dynamics in a mathematical model of early atherosclerosis[END_REF][START_REF] Hao | The ldl-hdl profile determines the risk of atherosclerosis: a mathematical model[END_REF][START_REF] Friedman | A mathematical model of atherosclerosis with reverse cholesterol transport and associated risk factors[END_REF][START_REF] Calvez | Mathematical modeling of the atherosclerotic plaque formation[END_REF][START_REF] Calvez | Mathematical and numerical modeling of early atherosclerotic lesions[END_REF] and other computational models [START_REF] Parton | Computational modelling of atherosclerosis[END_REF]. However, it is now clear that there is a continuous distribution of lipid loads in macrophage populations, ranging from monocytes with low lipid loads to macrophages that can be labeled as foam cells in atherosclerosis. The variation in lipid load within a macrophage population suggests that it may be instructive to develop structured models in which macrophages are characterized by their intracellular lipid load, similar to adipocytes [START_REF] Soula | Modelling adipocytes size distribution[END_REF], in the spirit of [START_REF] Meunier | Mathematical Study of an Inflammatory Model for Atherosclerosis: A Nonlinear Renewal Equation[END_REF][START_REF] Ford | A lipid-structured model for macrophage populations in atherosclerotic plaques[END_REF].

In this paper, we model a single macrophage that is moving randomly through tissue where lipid is distributed in a non-uniform way. When the macrophage is at the same position as the lipid, it ingests lipid. Ingested lipid reduces the macrophage's motility and this is modelled as a reduction in the diffusion coefficient that governs the macrophage's rate of motion. We assume that lipid is concentrated at a point which we denote x = 0 and whenever a macrophage passes x = 0 it ingests lipid. More precisely, to model the undergoing random movement of a single macrophage in an arterial wall, we identify the arterial wall with the real line and the lipidenriched region with the point 0 (see Fig. 1). Our model gives the position X t of a macrophage cell and its diffusion coefficient A t , at time t, as solution of the following system of stochastic differential equations (1)

dX t = √ 2A t dW t , dA t = -A α t dL 0 t (X) + λA t dt,
with initial condition (X 0 , A 0 ), X 0 ∈ IR, A 0 > 0, where λ and α are two nonnegative parameters, (W t , t ≥ 0) is a given one dimensional Brownian motion, and (L 0 t (X), t ≥ 0) is the local time process at 0 of the unknown process X = (X t , t ≥ 0). The aim of our analysis is to find the conditions which lead to a macrophage becoming unable to move due to lipid ingestion. The system (1) will provide answers depending on the parameters X 0 , A 0 , α, λ. They are presented in Proposition 2.3 and analyzed in section 6. We start the biological justification of (1) by explaining the presence of local time. The notion of local time has been already used to model ((X t , A t ), t ≥ 0) in [START_REF] Bénichou | Non-Gaussianity and dynamical trapping in loccaly activeted random walks[END_REF][START_REF] Meunier | Long time behavior in locally activated random walks[END_REF]. Indeed their model corresponds to the case λ = 0. For a rigorous introduction to the notion of local time one can see e.g. [START_REF] Revuz | ume 293 of Grundlehren der Mathematischen Wissenschaften[END_REF]. Intuitively, L 0 t (X) represents the time spent by X at point 0 during the time interval [0, t], and actually (L 0 t (X), t ≥ 0) is an increasing process which increases only at times when X equals 0. Hence according to (1) each crossing of the origin by the macrophage cell tends to decrease its diffusion coefficient (A t , t ≥ 0), or equivalently to decelerate the cell. Such a decrease in local mobility is due to the internalization and accumulation of lipid particles by macrophages. It can be related to several phenomena: increase in volume or change in cell fate depending on the amount of LDL ingested. The parameter α describes the intensity of the internalization: the higher the α, the lower the internalization.Indeed, an healthy macrophage (without lipids) moves slowly through tissue (A 0 is not large i.e. A 0 < 1) and A t decreases with each passage through the lipid-dense zone. In such a case, the positive parameter α describes the change in mobility when passing through the lipid-dense region. The smaller the α parameter, the greater the change. If α = 0, then the internalization is the same, and constant, at each macrophage passage through the origin. Therefore, this parameter is related to both the amount of lipid located at point 0 of the arterial wall and to the capacity of each macrophage to internalize lipid. The case α = 0 corresponds to a very large amount of lipids located at point 0 and a very high capacity of macrophages to internalize lipids. In summary, the α parameter is related to the inflammatory state of the lesion. In contrast, the term λA t dt describes the loss of lipids by the macrophage and its tendency to recover its natural state. Therefore, the parameter λ describes the natural immune defenses or medical treatment. Finally, the initial condition (X 0 , A 0 ) is related to the condition of the arterial wall. Indeed X 0 represents the distance between the macrophage penetrating the wall and a lipid-enriched region. A small value of X 0 can be interpreted as a situation where lipid-enriched regions are frequent. Similarly, a small value of A 0 means a very inflamed macrophage with a low ability to move. Note that (A t , t ≥ 0) is continuous and starts from A 0 > 0. The first time t when A t = 0, is the time of the cell arrest. It allows us to quantify the necrotic core. The couple ((X t , A t ), t ≥ 0) describes the dynamics of a monocyte -derived macrophages population structured by its lipid accumulation, see Fig. 2. It accounts for the impact of lipids ingested by macrophages on their size, dynamics and structure. At the core of our model is this state variable, which plays a similar role to age in age-structured population models and to size in size-structured models [START_REF] Werner | The ontogenetic niche and species interactions in sizestructured populations[END_REF][START_REF] Charlesworth | Evolution in age-structured populations[END_REF].

Although simple, the biophysical question that the model (1) raises is whether spatially localized mobility perturbations have consequences on diffusion at all scales, such as the emergence of a dynamical transition to an absorbing static state. In the context of atherosclerosis, this dynamical transition can be viewed as a minimal mechanism that causes macrophages to aggregate in lipid-enriched regions and thereby to the formation of atherosclerotic plaques (necrotic core). The model [START_REF] Bénichou | Non-Gaussianity and dynamical trapping in loccaly activeted random walks[END_REF] incorporates cellular processes which are known to play a significant role in the dynamics of macrophages, apoptotic cells and necrotic material within plaques.

The remainder of this paper is organised as follows. In section 2 we describe the main results: the explicit unique solution of (1) and its arrest time. To establish them, we first give in section 3 an equivalent system to [START_REF] Bénichou | Non-Gaussianity and dynamical trapping in loccaly activeted random walks[END_REF]. In section 4 we prove the main results. In section 5 we show that the solution of (1) in the case λ = 0 (corresponding to the model introduced in [START_REF] Bénichou | Non-Gaussianity and dynamical trapping in loccaly activeted random walks[END_REF][START_REF] Meunier | Long time behavior in locally activated random walks[END_REF]) can be obtained as the limit of the cases λ > 0 as λ tends to 0. In section 6 we discuss our results and draw conclusions. In particular, we underline the possibility when α < 1/2 (cases of high inflammatory intensity) for the initial condition (X 0 , A 0 ) to determine the fate of a macrophage governed by [START_REF] Bénichou | Non-Gaussianity and dynamical trapping in loccaly activeted random walks[END_REF]. This sensitivity to the initial condition, that is the initial inflammatory state, should strenghten the interest of the model (1).

Main results

Our objective is to study the competition between cell mobility and intracellular lipid accumulation to see whether, after several visits to the lipid-enriched region, the macrophage cell will eventually be trapped (corresponding to necrosis) and thus contribute to atherosclerotic plaque formation or whether the offload will be strong enough and prevent plaque formation. As introduced in section 1, the cell arrest is the first time A hits 0. Denote it by T 0 (A). More precisely: T 0 (A) = inf {t ≥ 0 : A t = 0}, with the usual convention inf ∅ = +∞. Note that T 0 (A) = +∞ is compatible with lim t→+∞ A t = 0, but does not imply it. The event {T 0 (A) < ∞} corresponds to the situation where the particle gets trapped at the origin in finite time.

In the present framework it makes no sense to go beyond the time T 0 (A), hence (3) after many visits to the lipid-enriched region, the macrophage eventually gets trapped, resulting in the formation of an atherosclerotic plaque.

a solution of ( 1) is actually a process ((X t , A t ), 0 ≤ t < T 0 (A)) satisfying (1) up to time T 0 (A). A solution of ( 1) is said to be strong if (T 0 (A), ((

X t , A t ), 0 ≤ t < T 0 (A))) is such that T 0 (A) is a stopping time of the Brownian motion (W t , t ≥ 0) and ((X t , A t ), 0 ≤ t < T 0 (A)) is a continuous process adapted to (W t , t ≥ 0) satisfying (1) until time T 0 (A).
Define (Y t , t ≥ 0) as the Ornstein-Uhlenbeck process solution of (2)

Y t = X 0 √ 2A 0 + W t - λ 2 t 0 Y s ds.
Denote by L 0 t (Y ), t ≥ 0 the local time process at 0 of Y . The following theorem gives the explicit solution of (1) for any α ≥ 0, λ ≥ 0, A 0 > 0 and any real X 0 . Theorem 2.1. For any couple (α, λ) of nonnegative real numbers, the system (1) with initial condition (X 0 , A 0 ), A 0 > 0, has a unique strong solution given by

(3) ((X t , A t ), 0 ≤ t < T 0 (A)) = ( 2A t Y t , A t ), 0 ≤ t < T 0 (A) . If α = 1/2, then a.s. T 0 (A) = +∞ and (A t , t ≥ 0) = A 0 exp{λt - √ 2L 0 t (Y )}, t ≥ 0 . If α = 1/2, then (A t , 0 ≤ t < T 0 (A)) = A 0 e λt 1 + √ 2(α - 1 2 )A α-1 2 0 t 0 e λ(α-1 2 )s dL 0 s (Y ) 1 1 2 -α , 0 ≤ t < T 0 (A) .
For α = 1/2, the nature of T 0 (A) depends on the value of the parameters α, λ, X 0 and A 0 . To reduce the general case (α, λ, X 0 , A 0 ) to the case (α, λ, 0, A 0 ), we will establish Proposition 2.2 below. To enunciate it, we make use of T 0 (X) the first time X hits 0 and define T 0 A T 0 (X)+ • on {T 0 (X) < ∞} by: T 0 A T 0 (X)+ • = inf t ≥ 0 : A T 0 (X)+t = 0 . Note that T 0 (X) ≤ T 0 (A) a.s. Proposition 2.2. For ((X t , A t ), 0 ≤ t < T 0 (A)) solution of (1) with initial condition (X 0 , A 0 ), one has: T 0 (X) < ∞ a.s. and the process (X T 0 (X)+t , A T 0 (X)+t ), 0 ≤ t < T 0 (A T 0 (X)+ • ) satisfies (1) with respect to another Brownian motion W , and initial condition (0, Ã0 ) where

Ã0 = A 0 + λ X 2 0 2N 2 ,
and N is a real standard Gaussian variable independent of W .

To describe the nature of T 0 (A) when α = 1/2, it is hence sufficient to consider the cases with initial condition (0, A 0 ), A 0 > 0. Note that when X 0 = 0, the Ornstein-Uhlenbeck process (Y t ) t≥0 is independent of A 0 . For α < 1/2 and λ > 0, set

V α,λ = ∞ 0 exp -λs( 1 2 -α) dL 0 s (Y ). (4)
Note that V α,λ is a random time decreasing with λ( 1 2 -α).

Proposition 2.3. Let ((X t , A t ), 0 ≤ t < T 0 (A)) be the solution of (1) with initial condition (0, A 0 ) (A 0 > 0).

• If α < 1/2 and λ > 0,

then IE[V α,λ ] < ∞ and the following dichotomy holds:

-on {V α,λ > A 1 2 -α 0 √ 2( 1 2 -α) }, one has T 0 (A) < ∞ and lim t→T 0 (A) (X t , A t ) = (0, 0), -on {V α,λ ≤ A 1 2 -α 0 √ 2( 1 2 -α) }, one has T 0 (A) = ∞, lim t→∞ A t = +∞,
and X is recurrent.

• If α < 1/2 and λ = 0 then a.s. T 0 (A) < ∞ and lim t→T 0 (A) (X t , A t ) = (0, 0).

• If α ≥ 1/2 and λ ≥ 0, then a.s. T 0 (A) = ∞, and 0 is a recurrent point for X.

Moreover one has:

-if α = 1/2, then a.s. lim sup t→∞ e -λt A t ≤ A 0 and lim sup t→∞

te -λt |Xt| log t ≤ 2 √ A 0 , -if 1/2 < α < 1, then a.s. lim inf t→∞ A t = 0,
-if α = 1, then a.s. lim inf t→∞ A t = 0, lim sup t→∞ A t = +∞, and X is recurrent.

-if α > 1, then a.s. lim sup t→∞ A t = +∞ and X is recurrent.

We will interpret in section 6 the results presented by Proposition 2. Nevertheless we already make the following comments.

For α ≥ 1/2: one has IP [T 0 (A) < ∞] = 0, whatever the value of the other parameters.

For α < 1/2:

IP [T 0 (A) < ∞] = IP [V α,λ > A 1 2 -α 0 √ 2( 1 2 -α) ]. Since for X 0 = 0, (Y t ) t≥0 is independent of A 0 , so is V α,λ . One sees hence immediately that IP [T 0 (A) < ∞] is a decreasing function of A 0 . Moreover, since V α,λ is a decreasing function of λ, IP [T 0 (A) < ∞] is also a decreasing function of λ.
As it has been mentioned in the introduction of the α parameter, to make sense with the modelisation with the macrophage dynamic, A 0 must be taken small. With elementary computations, one obtains when

A 0 ≤ e 2 that IP [T 0 (A) < ∞] is a decreasing function of α.
In Proposition 2.3, the results concerning λ = 0 have all been already established in [START_REF] Meunier | Long time behavior in locally activated random walks[END_REF]. Since when λ = 0, the Ornstein-Uhlenbeck process Y is the Brownian motion W , their asymptotics when α ≥ 1/2 are more precise.

An equivalent system to (1)

We give a system equivalent to (1), which will be simpler to handle. Let (X, A) be a solution of (1) on [0, T 0 (A)) with A 0 > 0. On [0, T 0 (A)) define Y by ( 5)

Y t = X t √ 2A t .
We obtain that (Y, A) is a solution of the following system (6)

Y t = Y 0 + W t -λ 2 t 0 Y s ds, dA t = - √ 2A α+ 1 2 t dL 0 t (Y ) + λA t dt, with Y 0 = X 0 √ 2A 0 .
Indeed, using Ito's formula together with the fact that the measure dL 0 s (X) is supported by the set {s : X s = 0}, one has

Y t = Y 0 + t 0 dX s √ 2A s - t 0 X s 2 √ 2 A -3/2 s dA s = Y 0 + W t + t 0 X s (2A s ) 3/2 A α s dL 0 s (X) -λ t 0 X s A s (2A s ) 3/2 ds = Y 0 + W t - λ 2 t 0 X s √ 2A s ds = Y 0 + W t - λ 2 t 0 Y s ds.
Besides, there is a connection between the local times at 0 of X and Y :

dL 0 t (X) = 2A t dL 0 t (Y )
. This is a consequence of Tanaka's formula (see e.g. [START_REF] Revuz | ume 293 of Grundlehren der Mathematischen Wissenschaften[END_REF]) according to which

L 0 t (X) = |X t | -|X 0 | - t 0 sgn(X s )dX s = 2A t |Y t | -2A 0 |Y 0 | - t 0 sgn(Y s )dX s = 2A t |Y t | -2A 0 |Y 0 | - t 0 sgn(Y s ) 2A s dY s - t 0 |Y s | dA s √ 2A s . Since 2A t |Y t | = 2A 0 |Y 0 | + t 0 2A s d(|Y s |) + t 0 |Y s | dA s √ 2A s , it yields L 0 t (X) = t 0 2A s dL 0 s (Y ),
which leads to [START_REF] Cohen | Athero-protective effects of high density lipoproteins (hdl): an ode model of the early stages of atherosclerosis[END_REF].

Conversely any couple (Y, A) solution of ( 6) on [0, T 0 (A)), corresponds to a solution (X, A) of ( 1) on [0, T 0 (A)), by setting X t = √ 2A t Y t . We now work with (Y, A) solution of ( 6) on [0, T 0 (A)). The first equation in ( 6) defines an Ornstein-Uhlenbeck process (Y t , t ≥ 0) which can be explicitly computed ( 7)

Y t = e -λt 2 Y 0 + e -λt 2 t 0 e λs 2 dW s ,
where Y 0 is the given initial value of Y . Alternatively, one can use the following expression of Y which is due to Doob [START_REF] Doob | The Brownian movement and stochastic equations[END_REF] (8)

(Y t , t ≥ 0) = e -λt 2 Z u(t) + Y 0 , t ≥ 0 ,
where Z is a real Brownian motion starting from 0 and ( 9)

u(t) = e λt -1 λ .
Actually, Z is defined by

(10) (Z t , t ≥ 0) = 1 λ log(1+λt) 0 e λs 2 dW s , t ≥ 0 .
We have to study the behavior of (A t ) t≥0 solution of the second equation of ( 6)

(11) A t = A 0 - √ 2 t 0 A α+ 1 2 s dL 0 s (Y ) + λ t 0 A s ds,
with A 0 > 0 and t ∈ [0, T 0 (A)).

In [START_REF] Eisenbaum | A Lifetime of Excursions Through Random Walks and Lévy Processes[END_REF], using [START_REF] Eisenbaum | A Lifetime of Excursions Through Random Walks and Lévy Processes[END_REF] allows to obtain the following expression for the local time process at 0 of Y

(12) (L 0 t (Y ), t ≥ 0) =   e λt -1 λ 0 1 √ 1 + λs dL -Y 0 s (Z), t ≥ 0   . 4 Proofs 4.1 Proof of Proposition 2.2
We use of Y and Z introduced in ( 7), ( 8) and [START_REF] Friedman | A mathematical model of atherosclerosis with reverse cholesterol transport and associated risk factors[END_REF]. The parameters α and λ are fixed. For X 0 = 0, define

T 0 (X) = inf {t ≥ 0 : X t = 0} .
Since A is increasing on [0, T 0 (X)), recalling the definition ( 5) of Y , one has T 0 (X) = T 0 (Y ). Using [START_REF] Doob | The Brownian movement and stochastic equations[END_REF] it follows that a.s. T 0 (X) < ∞.

The time T 0 (Y ) is a stopping time for the natural filtration of W . Hence the process ( Wt , t ≥ 0) defined by W = (W T 0 (X)+t -W T 0 (X) , t ≥ 0) is a Brownian motion starting from 0, independent of (W t , t ∈ [0, T 0 (X)]). The process (Y T 0 (X)+t , t ≥ 0 is the unique solution of

y t = Wt - λ 2 t 0 y s ds.
By [START_REF] Hao | The ldl-hdl profile determines the risk of atherosclerosis: a mathematical model[END_REF], we see that A t = A 0 exp(λt) for t ∈ [0, T 0 (X)] and in particular it holds that A T 0 (X) = A 0 exp(λT 0 (Y )). Hence, A T 0 (X) is independent of W . Moreover, (A T 0 (X)+t , t ≥ 0) is solution of

a t = A T 0 (X) - t 0 √ 2a α+ 1 2 s dL 0 s (Y T 0 (Y )+. ) + λ t 0 a s ds.
Consequently, ((Y T 0 (X)+t , A T 0 (X)+t ), t ≥ 0) is solution of ( 6) with respect to W with initial condition (0, A T 0 (X) ). This implies that ((X T 0 (X)+t , A T 0 (X)+t ), t ≥ 0) is solution of (1) with respect to W with initial condition (0, A T 0 (X) ). Set Ã0 = A T 0 (X) . Thanks to (8), we have

T 0 (Y ) = 1 λ log 1 + λT -X 0 √ 2A 0 (Z) ,
where for all a ∈ IR, T a (Z) is defined by

T a (Z) = inf {t ≥ 0 : Z t = a} .
Thanks to the reflexion principle, one sees that T a (Z)

(law) = ( a Z 1 ) 2 , which leads to Ã0 = A 0 + λ X 2 0 2N 2 ,
where N is a real standard Gaussian variable independent of W .

Proof of Theorem 2.1 and Proposition 2.3

To solve (1) we use the equivalent system [START_REF] Cohen | Athero-protective effects of high density lipoproteins (hdl): an ode model of the early stages of atherosclerosis[END_REF]. We start with the case α = 1/2 and then study the case α = 1/2.

4.2.1

Case α = 1/2

In this case [START_REF] Hao | The ldl-hdl profile determines the risk of atherosclerosis: a mathematical model[END_REF] reads

A t = A 0 + t 0 A s (λds - √ 2dL 0 s (Y )), t ∈ [0, T 0 (A)),
which admits as unique solution on [0, T 0 (A)) ( 13)

A t = A 0 exp λt - √ 2L 0 t (Y ) .
Consequently, A can not reach 0 at finite time, hence

IP [T 0 (A) = ∞] = 1.
Assume now that X 0 = 0. Since:

X t = √ 2A t Y t , using (13)
together with (8) one easily obtains:

|X t | = 2A 0 |Z u(t) | exp - 1 √ 2 L 0 t (Y ) ,
with u defined by [START_REF] Ford | A lipid-structured model for macrophage populations in atherosclerotic plaques[END_REF]. Thanks to [START_REF] Jenkins | Local macrophage proliferation, rather than recruitment from the blood, is a signature of th2 inflammation[END_REF], one obtains

|X t | ≤ 2A 0 |Z u(t) | exp - 1 √ 2
L -Y 0 u(t) (Z) 1 + λu(t) ,
where L -Y 0 u(t) (Z) stands for the local time at level -Y 0 and time u(t) of the process Z. From Kesten's asymptotic results on the supremum of the local times of Z, see [START_REF] Kesten | An iterated logarithm law for local time[END_REF], it follows that a.s.

lim inf t→∞ L -Y 0 t (Z) √ log log t √ t < ∞, hence a.s. ( 14
) lim inf t→∞ L -Y 0 t (Z) √ t = 0.
Using [START_REF] Llodrá | Emigration of monocyte-derived cells from atherosclerotic lesions characterizes regressive, but not progress ive, plaques[END_REF] together with the law of iterated logarithm for Z, we get that a.s.

lim sup t→∞ te -λt |X t | log t ≤ 2 A 0 .
Besides, one has

A t ≤ A 0 (1 + λu(t)) exp - √ 2
L -Y 0 u(t) (Z) 1 + λu(t) ,
which with [START_REF] Llodrá | Emigration of monocyte-derived cells from atherosclerotic lesions characterizes regressive, but not progress ive, plaques[END_REF] gives: lim sup t→∞ e -λt A t ≤ A 0 a.s.

4.2.2

Case α = 1/2

Let β be a non zero real number. By Ito's formula, one obtains on [0, T 0 (A))

A β t = A β 0 - √ 2β t 0 A β+α-1 2 s dL 0 s (Y ) + λβ t 0 A β s ds. Choose β = 1 2 -α and set B t = A β t . For t in [0, T 0 (A)), we have (15) B t = B 0 - √ 2βL 0 t (Y ) + λβ t 0 B s ds, with B 0 = A β 0 . By differentiation of B t e -λβt , one sees that d B t e -λβt = - √ 2βe -λβt dL 0 t (Y ),
which leads to ( 16)

B t = e λβt B 0 - √ 2β t 0 e -λβs dL 0 s (Y ) .
We now study T 0 (A) according to the value of α.

Case α < 1/2 : In this case β > 0 and (16) leads to (17)

A t = e λt B 0 - √ 2( 1 2 -α) t 0 e -λ( 1 2 -α)s dL 0 s (Y ) 1 ( 1 2 -α)
, t ∈ [0, T 0 (A)).

Moreover we have

T 0 (A) = T 0 (B) = inf t ≥ 0 : t 0 e -λβs dL 0 s (Y ) ≥ B 0 √ 2β .
Assume now that X 0 = 0. Under this assumption Y starts from 0 and is independent of A 0 . Note that

IE ∞ 0 e -λβs dL 0 s (Y ) = IE L 0 e (Y )
,

where e is an exponential variable with parameter λ( 1 2 -α), independent of Y . Using the fact that Y t has a Gaussian law with mean 0 and variance 1-e -λt λ , one has

IE[L 0 e (Y )] = λ 2π ∞ 0 e -λβs (1 -e -λs ) -1/2 ds = 1 √ 2πλ ∞ 0 e -βs (1 -e -s ) -1/2 ds < ∞.
Consequently a.s.

V λ,α := ∞ 0 exp -λs( 1 2 -α) dL 0 s (Y ) < ∞. On {V λ,α > A 1 2 -α 0 √ 2( 1 2 -α) }: T 0 (A) < ∞ and lim t→T 0 (A) X t = 0. On {V λ,α ≤ A 1 2 -α 0 √ 2( 1 2 -α) }: T 0 (A) = +∞ and lim t→∞ B t = lim t→∞ e λ( 1 2 -α)t A 1 2 -α 0 - √ 2( 1 2 -α)V λ,α =
(Y ) = L -Y 0 u(t) (Z) (1 + λu(t)) 1-α -λ(α -1) u(t) 0 L -Y 0 s (Z)(1 + λs) α-2 ds.
To finish proving Proposition 2.2, we assume now that X 0 = 0 and consider the three different cases: α = 1, α > 1 and 1/2 < α < 1. If α = 1, then

B t ∼ t→∞ √ 2γ L 0 u(t) (Z) (1 + λu(t)) 1/2 .
Using the law of iterated logarithm for local times of a real Brownian motion [START_REF] Perkins | On the iterated logarithm law for local time[END_REF], we obtain a.s. [START_REF] Parton | Computational modelling of atherosclerosis[END_REF] lim sup

t→∞ L 0 t (Z) 2t log 2 t = 1,
hence lim inf t→∞ A t = 0. Using [START_REF] Llodrá | Emigration of monocyte-derived cells from atherosclerotic lesions characterizes regressive, but not progress ive, plaques[END_REF], it follows that lim inf t→∞ B t = 0, thus lim sup t→∞ A t = +∞. Using the same argument as in [START_REF] Meunier | Long time behavior in locally activated random walks[END_REF] one obtains that X is recurrent on IR with lim inf t→∞ X t = -∞ and lim sup t→∞ X t = +∞. If α > 1, then from [START_REF] Nathan | Nonresolving inflammation[END_REF], we deduce that there exists a positive constant C such that

(21) B t ≤ C L 0 u(t) (Z) u(t) 1/2
and hence ( 14) leads similarly to lim sup t→∞ A t = +∞. Moreover using ( 21), for t big enough it holds

|X ] = 2A t |Y t | ≥ C |Z u(t) | u(t) ( u(t) L 0 u(t) (Z) ) 1 2γ .
Finally, using (20), a.s. for all t big enough we have

u(t) L 0 u(t) (Z) ≥ C 1 log 2 u(t
B t ≥ C L 0 u(t) (Z) u(t) 1/2 ,
and (20) gives: lim inf t→∞ A t = 0 a.s.

5

The case λ = 0 as a limit case

For α ≥ 0, λ ≥ 0, X 0 ∈ IR and A 0 > 0 fixed, there exists a unique strong solution of (1). Various continuity properties of this solution with respect to this parameters, can be studied. In particular for α ≥ 0 and A 0 > 0 fixed, denote by (X(λ), A(λ)) the solution of (1) with initial condition (0, A 0 ). Motivated by the fact that the case λ = 0 with initial condition (0, A 0 ), has been previously explored in [START_REF] Meunier | Long time behavior in locally activated random walks[END_REF], we show that the solution (X(0), A(0)) of ( 22)

dX t = √ 2A t dW t dA t = -A α t dL 0 t (X),
with initial conditions X 0 = 0 and A 0 > 0, can be be obtained as a limit solution of (1) when λ > 0 and decreases to 0.

Proposition 5.1. Let α > 0 be fixed. Denote by (X(λ), A(λ)) the solution of (1) for λ > 0 with initial condition (0, A 0 ). As λ decreases to 0, (X(λ), A(λ)) converges in probability uniformly on compact sets to the solution of [START_REF] Revuz | ume 293 of Grundlehren der Mathematischen Wissenschaften[END_REF].

Nevertheless this continuity property at λ = 0 of the solution of ( 1) is compatible with the fact that for α < 1/2, in view of the nature of the cell arrest, 0 is a critical value for λ (Proposition 2.3).

Proof. For λ > 0 denote by Y (λ) = (Y t (λ), t ≥ 0) the solution of (2) starting from 0. One easily checks that a.s. Y (λ) converges to W pointwise, uniformly on compacts, as λ tends to 0. Denote by (L 0 t (W ), t ≥ 0) (resp. (L 0 t (Y (λ)), t ≥ 0)) the local time process at 0 of W (resp. Y (λ)). Let τ be a stopping time for W with finite second moment. Using Tanaka's formula and Burkholder-Davis-Gundy inequalities (see e.g. [START_REF] Revuz | ume 293 of Grundlehren der Mathematischen Wissenschaften[END_REF]), one obtains

IE sup 0≤s≤τ L 0 s (Y (λ)) -L 0 s (W ) 2 → λ→0 0.
Consequently, on has:

(23) (L 0 t (Y (λ), Y t (λ)), 0 ≤ t ≤ τ → λ→0 (L 0 t (W ), W t ), 0 ≤ t ≤ τ
in probability for the infinite norm of IR 2 -valued continuous functions on compact sets.

If α ≥ 1/2, then we choose for stopping time τ any deterministic time t > 0 and using Theorem 2.1 we obtain the convergence for the infinite norm in probability of ((X s (λ), A s (λ)), 0 ≤ s ≤ t) to ((X s , A s ), 0 ≤ s ≤ t) . If α < 1/2, then we choose

τ = inf    t ≥ 0 : L 0 t (W ) > A 1 2 -α 0 √ 2( 1 2 -α)    .
In this case the solution of ( 22) satisfies T 0 (A) = τ (Theorem 2.1 for λ = 0). Similarly to the previous case, we obtain the convergence (23) but since on {T 0 (A(λ)) < τ }, A t (λ) is not defined for t ∈ (T 0 (A(λ)), τ ) we can not immediately conclude. One has:

(T 0 (A(λ)) ≥ τ ) =   τ 0 e -λ( 1 2 -α) dL 0 s (Y (λ)) ≤ A 1 2 -α 0 √ 2( 1 2 -α)   ,
hence (1 {T 0 (A(λ))≥τ } (X t (λ), A t (λ)), 0 ≤ t ≤ τ ) converges uniformly in probability to the solution of [START_REF] Revuz | ume 293 of Grundlehren der Mathematischen Wissenschaften[END_REF]. Moreover for any λ > 0, using [START_REF] Jenkins | Local macrophage proliferation, rather than recruitment from the blood, is a signature of th2 inflammation[END_REF] we can write

Conclusion

In this work, we propose and study a model for macrophage dynamics taking into account their lipid load. We incorporated essential parameters: intensity of inflammation α, immune defenses and medical treatments λ, and the initial inflammatory state (X 0 , A 0 ). Our mathematical study provides information on the dynamics of the macrophage and hence on the atherosclerotic plaque dynamics. In particular, depending on the value of the parameters, it can lead to the accumulation of macrophages in a particular region of the arterial wall. The key quantity to evaluate this risk is IP [T 0 (A) < ∞].

Figure 1 :

 1 Figure 1: A diagram showing coordinates at the scale of a single cell where x=0 denotes a point in the neighbourhood of the macrophage where lipids are abundant and can be ingested by the macrophage.

Figure 2 :

 2 Figure 2: Sketch of the different stages of atherosclerosis plaque formation: (1) diffusion of a "free" macrophage cell; (2) upon entering a localized lipid-enriched region, the macrophage accumulates lipids, and thereby becomes less mobile; and(3) after many visits to the lipid-enriched region, the macrophage eventually gets trapped, resulting in the formation of an atherosclerotic plaque.

0 e

 0 λu(t)) 1-α + λ(1 -α) λs) 2-α ds, t ≥ 0 where W is a real Brownian motion starting from 0. Thanks to[START_REF] Tall | Cholesterol, inflammation and innate immunity[END_REF] one has λs) 2-α ds. The scaling property of W yields λ(1 -α) )s dL 0 s (Y (λ)) = +∞ in probability, i.e. for every c > 0, lim λ→0 IP ∞ -λ( 1 2 -α)s dL 0 s (Y (λ)) < c = 0.

  Case α > 1/2 : In this case β < 0. Set γ = -β = α -1/2 > 0. We haveB t = e -λγt B 0 +

									+∞,
	which implies that lim t→∞ A t = +∞ and
		lim inf t→∞	X t = -∞, lim sup
					√	2γ	t	e λγs dL 0 s (Y ) ,
									0
	hence B t > 0 for every t ≥ 0.
	The process A reaches 0 in finite time if B takes an infinite value in finite time.
	Using (15) we see that	B t ≤ B 0 +	√	2γL 0 t (Y ).
	Obviously (L 0 t (Y )) t≥0 can not blow-up in finite time. Hence IP [T 0 (A) < ∞] = 0 and
	(18)	A t = e λt B 0 +	√	2γ	t	e λγs dL 0 s (Y )	-1/γ	, t ≥ 0.
					0
	We have: B t ∼ t→∞ To compute t 0 e λγs dL 0 √ 2γe -λγt t 0 e λγs dL 0 s (Y ). s (Y ), we use (12), make the time change v(s) = 1 λ log(1 + λs),
	and then an integration by part to finally obtain:
	t
	(19)	e λγs dL 0 s
	0

t→∞ X t = +∞.

  ) .

	Recalling Strassen's Law of the Iterated Logarithm
	lim sup t→∞	Z t 2t log 2 t	= 1 a.s.

we obtain: lim sup t→∞ |X t | = +∞, which leads by symmetry and continuity to lim inf t→∞ X t = -∞ and lim sup t→∞ X t = +∞. If 1/2 < α < 1, then

[START_REF] Nathan | Nonresolving inflammation[END_REF]

leads to

Here are the main conclusions that one can reach adopting this model.

1 -If α ≥ 1/2, which corresponds to a relatively low intensity of inflammation, no atherosclerotic plaque will form whatever λ, X 0 and A 0 are (equivalently: IP [T 0 (A) < ∞] = 0). Unsurprisingly, in this case the introduction of the parameter λ in the system [START_REF] Revuz | ume 293 of Grundlehren der Mathematischen Wissenschaften[END_REF] previously studied in [START_REF] Meunier | Long time behavior in locally activated random walks[END_REF], does not change IP [T 0 (A) < ∞] which remains null.

2 -If α < 1/2, which corresponds to a relatively high intensity of inflammation, then the introduction of the parameter λ affects the dynamics of the particle in several ways. Indeed for λ = 0, whatever the value of A 0 , one has IP [T 0 (A) < ∞] = 1, which means cell arrest or apoptosis in finite time a.s. While for λ > 0, if A 0 is large enough (which means low intracellular load in lipids at time t = 0) with X 0 = 0, then T 0 (A) = ∞ and lim t→∞ A t = +∞. This infinite limit can be interpreted as the emigration of the macrophage and corresponds to no contribution from this cell to the formation of an atherosclerotic plaque at the origin. Moreover one can find an upper bound the risk IP [T 0 (A) < ∞] by using Markov inequality:

Since V α,λ is independent of A 0 , we clearly see that by increasing A 0 , one increases the probability of avoiding arrest or apoptosis. But X 0 also plays a role in the case where α < 1/2 and λ > 0. Using Proposition 2.2, we see that cell arrest (or apoptosis) occurs in finite time for the solution of (1) with initial condition (X 0 , A 0 ) if and only if it is so for the solution of (1) with initial condition (0, Ã0 ). The variable Ã0 increases stochastically with A 0 and |X 0 |. Thus, given N , the state of the arterial wall can be improved by increasing A 0 or |X 0 |. This corresponds to a lower intracellular load in lipids at time t = 0 or a greater dispersion of lipid sources at time t = 0, which means a lower initial inflammatory state.

Hence the introduction of λ not only offers the possibility to lower IP [T 0 (A) < ∞] in the cases when α, the inflammatory intensity, is high but this possibility takes into account the initial inflammatory state (X 0 , A 0 ).

Model (1) could be enriched by considering a population of macrophages instead of a single one. The role played by the accumulation of trapped macrophages on the still moving macrophages should be well understood to be taken into account by the enriched model.

Competing Interests

The authors declare that no competing interests exist.