Jonathan Stéphano

Benjamin Mauroy
email: benjamin.mauroy@univ-cotedazur.fr

Wall shear stress distribution in a compliant airway tree

The airflow in the bronchi applies a shear stress on the bronchial mucus, which can move the mucus. The air-mucus interaction plays an important role in cough and in chest physiotherapy (CP). The conditions under which it induces a displacement of the mucus are still unclear. Yet, the air-mucus interaction justifies common technics of CP used to help the draining of the mucus in prevalent diseases. Hence, the determination of the distribution of the shear stress in the lung is crucial for understanding the effects of these therapies and, potentially, improve their efficiency.

We develop a mathematical model to study the distribution of the wall shear stress (WSS) induced by an air flow exiting an airway tree. This model accounts for the main physical processes that determines the WSS, more particularly the compliance of the airways, the air inertia and the tree structure.

We show that the WSS distribution in the tree depends on the dynamics of the airways deformation and on the air inertia. The WSS distribution in the tree exhibits a maximum whose amplitude and location depend on the amount of air flow and on the "tissue" pressure surrounding the airways. To characterize the behavior of the WSS at the tree bifurcations, we derive new analytical criteria related to the airway size reduction in the bifurcations.

Our results suggest that a tuning of the airflow and of the tissue pressure during a CP maneuver might allow to control, at least partially, the air-mucus interaction in the lung.

Introduction

The lung forms an interface between the ambiant air and the blood. It is an easy entry point in the organism for toxic or infectious particles. As a protection, the walls of the bronchi are covered by a mucus layer. The mucus captures the inhaled particles and is incessantly moved toward the oesopahryngeal region by the mucociliary clearance [START_REF] Ximena | Cilia and Mucociliary Clearance[END_REF]. Once in the oesopharyngeal region, the mucus is either swallowed or expelled by coughing. In the non-mature or pathological lung, the mucociliary clearance might not be able to drain the mucus correctly and the lung relies on the cough to drain the mucus [START_REF] Ho | The effect of aging on nasal mucociliary clearance, beat frequency, and ultrastructure of respiratory cilia[END_REF][START_REF] Kurian Baby | Effect of cigarette smoking on nasal mucociliary clearance: A comparative analysis using saccharin test[END_REF]. The mucus is a complex viscoelastic material that moves as a fluid only when its inner stress overcomes its yield stress [START_REF] Samuel | Micro-and macrorheology of mucus[END_REF]. During lung's ventilation [START_REF] Stephano | Modeling shear stress distribution in a deformable airway tree[END_REF], during cough [START_REF] King | The role of mucus viscoelasticity in cough clearance[END_REF] or during chest physiotherapy [START_REF] Mauroy | Toward the modeling of mucus draining from the human lung: role of the geometry of the airway tree[END_REF][START_REF] Mauroy | Toward the modeling of mucus draining from human lung: role of airways deformation on air-mucus interaction[END_REF][START_REF] Stephano | Modeling shear stress distribution in a deformable airway tree[END_REF], the airflow interacts with the mucus and applies a shear stress that can potentially overcome the mucus yield stress. Chest physiotherapy is commonly used to compensate the dysfunctions of the mucocilliary clearance and of the cough occurring in prevalent diseases such as asthma, bronchiolitis, COPD or cystic fibrosis [START_REF] Pryor | Physiotherapy for airway clearance in adults[END_REF]. However, the distribution of the shear stress and the region in the lung where the mucus might be liquified are not well identified as of today. Hence, we propose in this work to model the airflows in the bronchial tree to study the distribution of the wall shear stress.

Past studies have developed models of the air fluid dynamics in the lungs. However, the wide range of scales covered by the airways makes lung's modelling a challenge [START_REF] Crampin | Computational physiology and the physiome project[END_REF][START_REF] Longest | Use of computational fluid dynamics deposition modeling in respiratory drug delivery[END_REF]. There are two main approaches. The first approach is to model the fluid dynamics in a subpart of the lung using 3D geometries, either idealised or reconstructed from CT-scans [START_REF] Mauroy | Interplay between geometry and flow distribution in an airway tree[END_REF][START_REF] Mauroy | 3D Hydrodynamics in the Upper Human Bronchial Tree: Interplay between Geometry and Flow Distribution[END_REF][START_REF] Kabilan | Characteristics of airflow in a CTbased ovine lung: a numerical study[END_REF][START_REF] Moghadas | Numerical investigation of septal deviation effect on deposition of nano/microparticles in human nasal passage[END_REF][START_REF] Paz | Glottis effects on the cough clearance process simulated with a CFD dynamic mesh and Eulerian wall film model[END_REF][START_REF] Farnoud | Large eddy simulations of airflow and particle deposition in pulsating bidirectional nasal drug delivery[END_REF][START_REF] Shuai Ren | Numerical Analysis of Airway Mucus Clearance Effectiveness Using Assisted Coughing Techniques[END_REF]. Some studies go further and model the air-mucus interaction in 3D geometries using the volume of fluids method [START_REF] Kumar | Numerical model of air-mucus interactions in a bronchus[END_REF][START_REF] Rahul R Rajendran | Mucus transport and distribution by steady expiration in an idealized airway geometry[END_REF] or the thin layer theory [START_REF] Paz | Glottis effects on the cough clearance process simulated with a CFD dynamic mesh and Eulerian wall film model[END_REF][START_REF] Shuai Ren | Numerical Analysis of Airway Mucus Clearance Effectiveness Using Assisted Coughing Techniques[END_REF]. Such studies are mainly focused on the upper respiratory region and assume the geometries to be rigid. Micro scale models have also been developed to study mucociliary clearance [START_REF] Lubkin | A Viscoelastic Traction Layer Model of Muco-Ciliary Transport[END_REF][START_REF] Jayathilake | Three-dimensional numerical simulations of human pulmonary cilia in the periciliary liquid layer by the immersed boundary method[END_REF] or the respiratory gaz flows in the deep lung (acini) [START_REF] Felici | Diffusional screening in real 3D human acini-a theoretical study[END_REF][START_REF] Swan | Evidence for minimal oxygen heterogeneity in[END_REF]. The models are analysed with high-end computational fluid dynamics based on finite elements or finite volumes, which are not well adapted to large change of scales. In such models, good boundary conditions to mimic the response of the regions non accounted for are not trivially determined. Nevertheless, the resulting predictions are often The number on the airways corresponds to the generation index of that airway. The generation index corresponds to the number of airways between the airway considered and the root of the tree. From the generation i to the generation i`1, the diameter of the airways is multiplied by a ratio h i ă 1.

very detailed and give rich insights on the local dynamics. The second approach, which we will use in this work, is based on idealised representations of the lungs that reduce the behavior of groups of airways to a single behavior. Typically the airways are regrouped by generations, i.e. the number of bifurcations between the airway considered and the root of the tree [START_REF] Mauroy | An optimal bronchial tree may be dangerous[END_REF][START_REF] Mauroy | Toward the modeling of mucus draining from the human lung: role of the geometry of the airway tree[END_REF][START_REF] Mauroy | Toward the modeling of mucus draining from human lung: role of airways deformation on air-mucus interaction[END_REF][START_REF] Stephano | Modeling shear stress distribution in a deformable airway tree[END_REF]. These approaches allow to mimic the whole airway tree and have led to interesting fundamental insights on the lung [START_REF] Lambert | A computational model for expiratory flow[END_REF][START_REF] Mauroy | An optimal bronchial tree may be dangerous[END_REF][START_REF] Mauroy | Toward the modeling of mucus draining from the human lung: role of the geometry of the airway tree[END_REF][START_REF] Mauroy | Toward the modeling of mucus draining from human lung: role of airways deformation on air-mucus interaction[END_REF][START_REF] Noël | Interplay Between Optimal Ventilation and Gas Transport in a Model of the Human Lung[END_REF][START_REF] Noël | The origin of the allometric scaling of lung's ventilation in mammals[END_REF]. However, in either approach, no model did include both the role of the air inertia and of the airways compliance, although these are known to affect notably the flow properties in the lung [START_REF] Mauroy | Interplay between geometry and flow distribution in an airway tree[END_REF][START_REF] Mauroy | Toward the modeling of mucus draining from human lung: role of airways deformation on air-mucus interaction[END_REF]. Indeed, the airways diameters depend on the airways inner air pressure, the airways inner air pressure depends on the air fluid dynamics, potentially affected by inertia, and the air fluid dynamics depends on the airways diameters. These interactions occur in a tree-like structure which increases the complexity of the system [START_REF] Mauroy | An optimal bronchial tree may be dangerous[END_REF]. In this framework, we propose an original model to study for the first time the physics of an inertial flow in an airway tree with compliant airways. Because the mucus thickness is far smaller than the airway diameters [START_REF] Karamaoun | New insights into the mechanisms controlling the bronchial mucus balance[END_REF], we neglect the mucus layer and assume that the shear stress induced by the air on the airways wall is a good estimation of the shear stress felt by the mucus in the lung.

For the first time, a study is able to propose a detailed physical analysis of the interactions between the airflows and the deformation of the airways in the tree structure. Our model predicts the existence of a maximum of wall shear stress (WSS) within the tree. We uncover a new analytical criterion that allows to explain the existence of this maximum. This criterion is based on the physics occurring in the bifurcations and allows to determine whether the shear stress will increase or decrease when the air is going through a bifurcation. The criterion depends on the geometries of the deformed airways and on the amount of air inertia. We show that the location of the maximum of WSS is a function of the total airflow rate in the tree and of an homogeneous "tissue pressure" that surrounds the airways. Modulating these two quantities allows, to some extent, to control the location and amplitude of the maximum of wall shear stress.

This model can be considered as an idealised representation of the lung, with the tissue pressure being a representation of the pressure in the lung's parenchyma. In that framework and in the limit of our hypotheses, our model brings for the first time a better understanding of the shear stress induced by the air on the mucus in the lung. Moreover, our work highlights the underlying biomechanics involved in high ventilation regimes and in some of the most common chest physiotherapy technics. This study might pave the way for future improvements of the manipulations of chest physiotherapy.

Model

The airflow inside the airway induces an air pressure that opposes the mechanical pressure that surrounds the airways and that we will call from now on the tissue pressure. The balance between these two pressures affects the diameters of the deformable airways, which in turn affects the airflow inside the airway.

Model of the bronchial tree

We mimic the bronchial tree using the model of the human lung from [START_REF] Lambert | A computational model for expiratory flow[END_REF]. The bronchial tree is represented by a cascade of bifurcating cylinders with radius R and length L representing the bronchial airways. The generation of an airway is an index that counts the number of bifurcations on the path between the root of the tree, that mimics the trachea (index 0), and the airway under consideration. Each cylinder has a fixed length and a variable diameter that is adjusted depending on the transmural pressure that it is submitted to. The transmural pressure of an airway is the difference between its inner air pressure and the pressure surrounding the airways, the tissue pressure. In our model, the tissue pressure is assumed homogeneous and the air pressure results from the air flow inside the airway. The bifurcations are assumed symmetric, consequently all the cylinders with the same generation index have exactly the same physical and geometrical properties. Hence, it is sufficient to study one single airway per generation. The ratio between the diameter of a cylinder in the generation i `1 and a cylinder in the generation i is called the reduction ratio and is denoted h i .

In [START_REF] Lambert | A computational model for expiratory flow[END_REF], the authors propose a model for static compliant airways based on bronchi data. In a generation with index i, this model links a normalized section of the airways α i to the transmural pressure P i . The section is normalized relatively to a maximal surface area S max,i reached for an infinite transmural pressure. Since we hypothesize that the airways remain cylindrical whatever the transmural pressure, the model can be expressed in terms of the airways radii. Hence, the maximal surface area S max,i corresponds to a maximal radius R max,i and S max,i " πR 2 max,i . Finally, the radius R i of an airway in the generation i is related to the normalised surface area of the airway by R i " R max,i ? α i . The reformulated Lambert's model writes

$ ' ' & ' ' % α i " P i ď0
α 0,i p1 ´Pi {P í q ´ní α i "

P i ě0 1 ´p1 ´α0,i q p1 `Pi {P ì q ´nì R i " R max,i ? α i (1)
The quantities α 0,i , P ȋ , n ȋ and the maximal possible radius for the airway in generation i, R max,i , are positive fixed data given by the model. They depend on the generation index i only.

Since an healthy mucus layer is very thin, from about 3 micrometers to 10 micrometers [START_REF] Karamaoun | New insights into the mechanisms controlling the bronchial mucus balance[END_REF], its influence on the airways radii and on the air fluid dynamics can be neglected.

Fluid dynamics in the airways

We consider a cylindrical airway with a radius R and a length L. We will assume the velocity profile to be stationary, axisymmetric, invariant along the airway axis and parallel to the airway axis. A [START_REF] Lambert | A computational model for expiratory flow[END_REF]. The curves represent the normalised surface areas of the airways of a same generation as a function of the transmural pressure in these airways. The normalisation is made relatively to a maximal possible surface area of the airways. The transmural pressure is the difference between the inner (air) and outer pressures (surrounding tissue) of the airways. The boxed index indicates the generation of the corresponding airways. volumetric integration of the Navier-Stokes equations with the above hypotheses leads to a relationship between the wall shear stress σ and the pressure drop q through the airway,

σ " ´q R 2L (2)
Details of the computation are given in the Supplementary Materials (A). Since we are studying the outgoing flow that is positive, the pressure drop is negative. Hence, the wall shear stress is positive.

For each airway, the pressure drop q between its two extremities is linked to the flow rate Φ going through the airway and to the hydrodynamic resistance R of the airway,

q " R Φ (3)
The resistance depends on the airway's geometry and on the air fluid dynamics. In the case of a Poiseuille flow -i.e. slow, fully developed regime-, the hydrodynamic resistance R P depends only on the airway geometry,

R P " 8µL πR 4 (4)
The quantity µ is the viscosity of the air. However, in higher regimes, the air flow is affected by the inertial effects and the hydrodynamic resistance becomes larger than in Poiseuille's regime, R " Z R P where the prefactor Z is larger than 1.

Various models have been proposed to express the prefactor Z [START_REF] Pedley | The prediction of pressure drop and variation of resistance within the human bronchial airways[END_REF][START_REF] Florens | The optimal branching asymmetry of a bidirectional distribution tree[END_REF]. We choose to use the velocity profile proposed in [START_REF] Merryn | Modeling Water Vapor and Heat Transfer in the Normal and the Intubated Airways[END_REF]. It accounts for inertia using a power law. The axial velocity depends on the position r on the radius of the airway,

v z prq " v 0 ´1 ´´r R ¯α¯w ith α " max ˆ2, Re 150
ẇhere Re is the Reynolds number Re " 4ρΦ µπR . Then, we can compute the flow rate Φ through the airway, Φ " α α`2 πv 0 R 2 , and using the relationship between the wall shear stress σ " µ dvz dr pRq " ´µαv 0 {R and the pressure drop q, σ " ´qR{p2Lq, we can deduce that q " α`2 4 R P Φ. Hence, the expression of Z in our model is

Z " max ˆ1, 1 2 `Re 600
Ċonsequently, the resistance factor Z has a lower bound equal to 1 that is reached for Re " 300 and that corresponds to the Poiseuille's regime.

Finally, we set the density of the air to ρ " 1.225 kg.m ´3, and its viscosity to µ " 1.70 10 ´5 Pa.s. The reference pressure is set to 0 at the opening of the first generation, that corresponds to the trachea. Notice that, with this hypothesis, we do not account for the nasopharyngeal or buccopharyngeal pathway to the air pressure. Actually, the error made with this hypothesis may be of various importance, depending on the flow rate and on the sizes of the airways.

Resolution of the model's equations and validation.

The global model leads to a set of non linear equations for each generation of the airway tree. The generations are coupled by the conservation of the air flow and by the continuity of the air pressure. The equations are solved using a Newton algorithm, implemented in C++ with the library Eigen [33]. The set of equations, its properties and the numerical process are detailed in the Supplementary Materials (B, C and D).

To validate our model, we ran it with simple configurations whose behaviors were either straightforward or predicted by the analytical results presented in the section 3.1. Hence, we tested trees with fractal geometries and rigid airways as in [START_REF] Mauroy | An optimal bronchial tree may be dangerous[END_REF], trees at different scales reduced to a single deformable airway or to a single bifurcation. 3 Links between the tree geometry and the wall shear stress distribution

Inertial and viscous critical ratios

Our model predicts that the wall shear stress exhibits a bell distribution with a maximal value somewhere in the airway tree. On each side of the maximum, the wall shear stress varies exponentially, increasing in the proximal airway and decreasing in the distal airway, as shown in figure 3A. Actually, we will show that this behavior results from the geometry of the airway tree. In a bifurcation, the radius R i`1 of the daughter airways is a reduction of the radius R i of the mother airway, R i`1 " h i R i . The ratios h i allow to derive a scaling for the mean velocity in the airway v " Φ{pπR 2 q, with v i`1 " 1 2h 2 i v i , and for the Reynolds number Re " 4ρΦ{pµπRq, with Re i`1 " 1 2h i Re i . For Poiseuille's regime (low flow), we have shown that the change of the wall shear stress through a bifurcation determines the behavior of the wall shear stress at the scale of the airway tree, see [START_REF] Stephano | Modeling shear stress distribution in a deformable airway tree[END_REF]. Indeed, in Poiseuille's regime, the change of the wall shear stress σ p " 4µΦ{pπR 3 q through a bifurcation scales as

σ p i`1 " 1 2h 3 i σ p i .
Depending on the position of the ratio h i relatively to the viscous critical ratio `1 2 ˘1 3 » 0.793, the wall shear stress can either increase or decrease through the bifurcation. The different scenarios for the wall shear stress and for other physical quantities are given in table 2. For example, in the case of a constant ratio h i " h throughout the airway tree, the radii R i , the velocities v i and the mean wall shear stresses σ p i will either grow exponentially, remain constant or decrease exponentially.

In the case studied in this paper, we also account for the inertial effects by deforming the parabolic velocity profile of the linear regime into a plug-like velocity profile that depends on the Reynolds number [START_REF] Merryn | Modeling Water Vapor and Heat Transfer in the Normal and the Intubated Airways[END_REF]. The plug like profile induces higher wall velocity gradients than the parabolic profile. Consequently, inertial effects increase the pressure drops and the wall shear stresses in the airways. In the presence of inertia, the shear stress decomposes into σ i " σ p i {2 `σI i where σ I i reflects the influence of the inertia,

σ I i " $ & % σ p i 2 " 2µ π Φ i R 3 i if Re i ď 300 8ρ 600π 2 Φ 2 i R 4 i if Re i ą 300
When inertia is dominant (Re i ąą 1), the change in wall shear stress between the generations i and i`1 depends on the location of h i relatively to the inertial critical ratio `1 2 ˘1 2 , since

Φ 2 i`1 R 4 i`1 " ´1 2h 2 i ¯2 Φ 2 i R 4 i .
When inertia is small (Re ď 300 in our model), the change in wall shear stress between the gener- Distribution of the radii on the pressure-radius relationships of the airways [START_REF] Lambert | A computational model for expiratory flow[END_REF] Figure 3: A, B: Shear stress distribution along the airway tree as a function of the generation index, for different airflows Φ and tissue pressures P tissue (A: with inertia, B: without inertia). The dashed horizontal lines represent an estimation of the the yield stress for healthy mucus (1 Pa) and for cystic fibrosis mucus (15 Pa). The red dots indicates regions where inertia is large, i.e. Z ď 2 which is equivalent to σ i ą 2σ i P . C: Radius reduction ratios h i between two successive generations of the airway tree as a function of the generation index i, for different airflows Φ and tissue pressures P tissue . The ratios exhibit up to three different regimes. D: Radii of the airways located on their compliance curve as a function of their transmural pressure, for a subset of the cases tested in this study. One symbol represents a specific configuration of air flow and pressure. The horizontal shift between two generations corresponds to the air pressure drop between these geenrations. reduction ratio h i airway radius R i mean air velocity v i Reynolds number Re i wall shear stress (Poiseuille) σ p i wall shear stress (inertia) σ I i (Re i ą 300)

1 2 b 1 2 3 b 1 2 1 OE OE OE OE OE OE OE Ñ Õ Õ Õ Õ Ñ OE OE OE OE OE Õ Ñ OE OE OE OE OE OE OE Õ Õ Õ Õ Õ Ñ OE OE OE Õ Õ Õ Ñ OE OE OE OE OE
Table 2: Variations of the core quantities involved in our study from one generation i to the next generation i `1, depending on the value of the ratio h i between the radius of the airways in generation i `1 and the radius of the airways in generation i. The upwards arrows indicate that the quantity is increasing from the mother airway to the daughter airways, the horizontal arrow that it remains constant and the downward arrows that it is decreasing. ations i and i `1 is that of the Poiseuille regime and, as detailed above, it depends on the location of h i relatively to the viscous critical ratio `1 2 ˘1 3 . The Reynolds number is multiplied from the generation i to the next by the ratio 1 2h i . Hence, the influence of inertia on the shear stresses is also changing from one generation to the next with the same rate. Since h i is in general larger than 1 2 , the influence of inertia is decreasing along the generations of the tree and disappears in the distal airways.

A bell-shaped distribution of the wall shear stress

The bell shape exhibited by the shear stresses in the airway tree results from the interaction of the flow with the geometry of the airway tree. This interaction is driven by the reduction ratios h i that correspond to the relative change of the radii between one generation and the next. If we compute these ratios during idealised manipulations, shown in figure 3B, we can see that the bifurcations along which the shear stress is increasing correspond to the low values of h i and the bifurcations along which the shear stress is decreasing correspond to the high values of h i .

In the case of a low air flow rate (0.1 L.s -1 ), inertia is low and the behavior is controlled by the location of the ratios h i relatively to the viscous critical ratio `1 2 ˘1 3 . In this configuration, the generation at which the maximum of the shear stress is reached is exactly where the transition

h i " `1 2
˘1 3 occurs, see figure 3A and B (0.1 L.s -1 ), table 2 and [START_REF] Stephano | Modeling shear stress distribution in a deformable airway tree[END_REF]. In the presence of inertia, the shear stress in the proximal airways is increased relatively to the Poiseuille regime, see figure 3A andC. The shape of the wall shear stress distribution and the location of its maximum result from the relative contribution and evolution along the generations of the two terms in the decomposition of the shear stress, σ i " σ p i {2 `σI i . For the large flows, typically 4 L.s -1 , the shear stress in the proximal airways is dominated by inertia, i.e. σ I i ąą σ p i . Then, the maximum of the shear stress is located where h i exceeds the inertial critical ratio `1 2 ˘1 2 » 0.707, see figures 3A and B (4 L.s -1 ) and table 2.

Role of the fluid dynamics 4.1 Three regions, three behaviors

The distribution of the ratios h i along the generations of the tree can be decomposed into three regions, see figure 3B: a proximal region, where the diameters reduction ratios are low; A transitional region, where the diameters reduction ratios increase; A distal region, not present at low flow rates, where the diameters reduction ratios are high. The three regions arise from different resolutions of a trade-off between the air pressure in the airways and the air pressure drops in the airways bifurcations.

Proximal region. The air pressure drops in the large airways are low, even in the pres-ence of inertia. Hence, the airway diameters are mainly affected by the tissue pressure. As the daughter airways in a bifurcation are more compliant than the mother airway, they are more constricted. Consequently, their radii reduction ratios are smaller than the ratios at functional residual volume (FRC) shown in table 1. Hence, the ratios tend to be close to the inertial critical ratio `1 2 ˘1 2 » 0.707 and smaller than the viscous critical ratio `1 2 ˘1 3 » 0.793. Consequently, in the proximal region, the shear stress is in general increasing from one generation to the next.

Transitional region. The diameters of the airways are decreasing along the generations, see figure 3D. Hence, the pressure drops are increasing, at least in the first generations of the transitional region. The change of the transmural pressure through a bifurcation starts to compensate the lower compliance of the daughter airways. The daughter airways are actually less compressed relatively to the mother airway than in the proximal region. This results in higher ratios h i . This effect becomes stronger and stronger along the generations, and the ratios increase steadily. Eventually, they exceed `1 2 ˘1 2 and `1 2 ˘1 3 , and the maximum of the shear stress and of the pressure drop is reached.

Distal region. In the distal region, the transmural pressure is slightly positive, at least in the cases that we tested. The small transmural pressures in the distal region corresponds to near maximal compliances values, as they correspond to large slopes zones on the pressure-section curves shown on figure 3D. Moreover, since the smaller airways have larger compliances and higher inner air pressures, they are, relatively, more dilated than the larger ones. As a consequence, the reduction ratios are relatively high in the distal region, reaching values near 0.9 or higher.

Dependance on air flow rates and tissue pressures

The distribution and amplitude of the shear stresses in the airway tree depend on the flow rate Φ and on the tissue pressure P tissue , more particularly the position and amplitude of its maximum, see figure 3. The shear stress depends on the ratio between the air flow and the cube of the radius of the airway, σ " α`2

4 8µ π |Φ| R 3 .
An increase of the air flow rate induces higher pressures in the airways, which are more dilated. The variation of the shear stress is the result of the balance between the air flow rate increase and the increase of the cube of the airways' radii. Since our model predicts that the amplitude of the maximal shear stress increases with the air flow rate, this suggests that the increase of the air flow rate might be dominant over the increase of the airways' radii, at least in the range of the parameters tested. In this configuration, pressure drops in the airways are larger because of the higher flows and the air pressure in the airways increases more quickly along the generations. Hence, both the proximal and the transitional regions are covering less generations. The location of the maximum of the shear stress is then shifted toward the proximal generations of the airway tree.

When the tissue pressure increases at constant air flow rate, the amplitude of the shear stresses also increases, see figure 3. This phenomenon is due to the resulting decrease of the radii of the airways. Hence, the pressure drops between two successive generations are increased. In the proximal region, the reduction ratios are smaller due to the increase of the compliance with the generation index. However, the inner air pressures needed to increase the ratios are also higher. As long as these inner air pressures remain low, the transmural pressures change very little. As a consequence, the range of generations covered by the proximal region remains the same. The two other regions are however affected by the larger pressure drops. In the transitional region, the reduction ratios increase more quickly, hence reducing the range of generations covered by the region. In the distal region, the air pressure is high as it results from the addition of all the pressure drops. The airways are then closer to their maximal dilation, and for high tissue pressures the reduction ratios are near 0.9, see table 1.

Our model predicts that the location in the airway tree of the maximum of the shear stress depends only slightly on the tissue pressure. This phenomenon results from the variation of the reduction ratios in the transitional region. In that region, if the tissue pressure is increased, the reduction ratios are increasing from values that are lower, as the higher tissue pressure decreases the ratios in the proximal region. Nevertheless, in the transitional region, the increase of the ratios from one generation to the next occurs at a higher pace, due to the higher pressure drops. This results in an increase in reduction ratios that compensates the lower values in the proximal region. This allows to reach the critical value of the ratio, which corresponds to the maximum of the shear stress, in a generation index close to that of a lower tissue pressure.

Model limitations

Our predictions should be considered as qualitative only, since our model is based on several simplification hypotheses. First, both the fluid dynamics and the airways deformation [START_REF] Lambert | A computational model for expiratory flow[END_REF] are considered stationary. Moreover, we assume that the air velocity profiles in the airways are not affected by the bifurcations and are unchanged along the axis of the airways. Turbulence was neglected as it can scarcely be established in the airways, with the notable exception of the trachea [START_REF] Lin | Characteristics of the turbulent laryngeal jet and its effect on airflow in the human intra-thoracic airways[END_REF]. Actually, in most of the airways, the air flow reaches the bifurcation before the turbulence could really develop and have some influence on the pressure drops. Also, neglecting turbulence allows to get a more tractable model. Nevertheless, a more precise evaluation of the role of turbulence will be made in future evolutions of our model. In term of geometry, the airway tree mimicking the lung in our model has symmetric bifurcations, unlike the lung where asymmetric branching is known to affect the fluid dynamics [START_REF] Mauroy | An optimal bronchial tree may be dangerous[END_REF][START_REF] Buess | Energy dissipation in an asymmetric tracheobronchial tree: a CFD analysis[END_REF]. Also, our airways remain perfectly cylindrical whatever their transmural pressure. The typical star shape of the lumen of a constricted airway is not accounted for [START_REF] Mauroy | Role of bronchus morphology on mucus mobilization during chest physiotherapy[END_REF]. A total collapse of the airways is also not possible in our model since the air flow rate is constrained. The computation of the transmural pressure is an approximation, as we assume an homogeneous tissue pressure all throughout the lung. Finally, we neglect the pressure drops induced by the bifurcations and by the nasopharyngeal or buccopharyngeal pathway. Actually, only a full 3D analysis of the different regimes could bring adequate estimations of these pressure drops [START_REF] Paz | Glottis effects on the cough clearance process simulated with a CFD dynamic mesh and Eulerian wall film model[END_REF][START_REF] Longest | Use of computational fluid dynamics deposition modeling in respiratory drug delivery[END_REF]. Such a study is out of the scope of this paper but should be performed in the near future in order to reach better estimations of the pressure drops.

Nevertheless, our model includes the main physical and physiological phenomena affecting the wall shear stress in the lung: the airway tree structure with deformable airways and the resulting fluid dynamics, including some inertial effects. Hence, the predictions of our model should highlight characteristic behaviors for the wall shear stress distributions in airway trees such as the lung.

Discussion

Our work predicts the distribution of the wall shear stress induced by an air flow exiting an airway tree surrounded by an external pressure. The wall shear stress exhibits a bell-shaped distribution along the generations of the tree. This shape results from the complex, non-linear interactions between the fluid dynamics, the geometry of the tree and the compliance of the airways. Our analyses allow to uncover the dynamics of these interactions. We show that the change in diameter of the airways through a bifurcation determines if the wall shear stress increases or decreases through that bifurcation. We derived analytical criteria that characterize specific regimes. They are based on the ratio h i between the diameters of the daughter airways and the mother airway, see table 2. In Poiseuille's regime, the wall shear stress increases through a bifurcation only if h i is larger than the viscous critical ratio `1 2 ˘1 3 » 0.793 and decreases otherwise. In regimes where inertia is dominant, the inertial critical ratio is lower and equal to `1 2 ˘1 2 » 0.707. When both viscous and inertial influences are of the comparable order of magnitude, then the corresponding critical ratio stands between the viscous and the inertial ones.

Our results give insights about the air-mucus interactions at expiration, not only during normal ventilation but also during chest physiotherapy. Chest physiotherapy uses the air wall shear stress to help the draining of the mucus from the bronchial tree. But in order for mucus to be motioned, the mucus yield stress has to be overcome. Our model suggests that a maximum of the wall shear stress can be, to some extent, localized in the airway tree by tuning the tissue pressure and the air flow rate. The existence of an adjustable maximum of wall shear stress by accessible physical parameters is of great importance for validating scientifically the technics used by chest physiotherapists. Our results suggest that the chest physiotherapists might intuitively manipulate the location and the amplitude of the maximum of the wall shear stress in order to overcome the yield stress of the mucus in specific regions of the lung.

Our results and analyzes can be obtained only because our model is a "minimal model". Minimal is meant here in the sense that the model includes only the core physical properties of the lung's biomechanics: the tree structure, the deformation of the airways, the air fluid dynamics with inertia and an external tissue pressure. The first benefit in using a minimal model is to catch and to interpret more easily the dynamics of the system. The effects of each physical process can be isolated and confronted to the others. The second benefit in using a minimal model is that the model is numerically tractable. This allows the exploration of many scenarios as its computation time is cheap. Tractable models such as ours are crucial for analyzing future models with higher complexity.

Future evolutions of our model will include a more realistic tree structure that accounts for the asymmetry of the bifurcations [START_REF] Merryn | Modeling Water Vapor and Heat Transfer in the Normal and the Intubated Airways[END_REF]. Moreover, we plan to mimic both the entering (inspiration) and the exiting (expiration) airflows. We will also evaluate the role of the turbulence in the airways with large Reynolds numbers. Finally, in the lung, the tissue pressure and the amount of airflow are related. In order to link these two physical quantities, we plan to couple an idealised 0D model of the lung's tissue mechanics to the model presented in this paper. This will allow to mimic lung's configurations that reflect better the physiology of the lung, most particularly during chest physiotherapy.

Conclusion

This work aims to improve our understanding of how the wall shear stresses are distributed during expiration in an airway tree with deformable airways. The wall shear stress has a bell-shaped distribution along the generations of the tree. This shape results from the interactions between the geometry of the tree, the compliance of the airways and the air fluid dynamics. The proximal part of the bell and the location of the maximum are dependent on the presence of inertia. Our analyses show that the properties of the distribution of the wall shear stress, and more particularly the location of the maximum, are affected by the air flow rate and the tissue pressure. In the limit of our model, these results could allow to do physics-based analyses of several commonly used chest physiotherapy technics. Moreover, our model suggests that the control of these two quantities during chest physiotherapy manipulations might allow to focus the draining of the mucus on a specific region in the lung. Hence, our work might help to understand, or possibly improve, the current mucus draining strategies. Then regrouping all the previous transformed terms as in equations 5, we obtain:

σ " |q R 2L | (7)
The wall shear stress is directed in the same direction as the airflow and in the opposite direction of the pressure drop. The signed mean shear stress relationship is then σ " ´q R 2L .

pressure is lower and the airway is more constricted. A smaller airway radius induces an airway Reynolds number larger than 300. Consequently, we have to account for inertia in our model. The second case is the opposite and a similar reasoning leads to an airway Reynolds number for the solution of the equations that is smaller than 300. The air fluid dynamics in our model is then in Poiseuille regime.

Finally, we use a Newton algorithm to search for the solution of the equation F i pqq ´q " 0 for each generation i from the root of the tree to the leaves of the tree. Starting with the initial condition q 0,i " 0, the numerical scheme in the generation i writes: q j`1,i Ð q j,i ´Fi pq j,i q ´qj,i pB q F i q pq j,i q ´1 (16)

Here the quantity q j,i is the value computed for the pressure drop q at the j th step of the iterative algorithm.

Figure 1 :

 1 Figure 1: Schematics of the airway tree structure.The number on the airways corresponds to the generation index of that airway. The generation index corresponds to the number of airways between the airway considered and the root of the tree. From the generation i to the generation i`1, the diameter of the airways is multiplied by a ratio h i ă 1.

FRCFigure 2 :

 2 Figure2: Model of compliant airways from[START_REF] Lambert | A computational model for expiratory flow[END_REF]. The curves represent the normalised surface areas of the airways of a same generation as a function of the transmural pressure in these airways. The normalisation is made relatively to a maximal possible surface area of the airways. The transmural pressure is the difference between the inner (air) and outer pressures (surrounding tissue) of the airways. The boxed index indicates the generation of the corresponding airways.

  the healthy human pulmonary acinus. Journal of Applied Physiology, 110(2):528-537, February 2011.The viscous term µ∆v becomes: Br Under the assumption that the air velocity is invariant along the axis of the airway, the quantity Bv{Bz is zero, henceµ £ S b dS ¨∇ b v ¨ẑ " µ ij S Wall dS Bv z BrFinally, the hypotheses of axial and axisymmetric air velocities that are unchanging along Ẑ allow to write ||µBv{Br|| " |µBv z {Br|. Using our definition for the mean wall shear stress σ leads to µ

Table 1 :

 1 Radii and radii reduction ratios in the airway tree as a function of the generation index, from

	Generation	FRC Radius (mm) Ratio Radius (mm) Ratio Maximal dilation
	1	8.18	0.708	8.69	0.707
	2	5.79	0.687	6.14	0.769
	3	3.98	0.716	4.72	0.790
	4	2.85	0.740	3.73	0.802
	5	2.11	0.711	2.99	0.769
	6	1.50	0.733	2.30	0.783
	7	1.10	0.736	1.80	0.783
	8	0.810	0.756	1.41	0.801
	9	0.612	0.755	1.13	0.787
	10	0.462	0.751	0.889	0.792
	11	0.347	0.769	0.704	0.805
	12	0.267	0.798	0.567	0.834
	13	0.213	0.831	0.473	0.879
	14	0.177	0.847	0.416	0.882
	15	0.150	0.873	0.367	0.902
	16	0.131	0.878	0.331	0.894
	17	0.115		0.296

[START_REF] Lambert | A computational model for expiratory flow[END_REF]

. Two states are presented: functional residual volume (FRC) (no airflow, transmural pressure of 5 cmH 2 O) and maximal dilation (no air flow, infinite transmural pressure).

Acknowledgement

The authors would like to thank the physiotherapists Jean-Claude Jeulin, Christian Fausser and Dominique Pelca for fruitful discussions.

This work has been supported by the Centre National de la Recherche Scientifique (CNRS) and by the Agence National de la Recherche in the frame of the project VirtualChest (ANR-16-CE19-0014) and of the IDEX UCA JEDI (ANR-15-IDEX-01).

Data Availability Statement

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Appendix

A Derivation of the relationship between the wall shear stress and the pressure drop

We consider a cylindrical airway b with a radius R and a length L. We use cylindrical coordinates and the axis of the airway is aligned with the axis Ẑ. The airway is delimited by three surfaces: the airway wall S Wall where the air velocity is assumed to be zero (non slip boundary conditions), the normal vector to S Wall is denoted r; the airway inlet S In , where the air flow rate is Φ " ť S In v z dS with v z the component of the air velocity along the axis of the cylinder, the normal vector to S In is ´ẑ; the airway outlet S Out where the air flow rate equals the inlet flow rate, Φ " ť S Out v z dS, the normal vector to S Out is ẑ. The total surface of the airway is denoted S b " S In Y S Out Y S Wall . A schematics of the airway is shown in figure 4.

The fluid dynamics of the air in the airway is driven by the Navier-Stokes equations that relates the pressure p and the velocity v of the air,

The shear stress is the consequence of the friction between the air layers and between the air and the airway's wall. The friction is related to the viscosity of the air, i.e. the term µ∆v in the Navier-Stokes equations. Due to the stationary hypothesis, we have Bv{Bt " 0. We denote p In and p Out the mean pressures at the inlet and the outlet: p ˚" pπR 2 q ´1 ť S˚p dS with ˚" In or Out. We define the wall shear stress σ as the mean norm of the friction constraint on the airway wall: σ " µp2πRLq ´1 ť S Wall ||Bv{Br||dS. Then, by integrating the Navier-Stokes equations ( 5) on the volume of the airway V b , using the Stokes theorem and projecting the result on the Ẑ axis, we have

The term ρv 2 z is the local kinetic energy of the air in the direction of the axis of the airway. Due to the assumption that the velocity profile is invariant along the airway axis Ẑ, the axial kinetic energy is conserved all along the airway and the previous term is zero.

The term resulting from the ∇p in the Navier-Stokes equations can be expressed using the airway pressure drop q " p In ´pOut :

B Solving the model's equations

The section 2 describes the model used in this work. The system of equations that needs to be solved is based on one sub-system for each generation i of the airway tree that writes

The pressure drop q i in the airways in generation i depends on their radius R i which depends in turn on the pressure drop q i . The shear stress σ i can be computed only when the pressure drop q i is known. So, solving the system is equivalent to find a fixed point q i of the previous system of equations [START_REF] Mauroy | Toward the modeling of mucus draining from human lung: role of airways deformation on air-mucus interaction[END_REF]. Hence, the system can be reformulated into finding the fixed point q i of an application F i that computes a pressure drop F i pqq from a given pressure drop q in generation i.

For any generation index i, the pressure drops q k for k ă i are the fixed points of the systems associated to the generations k ă i. Hence, we can compute all the fixed points of the tree-wide system by starting with the generation 0 and progressing downward the tree, down to the final one. A Newton algorithm is ran to find the zero of the application F i pqq ´q to determine the fixed point q i . In the Supplementary Materials (C), we show that the system has only one solution, hence ensuring that the convergence of the Newton algorithm leads to the correct pressure drops distribution. In the Supplementary Materials (D), details about the numerical method are given.

C Existence and uniqueness of the solution of the model equations

In order to be sure that our algorithm converges towards the correct physical solution, we need to highlight the existence and unicity of the solution.

The first step is to study the behavior of the quantity F i pqq ´q. The derivative of F i pqq relatively to q writes pB q F i q pqq " pB R i F i ¨Bα i R i ¨BP i α i ¨Bq P i q pqq (9)

F i pqq is a function of R i that can be expressed as F i pqq " aR ´4 i `bR ´5 i , with a and b data that depends on the flow rate Φ i and on the hydrodynamic resistance. Since, we know that a and b are always positive, B R i F i pqq is always negative:

Now, for B α i R i , we have:

Next, for B P i α i , we have:

"

Finally, for B q i P i , we have:

This leads to @ q P R pB q F i q pqq ´1 ă ´1

So, F i pqq ´q is strictly decreasing and lim qÑ`8 F i pqq ´q " ´8.

The next step is to show that the application F i pqq ´q is bijective from R onto R, i.e. that it can reach any real value with a unique value of q. The model for airway compliance makes the values of the airway radius R i range unequivocally from 0, at a limit transmural pressure P i " ´8, to R max,i , at a limit transmural pressure P i " `8. The dependance of R i on q is identical, since the transmural pressure P i is linearly and positively dependant on the pressure drop q. Hence, the application F i pqq that depends on 1{R 4 i ranges from `8 when q goes to ´8 to a limit value q inf,i when q goes to `8. We can conclude that the application F i pqq ´q ranges from `8 when q goes to ´8 to ´8 when q goes to `8. Finally, adding the fact that the application F i pqq ´q is strictly decreasing, we can conclude that F i pqq ´q is bijective from R to R Hence, there is one, and only one, solution F i pqq ´q " 0 for the equations (8).

D Numerical method for solving the system

For each generation i, we start by testing if the inertia is affecting the flow profile. To get this information, we use the fact that F i pqq is decreasing. A simple graphical analysis shows that, for any pressure drop q a , the solution of F i pqq " q lies between the values of q a and F i pq a q. We define the pressure drop q| Re"300 that corresponds to a Reynolds number in the airway equal to 300. The pressure drop q| Re"300 at Re " 300 is computed with the Reynolds number formula to compute the corresponding airway radius R i | Re"300 . The radius R i | Re"300 is then injected into the formula that links the pressure drop to the air flow, see equation [START_REF] Mauroy | Toward the modeling of mucus draining from human lung: role of airways deformation on air-mucus interaction[END_REF]. As the solution q ‹ of F i pqq ´q " 0 stands between q| Re"300 and F i pq| Re"300 q, two cases are possible. The first case corresponds to q ‹ ď q| Re"300 .

As the solution of the equation has a pressure drop which is smaller than q| Re"300 , the transmural


Jabrane Karkouri

Hélène Ratiney

Magalie Viallon

Rémy Prost

Life Member IEEE Fabien Millioz

TIME UNDERSAMPLED ACQUISITION FOR MULTIDIMENSIONAL SPARSE SIGNALS WITH APPLICATION TO MAGNETIC RESONANCE SPECTROSCOPIC IMAGING

Keywords: Irregular sampling, undersampling, sample selection, sensor placement, least-square, magnetic resonance spectroscopic imaging, spiral spectroscopic imaging

This paper deals with acceleration of multidimensional signal acquisition. The signal is assumed to have multiple discrete spatial dimensions where each point is time varying. The one-dimensional Fourier transform of the time evolution of each point is assumed to have an a priori known bounded support. The Fourier transform of the spatial domain is divided into several partitions, each of which can be sequentially acquired over time. We propose a method for undersampling the time dimension that enables interleaving of the samples hence speeding up signal acquisition. This method, applied to realistic simulated magnetic resonance spectroscopic imaging data (MRSI), leads to a reduction in the acquisition time by a factor of three.

INTRODUCTION

Compressed Sensing (CS) [START_REF] Candès | Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[END_REF], [START_REF] Donoho | Compressed sensing[END_REF] introduced a new way to sample, with a lowered sampling rate compared to the Nyquist-Shannon criterion. The signal to be sampled has to be sparse in a given space so that CS can be applied. CS has already been used for several applications in medical imaging [START_REF] Lustig | Compressed sensing mri[END_REF], digital imaging [START_REF] Duarte | Single-pixel imaging via compressive sampling[END_REF] and radar imaging [START_REF] Herman | High-resolution radar via compressed sensing[END_REF]. It is based on three main steps: sparse representation, measurement encoding, and sparse signal reconstruction [START_REF] Candès | An introduction to compressive sampling[END_REF].

For a known sparse spectra support, it is possible to accelerate and simplify the reconstruction process. In [START_REF] Gao | Optimal k-space sampling in MRSI for images with a limited region of support[END_REF], [START_REF] Merhej | Fast multidimensional NMR spectroscopy for sparse spectra[END_REF] the l 2 norm was used to reconstruct the signal with a known sparse spectra support using the least square (LS) overdetermined solution. In [START_REF] Gao | Optimal k-space sampling in MRSI for images with a limited region of support[END_REF] optimal k-space sampling in MRSI for images with a limited region of support is investigated while in [START_REF] Merhej | Fast multidimensional NMR spectroscopy for sparse spectra[END_REF] a fast multidimensional NMR spectroscopy for sparse spectra is proposed. Compared to the orthogonal matching pursuit algorithm (OMP) [START_REF] Wang | Generalized orthogonal matching pursuit[END_REF], the LS solution corresponds to only the last step of this algorithm when all of the samples of the support are found [START_REF] Tropp | Signal recovery from random measurements via orthogonal matching pursuit[END_REF]. In a noisy scenario, the reconstruction error depends on the choice of the acquired time samples, but selecting the optimal choice is combinatorial. Computationally efficient methods leading to good suboptimal solutions are required. The first proposed algorithm was the Sequential Backward Selection (SBS) algorithm [START_REF] Reeves | Selection of observations in signal reconstruction[END_REF], iteratively removing samples up to the desired number, minimizing the Mean Squared Error (MSE) at each iteration. The same authors proposed an equivalent constructive algorithm called the Sequential Forward Algorithm (SFS), [START_REF] Reeves | Sequential algorithms for observation selection[END_REF], but it requires a full rank matrix to start with, and thus some already selected samples. This sample selection problem is equivalent to a sensor placement problem, and thus has several solving algorithms and methods, such as approximating the problem using a convex relaxation [START_REF] Joshi | Sensor selection via convex optimization[END_REF], [START_REF] Jamali-Rad | Distributed sparsity-aware sensor selection[END_REF], or greedy methods. The latter consists in selecting samples one by one by optimizing a proxy of ideal acquisition, such as the MSE [START_REF] Jiang | Sensor placement by maximal projection on minimum eigenspace for linear inverse problems[END_REF], the worst case error variance [START_REF] Jiang | Sensor and cfd data fusion for airflow field estimation[END_REF], the volume of the confidence ellipsoid [START_REF] Shamaiah | Greedy sensor selection: Leveraging submodularity[END_REF], or the frame potential [START_REF] Ranieri | Near-optimal sensor placement for linear inverse problems[END_REF]. An exhaustive bibliography is given in [START_REF] Jiang | Group greedy method for sensor placement[END_REF], with an extension to these greedy methods, by searching samples among a given number of groups. One can note that SBS is a backward version of the greedy method minimizing the MSE, while SFS is identical as soon as the number of selected samples leads to a full rank matrix.

In the present study, we tackle the problem of the design of a fast sampling strategy for multidimensional signal acquisition, with its application to simulated magnetic resonance spectroscopic imaging (MRSI) data [START_REF] Gao | Optimal k-space sampling in MRSI for images with a limited region of support[END_REF]. Here we use the sample selection algorithm in a different way, and we push further the method to design a new sampling strategy taking advantage of the multidimensional characteristic. The first results from this study were presented in [START_REF] Karkouri | Time samples selection in spiral acquisition for sparse magnetic resonance spectroscopic imaging[END_REF].

Magnetic resonance spectroscopic imaging (MRSI) signals consist of several spatially distributed one-dimensional spectra. As for magnetic resonance imaging, the acquisition is performed in a spatial two-dimensional spectrum domain called the k-space, for which (k x ,k y ) represents the wave numbers in each dimension. In MRSI, this k-space evolves in time according to an oscillatory function with an exponential decay. The oscillation frequencies are characteristic of the molecular content under analysis. A two-dimensional Fourier transform of the k-space is then necessary to spatially localize the biochemical information, and a one-dimensional Fourier transform along the time direction is needed to obtain the magnetic resonance spectra. The time domain and the k-space domain are separable. At each (k x ,k y ) point, the associated spectrum is bandlimited. In in vivo magnetic resonance spectroscopy, it is customary to consider that molecular resonance frequencies in living tissues [START_REF] Govindaraju | Proton nmr chemical shifts and coupling constants for brain metabolites[END_REF], [START_REF] De Graaf | vivo NMR spectroscopy: principles and techniques[END_REF] and the expected 'full width at half maximum' (FWHM) of the resonating peaks are known a priori. Moreover, frequencies and FWHM are used as starting values and prior knowledge in magnetic resonance signal quantification [START_REF] Ratiney | Time-domain semiparametric estimation based on a metabolite basis set[END_REF]; here, these are used to define this limited band a priori, i.e, the spectral support.

For the purpose of generality of this work, we deal with a spatio-temporal signal model y(U, n), where U is a Kmultidimensional space and n is the discrete time dimension. It is assumed that y(U, n) can be fully acquired by sampling it sequentially over different partitions of U .

For this particular MRSI application, U is the k-space. Due to physical constraints, it is not possible to scan the whole joint spatial and temporal dimensions of the signal all at once. Indeed, at each scan, data points are sequentially acquired over time and over limited regions of the k-space. Given that the signal is transient, several scans are required to obtain the full U space over time. As a result, the full and dense sampling of this signal is time consuming. This approach is similar to multiple measurement vector methods [START_REF] Cotter | Sparse solutions to linear inverse problems with multiple measurement vectors[END_REF], [START_REF] Chen | Theoretical results on sparse representations of multiple-measurement vectors[END_REF], which represent an extension of the single measurement vector methods, with the difference being that the acquisitions are sequential and the spectral support is known.

Several fast k-space scanning or CS methods have been proposed for MRSI [START_REF] Posse | MR spectroscopic imaging: principles and recent advances[END_REF], [START_REF] De Graaf | vivo NMR spectroscopy: principles and techniques[END_REF]. The undersampling of the temporal dimension has never been proposed except for a CS scenario [START_REF] Hu | Compressed sensing for resolution enhancement of hyperpolarized 13C flyback 3D-MRSI[END_REF]. In the special case of sparse spectra, we show that subsampling the temporal dimension enables interleaving of the acquisitions of the partitions of U of the signal y(U, n). Therefore, it leads to a decrease in the number of scans, and a reduction of the acquisition time of the signal. We propose an algorithm that specifies the signal acquisition scheme.

In this work, our main contribution is:

• a method to sample a multidimensional signal in a multidimensional framework in order to accelerate the acquisition while minimizing the MSE; Other contributions include:

• a group extension of SBS, and its comparison to the state-of-the-art of one-dimensonnal sensor placement; • theorems required for our proposal;

• an application of the proposed method to a realistic MRSI simulation with performances evaluation. This paper is organized as follows. In section 2, we first extend SBS to a group method and compare it with a generic state-of-the-art method, and then apply it to a one dimensional sparse-spectrum signal to define an acquisition strategy. Finally two theorems are derived from this special case. Section 3 uses these theorems to propose a multidimensional signal acquisition algorithm. The MRSI application is presented in section 4, followed by our conclusion and some perspectives.

SAMPLE SELECTION FOR MSE MINIMIZATION 2.1. Problem statement and state-of-the-art

We want to solve a linear inverse problem, expressed as

y = Ψx + ν, (1)
where y is a N -samples measurement vector, x a M -samples unknown vector, Ψ an N × M observation matrix, and ν a white Gaussian zero mean noise which is identically distributed, with variance σ 2 . We use N measurements to estimate the M -dimensional x vector, assuming N > M . In this case, the matrix system is overdetermined, and the minimum variance unbiased estimate of x is the Least Square (LS) solution

x = (Ψ * Ψ) -1 Ψ * y, (2)
where * denotes the conjugate transpose. Given that (1) is overdetermined, it is possible to subsample y, by selecting only P rows from Ψ, M ≤ P < N .

This sample selection problem is strictly equivalent to a sensor placement problem: from the N possible samples, which ones are the P best leading to a minimal error?

Here we arbitrarily fix P for targeted subsampling factor of N/P . The lower the subsampling, the higher the error.

The choice of P is a tradeoff between subsampling and increased MSE. The optimal solution can be found by an exhaustive search, but this is computationally unaffordable. Algorithms leading to good solutions in a reasonable time, even if not optimal, are thus required.

Among these algorithms, greedy methods iteratively select rows of Ψ one by one by optimizing a given criterion, arbitrarily chosen as the Mean Squared Error (MSE).

In this work, we use the MSE criterion as it is commonly used in MRSI.

The MSE of the estimate x of (2) is

M SE(x) = E||x -x|| 2 = σ 2 tr[(Ψ * Ψ) -1 ]. (3)
The idea of the greedy method is to iteratively select the row of Ψ that minimizes M SE(x) at each iteration. Note that when the number of selected rows is smaller than the number of columns M , Ψ * Ψ is singular. To avoid this problem, Jiang et al. [START_REF] Jiang | Group greedy method for sensor placement[END_REF] replace this singular matrix by Ψ * Ψ + I, where is a small number and I the identity matrix. In addition, Jiang et al. extended the greedy method by searching for L parallel groups at each iteration, with the intuition that it expands the solution search, and should be closer to optimal. They also compared their group greedy method [START_REF] Jiang | Group greedy method for sensor placement[END_REF] with five methods: FrameSense [START_REF] Ranieri | Near-optimal sensor placement for linear inverse problems[END_REF], convex relaxation, [START_REF] Joshi | Sensor selection via convex optimization[END_REF], SparSense [START_REF] Jamali-Rad | Distributed sparsity-aware sensor selection[END_REF], Greedy [START_REF] Jiang | Sensor and cfd data fusion for airflow field estimation[END_REF] [17] and even random selection, demonstrating that the group greedy method leads to the best results.

We propose here to expand the SBS method by using the group idea of [START_REF] Jiang | Group greedy method for sensor placement[END_REF], and compare it to the group greedy method using the same experiment.

The observation matrix Ψ is a 100 × 20 uniform random matrix with independent components following a uniform distribution in [0; 1]. The greedy method and SBS are applied to this matrix, for groups sizes of L ∈ [START_REF] Candès | Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[END_REF][START_REF] Herman | High-resolution radar via compressed sensing[END_REF][START_REF] Bogner | Accelerated mr spectroscopic imaging-a review of current and emerging techniques[END_REF], for a number of selected rows ranging from 20 to 30. The MSE is computed for each case and averaged over 100 realizations of Ψ.

Fig. 1 shows the results. Overall, SBS gives better results than the greedy method. Moreover, as expected, the bigger the group, the lower the MSE, leading to a better sample selection.

Application to sparse spectrum signals

In this section, we consider a sparse spectrum signal

y = F -1 x, ( 4
)
where y is an N -samples time-domain signal, with a sampling period T s , x its N -samples associated spectrum, and F the unitary discrete Fourier matrix of size N × N , F n,k = exp(-i2πnk/N )/ √ N . Let x m be the restriction of x to its M non-zero values. Equation (4) becomes

y = Ax m , ( 5
)
where A is the N × M matrix that is obtained by selecting from F -1 the columns that correspond to the M elements of x in its sparse support S. A is the special case of Ψ in [START_REF] Candès | Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[END_REF]. We need to determine the best P samples of y to be acquired using a MSE-minimization algorithm, leading to a final A matrix of size P × M .

The SBS algorithm and the greedy method have been compared for group sizes of L ∈ [START_REF] Candès | Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[END_REF][START_REF] Herman | High-resolution radar via compressed sensing[END_REF][START_REF] Bogner | Accelerated mr spectroscopic imaging-a review of current and emerging techniques[END_REF] with three types of support, which we define as compact, block, and random. Each of these support types has M = 128 nonzero values and the total sample number N = 1024. A compact support is defined as M consecutive nonzero frequency locations. A block support is defined as eight different compact supports of size M/8 that are randomly distributed. A random support is defined as M randomly selected nonzero frequency locations.

Figure 2 shows the MSE evolution that is averaged on 50 random realizations of each support, versus the number P of selected rows with N = 1024. Several conclusions arise.

i) The more "compact" the support, the lower the MSE.

ii) Except for random support, using groups is harmful. When looking at the evolution of the MSE along the iterations, using groups shows better MSE at the beginning. However, selecting these better groups finally leads to higher MSE at the end. This shows the non-linearity of the iterative search. However, when looking at individual realizations, sometimes using large groups improves the MSE, while on average it does not.

iii) The greedy algorithm is on average better than the SBS one. Note that for certain realizations, SBS is better than the greedy method.

Based on these observations and given the particularities of the Fourier matrix, in the following, we choose to use as an MSE-minimization algorithm the greedy method, with a group size of 1.

Properties of the subsampled sparse Fourier matrix

The following theorems give some properties of the matrix A in (5) with only P selected rows, which will be useful for our acquisition method described in section 3.

Theorem 1.

Let A be the P × M , M ≤ P ≤ N , matrix constructed with the discrete Fourier transform matrix F -1 of size N × N where N -M columns and N -P rows have been deleted. The entries of A are (A) ns,ks = e 2πinsks/N / √ N , where k s are the columns indexes corresponding to the support of x. The indexes n s are in the list p of the P selected rows in F -1 . Consider a circular shift d, with period N , of the indexes in p and select the corresponding rows in F -1 : if

A * A is nonsingular, then tr[(A * A) -1 ] is invariant
= e 2πi(ns+d)ks/N / √ N . with (n s + d) (mod N ), Then (A d
tr[(A d * A d ) -1 ] = tr[(D d * A * AD d ) -1 ] = tr[(D d ) -1 (A * A) -1 (D d ) -1 * ] = tr[(D d * ) -1 (D d ) -1 (A * A) -1 ] = tr[(A * A) -1 ].

Involvement of Theorem 1

For noisy data with the same hypothesis as in section 2.2, the selected time samples to be acquired can be circularly shifted, without increasing the resulting trace of the square matrix (A * A) -1 , and therefore without increasing the noise amplification. This enables more flexibility in the choice of time samples used to solve (5) during the signal acquisition and this is a key point for our method based on the interleaved acquisition described in section 3.

Theorem 2.

Given the matrix A defined in Theorem 1, consider a circular shift r of the support of x with period N , the columns of the matrix A. If A * A is nonsingular, then A * A is invariant from the circular shift, and so is tr[(A * A) -1 ].

Proof.

Consider a circular shift r of the M selected columns.

Write (A r ) ns,ks = (A) ns,ks+r = e 2πins(ks+r)/N / √ N .

Then (A r ) ns,ks+r = e 2πirns/N (A) ns,ks

It follows that A r = D r A
where D r = diag(e 2πirns/N ), n = 0, 1, . . . , P -1.

Note that D r * D r = I P , where I P is the P × P identity matrix.

Then A r * A r = A * D r * D r A = A * A Involvement of Theorem 2 tr[(A * A) -1
] is invariant from any circular shift of the spectrum support (mod N ), which does not require a new computation of the MSE minimization algorithm.

MULTI-SETS ACQUISITION

This section presents the signal model used in this study, and an algorithm whose goal is to minimize the acquisition time required to sample this signal.

Model of the temporal signal y(U, n)

Let y(U, n) be a multidimensional signal, where n is the discrete time dimension, and U a K-dimensional space. U is assumed to be divisable into N u partitions U u , such that

U = ∪ Nu-1 u=0 U u . (6)
Consequently, the signal y(U, n) over all partitions is the union of each signal acquired over partitions

y(U, n) = ∪ Nu-1 u=0 y u (n) ∀n ∈ [0, N -1], (7)
where y u (n) is the signal acquired only on partition U u .

At each n, it is assumed that only a single partition can be sensed, leading to a single y u (n). The acquisition of a single y u (n) over time is called time-scan in the following. N u time-scans of different signals y u (n) are required to collect the full signal y(U, n) over of the N u partitions. For a uniformly sampled temporal dimension, the acquisition time is

T f ull acq = T scan N f ull acq , (8)
where T scan = N T s is the duration of one time-scan, T s is the sampling time period of the signal y u (n) and N f ull acq = N u the total number of time-scans. To reduce the acquisition time, we subsample y u (n) to lower the number of data points along the temporal dimension. Each signal y u (n) is assumed to have an identical known sparse spectral support S. Thus, the signals y u (n) can be sub-sampled as described in section 2.2. Let {p} be the set of all the selected time samples by the MSE-minimization algorithm with Card({p}) = P < N , P ≥ M . As each y u (n) has the same spectral support S, it follows that the selected samples {p} are also identical for each y u (n). At each time sample p, N u partitions have to be acquired simultaneously and thus N u time-scans are still required. To reduce the acquisition time, the idea is to take advantage of the unused time zones. The time indices n / ∈ {p} will be used to acquire other partitions of U . From Theorem 1, these unused time zones can be filled by a circular shift of {p}. The ideal scenario would be to fill the N u P samples to be acquired without any empty time zones during the N -points time-scan. This results in N u P/N time-scans compared to N u time-scans in the full sampled acquisition. The best acceleration factor is thus N/P . In practice, this is very unlikely to happen.

Interleaving the acquisition in the time domain

The objective is to minimize the total number of time-scans N acq to acquire y(U, n). This objective requires that the maximum number of partitions to be acquired for any n should be minimal.

Let {p} u be the set of selected time samples for the signal of the partition u, y u (n). Let I u (n) be an index vector of length N that indicates the P -samples to be acquired in partition u, I u (n) = 1 if n ∈ {p} u , 0 otherwise. The total number of time-scans is then

N acq = max n ( Nu u =1 I u (n)).
The proposed algorithm minimizes the maximum number of partitions to be acquired for any n. It searches at each iteration j for a set {p} j that minimizes the number of acquisitions considering only the first j partitions,

N j acq = max n j u =1 I u (n) , (9)
where j is from 1 to N u . The algorithm selects the best circular shift of {p} that minimizes N j acq , to determine {p} j . In the case of multiple solutions, the smaller shift is selected. The first iteration sets {p} 1 = {p}, and the final iteration acquire for each partition. Let us define I 1 (n -d mod N ) as the d circular shift of I 1 (n).

Algorithm 1 Proposed algorithm -Run the MSE-minimization algorithm to get {p} -Set {p} 1 = {p} and derive

I 1 (n) for j = 2 to N u do -Search d j such that min dj (max n j u =1 I u (n) ), with I j (n) = I 1 (n -d j mod N ) -Set {p} j = {p, circularly shifted by d j }, derive I j (n) end for -N acq = N Nu acq = max n Nu u =1 I u (n)
When the algorithm finishes, the N u set of points {p} u needed to acquire each partition are determined, as well as the number of acquisition required to do so.

The new acquisition time of y(U, n) is then

T acq = T scan N acq . (10)
Compared to Equation ( 8), the time acceleration factor is N f ull acq /N acq = N u /N acq ≥ 1.

Illustration of the proposed algorithm

An illustrative example of the proposed algorithm is used here, as a signal of size N = 16 with a spectrum block support size M = 7. The number of partitions is set to N u = 3. The MSE-minimization algorithm is run with P = 8 time samples. Figure 3 shows the iterative steps of the algorithm and the sets {p} u and I u (n) generated at each step. In this example, the final number of acquisitions N acq is 2, to give a time acceleration factor of 3/2. For each iteration, the new time samples are represented in red in Figure 3. In the first step, {p} 1 = {p} and I 1 (n) = 1 if n ∈ {p} 1 , 0 otherwise. In the second step, there is no circular shift of {p} 1 , such that the maximum of I 1 (n) + I 2 (n) is 1. In the third step, the algorithm finds a solution such that the maximum of I 1 (n)+ I 2 (n) + I 3 (n) leads to two time-scans.

Figure 4 illustrates how the partitions are acquired at each time-scan. During the first time-scan, it is always possible to entirely acquire the signal over the first partition (in blue), given that by definition a time-scan is the duration required to acquire the signal. Given than the signal over the first partition is subsampled, some time bins are unused. The algorithm shifts the selected samples to fill these holes for the other partitions. Consequently, 6 time samples related to the second partition and 2 time samples related to the third one are sampled during the first time scan. As a consequence, the remaining points of the last two partitions can be sampled in

= NuM/N 1 1 1 1 1 1 1 2 1 1 2 4
Table 1. Given N u , M and different support types, the highest number of acquisitions N f ull acq (= N u ), the lowest N f ull acq M/N (rounded up) and the computed N acq by the proposed algorithm are given, averaged over 100 realisations of the randomly chosen support. The greedy algorithm was used with P = M and N = 1024. a second time scan, avoiding the third time scan that would have been required in the case of an acquisition without interleaving.

Study of the time acceleration factor

To evaluate the time acceleration factor reached by the algorithm proposed in the previous section, the algorithm is run to test different spectral supports.

The signal size is set to N = 1024; the spectrum has M nonzero values, with M ∈ {64, 128, 256}; the support type of the spectrum is compact, block or random (See section 2.2); the number of partitions is N u ∈ {2, 4, 8, 16}.

The MSE-minimization algorithm is run with P = M time points, which is the largest feasible subsampling with a given sparse spectrum.

For each set of parameters, the algorithm is run on 100 realizations of the randomly chosen support. For each run, the following are computed: the highest number of acquisitions required N f ull acq = N u ; the ideal lowest number of acquisitions N f ull acq M/N (rounded up); and the number N acq of acquisitions given by the proposed algorithm. The results are averaged over the 100 realizations, as in Table 1.

We can see that the acquisition speed-up, which is defined as N f ull acq /N acq , increases with the number of partitions and decreases with the number of samples in the spectrum support. As expected, the sparser the spectrum, i.e., the lower the M , the lower the number of acquisitions. In comparison, the number of partitions N u has less influence. The support type has strong influence: the more "compact" it is, the better acquisition gain it gives. Interestingly, for a sparse spectrum of M contiguous frequency samples, the algorithm always gives the lowest possible number of acquisitions.

Performance and complexity of the algorithm

The algorithm proposed in section 3.2 searches sequentially in N u × N circular shifts for the combination of shifts that leads to the lowest number of time-scans. The total number of possible combinations is N N u . As seen in the previous section, this algorithm actually gives the lowest possible number of time-scans in the case of a compact support.

For the two other cases, an extensive search in the N N u possibilities is run over 100 iterations of both the block and random supports, with a signal size N = 128, 3. The first partition is in blue with a square symbol, the second is in red with a circle symbol, the third is in green with a cross symbol. Two time-scans are enough to sample all of the time samples of all of the partitions.

N u = 4 partitions and spectra of M nonzero values, M ∈ {8, 16, 32, 64}. When the spectrum is very sparse, M = 8, many combinations lead to the best possible N acq , and the proposed algorithm always finds this. On the other hand, when half the spectrum points are used, M = 64, the extensive search gets a better result than our algorithm in each run for a random support, and in 44 runs for a block support. For M = 16 and M = 32, the extensive search gets a better N acq in 3% of the runs for a block support, and in 27.5% of the runs for a random support.

The proposed algorithm first runs the MSE-minimization algorithm, then performs circular shifts as detailed in section 3.2. The greedy MSE-minimization algorithm has a complexity of O(N P M 2 ) [START_REF] Jiang | Group greedy method for sensor placement[END_REF]. The proposed algorithm in section 3.2 (excluding the MSE-minimization) has a complexity of O(N 2 N u ). The LS reconstruction has a complexity of O(M 2 P ) for each partition. Note that the selection of the temporal samples to acquire (MSE-minimization) is used only once before the acquisitions.

APPLICATION

Magnetic resonance spectroscopic imaging (MRSI)

Nuclear magnetic resonance relies upon in a physical phenomenon that induces a perturbation of the macroscopic magnetization of an ensemble of nuclear spins that interact with a surrounding static magnetic field. This perturbation is caused by radiofrequency pulses, which are more commonly called excitation pulses, and lead to an electromagnetic signal response with frequencies related to the atoms considered and their chemical environment [START_REF] De Graaf | vivo NMR spectroscopy: principles and techniques[END_REF]. This signal is called the free induction decay and it is encoded in the k-space with the help of magnetic field gradients. It is complex, of the form y(k x , k y , t), depending on the time t and the spatial two-dimensional excited environment. In magnetic resonance imaging, the signal y(k x , k y , t) is acquired in the Fourier domain, which is called the k-space excitation, at a given t. To obtain the corresponding anatomic image, an inverse Fourier transform is performed. In MRSI, we are also interested in the temporal evolution of the signal, which is why it is acquired for t = nT s , n ∈ [0, (N -1)]. For each (k x , k y ), the signal is modeled by

y kx,ky (n) = j a j (k x , k y )exp(i2πf j nT s -nT s /T * 2 ), ( 11
)
where the time constant T * 2 is the transverse relaxation that accounts for the spin-spin relaxation and the field inhomogeneity. Due to the term exp(-nT s /T * 2 ), its spectrum is not strictly bandlimited because each spectral line becomes Lorentzian shaped; however, it can be approximately considered as sparse [START_REF] Merhej | Fast multidimensional NMR spectroscopy for sparse spectra[END_REF]. The coefficients a j (k x , k y ) are complex and are the Fourier transform of the metabolic image characterized by the frequency f j .

Each radiofrequency pulse excitation results in the same free induction decay signal. Excitations can be repeated as many times as necessary to get the joint spatial spectral information. At the end of the scan, a delay can be added to ensure that the macroscopic magnetization returns to its original state. The time between two excitations is called the repetition time, T R ≥ N T s .

In routine clinical work, the sequences that are commonly used acquire the k-t space, (k x , k y )-point by (k x , k y )-point, with the phase encoding method [START_REF] Posse | MR spectroscopic imaging: principles and recent advances[END_REF]. A single point of the k-space and the whole temporal dimension is acquired at each excitation. This results indeed in a long acquisition time because, in that case, N u , the number of partitions of the kspace, will be very large (see equation (8)). For example, for a 32 × 32 matrix, the number of (k x , k y ) points is N u = 32 × 32 = 1024. To speed up the acquisition, it is possible to group the k-space points into partitions, that have to be acquired after each excitation, which leads to the model of [START_REF] Candès | An introduction to compressive sampling[END_REF]. However, due to technical constraints, the acquisition of one partition is not instantaneous. The larger the partition is, the longer the acquisition is. There is an issue if the partition acquisition time is longer than the sampling time period T s , because this means that the next time sample is to be acquired while the previous one is still in the process of being acquired. Let T part be the acquisition time of each partition. The acquisition of each partition lasts for N p = Tpart Ts time sampling periods and forbids the acquisition of the next N p -1 time sample. N p -1 additional excitations are required to acquire the full signal. This is called temporal interleaving. Moreover, as the repetition time T R is longer than N T s , the total time acquisition time can be rewritten as

T f ull acq = T R N u N p . (12)
Acquiring the whole signal requires N u × N p excitations.

) Fig. 7. Partitions of the k-space with N u = 4 spirals Choosing small partitions, which are faster to sample, decreases N p but increases the number of partitions N u and vice versa.

Time interleaved application to MRSI acquisition

The additional constraint of the noninstantaneous acquisition of a partition requires an update of the proposed algorithm in section 3.2. In the previously proposed algorithm, the assumption of instantaneous acquisition of a partition is equivalent to N p = 1. I u (n) is the index vector that indicates whether the time sample n of the partition u is to be acquired, consequently forbidding any other acquisition at time index n. With the constraint of the duration of one partition, a time sample at position n forbids the N p -1 following sample to be acquired. As an example for N p = 2, if there are two consecutive points to be acquired, the acquisition of the first one blocks out the acquisition of the second. This leads to at least two excitations to obtain the two points. By setting

I + u (n) = N p-1 d=0 I u (n-d)
, the number of excitations to acquire the partition u is max(I + u (n)), as illustrated in Figure 5. It follows then that the algorithm in section 3.2 requires the replacement of I u (n) by I + u (n). To promote the scan of the beginning of the free induction decay where the signal magnitude is maximal, the {p} u are circularly left shifted, if necessary, using Theorem 1.

The steps to fully reconstruct the spectroscopic image y(x, y, f ) are the following:

• determine the time samples {p} u for each partition, through the algorithm in section 3.2 with I + u (n) instead of I u (n);

• acquire each partition y u (n) for n ∈ {p} u of the kspace; • apply LS reconstruction on each partition to obtain y u (f ) (Equation ( 3)); • merge all partitions to get the whole k -f space: y(k x , k y , f ) = ∪ Nu-1 u=0 y u (f ); • apply a two-dimensional nonuniform discrete Fourier transform [START_REF] Fessler | Image reconstruction toolbox[END_REF], [START_REF] Fessler | Nonuniform fast fourier transforms using min-max interpolation[END_REF] to reconstruct the spectroscopic image y(x, y, f ). Note that reproduction of the experiments of section 3.4 with the additional parameter N p shows that N p has similar influence to N u on the speed-up gain.

Simulation and Results

To compare the fully sampled method with the proposed interleaved time undersampling method, we use a numerical phantom of 31 P spiral MRSI data, mimicking acquisition of the thigh muscle, with three compartments as defined in Figure 6-a and in which different metabolite ratios are used. The acquisition is simulated with a field of view of 25 cm, and an image matrix size of 32 × 32. The k-space is partitioned in both methods with four spirals, as illustrated in Figure 6. The simulation parameters are N u = 4 partitions, time sampling period T s = 0.25 ms, N = 1024 time samples, repetition time T R = 2s, and spiral of N p = 10 discrete time samples which corresponds to a duration of T part = 2.5 ms. Each partition is a discrete spiral of 2048 points, which leads to 8192 points to sample the whole kspace. The support of the phosphorus spectrum contains two blocks, of total size M = 217 samples. Here, the support is chosen such that the truncation of the spectral peaks is less than 1% of the total peak energy. In practice, the resonance frequency of the molecule is well known and the full width at half maximum depends on the static magnetic field homogeneity which is rather reproducible, so a support can always be determined a priori. For the proposed interleaved method, an objective of a best acceleration factor of 4 is used by setting P = 256. As the spiral is not instantaneous in practice, every sample of the spiral is phase-corrected according to the Fourier shift theorem.

For the fully sampled spiral method, the acquisition time is T f ull acq = 80s. When applying the algorithm of section 4.2, this results in an acquisition time T acq = 26s, corresponding to a reduction by a factor of 3.08.

The two methods are also compared in terms of the signalto-reconstruction error ratio (SRER), defined as

SRER = 10 × log E||x|| 2 E||x -x m || 2 , (13)
with x the original spectrum and x m the reconstructed one, both restricted to the spectral support. Note that x, the original spectrum, is unknown. Incorporating (3) into [START_REF] Joshi | Sensor selection via convex optimization[END_REF], it follows

SRER = SNR -10 × log(tr[(A * A) -1 ]), (14)
where SNR is the unknown signal-to-noise ratio

SNR = 10 × log E||x|| 2 σ 2 . ( 15
)
Only the full sampled spectrum x f can be known. Its LS reconstruction is given by all the N rows of the matrix A.

The reconstruction error follows

E||x -x f || 2 = σ 2 M ( 16
)
leading to the SRER for full sampling results

SRER f = SNR -10 × log(σ 2 M ). (17)
The difference between the SRER f and the SRER cancels the unknown SNR

SRER f -SRER = 10 × log( tr[(A * A) -1 ] M ). (18)
In our experiments, the obtained undersampling scheme leads to SRER f -SRER = 7.5 dB.

SRER performance for a realistic simulation

Noisy spiral spectroscopic imaging data sets were generated, for both full sampling and subsampling with the proposed time interleaved acquisition, by adding complex white Gaussian noise samples to y(U, n) (in the k-space, for all the partitions). The noise standard deviations were chosen to get two SN Rs: high SNR (≈ 27dB on average over all the samples containing signal energy), and low SNR (≈18dB, on average). 50 noisy data sets per SNR were drawn and SRER as well as SRER f were computed for each set, and SRER performance were analyzed on average over all 50 sets. The T * 2 used for the exponential decay of the free induction decay signal was 100 ms.

The results for the generated phosphorus numerical phantom of the thigh are shown in figures 6 (b,c,d) and 8 and table 2. Figure 6 shows the reconstructed phantom with the use of spiral spatial sampling, full time sampling, SN R =27.38 dB for different molecular frequencies (i.e. (b) phosphocreatine (PCr) for f = 0 Hz, (c) inorganic phosphate (Pi) for f = 258 Hz and (d) γ-adenosine triphosphate (γ-ATP) for f = -129 Hz). Figure 8 shows the reconstructed spectra at the location indicated by the red cross on Fig. 6, with a noisy acquisition for the fully sampling method and for the sub-sampling proposed method, with their respective errors.

In table 2, results over all the voxels containing signals and results restricted to the left ellipse of the signal from Fig. 7a) (where all molecules are present, and thus with a higher SNR with respect to the other regions) are given. The average of the SRER f and SRER are computed over the 50 noisy sets and are given over two selected region: "Ellipse", the voxels in the left ellipse, and "All voxels", all the voxels with non-zero signal, plus or minus their associated mean standard deviation. One can see that in the noisiest scenario (SNR=18.67±2.61 dB), the theoretical SRER difference between full and undersampled case (Equation 18) is recovered. With a higher SNR, the difference is higher than expected: the choice of the spectral support, which truncates the spectra that are not perfectly sparse, leads to an additional error reconstruction compared to the fully sampled case.

CONCLUSION & PERSPECTIVES

A new acquisition method of interleaved time samples in undersampled signals is proposed here to speed up acquisition times, assuming a sparse known support.

Undersampling such a monodimensional signal can be done by minimizing the MSE of the LS reconstruction. We show that a circular shift of the selected samples does not change the reconstruction error. For multidimensional Table 2. SRER performances for noisy phantoms of phosphorus spiral spectroscopic imaging, described in Fig 6 . Two standard deviations of noise were chosen, resulting in different averaged SNR according to the spatial regions analyzed ("Ellipse", voxels containing all molecules responses, "All voxels", any voxel containing nonzero signal).

signals, it means that different equivalent sets of samples can be selected for each additional dimension. By interleaving the partition sample acquisitions in the time domain, we can then reduce the total acquisition time. The proposed algorithm sequentially searches for the circular shift of the undersampling for each partition that will minimize the acquisition time. The proposed algorithm is applied here to MRSI. Our MRSI acquisition and reconstruction solution was compared to the results of a standard fully sampled spiral acquisition method which is here taken as the state-of-the art method. Note that among the accelerated sampling schemes [START_REF] Bogner | Accelerated mr spectroscopic imaging-a review of current and emerging techniques[END_REF], the spiral ones are the most efficient trajectories, offering the best acquisition time and spectral bandwidth benefit with the smallest SNR loss, but they are limited by a susceptibility to gradient infidelities. The results demonstrated that the acquisition time is reduced by a factor of 3.08 but the signal error SRER decreases by 7.5 dB in a realistic noisy acquisition.

The proposed interleaved acquisition strategy can, in theory, be generalized to three dimensional data (3D), with 3D partitioning, and (any) 3D trajectories that would regularly return to the center of k-space but that would also require to be interleaved in the spectroscopic time dimension to reach a desired spectral bandwidth. Other reconstruction algorithms that require knowledge of the support have been proposed in the literature and are interesting in terms of reconstruction performance compared to the L 2 reconstruction when the SNR deteriorates. We mention in particular the Approximate Message Passing for multiple measurement vector method [START_REF] Ziniel | Efficient high-dimensional inference in the multiple measurement vector problem[END_REF] or the modified CS-based method for problems with partially known support proposed in [START_REF] Vaswani | Modified-cs: Modifying compressive sensing for problems with partially known support[END_REF].

A future goal will be to apply this method to exercise dynamics 31 P applications, accurate fat/water quantification, and 1 H MRSI or hyperpolarized 13 C studies. Magnetoencephalography can also be considered here, as an inverse problem used for brain imaging [START_REF] Phillips | Megbased imaging of focal neuronal current sources[END_REF], [START_REF] Gorodnitsky | Neuromagnetic source imaging with focuss: a recursive weighted minimum norm algorithm[END_REF]. Magnetoencephalography signals come from activated regions in the brain. Taking several signals over time enables the magnetic response of the brain to be followed over time. The activated regions could be known a priori based on the knowledge of brain functional structure, and our method can be applied.

ACKNOWLEDGEMENTS

This study was partly supported by the LABEX PRIMES (ANR-11-LABX-0063), program "Investissements d'Avenir" (ANR-11-IDEX-0007) and Siemens Healthineers.

Fig. 1 .

 1 Fig.1. MSE versus the sensor (or sample) number, for a 100 × 20 matrix Ψ, whose components follows a uniform distribution in [0, 1], for the group greedy method and the group SBS of size L ∈[START_REF] Candès | Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[END_REF][START_REF] Herman | High-resolution radar via compressed sensing[END_REF][START_REF] Bogner | Accelerated mr spectroscopic imaging-a review of current and emerging techniques[END_REF]. MSE is averaged over 100 realizations.

  ) ns,ks = e 2πidks/N (A) ns,ks . Then (A d ) ns,ks = e 2πidks/N (A) ns,ks . It follows that A d = AD d , where D d = diag(e 2πidsks/N ), k s = 0, 1, . . . , M -1. Note that D d * D d = I M , with I M the M ×M identity matrix. Finally,

2 IFig. 3 .

 23 Fig. 3. Illustration of the proposed algorithm for N = 16, M = 7, and P = 8. At each iteration, {p} u is a circular shift of {p} 1 (in this case, 1 shift for {p} 2 and 2 shift for {p} 3 ) which minimizes the maximum value of u 1 I u (n). The maximum corresponds to the number of time-scans needed to sample all of the partitions.

Fig. 4 .

 4 Fig.[START_REF] Duarte | Single-pixel imaging via compressive sampling[END_REF]. Each row corresponds to a time-scan of the algorithm illustrated in Figure3. The first partition is in blue with a square symbol, the second is in red with a circle symbol, the third is in green with a cross symbol. Two time-scans are enough to sample all of the time samples of all of the partitions.

Fig. 5 .

 5 Fig. 5. Illustration of the difference between I 1 (n) and I + 1 (n) = N p-1 d=0 I 1 (n -d) with N p = 4. The same parameters as in Figure 3 are used here. I 1 (n) is shown with a cross symbol on these two panels, and I + 1 (n) with a circle symbol in the right panel. I + 1 (n) shows the number of excitations required for each n.

Fig. 6 .

 6 Fig. 6. a) Referenced phantom at t=0 (based on 31 P in vivo Spectral images for (b) phosphocreatine, frequency = 0 Hz, (c) inorganic phosphate, f = 258 Hz, and (d) γadenosine triphosphate, f = -129 Hz. The red cross is used to locate the spectrum in Figure 8.

Fig. 8 .

 8 Fig.8. Noisy reconstructed modulus spectrum from a voxel in the left ellipse (red cross in Figure6-b) for the fully sampled spiral method (left) and for the new proposed approach (right). The reconstruction error (shifted by -10) is given at the bottom.

  from this MSE versus the sensor (or sample) number, for a Fourier matrix restricted to a random spectral support, for the group greedy method and the group SBS of size L ∈[START_REF] Candès | Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[END_REF][START_REF] Herman | High-resolution radar via compressed sensing[END_REF][START_REF] Bogner | Accelerated mr spectroscopic imaging-a review of current and emerging techniques[END_REF]. Left: compact support. Center: block support. Right: random support. MSE is averaged over 50 realizations of each support.

	Mean MSE	128 990 1000 1010 1020 1030 1040	130 Number of samples 132 134 Support type : compact SBS, L=1 136 SBS, L=5 SBS, L=30 Greedy, L=1 Greedy, L=5 Greedy, L=30	Mean MSE	128 1800 2000 2200 2400 2600 2800 3000 3200	130 Number of samples 132 134 Support type : block SBS, L=1 136 SBS, L=5 SBS, L=30 Greedy, L=1 Greedy, L=5 Greedy, L=30	Mean MSE	128 4000 5000 6000 7000 8000 9000	130 Number of samples 132 134 Support type : random SBS, L=1 136 SBS, L=5 SBS, L=30 Greedy, L=1 Greedy, L=5 Greedy, L=30
	Fig. 2.

circular shift. Proof. Consider a circular shift d of the selected P rows. Write (A d ) ns,ks = (A) ns+d,ks