Mikaël Jacquemont
email: jacquemont@lapp.in2p3.fr

Thomas Vuillaume
email: vuillaume@lapp.in2p3.fr

Alexandre Benoit
email: alexandre.benoit@univ-smb.fr

Gilles Maurin
email: maurin@lapp.in2p3.fr

Patrick Lambert
email: patrick.lambert@univ-smb.fr

Deep Learning for Astrophysics, Understanding the Impact of Attention on Variability Induced by Parameter Initialization

Keywords: multitasking, artificial neural networks, gamma rays, attention, visual explanation

In the astrophysics domain, the detection and description of gamma rays is a research direction for our understanding of the universe. Gamma-ray reconstruction from Cherenkov telescope data is multi-task by nature. The image recorded in the Cherenkov camera pixels relates to the type, energy, incoming direction and distance of a particle from a telescope observation. We propose γ-PhysNet, a physically inspired multi-task deep neural network for gamma/proton particle classification, and gamma energy and direction reconstruction. As ground truth does not exist for real data, γ-PhysNet is trained and evaluated on large-scale Monte Carlo simulations. Robustness is then crucial for the transfer of the performance to real data. Relying on a visual explanation method, we evaluate the influence of attention on the variability due to weight initialization, and how it helps improve the robustness of the model. All the experiments are conducted in the context of single telescope analysis for the Cherenkov Telescope Array simulated data analysis.

Introduction

Gamma-ray astronomy aims to study astronomical phenomena (supernova remnants, dark matter annihilation...) based on the gamma radiation generated by these phenomena. The analysis of this radiation is performed through the observation of telescope images of the particle shower resulting from the penetration of high-energy particles in the atmosphere (Cherenkov effect [START_REF] Hillas | Cerenkov light images of eas produced by primary gamma[END_REF], see Fig. 1). The purpose of the image analysis is twofold:

1. identifying the gamma rays in the cosmic ray background mainly composed of protons (with a signal-to-noise ratio typically lower than 1/1000), which is a classification problem, 2. estimating the energy and direction of the identified gamma rays, which is a regression problem.

The Cherenkov Telescope Array (CTA) 1 is the next generation of Imaging Atmospheric Cherenkov Telescopes (IACTs). Composed of ∼ 100 telescopes of different sizes, it will improve sensitivity and accuracy in gamma-ray analysis. However, the huge amount of data (210 PB of raw data per year when in full operation) requires a shift towards new methods of analysis, in particular deep neural network approaches. The work presented in this paper is carried out on the large-scale simulation data shared within the CTA international collaboration. As ground truth is not available in the field of astrophysics, the analysis toolchains are prepared with very high-quality simulation relying on a wellknown physical model of the phenomenon and on the detector simulation [START_REF] Bernlöhr | Monte carlo design studies for the cherenkov telescope array[END_REF]. Besides, the first real data are just available.

In this paper, we first propose a new deep multi-task architecture, named γ-PhysNet, taking into account physics considerations and designed for single telescope gamma event analysis. We evaluate this model on the simulation data of the Large Size Telescopes 1 (LST1 [START_REF] Ambrosi | for the CTA Consortium: The Cherenkov Telescope Array Large Size Telescope[END_REF]), the first prototype installed at the Northern CTA site in La Palma. In a second step, with the help of a visual explanation method for neural networks, we analyze the impact of attention (mechanism that reinforces relevant features) on the variability introduced by model weights initialization, and show that augmenting the network with attention improves the robustness of the model.

The rest of the paper is organized as follows. In Section 2, we give a short state of the art related to this work. Section 3 is a presentation of γ-PhysNet. Understanding the impact of attention mechanism and weight initialization are discussed in Section 4. Finally, Section 5 gives conclusions and some perspectives for future works.

Related Work

Over the past decade, deep learning has emerged as the leading approach in many computer vision tasks. IACT data has not escaped this trend [START_REF] Mangano | Extracting gamma-ray information from images with convolutional neural network methods on simulated cherenkov telescope array data[END_REF][START_REF] Nieto Castaño | Exploring deep learning as an event classification method for the Cherenkov Telescope Array[END_REF][START_REF] Parsons | Background rejection in atmospheric Cherenkov telescopes using recurrent convolutional neural networks[END_REF][START_REF] Shilon | Application of deep learning methods to analysis of imaging atmospheric cherenkov telescopes data[END_REF].

Fig. 1. Imaging Atmospheric Cherenkov Telescope. When a gamma particle enters the atmosphere, it generates an electromagnetic shower of secondary particles emitting Cherenkov light. This blue light is captured by the telescope on the ground and recorded as the signal by the camera. The event reconstruction consists in analyzing this signal to retrieve the physical parameters (type, energy, altitude, azimuth) of the primary particle.

These papers present promising results, especially for gamma/proton classification. However, they have handled the different reconstruction problems as single tasks, without considering their strong interdependence. Further, multiple single task models can increase the computational cost, thus making their use complicated in our massive data context.

Multi-Task Models. When dealing with a multi-task problem, it has been shown [START_REF] Thrun | Is learning the n-th thing any easier than learning the first?[END_REF] that transferring knowledge across related tasks improves the generalization of the deep models with fewer data. Thus, in Multi-Task Learning (MTL), the different tasks to address are trained simultaneously, using a partially shared modeling. In hard parameter sharing architectures, the most frequently used, a whole part of the network (generally the encoder or its first layers) is shared between all tasks [START_REF] Ruder | An overview of multi-task learning in deep neural networks[END_REF]. In soft parameter sharing architectures [START_REF] Cao | Partially shared multi-task convolutional neural network with local constraint for face attribute learning[END_REF], each task is learned with its own network. Some additional layers are shared and constrained in order to encourage their weights to be similar. In MTL, balancing the tasks is the critical point. For most of the MTL related papers [START_REF] Luvizon | 2d/3d pose estimation and action recognition using multitask deep learning[END_REF][START_REF] Ren | Cross-domain self-supervised multi-task feature learning using synthetic imagery[END_REF], this is done manually, requiring an extensive optimization process. However, adaptive methods have been proposed, with different balancing strategies such as modeling homoscedastic uncertainty for each task [START_REF] Kendall | Multi-task learning using uncertainty to weigh losses for scene geometry and semantics[END_REF], using the task loss gradients [START_REF] Chen | GradNorm: Gradient normalization for adaptive loss balancing in deep multitask networks[END_REF], using learning progress signals as key performance indicators [START_REF] Guo | Dynamic task prioritization for multitask learning[END_REF] or regarding the problem as a multi-objective optimization [START_REF] Sener | Multi-task learning as multi-objective optimization[END_REF].

Complementary to MTL, the development of deep learning models integrates attention mechanisms.

Attention Mechanisms. Attention is a mechanism that helps deep learning model focus on relevant features based on a defined context through trainable weights. A distinction can be made between restricted self-attention, focusing on local spatial neighborhoods and global self-attention as a non-local operation [START_REF] Wang | Non-local neural networks[END_REF].

Local and global self-attention can be considered as spatial attention mechanisms, as they capture short-and long-range dependencies in data, by weighting each pixel. On the contrary, Hu et al. [START_REF] Hu | Squeeze-and-excitation networks[END_REF] introduce a lightweight channel-wise attention denoted Squeeze-and-Excitation. The squeeze operation produces descriptors for each input channel, and is followed by an adaptive recalibration, the excitation, and a scale operation that weights the input channels. Dual Attention [START_REF] Sun | Saunet: Shape attentive u-net for interpretable medical image segmentation[END_REF] has been proposed to improve model interpretability and robustness of U-Net models for a semantic segmentation task. It combines Squeeze-and-Excitation and a simple spatial attention path. The latter first compresses the number of input channels to one, and then applies a sigmoid function to the resulting pixel values, and adds one to them to produce an attention map. This way, the spatial attention map can only increase pixel values. This map is then used to rescale the output of the Squeeze-and-Excitation. Such attention strategy then allows for the weighting of both spatial and channel-wise information in a low-cost operator.

Model Explanation. While traditional expert systems are highly interpretable and explainable, deep neural networks are often perceived as black boxes. Recently, the deep learning community has put an increasing effort on opening the black box. Some methods explore the role of individual neurons or linear combination of units through ablation [START_REF] Morcos | On the importance of single directions for generalization[END_REF][START_REF] Zhou | Revisiting the importance of individual units in cnns via ablation[END_REF] or optimization [START_REF] Olah | Feature visualization[END_REF]. Others estimate the importance of input features for a particular output activation. They produce saliency maps [START_REF] Cao | Look and think twice: Capturing top-down visual attention with feedback convolutional neural networks[END_REF][START_REF] Springenberg | Striving for simplicity: The all convolutional net[END_REF][START_REF] Srinivas | Full-gradient representation for neural network visualization[END_REF] also called heatmaps in [START_REF] Selvaraju | Gradcam: Visual explanations from deep networks via gradient-based localization[END_REF] for network visualization. In the rest of this paper we rely on Grad-CAM [START_REF] Selvaraju | Gradcam: Visual explanations from deep networks via gradient-based localization[END_REF], a highly class-discriminative localization method that produces a heatmap of the region that is the most relevant for the model's decision. It is applicable to any network structure. The heatmaps are computed based on the last convolutional layer. In the context of multi-task learning, Grad-CAM is then especially adapted to hard parameter sharing networks that share a convolutional encoder between all tasks. In addition, when available, attention maps, produced by attention mechanisms (see previous section), can provide insights on the model regions of interest that can explain the predictions.

Proposed Architecture and Performance

To achieve full event reconstruction from IACT data, we propose a MTL architecture, denoted γ-PhysNet. We then compare the vanilla model to a version augmented with attention modules. Our aim is to understand model robustness against weight initialization, and the regularization effect provided by attention mechanisms. Model performances are assessed relying on different random seeds.

In this paper we restrict our analysis to a specific model with optional attention modules. This model is a performing baseline that originates from preliminary extensive comparative studies, and the scope of the paper is to study its robustness towards weight initialization variability.

γ-PhysNet Architecture

As illustrated in Fig. 2, γ-PhysNet is a hard parameter sharing architecture composed of a backbone encoder and a multi-task block inspired by the physics of the reconstruction. The network is fed with two-channel IACT data discussed in Section 3.2. It performs gamma rays from background noise separation, and primary particle energy and arrival direction regression as the altitude and azimuth. The regression of the virtual impact point on the ground of the particle is an auxiliary task (see 1 for the meaning of these characteristics). Even though it is not required for higher-level analysis, physics shows that this parameter provides meaningful information to solve the other tasks. The baseline backbone of γ-PhysNet is the convolutional part of a ResNet-56 [START_REF] He | Identity mappings in deep residual networks[END_REF][START_REF] He | Deep residual learning for image recognition[END_REF], CIFAR-10 version, with full pre-activation implemented with IndexedConv [START_REF] Jacquemont | Indexed operations for non-rectangular lattices applied to convolutional neural networks[END_REF]. The latter allows the direct processing of the hexagonal pixel images of the LST1 data used in this paper without transforming them to square pixel ones. We also propose a refined version of the model, named γ-PhysNet DA (DA for Dual Attention), that makes use of attention. Dual Attention modules are inserted into the backbone after every stage, i.e., processing scale, to benefit from attention at each feature scale. As mentioned in Section 2, Dual Attention is composed of a spatial attention path and a channel-wise attention path. The latter consists of a Squeeze-and-Excitation module that has a ratio parameter to control its bottleneck. An extensive study of this hyperparameter has shown that, in the context of the experiments carried out for this paper, the default ratio of 16 allows obtaining the best results. The physically inspired multi-task block finalizes the model. It is composed of a global feature network and a local feature one, both based on fully connected layers. The global feature part, dedicated to energy regression, starts with a global average pooling. This strategy follows the physics of the phenomenon, as for a given arrival direction and impact point, the amplitude of the image is roughly proportional to the primary gamma ray energy [START_REF] Völk | Imaging very high energy gamma-ray telescopes[END_REF]. The local feature part is devoted to gamma/proton classification, and regression of the arrival direction and impact point. It is fed with the flattened feature maps produced by the backbone. The aim is to exploit local (the shape of the signal in the image) and spatial (the position and orientation of the signal) information that is more deeply related to the particle type, its arrival direction and virtual impact point.

Experiments

Dataset. As ground truth is impossible to obtain from real data, γ-PhysNet and γ-PhysNet DA are evaluated on the LST4 mono-trigger production (from 2019/04/15), the reference large-scale Monte Carlo production generated by the LST international consortium for the LST1 commissioning. This dataset has been calibrated and integrated with DL1DataHandler [START_REF] Kim | DL1-Data-Handler: DL1 HDF5 writer, reader, and processor for IACT data[END_REF]. Each sample corresponds to a single event (particle) and is composed of two-channel images: one with pixel intensities in number of photoelectrons and the other one containing per-pixel temporal information in nanoseconds. Data amplitude is not normalized since it is related to the energy of the detected particles. A data selection step is applied, following the standards in the domain and the project collaboration. It consists of a series of relatively loose selection cuts on image amplitude, shower size and truncated showers applied to the data in order to discard bad quality events that would not be reconstructed by standard methods either. The training set is finally composed of 874k gamma events and 506k proton events, the validating set 201k and 136k, and the test set 209k and 38k.

Training. In order to analyze the robustness of both models to parameter initialization, we repeat the experiments with six different random seeds. We define the following optimization criteria: cross-entropy loss for the classification task and the L1 loss for each regression task. Models are trained for 25 epochs with Adam [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF] as the optimizer, and a weight decay of 10 -4 for regularization purpose. The learning rate is set to 10 -3 , decayed by a factor of 10 every 10 epochs. The different tasks are balanced with the uncertainty estimation method presented in [START_REF] Kendall | Multi-task learning using uncertainty to weigh losses for scene geometry and semantics[END_REF]. In this setup, both models reach their best performance plateau on the validation set. In gamma-ray astronomy, proton events are considered as background noise. To prevent them from penalizing the learning of energy and direction regression for gamma events, we rely on a masked loss method, setting to zero the loss of the regression parameters (energy, arrival direction and impact point) when particles are protons.

Evaluation Metrics. The performance on energy and direction reconstruction tasks is measured through resolution curves. The angular resolution is defined as the 68% containment radius of the point-spread (distribution) function and the energy resolution as the 68% containment radius of the relative absolute deviation. Lower resolutions are better.

Weight initialization plays an important role in neural network performance. We then repeat the experiment six times for both models, and we illustrate the variability of these different runs by drawing the average resolution curve per energy bin, the surfaces representing the standard deviation. The latter, referred to as dispersion in this paper, serves as measure of the robustness of the models.

For the gamma/proton classification task, the overall performance of the network is given by the area under the ROC curve (AUC), the precision and the recall.

Results. As shown in Table 1, both models with and without attention have comparable results on the classification task, within the standard deviation range. However, for the energy and direction regression, as illustrated in Fig. 3, the model with Dual Attention (γ-PhysNet DA) obtains slightly better average results above 100 GeV. Furthermore, we observe that the network with attention has significantly less spread results. In particular, γ-PhysNet DA has a constantly lower dispersion on the direction reconstruction task. On the energy reconstruction one, at energies above 200 GeV, the model without attention has dispersion up to three times higher. This lower dispersion of the results of γ-PhysNet DA denotes a better robustness to parameter initialization. This will lead to a more reliable estimation of the particle parameters when we analyze real data.

Understanding the Impact of Dual Attention

As we have seen in Section 3.2, the addition of Dual Attention modules to γ-PhysNet improves the robustness of the model, especially for the energy and arrival direction reconstruction tasks, by reducing the dispersion due to parameter initialization during the training step. However, it has no clear impact on the classification task. To understand how and why attention acts on the model predictions, we carefully observe the Grad-CAM heatmaps produced by γ-PhysNet Table 1. Classification performance of both models with and without attention. The AUC represents the overall performance, while the precision shows the ability of the model to discard protons (the background noise), and the recall highlights the ability of the model to retrieve gammas. and γ-PhysNet DA for 25 well and badly reconstructed events of the test set. We also observe the spatial attention maps of the three Dual Attention modules.

Additionally, for each input data and model seed, we combine these three spatial attention maps with the Hadamard product to obtain a global representation of the spatial attention in the model. Finally, to highlight the variability brought by the different initialization seeds for a particular event, we compute the mean and standard deviation of Grad-CAM heatmaps and spatial attention maps at the pixel level (to lighten the paper, only these means and standard deviations will be presented in the following figures).

Observation of the Grad-CAM Heatmaps and the Spatial Attention Maps. For all the 25 events analyzed, we observe a common trend in the averaged heatmaps produced by the Grad-CAM. Fig. 4, Fig. 5 and Fig. 6 illustrate 3 typical examples respectively, a well-reconstructed but partially truncated gamma shower, a well-reconstructed and centered gamma shower and a badly reconstructed one.

It is worth noticing that the different maps are represented using different color scales. For the classification task (denoted "class" in the figures), with or attention, the most relevant pixels highlighted by Grad-CAM are located in the signal area. On another hand, for the regression tasks (denoted "Energy", "Altitude" and "Azimuth" in the figures) we observe different behaviors. Without attention, the most relevant ones are situated on the border of the camera, while the signal pixels are of less importance. In our understanding, this phenomenon serves as an evaluation of the signal outside of the camera, and thus not acquired, that is important for regression tasks. With attention, all the relevant pixels are located in the shower area, thus better taking into account the signal pixels and relevant information. Besides, for all tasks, the model with attention focuses on a larger part of the image. It is worth noticing that, in addition to the signal pixels themselves, pixels situated in the signal neighborhood contain useful information about the shower shape. Then, the deviation measures (denoted "std" in the figures) of Grad-CAM heatmaps highlight the same general trend related to the robustness against the initialization. Although the variability of pixel relevance is quite important in both cases, for γ-PhysNet with attention the relevant pixels fluctuate in the signal area, while without attention they vary on a larger extent between the shower area and the image boundaries.

Then, the observation of the spatial attention maps of γ-PhysNet DA shows that the action of Dual Attention is different depending on the feature scale. The output of stage 1 has the same resolution as the input data. At this scale, the attention module mainly focuses on the shower pixels, and strongly rescales their value. The value of the signal pixels is multiplied by up to 1.8 in average, while, by construction, the rescaling factor computed by the spatial attention path of Dual Attention modules ranges in [START_REF] Ambrosi | for the CTA Consortium: The Cherenkov Telescope Array Large Size Telescope[END_REF][START_REF] Bernlöhr | Monte carlo design studies for the cherenkov telescope array[END_REF], as explained in Section 2. However, the attention maps are quite noisy at this scale, as also highlighted by the deviation measures. After stage 2, the attention modules also strongly favor the pixels in the signal area, and the attention maps are less noisy. Then, the last attention modules, after stage 3, have a lighter impact on the feature maps values. Finally, the observation of the combined spatial attention maps highlights that attention strongly helps the model focus on the signal pixels, which is consistent with the observation of the Grad-CAM heatmaps.

Impact of the Attention on the Classification Task. Table 1 shows no significant effect of the attention of the classification task. This can be explained by the attention and Grad-Cam heatmaps that always focus on the event shower area. More into the details, attention allows for a larger extent of the heatmaps around the shower and also avoids interest on the noisy boundaries of the images. However since both models report similar predictions performances, one can expect that they rely on the same most contributing, shower centered, features.

Impact on the Regression Tasks. For the energy and direction reconstruction tasks, Dual Attention forces the model to focus on the shower area instead of the image border, by strongly rescaling the signal pixels and their close neigh-bors. This is significant enough to improve the results of both tasks presented in Section 3.2, and to reduce the dispersion introduced by the parameter initialization.

Besides, it is worth noticing that for the event represented in Fig. 4, the most important pixels are still on the border of the image. Indeed, the signal is truncated, and the network has learned that it is an important information to take into account for the estimation of the energy and the direction of the gamma ray. However, focusing on the signal part of the image does not automatically imply a good reconstruction, as exemplified by the badly reconstructed event 24 shown in Fig. 6.

Conclusion

In this paper we have presented γ-PhysNet, a deep multi-task architecture for gamma-ray full-event reconstruction for IACT single telescope images. We have shown that augmenting the model with attention allows for a reduction of the performance variability induced by parameter initialization. Relying on a visual explanation method, we have then realized the first steps to understand how attention modifies the behavior of the model. In a future work, in order to deepen this analysis, we will define statistical criteria to quantify this effect. Correlating the pixels highlighted by Grad-CAM and the true signal over the whole test set is an option. We also consider analyzing statistically the distance of the relevant pixels to the shower centroid. Finally, it would also be interesting to study if the robustness brought by attention mechanisms also holds with slightly different datasets (different levels of noise, altered images, etc.).

Fig. 2 .

 2 Fig. 2. γ-PhysNet. Left: the model architecture composed of a convolutional backbone (ResNet-56) and a physically inspired Multi-task head block based on fully connected layers (FC). The latter is divided into two paths: a global path for the energy regression and a local path for the particle type classification, the direction and the virtual impact point regression. Right: refined model backbone with Dual Attention modules inserted after each stage.

Fig. 3 .

 3 Fig. 3. Angular (left) and energy (right) resolution as a function of the energy in the LST1 energy range (lower is better). Both resolution curves represent the error of the model for the regression of respectively the direction and the energy of the detected gamma particle. The dispersion, representing the induced by weight initialization, is a measure of the robustness of the models.

Fig. 4 .

 4 Fig. 4. Well-reconstructed gamma, event 16. The two upper rows represent the Grad-CAM heatmaps (mean and standard deviation) obtained for the vanilla γ-PhysNet, the next two rows correspond to the Grad-CAM heatmaps of γ-PhysNet (with attention), and the last two rows show the corresponding spatial attention maps.

Fig. 5 .

 5 Fig. 5. Well-reconstructed gamma, event 41. The two upper rows represent the Grad-CAM heatmaps (mean and standard deviation) obtained for the vanilla γ-PhysNet, the next two rows correspond to the Grad-CAM heatmaps of γ-PhysNet DA (with attention), and the last two rows show the corresponding spatial attention maps.

Fig. 6 .

 6 Fig. 6. Badly reconstructed gamma, event 24. The two upper rows represent the Grad-CAM heatmaps (mean and standard deviation) obtained for the vanilla γ-PhysNet, the next two rows correspond to the Grad-CAM heatmaps of γ-PhysNet DA (with attention), and the last two rows show the corresponding spatial attention maps.

https://www.cta-observatory.org/

We gratefully acknowledge financial support from the agencies and organizations listed here: www.cta-observatory.org/consortium acknowledgment. This project has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 653477, and from the Fondation Université Savoie Mont Blanc. This work has been done thanks to the facilities offered by the Univ. Savoie Mont Blanc -CNRS/IN2P3 MUST computing center and HPC resources from GENCI-IDRIS (Grant 2020-AD011011577) and computing and data processing ressources from the CNRS/IN2P3 Computing Center (Lyon -France). We gratefully acknowledge the support of the NVIDIA Corporation with the donation of one NVIDIA P6000 GPU for this research.


Jérôme Gilleron
email: gilleron@unice.fr

Cindy Meziat
email: cindy.meziat@univ-cotedazur.fr

André Sulen
email: andre.sulen@uib.no

Stoyan Ivanov
email: stoyan.ivanov@univ-cotedazur.fr

Jennifer Jager
email: jennifer.jager@univ-cotedazur.fr

David Estève
email: david.esteve@ipbs.fr

Catherine Muller
email: catherine.muller@ipbs.fr

Jean-François Tanti

Mireille Cormont
email: mireille.cormont@unice.fr

Exploring Adipose Tissue Structure by Methylsalicylate Clearing and 3D Imaging

Keywords: Adipose tissue, clearing, light microscopy, 3-D whole tissue, obesity, morphology

Here, we describe a simple, inexpensive and fast clearing method to resolve the 3D structure of both mouse and human white adipose tissue using a combination of markers to visualize vasculature, nuclei, immune cells, neurons, and lipid-droplet coat proteins by fluorescent imaging.

to expand during obesity is pivotal to the development of obesity-associated pathologies. This AT expansion is an important homeostatic mechanism to enable adaptation to an excess of energy intake and to avoid deleterious lipid spillover to other metabolic organs, such as muscle and liver. Therefore, understanding the structural remodeling that leads to the failure of AT expansion is a fundamental question with high clinical applicability. In this article, we describe a simple and fast clearing method that is routinely used in the laboratory to explore the morphology of mouse and human white adipose tissue by fluorescent imaging. This optimized AT clearing method is easily performed in any standard laboratory equipped with a chemical hood, a temperature-controlled orbital shaker and a fluorescent microscope. Moreover, the chemical compounds used are readily available. Importantly, this method allows one to resolve the 3D AT structure by staining various markers to specifically visualize the adipocytes, the neuronal and vascular networks, and the innate and adaptive immune cells distribution.

INTRODUCTION:

Obesity is characterized by an increase in adipose tissue mass and has become a major worldwide public health issue, given that people with obesity have increased risk of developing cardiovascular disease, type-2 diabetes, liver diseases and some cancers.

A fundamental physiological function of adipose tissue is to modulate whole-body glucose and lipid homeostasis 1,2 . During the feeding period, the adipocytes (i.e., the main cells of the adipose tissue) stores the excess glucose and lipids provided by a meal into triglycerides. During fasting, the adipocytes break down the triglycerides into non-esterified fatty acids and glycerol to sustain the energy demand of the body. During the development of obesity, adipose tissue expand by increasing the size (hypertrophia) and/or the number (hyperplasia) of adipocytes 1 , to increase their storage capacity. When the expansion of adipose tissue reaches its limit, a constant highly variable among patients, the remaining lipids accumulates into other metabolic organs including muscles and liver 3,4 , leading to their functional failure and initiating obesityrelated cardio-metabolic complications 1,[START_REF] Moreno-Indias | Impaired adipose tissue expandability and lipogenic capacities as ones of the main causes of metabolic disorders[END_REF] . Therefore, identifying the mechanisms that govern adipose tissue expansion is a key clinical challenge.

The morphological modifications documented within adipose tissues during obesity are linked to its pathological dysfunction. Several staining procedures have been used to describe the tissue organization of the adipose tissue, including actin [START_REF] Vergoni | DNA Damage and the Activation of the p53 Pathway Mediate Alterations in Metabolic and Secretory Functions of Adipocytes[END_REF] , vascular markers [START_REF] Xue | Preventing diet-induced obesity in mice by adipose tissue transformation and angiogenesis using targeted nanoparticles[END_REF] , lipid-droplet markers [START_REF] Zwick | Adipocyte hypertrophy and lipid dynamics underlie mammary gland remodeling after lactation[END_REF] , and specific immune cell markers [START_REF] Zhang | The inflammatory changes of adipose tissue in late pregnant mice[END_REF][START_REF] Yang | Obesity increases the production of proinflammatory mediators from adipose tissue T cells and compromises TCR repertoire diversity: implications for systemic inflammation and insulin resistance[END_REF] . However, because of the huge diameter of adipocytes (50 to 200 µm) [START_REF] Laforest | Comparative analysis of three human adipocyte size measurement methods and their relevance for cardiometabolic risk[END_REF] , it is essential to analyze a large portion of the whole tissue in three dimensions in order to accurately analyze the dramatic structural AT changes observed during obesity. However, because the light does not penetrate an opaque tissue, imaging in 3D within a large tissue samples using fluorescence microscopy is not possible. Methods of tissue clearing to make them transparent have been reported in the literature (for a review, see [START_REF] Azaripour | A survey of clearing techniques for 3D imaging of tissues with special reference to connective tissue[END_REF] ) allowing one to clear tissues and to perform in-depth, whole tissue fluorescence microscopy. These methods offer unprecedented opportunities to assess the 3D cellular organization in healthy and diseased tissue. Each of the described methods have advantages and drawbacks, and therefore need to be carefully selected depending on the studied tissue (for a review, see [START_REF] Richardson | Clarifying Tissue Clearing[END_REF] ). Indeed, some approaches require a long incubation period and/or the use of materials or compounds that are either expensive, toxic or difficult to obtain [START_REF] Chi | Adipo-Clear: A Tissue Clearing Method for Three-Dimensional Imaging of Adipose Tissue[END_REF][START_REF] Roberts | Optical Clearing of the Mouse Central Nervous System Using Passive CLARITY[END_REF][START_REF] Woo | Optimization of the optical transparency of rodent tissues by modified PACT-based passive clearing[END_REF][START_REF] Zhang | 3D imaging of optically cleared tissue using a simplified CLARITY method and on-chip microscopy[END_REF][START_REF] Ke | Optical clearing of fixed brain samples using SeeDB[END_REF][START_REF] Hahn | High-resolution imaging of fluorescent whole mouse brains using stabilised organic media (sDISCO)[END_REF] . Taking advantage of one of the first compounds used a century ago by Werner Spalteholz to clear tissues [START_REF] Spalteholz | Über das Durchsichtigmachen von Menchlichen und Tierichen Präparaten und Seine Theoretischen Bedingungen[END_REF] , we set up a user-friendly and inexpensive protocol that is very well adapted for the clearing of all mouse and human adipose tissue depots in any laboratory with typical equipment including a chemical hood, a temperature-controlled orbital shaker and a confocal microscope.

PROTOCOL:

This protocol was tested and is validated for all mouse and human white adipose tissue depots. Human and mouse adipose tissues were collected accordingly to European laws and approved by French and Swedish Ethical committees.

1.

Fixation of mouse and human white adipose tissue 1.1. Immerse the harvested mouse or human white adipose tissues in at least 10 mL of PBS containing 4% paraformaldehyde (PFA) in a 15 mL plastic tube.

1.2. Shake the plastic tube at room temperature on a rolling plate for 1 h.

1.3. Leave the plastic tube at 4 °C on a rolling plate overnight, to complete the fixation.

NOTE: This protocol is adapted for large adipose tissue samples, such as whole epididymal fat pads obtained from mice fed a normal diet (≈250 mg -≈0.6 cm 3 ). For larger samples like epididymal fat pads obtained from mice fed a high fat diet (≈1.5 g -≈4 cm 3 ) or human adipose tissue samples, although the clearing protocol should perfectly work for the whole sample (by scaling-up the PFA and the antibody mixtures), in general, we cut tissue pieces around 1 g (≈2.5 cm 3 ) to reserve the remaining samples either for additional staining or applications. For the PFA (step 1.1.) we recommend using roughly 10 times the volume of the tissue. For the antibody mixtures (steps 3.1. and 3.4.), increase the volume to completely immerse the tissue, and use a 15 mL plastic tube (see

3.

Staining procedure for mouse and human white adipose tissue 3.1. Transfer the tissue in a 1.5 mL plastic microtube containing 300 µL of PBS supplemented with 0.2% Triton X-100, 10% DMSO, 3% BSA and the primary antibodies (10x more concentrated than for cryosection staining but optimal antibody concentration should be evaluated for each antibody). Protect the tube from light by covering with aluminium foil and shake the tube at 37 °C in a temperature-controlled orbital shaker at 100 rpm for at least 2 days (see the note in step 1.1 and step 3.7.1).

3.2. Rinse the tissue in a 15 mL plastic tube containing 10 mL of PBS supplemented with 0.2% Triton X-100, 10% DMSO and 3% BSA and shake the tube protected from light at 37 °C in a temperature-controlled orbital shaker at 100 rpm for 5 h. Perform this step twice.

3.3. Rinse the tissue in a 15 mL plastic tube containing 10 mL of PBS supplemented with 0.2% Triton X-100, 10% DMSO and 3% BSA and shake the tube, protected from light, at 37 °C in a temperature-controlled orbital shaker at 100 rpm for one night to two days.

3.4. Transfer the tissue to a 1.5 mL plastic microtube containing 300 µL of PBS supplemented with 0.2% Triton X-100, 10% DMSO, 3% BSA and the secondary antibodies (10x more concentrated than for cryosection staining). Protect the tube from light by covering with aluminium foil and shake the tube at 37 °C in a temperature-controlled orbital shaker at 100 rpm for at least 2 days (see the note in step 1.1 and step 3.7.1).

3.5. Rinse the tissue in a 15 mL plastic tube containing 10 mL of PBS supplemented with 0.2% Triton X-100, 10% DMSO, and 3% BSA and shake the tube, protected from light, at 37 °C in a temperature-controlled orbital shaker at 100 rpm for 5 h. Perform this step twice.

3.6. Rinse the tissue in a 15 mL plastic tube containing 10 mL of PBS supplemented with 0.2% Triton X-100, 10% DMSO, and 3% BSA and shake the tube, protected from light, at 37 °C in a temperature-controlled orbital shaker at 100 rpm for one night to two days.

3.7.

Rinse the tissue in a 15 mL plastic tube containing 10 mL of PBS and shake the tube, protected from light, at 37 °C in a temperature-controlled orbital shaker at 100 rpm for 2 h.

3.7.1. For the labelling of specific structures such as nuclei or actin, add 4′,6-diamidino-2phenylindole (DAPI or equivalent) and/or the fluorescently labelled-phalloidin at step 3.1. when only fluorescently labelled primary antibodies are used, or at step 3.4. when secondary antibodies are used.

NOTE: For the staining procedure, primary antibodies already conjugated with fluorochromes (like those used for flow cytometry) can be used and we recommend it for the gain of time and specificity. Indeed, even if mouse primary antibodies can be used followed by secondary antimouse antibodies, as we demonstrated it here, we have often observed non-specific staining of blood vessels due to circulating immunoglobulin. Therefore, to avoid this issue, primary antibodies that are already labelled are highly recommended and here we provide evidence that the antibodies used in flow cytometry are compatible with our procedure. However, in the specific case of an antigen that is expressed at low levels, a signal amplification step via secondary antibodies is mandatory; when the only available primary antibody is made in mouse, a cardiac-perfusion of the mice with PBS for at least 5 min at the sacrifice can remove a large proportion of circulating immunoglobulin and thus the non-specific staining.

4.

Clearing procedure for mouse and human white adipose tissue 4.1. Immerse the tissue in a 15 mL plastic tube containing 10 mL of 50% ethanol and shake the tube, protected from light, at room temperature in an orbital shaker at 100 rpm for 2 h. 4.2. Immerse the tissue in a 15 mL plastic tube containing 10 mL of 70% ethanol and shake the tube, protected from light, at room temperature in an orbital shaker at 100 rpm for 2 h. 4.3. Immerse the tissue in a 15 mL plastic tube containing 10 mL of 95% ethanol and shake the tube, protected from light, at room temperature in an orbital shaker at 100 rpm for 2 h.

4.4. Immerse the tissue in a 15 mL plastic tube containing 10 mL of 100% ethanol and shake the tube, protected from light, at room temperature in an orbital shaker at 100 rpm for 2 h. 4.5. Immerse the tissue in a 15 mL plastic tube containing 10 mL of 100% ethanol and shake the tube, protected from light, at room temperature in an orbital shaker at 100 rpm overnight.

4.6. Immerse the tissue in a 20 mL glass bottle with a plastic cap (see table of materials) containing 5 mL of methyl salicylate (see table of materials) under a chemical hood and shake the glass container, protected from light, at room temperature in an orbital shaker at 100 rpm for at least 2 h.

NOTE:

The clearing procedure could be significantly accelerated (if necessary) by avoiding steps 4.3. and 4.5., although the final clearing quality would be slightly reduced.

3D-confocal imaging of cleared white adipose tissue

5.1. Transfer the tissue to a metallic imaging chamber equipped with a glass bottom (see 5.2. Place the imaging chamber on an inverted confocal microscope.

5.3. Image the tissue using a low magnification objective (e.g., 4x objective) to generate a few cm 3 3D maps of the whole adipose tissue or of the human tissue sample.

5.4. Select several areas for the tissue sampling at higher magnification. Typically, use a 20x long distance air objective that provides a good ratio between resolution and depth. We acquire large mosaic images with z-stacks between 600 and 2000 µm depth.

5.4.1. Use a long-distance objective to image deeper into the tissue.

Extraction of quantitative results from the 3D adipose tissue images

NOTE: The segmentation of the different structures and the subsequent extraction of the quantitative information from the 3D-image stack generated in point 5 can be performed using any of the many existing image analysis software options, either commercial or freeware. In the following points, we describe a strategy that is routinely used in our laboratory to extract quantitative information from 3D adipose tissue images using commercial software (see Table of Materials).

6.1. Convert the 3D stacks onto the software format to free-up memory space.

6.2. Segment the cells.

6.2.1. Open the Cell module of the software.

6.2.2. Change the Cell Detection setting to plasma membrane staining.

6.2.3. Choose the fluorescent channel of the marker used to delineate the cell periphery (phalloidin which labels cortical actin or plasma membrane markers such as F4/80 for macrophage, TCR-β for T cells, or cluster of differentiation CD proteins specific for immune cell subtypes).

6.2.4. Set up the thresholds and provide a range of expected cell size (i.e., 1 to 200 µm for cell detection).

6.2.5. Run the segmentation. A volume corresponding to the z-stack will be generated with the segmented cells color-coded with a different color for each of the neighbouring cells.

6.2.6. Apply statistical filters (size, roundness, circularity, etc.) in the statistic tab to exclude segmentation artefacts and/or to refine the segmented cells to the cell of interest based on their size (i.e., 20-200 µm for adipocytes; 5-25 µm for macrophages; 1-3 µm for lymphocytes).

6.2.7. Extract measurements and quantitative data (volume, number, size, diameter, etc.) from the statistics tab.

6.3. Segment cellular components (nuclei, vesicles, etc.) or vessels as a structure.

NOTE: Although the filament segmentation is very useful to study the connectivity of the vessels, we use the surface module segmentation to extract data on the size, volume and diameters of the vessels. Indeed, the quantification of the sizes of the vessels is lost when using the filaments module.

6.3.1. Open the Surface module of the software.

6.3.2. Select the fluorescent channel of the marker used to specifically label the subcellular component to reconstruct it in 3D.

6.3.3. Set up the thresholds and provide a range of the expected diameter size of the structure (i.e., 0.5-5 µm for nuclei; 0-100 µm for vessels; etc.).

6.3.4. Run the segmentation. A volume with the segmented cellular structure will be generated. Measurements and quantitative data (volume, number, size, diameter, etc.) can be extracted from the statistics tab.

6.4. Segment the vessels as a tubular continuum.

NOTE: Although the filament segmentation is very useful to study the connectivity of the vessels, we use the surface module segmentation to extract data on the size, volume and diameters of the vessels. Indeed, the quantification of the sizes of the vessels is lost when using the filaments module.

6.4.1. Open the Filaments module of the software.

6.4.2. Select the fluorescent channel corresponding to the vessel staining.

6.4.3. Set up the thresholds for the fluorescence intensity and select the appropriate number of expected connecting nodes.

6.4.4. Run the segmentation. Filaments representing vessel network will be generated in 3D.

Measurements can be extracted from the statistics tab, allowing to obtain quantitative data including vessel length, number of vessel branching, etc.

REPRESENTATIVE RESULTS:

Using the procedure described here and summarized in Figure 1, we were able to stain and optically clear human and mouse white adipose tissue as presented in Figure 2A and Figure 2B, respectively. The cleared tissue was transferred to the metallic imaging chamber to perform confocal imaging (Figure 3A). The clearing drastically improved the depth of the tissue images that we were able to acquire (Figure 3B and Figure 3C). The whole adipose tissue can be acquired in 3D at low magnification using a 4x objective (see supplementary movie 1 and supplementary movie 2), to generate a 3D map, enabling to select the different areas to be acquired at high magnification using a 20x objective (Figure 3D). The high magnification images are acquired in 3D with a depth of 2 mm (Figure 3E).

This procedure enables the labelling of numerous general cell markers, including nuclei by using DAPI and actin by using fluorescently-labelled Phalloidin. Specific staining of lipid droplets and adipocytes can be achieved by using perilipin antibody (see 4E). The peripheral nerve network can be detected by using the anti-Tyrosine Hydroxylase (TH) antibody (see table of materials; Figure 4F-G).

Additionally, this protocol allows the clearing of mouse epididymal adipose tissue (Figure 4A-B, E), mouse subcutaneous adipose tissue (Figure 4D, G), mouse brown adipose tissue (Figure 4F), and human adipose tissue (Figure 4C). Using a combination of these labelling and a commercial software (see table of materials) to segment these markers, we can determine within the adipose tissue i) adipocytes mean size and size distribution (Figure 5A-B) and ii) blood vessel network density (Figure 5C).

FIGURE LEGENDS:

Figure 1. Summary of the clearing procedure for white adipose tissue. Mouse or human white adipose tissues are fixed, permeabilized, saturated, stained and cleared. Images are then acquired in 3D and analyzed using a commercial software. (G) Z-projection of 600 µm z-stack of mouse subcutaneous white adipose tissue stained for nuclei and neuronal network using DAPI, anti-TH antibody and anti-rabbit-alexa647-conjugated antibody.

Figure 5. Analysis of 3D images using commercially available software (see table of materials).

(A) 3D volume rendering of human adipocytes segmented in 3D by commercially available software. Every adipocyte has a different color to its contacting adjacent adipocytes. Vessels are represented in red. (B) Size quantification and distribution of the adipocytes from the 3D volume rendering of panel A, which contains around 20,000 adipocytes. (C) 3D volume rendering of blood vessels segmented in 3D using commercially available software and represented in red or yellow for the large and small vessels, respectively.

Supplementary Movie 1. 3D volume rendering of the whole subcutaneous white adipose tissue shown in Figure 3D. The white box represents a 2 cm 3 cube.

DISCUSSION:

The modifications that occur within the adipose tissue over the course of pathological progression, such as that of obesity, is fundamental to the understanding of the mechanisms behind the pathology. Pioneering studies that revealed such mechanisms in adipose tissue have been based on global approaches such as whole adipose tissue proteomics [START_REF] Shields | Differential Adipose Tissue Proteomics[END_REF] , flow cytometry [START_REF] Bourlier | Remodeling phenotype of human subcutaneous adipose tissue macrophages[END_REF][START_REF] Hagberg | Flow Cytometry of Mouse and Human Adipocytes for the Analysis of Browning and Cellular Heterogeneity[END_REF] , and transcriptomics [START_REF] Hill | Distinct macrophage populations direct inflammatory versus physiological changes in adipose tissue[END_REF][START_REF] Acosta | Single cell transcriptomics suggest that human adipocyte progenitor cells constitute a homogeneous cell population[END_REF] . In addition, efforts have been made to explore the structural changes occurring in adipose tissue using histological analyses [START_REF] Coppack | Adipose tissue changes in obesity[END_REF][START_REF] Cancello | Reduction of macrophage infiltration and chemoattractant gene expression changes in white adipose tissue of morbidly obese subjects after surgery-induced weight loss[END_REF][START_REF] Cinti | Adipocyte differentiation and transdifferentiation: plasticity of the adipose organ[END_REF][START_REF] Wellen | Obesity-induced inflammatory changes in adipose tissue[END_REF] . However, the analysis of 5-10 µm adipose tissue sections, because of the limited size of the section, restricts the ability to appreciate important structural features. First, only a few adipocyte nuclei are detectable per section since each section represents only a small portion of the adipocytes (1/10 th ). Second, the size of the adipocytes estimated from adipose tissue sections is inaccurate because most adipocytes are cut below or behind their equator, leading to a biased (under)measurement of their mean size. Third, the blood vessel and nerve fiber continuum are lost during physical sectioning. Fourth, the distribution of each subtypes of immune cells within the tissue is difficult to establish because the three-dimensional coordinates are lost. In this context, the methodology we propose here allows any standard laboratory to image human and mouse white adipose tissue in three-dimensions, in a simple and inexpensive manner.

This method does have some limitations. First, the time required to prepare the tissue (fixation, permeabilization, staining, clearing) is long (more than a week). This would be difficult to reduce because the majority of this time is required to stain the tissue. The penetration of antibodies within such large tissue samples requires long incubation times. Decreasing the incubation time would increase the risk of having non-homogenous antibody penetration, which would lead to staining artifacts. Therefore, the only possible window to gain time and shorten the procedure is the time dedicated to clear the tissue, which takes around one day when optimal results are required, but this can be decreased to half a day without a dramatic drop in quality. The second limitation is the probable shrinkage of the tissue. Indeed, in any protocol that dehydrates tissue using solvents, it is likely that the tissues shrink due to water removal. Estimating this shrinkage is always challenging, and it is nearly impossible here because the edge of the tissue becomes transparent when lipids start to be extracted during the dehydration process. Hence, the size estimations of the cleared tissues need to be interpreted with caution. However, comparing changes in adipocyte sizes or vessel length between tissues derived from mouse with different genotypes and/or submitted to different environmental conditions remain informative [START_REF] Laforest | Comparative analysis of three human adipocyte size measurement methods and their relevance for cardiometabolic risk[END_REF] . The last limitation is linked to the large volume of the whole mouse adipose tissue or of a large human adipose tissue surgical sample. Indeed, although the protocol is perfectly adapted to clear these large samples, imaging a whole organ is challenging with a confocal microscope. The strategy to overcome this issue and to exploit the full potential of the whole organ clearing procedure, is to acquire a 3D low magnification map of the whole organ by using a 4x objective, and to select specific areas randomly dispersed within the tissue where higher magnification 3D images using a 20x long-distance objective can be acquired. The "high magnification" setup allows imaging the tissue volume of around 45 mm 3 (4.7 x 4.7 x 2 mm). Although this volume is limited compared to the volume of the whole organ (0.5 to more than 4cm 3 ), repeating the acquisitions at several positions within the tissue allows us to obtain a good sampling of the tissue at cellular to sub-cellular resolution. Note that imaging large tissues will generate a significant amount of imaging data that will need to be stored.

The procedure has significant differences when compared to methods that have previously been described to clear adipose tissue [START_REF] Chi | Adipo-Clear: A Tissue Clearing Method for Three-Dimensional Imaging of Adipose Tissue[END_REF][START_REF] Li | Co-staining Blood Vessels and Nerve Fibers in Adipose Tissue[END_REF] . First of all, unlike AdipoClear, that uses methanol, dichloromethane, and dibenzylether [START_REF] Chi | Adipo-Clear: A Tissue Clearing Method for Three-Dimensional Imaging of Adipose Tissue[END_REF] , the method presented here uses ethanol for dehydration and methyl salicylate for clearing. Therefore, the procedure is safer because it takes advantage of less toxic clearing solvents that are often used by food/drink-processing industries (ethanol and methylsalicylate), in comparison to the high toxicity of dichloromethane, methanol and dibenzylether. Moreover, the preparation of the tissue according to the protocol is two days faster than the AdipoClear protocol [START_REF] Chi | Adipo-Clear: A Tissue Clearing Method for Three-Dimensional Imaging of Adipose Tissue[END_REF] . More recently, another method was proposed by Li and colleagues to clear adipose tissue pieces by immersing them into glycerol [START_REF] Li | Co-staining Blood Vessels and Nerve Fibers in Adipose Tissue[END_REF] . This protocol is very simple even compared to thid procedure, and the quality of the images is good even at high magnification. However, this clearing protocol only works efficiently when using 2 mm 3 adipose tissue pieces. The sectioning of the tissue into these tiny samples prior to the clearing procedure prevents one from imaging tissue volumes larger than 2 mm 3 even at low magnification. The procedure presented here allows mapping of the whole tissue in 3D at low magnification and acquisition of 3D stacks of images around 45 mm 3 at several known positions within the whole cleared adipose tissue at high magnification.

This procedure may have several future applications. As with any method, the procedure described here can be simplified and/or further optimized. An interesting advance for the procedure could come from the combination of the clearing process, we described, with additional imaging approaches. For example, specific imaging setups, such as Optical Projection Tomography (OPT), Total Internal Reflection Fluorescence microscopy (TIRF), Selective Plan Illumination Microscopy (SPIM) or Stimulated Emission Depletion microscopy (STED) are in principle compatible with the procedure, although this will limit its applicability to those laboratories that are equipped with these systems. Alternatively, using an objective with a high numerical aperture index and an acquisition system that is able to acquire hundreds of images per second, which is found in many laboratories, one could perform super-resolution analyses by using the Super Resolution Radial Fluctuations (SRRF) post-processing image analysis freeware [START_REF] Gustafsson | Fast live-cell conventional fluorophore nanoscopy with ImageJ through super-resolution radial fluctuations[END_REF] . Interestingly, classical fluorochromes are adapted for the SRRF analysis, which would enable 3D super-resolution analyses of whole human and mouse white adipose tissues.

Many critical steps in the protocol are related to the tissue processing prior to the clearing itself. First of all, and as for classical histological analyses, the fixation step is fundamental. In the case of weak fixation, structural features are not preserved, whereas extended fixation leads to crosslinking and blocking of specific antigen sites and consequent non-homogeneous immuno-staining. Another sensitive step is the staining of the tissue, which presents several critical points that we will describe below. First, the antibodies need to be validated by testing them on cryostat sections to confirm that they are functional and to determine their optimal concentration. The final concentration that is used for the clearing procedure is ten times the optimal concentration for the cryostat sections. Second, the time of incubation is also critical because of the size of the tissue. An incubation time that is too short will produce weak or nonhomogeneous immunostaining (e.g., correct staining on the periphery of the tissue but no internal staining). Similarly, washing steps that are too short will prevent non-specifically-bound antibodies from being removed from the tissue leading to non-specific staining. To our knowledge, extended incubation of the tissue with antibodies does not impair the staining. Therefore, we strongly recommend long incubation and washing periods. Third, the choice of the fluorochrome must be adapted to the signal of interest. In tissue samples in general, and especially white adipose tissue, non-negligible auto-fluorescence is detectable at 488 nm and 555 nm. For highly expressed proteins this is not an issue and the staining will work well in these channels. However, for proteins with low expression we strongly recommend the use of far-red fluorochromes. For the clearing, there are no critical steps since the methodology is very simple. Keep in mind that methyl salicylate is not compatible with plastic materials, therefore we recommend glass bottles for the clearing steps and a metallic chamber equipped with a glass bottom to perform the imaging. Alternatives for this mounting procedure exist using i) a normal glass slide, dental resin and glass coverslip (to generate a small glass chamber), or ii) a glass petri-dish (for upright microscope setups).

Figure 2 .

 2 Figure 2. Photographs of white adipose tissue. (A) Photograph of human white abdominopelvic adipose tissue before and after the clearing procedure. (B) Photograph of mouse epididymal white adipose tissue before and after the clearing procedure.

Figure 3 .

 3 Figure 3. Improved imaging depth with adipose tissue clearing. (A) Mounting of the metallic imaging chamber equipped with glass bottom. (B) Z-series images of mouse epididymal white adipose tissue stained with phalloidin-Alexa488 (green) before and after clearing. (C) XZ projections corresponding to the white dotted lines across the panel B z-stack. (D) Z-projection of 0.4 cm z-stack of mouse subcutaneous white adipose tissue stained for actin using Phalloidin-Alexa488 and acquired at low magnification with a 4x objective. The small images on the right were acquired at the selected tissue position using a 20x objective. (E) Side view of the 3D volume rendering of image z-stack from cleared mouse white adipose tissue stained with Phalloidin-Alexa488 acquired similarly to the 20x images from (D). The 3D imaging depth achieved by our high magnification procedure is demonstrated here and underlined by a zdepth color coding (dark blue for z=0 mm to dark red for z=2 mm).

Figure 4 .

 4 Figure 4. Mouse and human white adipose tissues stained for several markers. (A)Single plane image of mouse epididymal white adipose tissue stained for lipid droplets using anti-perilipin1 antibody and anti-mouse-alexa647-conjugated antibody. (B) Mosaic single plane image of mouse epididymal white adipose tissue stained for nuclei and adipocytes using DAPI, anti-Glut4 antibody and anti-mouse-alexa647-conjugated antibody. (C) Z-projection of 600 µm z-stack of human white abdominopelvic adipose tissue stained for vessels using CD31-alexa647conjugated antibody. (D) Z-projection of 600 µm z-stack of mouse subcutaneous white adipose tissue stained for vessels and macrophages using lectin-DyLight649 intravenous injections and anti-CD301-alexa555 antibody. The lower right inset is a zoomed image of the white dotted box. (E) Single plane image of mouse epididymal white adipose tissue stained for actin and T-cells using phalloidin-Alexa488 and anti-TCRβ-Pacific blue-conjugated. The lower right inset is a zoomed image of the white dotted box. (F) Z-projection of 50 µm z-stack of mouse brown adipose tissue stained for nuclei and neuronal network using DAPI, anti-TH antibody and antirabbit-alexa647-conjugated antibody. (G) Z-projection of 600 µm z-stack of mouse subcutaneous white adipose tissue stained for nuclei and neuronal network using DAPI, anti-TH antibody and anti-rabbit-alexa647-conjugated antibody.

  Figure 1

Table of Materials) if the tissue is too large for a 1.5 mL plastic microtube (see Table of Materials). 2. Permeabilization and saturation of mouse and human white adipose tissue

	controlled orbital shaker for 2 h at 100 rpm.
	2.4. Immerse the tissue in a 15 mL plastic tube containing 10 mL of PBS supplemented with
	0.2% Triton X-100 and 20% DMSO (see Table of Materials) and shake the tube at 37 °C in a
	temperature-controlled orbital shaker at 100 rpm overnight.
	2.5. Immerse the tissue in a 15 mL plastic tube containing 10 mL of PBS supplemented with
	0.1% Tween-20 (see Table of Materials), 0.1% Triton X-100, 0.1% deoxycholate (see Table of
	Materials) and 20% DMSO and shake the tube at 37 °C in a temperature-controlled orbital
	shaker at 100 rpm for at least 24 h.
	2.6. Rinse the tissue in a 15 mL plastic tube containing 10 mL of PBS supplemented with 0.2%
	Triton X-100 and shake the tube at room temperature in an orbital shaker at 100 rpm for 1 h.
	2.7. Immerse the tissue in a 15 mL plastic tube containing 10 mL of PBS supplemented with
	0.2% Triton X-100, 10% DMSO and 3% BSA (see Table of Materials) and shake the tube at 37 °C
	in a temperature-controlled orbital shaker at 100 rpm for 12 h, to saturate any sites that could
	non-specifically bind antibodies.
	NOTE: In step 2.7, the BSA can be substituted by blood serum of the species of the secondary
	antibody used.
	2.1. Rinse the fixed white adipose tissue in 10 mL of PBS for 5 min at room temperature to
	remove all traces of PFA.
	2.2. Immerse the tissue in a 15 mL plastic tube containing 10 mL of PBS supplemented with
	0.3% glycine (see Table of Materials) and shake the tube at room temperature in an orbital
	shaker for 1 h at 100 revolutions per minute (rpm) to quench the remaining free aldehyde
	groups.
	2.3. Immerse the tissue in a 15 mL plastic tube containing 10 mL of PBS supplemented with
	0.2% Triton X-100 (see Table of Materials) and shake the tube at 37 °C in a temperature-

ACKNOWLEDGMENTS:

This work was supported by INSERM, Université Côte d'Azur, and by grants from the French National Research Agency (ANR) through the Investments for the Future Labex SIGNALIFE (ANR-11-LABX-0028-01), the program UCA JEDI (ANR-15-IDEX-01) via Academy 2 "Systèmes Complexes" and Academy 4 "Complexité et diversité du vivant", and the Young Investigator Program to J.G. (ANR18-CE14-0035-01-GILLERON). We also thank the Imaging Core Facility of C3M funded by the Conseil Départemental des Alpes-Maritimes and the Région PACA, and which is also supported by the IBISA Microscopy and Imaging Platform Côte d'Azur (MICA). We thank Marion Dussot for technical help in tissue preparation. We thank Abby Cuttriss, UCA International Scientific Visibility, for proof reading of the manuscript.

DISCLOSURES:

The authors have no conflicts to disclose.

REFERENCES:

1 Pellegrinelli, V., Carobbio, S., Vidal-Puig, A. Adipose tissue plasticity: how fat depots respond differently to pathophysiological cues. Diabetologia. 59 (6)


Liyun Cao

Zhichun Yang
email: yangzc@nwpu.edu.cn

Yanlong Xu

Shi-Wang Fan

Yifan Zhu

Zhaolin Chen

Brice Vincent

Badreddine Assouar
email: badreddine.assouar@univ-lorraine.fr

Disordered Elastic Metasurfaces

Keywords: Elastic metasurfaces, metamaterials, disorder, wave manipulation

The discovery of the disorder effect in traditional metamaterials has opened the possibilities in the search for the disordered metasurfaces. Photonic, dielectric and elastic metamaterials exhibiting the added value of the disorder effect on wave propagation physics, have been reported. Despite this extensive attention and progress in disordered metamaterials, the elastic metasurfaces however, involving disorder have not yet been reported. Here, we introduce the concept of disordered elastic metasurface composed of identical pillared resonators with a random arrangement in the subwavelength range. Based on theoretical formalism and direct acoustic measurement, we observe anomalous deflection and focusing effects of flexural waves in a plate, and elucidate in details the related physics. This research extends the disorder effect to metasurfaces and may lead to innovative acoustoelastic devices.

I. INTRODUCTION

Since the emergence of metamaterials, the study of the disorder has become an important research direction [START_REF] Baboux | Bosonic Condensation and Disorder-Induced Localization in a Flat Band[END_REF][START_REF] Conley | Light Transport and Localization in Two-Dimensional Correlated Disorder[END_REF][START_REF] Nau | Correlation Effects in Disordered Metallic Photonic Crystal Slabs[END_REF]. A deep understanding of disorder effect can help to improve metamaterials properties [START_REF] Hughes | Extrinsic Optical Scattering Loss in Photonic CrystalWaveguides: Role of Fabrication Disorder and Photon Group Velocity[END_REF][START_REF] Celli | Bandgap widening by disorder in rainbow metamaterials[END_REF] and enrich their physics [START_REF] Asatryan | Suppression of Anderson Localization in Disordered Meta-materials[END_REF][START_REF] Bellani | Experimental Evidence of Delocalized States in Random Dimer Superlattices[END_REF][START_REF] Liu | Disorder-Induced Topological State Transition in Photonic Metamaterials[END_REF][START_REF] Rahimzadegan | Disorder-Induced Phase Transitions in the Transmission of Dielectric Metasurfaces[END_REF]. For instance, the disorder can widen bandgap in elastic metamaterials [START_REF] Celli | Bandgap widening by disorder in rainbow metamaterials[END_REF] and induce topological state transition in photonic metamaterials [START_REF] Liu | Disorder-Induced Topological State Transition in Photonic Metamaterials[END_REF] or phase transitions in dielectric metamaterials [START_REF] Rahimzadegan | Disorder-Induced Phase Transitions in the Transmission of Dielectric Metasurfaces[END_REF]. In contrast, independent of the spatial organization of the material, the invariant band gaps of metamaterials can also be observed by taking advantage of strong localization effects [START_REF] Rupin | Experimental Demonstration of Ordered and Disordered Multiresonant Metamaterials for Lamb Waves[END_REF]. As a consequence, a detailed analysis of wave properties in disordered metamaterial is of interest, not only for physics, but also for pragmatic applications. Recently, as one kind of ultrathin metamaterials [START_REF] Catrysse | Routing of Deep -Subwavelength Optical Beams and Images without Reflection and Diffraction Using Infinitely Anisotropic Metamaterials[END_REF][START_REF] Zhu | Negative refraction of elastic waves at the deep-subwavelength scale in a single-phase metamaterial[END_REF][START_REF] Kweun | Transmodal Fabry-Perot Resonance: Theory and Realization with Elastic Metamaterials[END_REF], optic metasurfaces [START_REF] Yu | Light Propagation with Phase Discontinuities: Generalized Laws of Reflection and Refraction[END_REF][START_REF] Bok | Metasurface for Water-to-Air Sound Transmission[END_REF] and acoustic metasurfaces [START_REF] Assouar | Acoustic metasurfaces[END_REF][START_REF] Li | Tunable asymmetric transmission via Lossy acoustic metasurfaces[END_REF][START_REF] Zhu | Fine manipulation of sound via lossy metamaterials with independent and arbitrary reflection amplitude and phase[END_REF][START_REF] Fan | Tunable Broadband Reflective Acoustic Metasurface[END_REF] have been proposed to realize extraordinary wave manipulation. For elastic waves, due to complex wave coupling, the elastic metasurface design [START_REF] Zhu | Anomalous Refraction of Acoustic Guided Waves in Solids with Geometrically Tapered Metasurfaces[END_REF][START_REF] Liu | Source Illusion Devices for Flexural Lamb Waves Using Elastic Metasurfaces[END_REF] is of paramount significance to the current research efforts in a large variety of applications ranging from ultrasonic detection to vibration control. Very recently, disorder had been introduced to the optic metasurface [START_REF] Jang | Wavefront shaping with disorder-engineered metasurfaces[END_REF], which exhibits additional unique features for wavefront shaping such as a tailorable angular scattering profile. However, the acoustic/elastic metasurfaces involving disorder have not yet been reported.

In this research, we introduce the concept of disordered elastic metasurface (DEM) composed of identical pillared resonators with a random arrangement in the subwavelength range. We analytically and numerically reveal the decoupled region of the pillared resonators. The phase shifts and transmission amplitudes of the subunits are only linearly related to the number of resonators and independent of the spatial arrangement of resonators. The disorder releases the degree of freedom of pillar position in the subunit, which can greatly enhance the flexibility for the reconfigurable subunit. Both numerical simulations and experiments are carried out to demonstrate the extraordinary functionalities of DEM, including anomalous deflection and focusing.

The physics of the disorder in the proposed metasurface can be regarded as a universal method for transmissive metasurfaces and metastructure designs. The DEM may also have a direct pragmatic variety of applications in the field of controllable wave manipulation and related acoustoelastic devices and wave engineering.

II. MODEL AND DESIGN

Fig. 1(a) shows the schematic of DEM with predesigned geometries comprising an array of subunits, which is capable of manipulating flexural wave in a plate. The three-dimensional (3D) subunit is constructed based on pillared resonators in Fig. 1(b).

The corresponding two-dimensional (2D) model of the subunit is depicted in Fig. 1(c).

All the attached resonators have the same geometrical parameter: thickness d2 and height h are 3 mm and 30 mm, respectively. The resonators are attached on the plate with a disordered arrangement, i.e., the spacing sj between the jth resonator and (j+1)th resonator is arbitrary. Considering both flexural and longitudinal waves propagate in the hosting plate and resonators, the coefficient vectors of the wave fields at the point P of the incident region and the point Q of the transmitted region (marked in Fig. 1(c where n is the number of the pillared resonator, N1 is the transfer matrix for incident wave propagating from the left to the right of the pillared resonator. N (1) 2 , N (2) 2 and are the transfer matrixes for the flexural and longitudinal waves propagating in the plate regions with pillar spacing s 1 , s 2 and sn-1, respectively. and are the transfer matrixes for waves propagating in Region first and Region last, respectively. For the subunit with different number of resonators, we randomly adjust the pillar spacing sj between the jth resonator and (j+1)th resonator, which meets the relation of n‧d2< to make the total multiresonator thickness not exceed one wavelength λ, i.e., sub-wavelength. In this way, for the different number of attached resonators, 50 disordered subunits are randomly built without loss of generality. To calculate the phase shift and transmission amplitude analytically, we assume that the pillared resonators and plate are made of aluminum alloy. The excitation frequency is selected as 1 kHz and the corresponding flexural wave wavelength for the plate with a thickness of 3 mm is 57.1 times the pillar thickness (d2=3 mm). The transmission amplitude |tb| at the point Q is calculated analytically from Eq. (A6), and the corresponding phase shift is extracted from the complex transmission coefficient tb (see the details in Appendix A).

In Figs. 1(d) and 1(e), the analytical results of the phase shift and transmission amplitude of the flexural waves across the disordered subunits are plotted as blue circles.

(last ) ( 1) (2) (1) (first) out 2 1 2 1 2 1 2 1 2 in n N N N N N N N N N - = k k ! ! ( -1) 2 n N ! (first) 2 N (last ) 2 N -1 1 n j j s l = < å
Several interesting phenomena can be observed as follows. First, randomly disrupting the arrangement of resonators almost does not affect the phase shift and transmission amplitude. In other words, the phase shift and transmission amplitude are only related to the number of resonators in the subunits, but independent of the arrangement. It indicates that the designed metasurface has stronger robustness than conventional designs. Second, the phase shifts by varying the number of resonators can span over a full phase range of 2π with high transmission amplitude (bigger than 0.9). This ensures the high-efficiency of the metasurface for fine field manipulations [START_REF] Xu | Deflecting incident flexural waves by nonresonant single-phase meta-slab with subunits of graded thicknesses[END_REF][START_REF] Zhang | Metasurface constituted by thin composite beams to steer flexural waves in thin plates[END_REF]. Third, the phase shift increases linearly with the number of resonators. This allows us to simply predict the phase shift of these subunits with the different number of resonators.

Furthermore, the ordered arrangements of resonators are also investigated for comparison, with all the pillars spacing keep uniformly 5 mm without loss of generality.

The transmission coefficient tb can be calculated analytically from the modified transfer equation (see the details in Appendix A), where is the transfer matrix for the wave propagating in the plate region of the identical pillars spacing. The phase shift and the transmission amplitude by varying the number of resonators are solved analytically and numerically (plotted as black solid lines and red circles ) in Figs. 1(d) and 1(e). The simulated results are in good agreement with the analytical ones, which confirms the accuracy of the analytical method. Meanwhile, the results for the ordered resonators are very close to the disordered ones, showing the introduction of the disorder will not affect the phase shift and transmission amplitude.

( )

-1 (last ) (first) out 2 1 2 1 2 in n N N N N N * = k k 2 N *

III. THE PHYSICS OF DECOUPLED DISORDERED SUBUNITS

In order to explain these phenomena, we have investigated the model for two adjacent resonators in the subunit, as marked with the red dotted box in Fig. 1(c).

Specifically, we have analyzed the influence of the excitation frequency f, the typical structural parameter h (pillar height), and the pillars spacing s1 on the phase ϕ1 and transmission amplitude |t1| of the point Q in Fig. 1(c). The remaining parameters are the same as those discussed above. The phase ϕ1 and transmission amplitude |t1| can be expressed as ϕ1=Gϕ1(f,h,s1) and |t1|= G|t1|(f,h,s1), respectively. From the above analysis of the disorder unit, it can be seen that under certain conditions, the phase and transmission amplitude are independent of the pillars spacing. In order to give a quantitative evaluation of the dependence of phase and transmission amplitude on the pillars spacing, we define the coupling strengths as:

(A6) (A6)
We further obtain the coupling coefficients and by integrating the coupling strengths for all different pillars spacing (0.003λ, λ) and To further reveal the underlying physics, the reflection coefficient of a single resonator by varying the excitation frequency f and the resonator height h is analytically obtained (see the details in Appendix A), as shown in Fig. 2(c). We label the low reflection boundary line and observe that the decoupled boundary lines for the phase and transmission amplitude are consistent with the low reflection boundary line. This

( ) ( ) 1 1 1 1 , 1 1 , , , , s G f h s f h s s f f ¶ L = ¶ ( ) ( ) 1 1 1 1 1 , 1 , , , , t t s G f h s f h s s ¶ L = ¶ ( ) 1 1 , , s f h f L ( ) 1 1 , , t s f h L 1 1 ,s f L 1 1 , t s L 1 1 , =0 s f L 1 1 , =0 t s L decoupled
1 1 ,s f L 1 1 , t s L
indicates that the decoupling of the phase and transmission amplitude with respect to the pillars spacing depends on the low reflection intensity of the incident wave pass through the resonator. Therefore, the low reflection coupling between resonators induced from their low reflection intensity will lead to decoupling, releasing the degree of freedom of pillar position in the subunit. It should be pointed out that Fabry-Perot resonance, where the phases equal to an integer multiple of π, also has a low reflection intensity. However, it does not lead to the above decoupling. The proposed physics can broaden the understanding of metamaterials and be applied to other transmissive metasurfaces for optics, acoustics, and elastics.

For intuitively displaying the above analysis conclusions, predictive decoupled

IV. DEM MANIPULATING THE FLEXURAL WAVE

The above designed 2D disordered subunits can be extended directly to the corresponding three-dimensional (3D) ones with slots (see in Appendix B), while the

A. Simulations and experiments of the anomalous deflections

The DEM with the functions of anomalous deflection is designed. The

B. Simulation and experiment of the focusing

By employing the disordered subunits, our metasurface design can also be exploited to achieve flexural waves focusing, as shown in Fig. 4 The governing equation of flexural waves in Region j is a fourth-order partial derivative equation. The wave number has four solutions, i.e., two real wave numbers and two complex wave numbers. The real and complex wave numbers represent propagating flexural and evanescent flexural waves, respectively. Therefore, the general solution of displacement for the governing equation is:

(A2)
where the superscript (j) indicates Region j. , , and are complex coefficients. and correspond to the positive-going and negativegoing propagating flexural waves, respectively. and correspond to

( ) ( ) ( ) ( ) 4 2 1 1 1 1 4 2 2 2 1 1 2 2 2 1 , , 0 , , 0 -
w x t w x t E I d x t u x t u x t E x t r r n ¶ ¶ + = ¶ ¶ ¶ ¶ - = ¶ ¶ (1 ) 1 E 1 r 1 I 3 1 1 2 1 12(1 ) I d n - = 1 d 1 n ( ) (
) The relationships between slope and the displacement w, between shear force

1 1 1 1 - - ( ) ( ) ( ) ( ) ( ) , b b b b ik x ik x k x k x j j j j j i t w x t A e B e C e D e e w = + + + ( ) j A ( ) j B ( ) j C ( ) j D 1 ( ) b ik x j A e -
1/4 2 1 1 1 1 1 b d k E I r w ae ö = ç ÷ è ø 2 f w p = ( ) ( ) 1 1 - ( ) ( ) ( ) , l l ik x ik x j j j i t u x t P e Q e e w = + 1/2 2 2 1 1 1 1
(1 )

l k E r w n é ù - = ê ú ë û 1 ( ) l ik x j P e - 1 ( ) l ik x j Q e ( ) ( ) ( ) ( ) 2
+ + + ï ï í ï = + ï î (P ) j 1/4 2 2 2 2 2 2 b d k E I r w ae ö = ç ÷ è ø 1/2 2 2 2 2 2 2
(1 )

l k E r w n é ù - = ê ú ë û 2 E 2 r 2 I 3 2 2 2 2 12(1 ) I d n - = 2 d 2 n j
V and the displacement w, between bending moment M and the displacement w, between axial force F and the displacement u are: (A5) , respectively. Among the right interface of Region j-1, the left interface of Region j, and the lower interface of the resonator j, the following boundary conditions must be satisfied:

1. Displacement continuity in x and z directions:

= = - w V EIw M EIw Ed F u j n ¢ = ì ï ¢¢¢ = - ï ï ¢¢ í ï ï ¢ ï î (1 ) ( ) ( ) ( ) ( ) ( ) ( )
P -1 P -1 j j j j j j u u w w w u ì = =- ï í ï = = î ( ) ( ) ( ) P -1 j j j w w w x x z ¶ ¶ ¶ = = ¶ ¶ ¶ ( ) ( ) ( ) ( ) ( ) ( ) P -1 3 1 1 1 1 2 2 2 2 3 1 1 P -1 3 3 2 2 1 1 1 1 3 3 2 2 0 - - 0 - j j j j j j E d E d u u w E I x x z E d w w u E I E I x x z n n n ì ¶ ¶ ¶ ï - - = ¶ ¶ ¶ ï ï í ï ¶ ¶ ¶ ï - + - = ï ¶ ¶ ¶ î (1 ) (1 )
(1 )

( ) ( ) ( ) P -1 2 2 2 1 1 1 1 2 2 2 2 2 0 j j j w w w E I E I E I x x z ¶ ¶ ¶ - - = ¶ ¶ ¶ , , , , , j j j j j T j j A B C D P Q = k 1 ( ) ( 1) j j L R N - = k k ( 1) j R - k ( ) j L k 1 N ( ) j R k ( ) j L k ( ) ( ) ( ) ( ) 2 2 | = j j j j R L L j s s N N = = k k k 1 1 1
× -× × ae ö ç ÷ ç ÷ ç ÷ = ç ÷ ç ÷ ç ÷ ç ÷ ç ÷ è ø 1(c).
For Region first with length , the vector at the right interface and the vector in the left interface of Region first can be expressed as:

(A14)
The examined point for the transmission phase and amplitude is fixed at the point Q, as shown in Fig. 1(c). For Region last with the length , the vector at the right end of Region last and the vector in the left end of Region last can be expressed as:

(A15)
where, is the modified length of Region last.

According to Eqs. (A12)-(A15), for waves propagating from the left interface of

Region first to the right interface of Region last, the transfer relationship can be expressed as:

(A16)
From this relation, we can solve the transmission amplitude for the flexural wave transmitting through these n resonators. To do this, we make an incident positive-going flexural wave of at the left interface of Region first, so the wave fields at the left interface of Region first and the right interface of Region last can be expressed as follows:

(A17)

where , , , , , and are the amplitude ratios of the reflected propagating flexural wave, transmitted propagating flexural wave, reflected evanescent flexural

f =2 s l (first) R k in k 2 2 (first) (first ) in 2 in = R s N N l = = × × k k k l s out k (last) L k (last) (last) (last) out 2 2 = l L s s L N N * = = × × k k k ( ) 1 2 2 1 0.08 3 0.08 n l l j j s s nd s d l - * = = - » - + å (last ) ( 1) (2) (1)
(first) out 2 1 2 1 2 1 2 1 2 in n N N N N N N N N N - = k k ! 1 i t e w × ( ) ( ) ( ) ( ) ( ) ( ) (first) (first) (last ) (las ) * t * , 1 , , , i t b
b i t l i t b b i t l L L R R w x t r
ì = + + ï ï = ï í = + ï ï = ï î b r b t b r * b t * l r l t
wave, transmitted evanescent flexural wave, reflected longitudinal wave, and transmitted longitudinal wave to the incident positive-going propagating flexural wave, respectively. Therefore, the vectors and can be written as:

(A18)
The transfer equation for the single resonator can be expressed as:

(A19)

The transfer equation for the ordered multiresonator

For the ordered subunits, all the pillar spacings keep uniformly sj = 5 mm. The corresponding transfer matrix for wave propagating from the left interface to the right interface of the uniform Region j is written as For waves propagating from the left interface of Region first to the right interface of Region last in the ordered subunit, the transfer equation can be expressed as:

(A20)

Analytical solution for transmission phase shift and amplitude

According to Eqs. (A16), (A19) and (A20), , , , , , and can be solved with the aid of MATLAB code. is the amplitude ratio between the transmitted propagating flexural wave and the incident propagating flexural wave. The phase can be obtained from the complex amplitude ratio . The phase shift is obtained by subtracting the phase at the point Q for the plate with the resonators from that for the plate without resonators. It can be calculated by the following formula:

in k out k in out T T = 1, ,0, ,0, = , 0, ,0, , 0 b b l b b l r r r t t t * * ì é ù ï ë û í é ù ï ë û î k k (last ) (first) out 2 1 2 in N N N = k k 2 2 0.005 = | s N N * = ( ) -1 (last ) (first) out 2 1 2 1 2 in n N N N N N * = k k b r b r * l r b t b t * l t b t b t (A21)
This formula will make the value of the phase shift from 0 to .

APPENDIX B: DESIGN FROM THE 2D SUBUNIT TO THE 3D ONE

The above designed subunit is embodied in its 2D model. To manipulate flexural waves, the subunit should be designed as the corresponding 3D one in Fig. 5(a), which can be formed by simply stretching the 2D model in Fig. 1(c). The stretched width is denoted as . Because the metasurface is formed by periodically integrating different 3D subunits, the 3D subunit should be analyzed by applying the periodic boundary conditions on the two long boundaries in the 3D strip-like model in Fig. 5(a). Therefore, the formed 3D subunit is equivalent to the corresponding 2D one in Fig. 1(c). However, slots should be cut from the 3D subunit to separate adjacent different 3D subunits in the integrated metasurface, as shown in Fig. 5(b). The aim is to avoid the Fano interference [START_REF] Jin | Tunable Fano resonances of Lamb modes in a pillared metasurface[END_REF] between the adjacent different 3D subunits to independently design the subunits. The filling ratio of the slots in one subunit is denoted as , which is defined as dividing the slot width by the stretched width w1 in a subunit. In our designs, the widths of the processed pillared resonator w0 and the slot are 20.4 mm and 1 mm, respectively, i.e., w1=21.4 and =0.047. The slots are very small relative to the stretched width of the subunit. The little slot will lead to a very small flexural wave scattering, which will not affect the transmission amplitude and phase shift [START_REF] Cao | Deflecting flexural wave with high transmission by using pillared elastic metasurface[END_REF]. Therefore, the 3D subunit structure with the slot still can be equivalent to the corresponding 2D one.

APPENDIX C: GEOMETRY OF FABRICATED DEM AND TEST SET-UP

To accurately fabricate the metasurface, the Computer Numerical Control (CNC) milling machine with a manufacturing precision of 0.01 mm was adopted. The DEMs in the experimental measurements are composed of the identical processed resonator.

For the processed plate, slots should be cut to separate adjacent different 3D subunits in a metasurface, as shown in Figs. 678. The adjacent slots make up many strip-like models. Different numbers of the processed resonators are attached to the strip-like model to make up the different subunits, which consist of the metasurface. The numbers in yellow fonts represent the different number of attached resonators in the strip-like model, as shown in Figs. 678. These resonators can be conveniently disassembled and reattached, because the glue layer is very thin, in addition, there is no need for the ordered arrangement of these resonators. This process enables to design the disordered metasurface to achieve tunable multifunction by using one type of pillared resonator and the same hosting plate. For example, the fabricated metasurface in Fig. 6 An array of PZT-5A patches (20 mm × 20 mm × 0.3 mm) were bonded on the surface of the plate as actuators. The spacing of adjacent piezoelectric patches is the same as the piezoelectric patch width according to the phased array theory [START_REF] Zhu | Negative refraction of elastic waves at the deep-subwavelength scale in a single-phase metamaterial[END_REF]. The widths of the DEMs only need to be larger than the width of the excited flexural beams.

For simplicity, in the deflection experiments shown Fig. 6 and Fig. 7, we adopted 6 piezoelectric patches (the total width is 220 mm) to excite flexural waves and 16 subunits to build the DEMs (the total width is 342.4 m). In the focusing experiments shown Fig. 8, we adopted 12 piezoelectric patches (the total width is 460 mm) to excite flexural waves and 24 subunits to build the DEM (the total width is 513.6 mm). With the controlling of the waveform editor (based on Labview Software), a 5-cycle tone burst was generated by the signal generator (Agilent 33220A) and then amplified by a power amplifier (HVPA05). The wave fields in the incident and transmitted areas were measured by the measurement mode "time" of the PSV-400 scanning laser Doppler l vibrometer [START_REF] Jin | Tunable Fano resonances of Lamb modes in a pillared metasurface[END_REF].

APPENDIX D: THE DESIGN OF THE RECONFIGURABLE SUBUNIT

The schematic diagram of the reconfigurable subunit is shown in Fig. 9. It can be seen from Fig. 9 that the reconfigurable subunit has some improvements compared with the disordered subunit in Fig. 1 The working principle is very simple. Magnetic fields in different directions are formed under the three energized electromagnets. For example, magnetic fields are N, N and S poles for lower surfaces of # 1, # 2 and # 3 electromagnets, respectively. We can move the pillar through the plastic guide rod to select the repulsion or attraction applied by the magnetic field. Under the # 1 and # 2 electromagnet, the pillar is repulsed and attached to the host plate. Under the # 3 electromagnet, the pillar is separated from the host plate due to the magnetic attraction. Therefore, the number of pillars attached on the host plate can be reconstructed by moving the pillars to different magnetic fields.

In this way, the disorder has great advantages over the order. The reason is the disorder releases the degree of freedom of pillar position, which greatly enhance reconfigurability.

FIG. 10. The differential element of the pillar resonator.

We also have conducted the feasibility analysis for the designed structure. First, in the experiment, we can apply a small amount of oil on the lower surface of the pillar to fill the gap between the pillar and the host plate. This super subwavelength oil layer not only ensures a good connection between the pillar and the host plate, but also ensures that it does not affect the waves propagation at the joint interface. Then, the impedance of plastic and the thin iron layer is greatly mismatched, which will ensure that the lightweight plastic guide rod will not affect the wave propagation in the pillar. Finally, applying a repulsive force to the pillar is equivalent to giving the pillar an axial prestress T. The differential element of the pillar is taken for analysis, as shown in Fig. 10.

The equation of motion in the x direction, following Newton's second law ( ),

F ma S = It should be point that for the x-direction boundary condition at the junction of the pillar and the host plate, the non-glued pillar in the above reconfigurable subunit may weaken it compared with the glued pillar. We can reduce the height of the pillared resonator or decrease the rigidity of the pillared resonator, so that the vibration of the pillared resonator can couple into the host plate with the dominant effect coming from the compressional modes of the resonator [START_REF] Rupin | Experimental Demonstration of Ordered and Disordered Multiresonant Metamaterials for Lamb Waves[END_REF][START_REF] Colquitt | Seismic metasurfaces: Subwavelength resonators and Rayleigh wave interaction[END_REF], in other word, the y-direction

FIG. 1 .

 1 FIG. 1. Disordered elastic metasurface (DEM). (a) Schematic diagram of DEM. (b) Three-

  conducting normalization with respect to their maxima. The effects of changing f and h on the coupling coefficients and are shown analytically in Figs. 2(a) and 2(b), respectively. It clearly shows the existence of completely decoupled regions characterized by and , respectively, which means that the phase shift and the transmission amplitude are independent of the pillars spacing. The marked

  FIG. 2. The physics of decoupled disordered resonators. (a), (b) Effect of the changes of f and h on

  Figs. 2(d) and 2(e). One can observe that the phase and transmission amplitude at the

  phase shift and transmission amplitude almost remain unchanged. The elastic metasurfaces are designed by the 3D disordered subunits. The different number of the pillared resonators are attached on the hosting plate with slots (see in Appendix C) to achieve different function, including anomalous deflection and focusing. Importantly, these resonators can be conveniently disassembled and reattached, since the disordered and ordered resonators have the same effect on wave manipulation. It is convenient to continuously change the phase shift for every subunit by changing their number of resonators. It should be pointed out that the reconfigurable subunit can be designed (the detailed design can be found in Appendix D) based on the disordered subunit. The disorder can greatly enhance the flexibility for the reconfigurable subunit, because it releases the degree of freedom of pillar position in the subunit.

  FIG. 3. Deflecting flexural wave with DEM. (a), (b) The wave fields for vertically and obliquely

  (e). The DEM can be readily constructed by selecting units that match the requested hyperbolic phase profile , where F=1.5λ is the desired focal distance (detailed geometric parameters of the DEM are provided in the Appendixes C). The transmitted intensity and displacement fields with focus at 1.5λ are simulated and shown in Figs. 4(a) and 4(c), respectively. One clearly observes a focal spot after the incident waves pass through the DEM. The full field experimental measurement is shown in Fig. 4(d). The detailed propagation phenomenon is presented in Movie S3 of the Supplemental Material [25]. For a more quantitative comparison, the normalized intensity fields of the transmitted wave at the focal spot along the y-direction from the simulations and experiments are shown in Fig. 4(b). The sharp peak values are observed, which are about more than 3 times the intensity of the incident wave. The experimental measurement is in good agreement with the simulated result and demonstrates a fine

FIG. 4 .

 4 FIG. 4. Focusing flexural wave with DEM. (a), (c) The simulated results of the transmitted intensity

  positive-going and negative-going evanescent flexural waves, respectively. The flexural wavenumber is , in which is the circular frequency. The general solution of the displacement for the second-order partial derivative governing equation of longitudinal waves in the Region j is:(A3)where the longitudinal wavenumber is , and correspond to the positive-going and negative-going propagating longitudinal waves, respectively.The governing equations for the flexural wave and the longitudinal wave in the pillared resonator is similar to that in the host plate. The general solutions of the displacements for the flexural wave and longitudinal wave equations in the pillared resonator j are (A4), respectively. The superscript indicates the pillared resonator j. 's modulus, density, and area moment of inertia of the pillared resonator, respectively. , in which and are the thickness of the resonator and Poisson's ratio of the material, respectively.

3 . 5 .

 35 Force balance in x and z directions: Forces and moment are 0 at the free boundary of the resonator j:

2

 2

FIG. 5 .

 5 FIG. 5. (a) Model diagram of the 3D subunit without slots. (b) Model diagram of the 3D subunit

FIG. 8 .

 8 Fig. 7 and Fig. 8 can achieve negative deflection (the experimental result shows in Figs.

FIG. 9 .

 9 FIG. 9. Schematic diagram of the reconfigurable subunit

  w is the displacement of the particle of the pillar in the x direction. Substituting the relationships and into Eq. (D1), we obtain (D2) From Eq. (D2), the flexural wavenumber can be obtained (D3) where, 2EI= 353.5. When T/(2EI) < 0.06 << 1, we can get . It indicates that the pre-pressure [T/(2EI) < 0.06] does not affect the propagation of flexural waves in the pillar. Through the same analysis, we find the pre-stress [T/(2EI) < 0.06] hardly affect the propagation of the longitudinal wave in the pillar and the propagations of the flexural wave and the longitudinal wave in the plate. Therefore, the pre-stress [T/(2EI) < 0.06] does not affect the propagations of flexural waves and longitudinal wave, at the same time, it can make the pillar stably attached to the host plate.

FIG. 11

 11 FIG. 11. (a) and (b) show wave fields for vertical and oblique incident flexural wave Gaussian

FIG. 12 .

 12 FIG. 12. (a) and (c) numerically show the transmitted intensity and displacement fields with focus

(A10)

The transfer equation for the disordered multiresonator

The complex coefficients , , , , , for Region (j)

in the plate can be organized as a vector The relationship between the vector at the right interface of Region j and the vector at the left interface of Region j can be expressed as follows:

(A13) where is the transfer matrix for waves propagating between the two ends of Region j.

The exciting point is assumed as the point P on the plate surface as shown in Fig.

( ) ( ) ( )

boundary condition will be dominant at the junction of the resonator and the host plate.

In this way, it can eliminate the impact of weakened x-direction boundary conditions.

APPENDIX E: METHOD FOR NUMERICAL SIMULATIONS

All the pillared resonators and plate are made of aluminum alloy, with the elastic