Julien Claron

Vincent Hingot

Isabelle Rivals

Line Rahal

Olivier Couture

Thomas Deffieux

Mickael Tanter

Sophie Pezet
email: sophie.pezet@espci.fr

Large-scale functional ultrasound imaging of the spinal cord reveals in-depth spatiotemporal responses of spinal nociceptive circuits in both normal and inflammatory states

Keywords: Pain, fUS, Neurovascular, Plasticity, NMDA

teaching and research institutions in France or abroad, or from public or private research centers.

Large-scale functional ultrasound imaging of the spinal cord reveals in-depth spatiotemporal responses of spinal nociceptive circuits in both normal and inflammatory states

Introduction

The spinal dorsal horn is the first site of integration of multiple sensory modalities and their modulation. Most of our knowledge on the physiology and physiopathology of the spinal cord comes from electrophysiological and histological studies. Extracellular recordings produced invaluable information regarding the neuronal electrophysiological properties in health and disease, through peripheral stimulations of the neurons' receptive field. They evidenced the increased neuronal firing, as well as the involvement of N-methyl-D-aspartate (NMDA) receptordependent mechanisms underlying central sensitization associated with persistent pain, [START_REF] Latremoliere | Central sensitization: a generator of pain hypersensitivity by central neural plasticity[END_REF] shedding light on the mechanisms of action of ketamine in the clinical manifestations of central sensitization. [START_REF] Davies | Evidence for involvement ofN-methylaspartate receptors in wind-up of class 2 neurones in the dorsal horn of the rat[END_REF][START_REF] Dickenson | Differential effects of excitatory amino acid antagonists on dorsal horn nociceptive neurones in the rat[END_REF] Neuroanatomical studies, however, helped understand the role of intracellular signaling pathways in central sensitization and the nature and role of subtypes of excitatory and inhibitory neurons in the various laminae of the dorsal horn. [START_REF] Todd | Neuronal circuitry for pain processing in the dorsal horn[END_REF] All these techniques have their strengths and weaknesses. Ultimately, a technique that would have excellent temporal and spatial resolutions, and high sensitivity, enabling measurements at the cellular level and in a large field of view, would be perfect. Electrophysiological and optical imaging [START_REF] He | Hemodynamic and light-scattering changes of rat spinal cord and primary somatosensory cortex in response to innocuous and noxious stimuli[END_REF] studies only allow recordings in a very narrow space. Neuroimaging techniques, such as functional MRI, by contrast, image with a large field of view, but due to the small size of the spinal cord in rodents, it requires high magnetic fields and long exposure times to visualize the spinal cord, with sometimes a poor spatiotemporal resolution. For all these reasons, and due to the breathing-related motion artefacts, [START_REF] Cohen-Adad | Functional magnetic resonance imaging of the spinal cord: current status and future developments[END_REF] there are significantly fewer neuroimaging studies performed in humans' [START_REF] Brooks | Stimulus site and modality dependence of functional activity within the human spinal cord[END_REF][START_REF] Cahill | Mapping of neural activity produced by thermal pain in the healthy human spinal cord and brain stem: a functional magnetic resonance imaging study[END_REF][START_REF] Cohen-Adad | BOLD signal responses to controlled hypercapnia in human spinal cord[END_REF][START_REF] Nash | Functional magnetic resonance imaging identifies somatotopic organization of nociception in the human spinal cord[END_REF][START_REF] Sprenger | Spinal cord-midbrain functional connectivity is related to perceived pain intensity: a combined spinocortical fMRI study[END_REF][START_REF] Stroman | Mapping of neuronal function in the healthy and injured human spinal cord with spinal fMRI[END_REF][START_REF] Stroman | Somatotopic arrangement of thermal sensory regions in the healthy human spinal cord determined by means of spinal cord functional MRI[END_REF][START_REF] Stroman | The current state-of-the-art of spinal cord imaging: Methods[END_REF][START_REF] Summers | A quantitative comparison of BOLD fMRI responses to noxious and innocuous stimuli in the human spinal cord[END_REF] and rodents' spinal cord [START_REF] Lawrence | Correlation of functional activation in the rat spinal cord with neuronal activation detected by immunohistochemistry[END_REF][START_REF] Lilja | Blood oxygenation level-dependent visualization of synaptic relay stations of sensory pathways along the neuroaxis in response to graded sensory stimulation of a limb[END_REF][START_REF] Majcher | Simultaneous functional magnetic resonance imaging in the rat spinal cord and brain[END_REF][START_REF] Malisza | Functional magnetic resonance imaging of the spinal cord during sensory stimulation in diabetic rats[END_REF][START_REF] Malisza | Functional MRI of the rat lumbar spinal cord involving painful stimulation and the effect of peripheral joint mobilization[END_REF][START_REF] Pich É | Tight neurovascular coupling in the spinal cord during nociceptive stimulation in intact and spinal rats[END_REF][START_REF] Zhao | BOLD and blood volumeweighted fMRI of rat lumbar spinal cord during non-noxious and noxious electrical hindpaw stimulation[END_REF][START_REF] Zhao | fMRI investigation of the effect of local and systemic lidocaine on noxious electrical stimulation-induced activation in spinal cord[END_REF] than in the brain.

Functional ultrasound (fUS) imaging is a relatively new versatile neuroimaging approach that allows imaging and measurement of cerebral blood volume in human neonates, [START_REF] Demene | Functional ultrasound imaging of brain activity in human newborns[END_REF] in adult intraoperative settings, [START_REF] Imbault | Intraoperative functional ultrasound imaging of human brain activity[END_REF] nonhuman primates, [START_REF] Dizeux | Functional ultrasound imaging of the brain reveals propagation of task-related brain activity in behaving primates[END_REF] anesthetized [START_REF] Demene | 4D microvascular imaging based on ultrafast Doppler tomography[END_REF][START_REF] Mace | Functional ultrasound imaging of the brain[END_REF] and freely moving rodents [START_REF] Bergel | Local hippocampal fast gamma rhythms precede brain-wide hyperemic patterns during spontaneous rodent REM sleep[END_REF][START_REF] Sieu | EEG and functional ultrasound imaging in mobile rats[END_REF][START_REF] Urban | Realtime imaging of brain activity in freely moving rats using functional ultrasound[END_REF] with excellent spatial and temporal resolutions.

Due to the existence of a neurovascular coupling, an increased local neuronal activity induces a hemodynamic response. [START_REF] Iadecola | The neurovascular unit coming of age: a journey through neurovascular coupling in health and disease[END_REF] Functional ultrasound imaging is sensitive enough to measure the cortical hemodynamic changes induced by olfactory, [START_REF] Osmanski | Functional ultrasound imaging reveals different odor-evoked patterns of vascular activity in the main olfactory bulb and the anterior piriform cortex[END_REF] visual, [START_REF] Gesnik | 3D functional ultrasound imaging of the cerebral visual system in rodents[END_REF] and auditory [START_REF] Bimbard | Multi-scale mapping along the auditory hierarchy using high-resolution functional UltraSound in the awake ferret[END_REF] stimuli in anesthetized animals. The use of such a technique in the spinal cord, in the context of a better understanding of the physiology of the spinal cord in sensory integration, would be very valuable. Recent studies have shown as first proofs of concept that fUS imaging is an emerging modality of neuroimaging in the rodent spinal cord, in particular in the context of spinal cord lesions and its restoration. [START_REF] Khaing | Contrast-enhanced ultrasound to visualize hemodynamic changes after rodent spinal cord injury[END_REF][START_REF] Song | Functional ultrasound imaging of spinal cord hemodynamic responses to epidural electrical stimulation: a feasibility study[END_REF][START_REF] Soubeyrand | Very high resolution ultrasound imaging for real-time quantitative visualization of vascular disruption after spinal cord injury[END_REF] Our study aimed at (1) a detailed imaging of the vascular structures of the rat spinal cord by combining both fUS imaging and ultrasound localization microscopy (ULM), which increases even further the spatial resolution (10 mm), [START_REF] Errico | Ultrafast ultrasound localization microscopy for deep super-resolution vascular imaging[END_REF] (2) evaluating the sensitivity of this fUS technique to image the spinal neurovascular response to natural (noxious or innocuous) peripheral stimuli, (3) determining the typical hemodynamic response function (HRF), and (4) imaging their alterations due to the central sensitization induced by short-term inflammatory pain.

Material and method

Animals

All experiments were performed in agreement with the European Community Council Directive of 22 September 2010 (010/63/UE) and the local ethics committee (Comit é d' éthique en mati ère d'exp érimentation animale number 59, "Paris Centre et Sud," project 2018-05). Accordingly, the number of animals in our study was kept to the necessary minimum. Experiments were performed on N 5 51 male Sprague-Dawley rats (Janvier Labs; Le Genest St Isle, France), weighing 325 to 350 g at the beginning of the experiments. Animals (2 per cage) arrived in the laboratory 1 week before the beginning of the experiment and were kept at a constant temperature of 22˚C, with a 12-hour alternating light/ dark cycle. Food and water were available ad libitum.

In toto coloration of spinal cord vascular architecture

After an imaging session, N 5 3 Sprague-Dawley rats were deeply anesthetized with an intraperitoneal (IP) injection of sodium pentobarbital (150 mg•kg [START_REF] Eippert | Direct evidence for spinal cord involvement in placebo analgesia[END_REF] ). Then, a thoracotomy was performed and an incision in the right atrium was made. Animals were perfused with 2 mL of saline solution (NaCl 0.9%), followed immediately by 15 mL of DiI (0.012 mg•mL [START_REF] Eippert | Direct evidence for spinal cord involvement in placebo analgesia[END_REF] in 2% ethanol, 98% PBS [phosphate buffered saline] solution, 0.02 mol•L 21 ), as previously described [START_REF] Li | Direct labeling and visualization of blood vessels with lipophilic carbocyanine dye DiI[END_REF] and then 10 mL of paraformaldehyde 4% at the rate of 7 mL•min [START_REF] Eippert | Direct evidence for spinal cord involvement in placebo analgesia[END_REF] . Spinal cords were extracted and fixed for 2 days in 4% paraformaldehyde at 4˚C before cryoprotection in a 30% sucrose solution for 2 days. Spinal cords were frozen in an OCT matrix (optimal cutting temperature compound) in cooled isopentane (240˚C) on dry ice, cut in 50 mm slices using "freefloating" method in PBS, and positioned on Superfrost slides. Mosaic images were immediately acquired using Axio Imager M1 (Zeiss, Jena, Germany).

Laminectomy and spinal functional ultrasound imaging

Under deep anesthesia (intraperitoneal [IP] bolus of medetomidine [Domitor, 0.4 mg•kg 21 ] and ketamine [Imalg ène, 40 mg•kg 21 ]), a laminectomy was performed between the 12th thoracic vertebra and the second lumbar vertebra. This window allowed the positioning of the whole ultrasonic probe (14 mm) in a sagittal plane.

Once the laminectomy performed, the animal was placed on a "spinal cord" stereotaxic frame, where the lumbar rachis was suspended. Lidocaine (Xylovet, Ceva, France, 2 mg/mL) was topically applied on the muscles that were exposed, as well as the insertion points of the spinal holders. The anesthesia was maintained but reduced, using subcutaneous perfusion at the rate of medetomidine (0.1 mg/kg/h) and ketamine (12.5 mg/kg/h) using a syringe pump. To obtain reproducible results, it was necessary to wait approximately 1h30 after induction before imaging to obtain a stable and light level of anesthesia (respiratory frequency around 120 bpm). During the surgical procedure and the imaging session, the animal's body temperature was kept at 37˚C using a heating blanket and an intrarectal probe (Physitemp, Clifton, NJ), and the heart and respiratory frequencies were monitored (MouseOxPlus, Ugo Basile, Italy). Each imaging session lasted from 4 to 6 hours.

Two milliliters of saline was gently dropped on the spinal cord (the dura mater was kept intact), and the cavity made by the laminectomy was then filled with echographic gel. The ultrasonic probe (f 5 15 MHz, 128 elements, Vermon, Tours, France) was then positioned just above the window using a 3-axis motor system on which the ultrasound probe was fixed. The probe connected to an ultrasonic ultrafast scanner (Verasonics, Kirkland, WA: 128 channels, 62.5 MHz sampling rate) was driven with Neuroscan live acquisition software (Iconeus, Paris, France; and Inserm Accelerator of Technological Research in Biomedical Ultrasound, Paris, France).

Peripheral stimulations

Mechanical stimulations

Twenty-four hours before the imaging session, the threshold of nociceptive mechanical sensitivity was determined using calibrated forceps (Bioseb, France), as previously described, [START_REF] Thibault | Characterisation of sensory abnormalities observed in an animal model of multiple sclerosis: a behavioural and pharmacological study[END_REF][START_REF] Thibault | Molecular mechanisms underlying the enhanced analgesic effect of oxycodone compared to morphine in chemotherapy-induced neuropathic pain[END_REF][START_REF] Thibault | Activation-dependent subcellular distribution patterns of CB1 cannabinoid receptors in the rat forebrain[END_REF] by applying an increasing pressure on the foot pad of the hind paw in awake and habituated animals. This nociceptive threshold (withdrawal threshold .4.5 N) was further used in imaging sessions to apply the necessary force to elicit a nociceptive mechanical stimulation.

Electrical stimulations

A low-frequency generator (A310 Accupulser, World Precision Instruments, Hitchin, United Kingdom) was used with a current generator (Digitimer, DS3; Letchworth Garden City, United Kingdom) to generate a waveform (2 msec width, 40 Hz stimulation for 5 seconds) for subcutaneous hind paw stimulation using acupuncture needles. To stimulate different groups of fibers, we chose different intensities: supramaximal stimulation of C fibers required stimulation at 3 mA for C fibers, 1.5 mA for of Ad fibers, and 0.5 mA for Ab fibers. These parameters were chosen according to previous literature. [START_REF] Chang | A fMRI study of brain activations during non-noxious and noxious electrical stimulation of the sciatic nerve of rats[END_REF][START_REF] Bars | Diffuse noxious inhibitory controls (DNIC). I. Effects on dorsal horn convergent neurones in the rat[END_REF][START_REF] Silva | Early temporal characteristics of cerebral blood flow and deoxyhemoglobin changes during somatosensory stimulation[END_REF] During the sets of experiments, the stimulation patterns were randomized. Also, to avoid any spinal hyperexcitability induced by repetitive stimulation of C fibers, [START_REF] Mcmahon | Central hyperexcitability triggered by noxious inputs[END_REF] only a maximal number of 8 acquisitions at 3 mA were performed on each animal.

For the determination of the spinal HRF, a high-intensity and short-duration [START_REF] Rosa | On the distinguishability of HRF models in fMRI[END_REF] (5 mA, 2 seconds) waveform at 40 Hz was used. The waveform was triggered by the ultrafast imager. We chose the shortest duration inducing a strong and reproducible response.

2.2. Functional ultrasound imaging: sequences of imaging and signal processing

Sequences of imaging

Functional ultrasound imaging was performed using a linear ultrasound probe (128 elements, 15 MHz, 110 mm pitch, 8 mm elevation focus, Vermon) driven by an ultrafast ultrasound scanner (Verasonics, 128 channels, 62.5 MHz sampling rate), which was driven with Neuroscan live acquisition software (Iconeus, Paris, France; and Inserm Accelerator of Technological Research in Biomedical Ultrasound, Paris, France). The fUS imaging sequence operated as follows: the spinal cord was insonified by 10 successive tilted plane waves with angles varying from 210˚to 10˚with a 5.5-kHz pulse repetition frequency. The backscattered echoes were recorded by the transducer array and beamformed to produce a block of 200 consecutive ultrafast images with a framerate of 500 Hz. Each block acquisition was followed by a pause up to 600 ms for data treatment. To filter the spinal blood volume (SBV) and remove the tissue signal, we used a clutter filter based on singular value decomposition (SVD) algorithm applied to 200 successive frames [START_REF] Demen É | Spatiotemporal clutter filtering of ultrafast ultrasound data highly increases Doppler and fUltrasound sensitivity[END_REF] by removing the 60 first singular vectors that correspond mainly to the tissue space. Finally, a Power Doppler image was obtained by integrating the energy of the filtered frames, resulting in a Power Doppler image every 400 ms or 1 second depending on the pause.

Note that no additional motion correction was necessary because the lumbar spinal cord was tightened and suspended using spinal holders and the aforementioned SVD clutter filter can remove residual tissue motion.

Two planes were chosen for imaging

A sagittal plane and a transverse one were chosen. For the placement of the sagittal plane, we made several scans in the transverse direction and we chose the middle of the spinal cord by recognizing specific structures. Once at the midline, the probe was moved 500 mm towards the ipsilateral side of the simulation. Therefore, the functional studies in the sagittal plane were performed at 500 mm laterality in the ipsilateral side. For the top-down positioning, the image was centered on the lumbar enlargement.

For the transverse positioning, a preliminary positive stimulation test was made on the sagittal plane before rotating the probe with the 3 axis motors on the area highlighted by preliminary acquisition.

Experiments of functional activation in the transverse plane were difficult to perform because it is easy to miss the responsive area in the cord (slice thickness is about 400 mm). Also, within this plane, the position of the vertebra and the ribs induce artifactual reverberations, causing a high level of background noise compared to the sagittal plane. This is probably due to a difference in wave scattering in the intervertebral space that is different between the sagittal plane and the transverse plane. These are the reasons why we did not perform many experiments in transverse plane.

Doppler signal analysis

Experimental determination of the spinal hemodynamic response function

Because of the short duration of noxious electrical stimulations used and the lack of previous descriptions of hemodynamic response in the rat's spinal cord, we first sought to determine the rat spinal HRF in our experimental conditions. This was necessary to correlate the variations of Doppler signal with the stimulation pattern accurately, and to produce activation maps. A strong (supramaximal 5 mA) stimulation of C fibers was applied for the shortest duration that produced a reproducible response (2 seconds). To stimulate accurately and reproducibly, the stimulations were triggered by the ultrasound scanner using an Arduino board. Five acquisitions per rat in 2 rats were used to determine the HRF as measured by fUS imaging in the spinal cord of the rat experimentally. All acquisitions were averaged, and the analytical equations of the HRF were determined by fitting an inverse gamma distribution, as previously described by other authors [START_REF] Lindquist | Modeling the hemodynamic response function in fMRI: efficiency, bias and mismodeling[END_REF] using Matlab (lsqcurvefit, Optimization Toolbox).

Doppler signal analysis and activation maps

Doppler data were analyzed using a generalized linear model approach implemented in Matlab. In the design matrix, the stimulation pattern was convoluted with our fUS-determined HRF for shortest electrical stimulation and we obtained Z-score and Pvalue map. The activation maps show the Z-score of all the pixels in the image with P-value , 0.05 (before Bonferroni correction).

We chose our region of interest (ROI) by calculating a 7-by-7 pixelsquare centered on the maximum of the thresholded Z-score map, and the signal was averaged along the 2 spatial dimensions to obtain a single temporal signal. The spatially averaged signal was then expressed as an SBV increase (in percent) by subtracting the SBV baseline (calculated for each acquisition by averaging all the temporal data within the spinal cord where the stimulation pattern was strictly equal to 0) and by dividing by the SBV baseline. Spinal blood volume increases were calculated by isolating each peak and by averaging the temporal signal between t 1 and t 2 , where t 1 and t 2 were determined using the width at half maximum for the considered peak. The baselines between each stimulation peak couples were estimated separately.

Statistical analysis of the evoked spinal hemodynamic responses

The statistical analysis was performed using a linear mixed model (LMM). Such approach is particularly interesting in the case of nonindependent, hierarchical data, which is the case of the results produced (DSBV or Z-scores from individual or series of stimulations, in series of acquisitions, in several animals). In addition, both fixed and random effect factors can be introduced. In our set of experiments, the fixed effect factor was the condition (eg, location of hemodynamic response [dorsal horn, ventral horn]) or the intensity of fiber stimulation, whereas the animals and acquisitions are random effect factors.

Data were first transformed by taking their square root to homogenize the variance across groups, as assessed with Bartlett test. The LMM was fitted with the restricted maximum likelihood method using Matlab's function fitlme. Residual normality was tested with Shapiro-Wilk test. In the case of more than 2 conditions, post hoc tests of significance of the fixed factor effect between couples of conditions were performed with a likelihood ratio test using Matlab's function coefTest, using Bonferroni correction for multiple testing.

Results of the LMM are shown in detail in Supplementary table 1, whereas results of the post hoc tests are presented in the Supplementary table 2 (available at http://links.lww.com/PAIN/ B176).

(1) Condition: evaluation of the significance of the changes of SBV due to the condition, such as comparison in the dorsal horn vs ventral horn, inflamed vs naive, etc. (2) Animal: evaluation of the significance of the changes due to differences in animals.

(3) Acquisition: evaluation of the significance of the changes due to differences in acquisitions, ie, 4 trains of stimulations, each one given rise to a value of evoked SBV (indicated as "Total number of stimulations" in the table). ( 4) Homoscedasticity: homogeneity of variance is assumed by the LMM. P values .0.05 suggest that the variance was indeed homogenous after the square root transformation of the data. ( 5) Residual normality: error normality is assumed by the LMM.

Therefore, analysis of the residual normality is a good indicator of how well this hypothesis is fulfilled (it is the case with a few exceptions).

Analysis of the somatotopy

Spinal somatotopic activations were performed by electrically stimulating different dermatomes corresponding to lumbar levels L3, L4, and L5 using a standard stimulation protocol as described above (40 Hz, 5 seconds, 3 mA). Activation maps were computed as previously described. Once filtered, the centroid was calculated as the Z-score weighted barycenter. We only kept the horizontal (anteroposterior) coordinates. We measured the L5-L4 and L5-L3 distances; L5 was used as the horizontal axis origin for this somatotopic study. The averaged position of L4 and L3 was calculated by averaging all the corresponding distances.

Results were presented as level contour lines, those of the first, second, and third quartiles of the Z-score.

Ultrasound localization microscopy

For ULM, a catheter filled with saline was inserted in the rat jugular vein before the positioning of the animal and the laminectomy. Ultrasound localization microscopy was performed similarly to the methods described in. [START_REF] Couture | Ultrasound localization microscopy and super-resolution: a state of the art[END_REF][START_REF] Errico | Ultrafast ultrasound localization microscopy for deep super-resolution vascular imaging[END_REF][START_REF] Hingot | Microvascular flow dictates the compromise between spatial resolution and acquisition time in Ultrasound Localization Microscopy[END_REF] Briefly, 150 mL of SonoVue microbubbles (reconstituted in 5 mL of saline) were injected through the catheter in the jugular vein, followed by 100 mL of saline. Two hundred blocks of 600 successive compounded frames (with angles at 25˚, 0˚, 15˚, at 3000 Hz pulse repetition frequency with a 10-mm imaging depth) were acquired every second using the same probe (f 5 15 MHz, 128 elements, Vermon). The elevation focusing was achieved by a plastic lens to reach 500 mm at 8-mm elevation focal depth. The acquisitions were performed using the same setup and software as above. Small rigid motions were corrected through phase-correlation. [START_REF] Hingot | Subwavelength motioncorrection for ultrafast ultrasound localization microscopy[END_REF] Singular value decomposition filters [START_REF] Demen É | Spatiotemporal clutter filtering of ultrafast ultrasound data highly increases Doppler and fUltrasound sensitivity[END_REF] were used to extract bubble signals from surrounding tissues by removing only the 2 first singular values and by applying a second-order Butterworth highpass filter with a cutoff frequency of 20 Hz to remove remaining tissue signal. Microbubbles were localized using a weighted average method. [START_REF] Henriques | QuickPALM: 3D real-time photoactivation nanoscopy image processing in ImageJ[END_REF] Microbubbles were detected as the brightest local maxima on the image. Then, the position of the microbubbles within the pixel was determined using a weighted average on the intensity in the neighboring pixels. Microbubble tracks were formed by pairing microbubbles between frames using Matlab's ("simpletracker") tracking algorithm, which is based on the Hungarian method for assignment. [START_REF] Kuhn | The Hungarian method for the assignment problem[END_REF] Distances to all the microbubbles in the following frame were calculated. The Kuhn-Munkres algorithm then minimizes the total distance, thus connecting all the positions of a microbubble to form a collection of its position in the different frames. Velocities were calculated as the mean displacement between 2 consecutive frames. The microbubble track was completed using linear interpolation to fill missing points. The track was then smoothed using a sliding average along 5 positions, and the positions were rounded to the chosen pixel sizes.

Results

Vascular anatomy of the spinal cord

In a preliminary study, we tested the possibility to image the rat spinal cord through both the vertebral muscles and the vertebra. As expected, due to the reflection and attenuation of the wave propagation through the bone, the presence of the vertebra prevented the imaging of spinal blood vessels completely. As a consequence, as previously shown by others in the context of shear wave elastography and functional ultrasound, [START_REF] Song | Functional ultrasound imaging of spinal cord hemodynamic responses to epidural electrical stimulation: a feasibility study[END_REF][START_REF] Vergari | Intervertebral disc characterization by shear wave elastography: an in vitro preliminary study[END_REF] only small portions of spinal cord were visible through the intervertebral space. To study spinal cord physiology using fUS, we chose to use a large laminectomy at the level of the lumbar enlargement.

Imaging the baseline power Doppler signal (without external stimulation) at the rate of 1 Hz was sufficient to study the anatomical vascularization of the rat lumbar spinal cord.

Imaging first in a transverse view, we observed a strong signal in the posterior spinal vein (PSV), a vein that runs dorsally between both right and left dorsal horns; the anterospinal artery (ASA) that, on the contrary, runs ventrally between both ventral horns with fast blood flow speeds (.15 mm/s); and the dorsal extension of the ASA from the ventral to the dorsal horn (Figs. 1A andB). Surprisingly, no signal could be detected in some parts of the dorsal horn (laminae I-II for instance, Fig. 1B) and in the white matter of the ventral horn (Figs. 1B andD), apart from initial bifurcations of the ASA.

To identify these vascular structures, an in toto vascular staining of the vasculature using intracardiac perfusion of the fluorescent dye DiI was performed in 3 animals, as previously described. [START_REF] Li | Direct labeling and visualization of blood vessels with lipophilic carbocyanine dye DiI[END_REF] The staining produced (Fig. 1C) confirmed the presence of staining in the PSV and ASA and allowed the visualization of the vasculature in the parenchyma of the dorsal and ventral horns. Measurement in fixed tissues leads to a diameter value of 200 mm for the ASA, 20 mm for the CSA (central sulcal artery) and PSA (posterior spinal artery), and 5 mm for the vasculature in the parenchyma, as in the brain. The small size of the latter is the obvious reason why the vasculature of the parenchyma cannot be imaged in fUS imaging.

To circumvent this problem, we performed ultrafast ULM using intravenous injection of biocompatible microbubbles. This approach relies on the separation and localization of echoes of intravascular microbubbles to reconstruct super-resolved images. It is similar to the localization of blinking fluorescent sources in photoactivated localization microscopy, which improves the optical resolution by more than 10-fold. [START_REF] Betzig | Imaging intracellular fluorescent proteins at nanometer resolution[END_REF] We previously showed that this approach allows the imaging of small vessels (approx. 10 mm) in the rat cortex. [START_REF] Errico | Ultrafast ultrasound localization microscopy for deep super-resolution vascular imaging[END_REF] This in vivo microscopy localization confirmed that the spinal cord is a highly vascularized organ with a dense network of capillaries in both dorsal and ventral horns (Figs. 1D andF). The size of the blood vessels imaged using this approach is consistent with the measurements performed in fixed tissues (Fig. 1G, ie, a diameter inferior to 5 mm in capillaries). Blood flow speeds quantified using microbubble tracking in ULM showed fast flows at 15 to 20 mm/s and higher in the ASA, moderate flows at 5 to 8 mm/s in CSA and PSA arteries, and even lower ones in small-diameter vessels of the dorsal horn (1-5 mm•s 21 for 10-50 mm arteriole diameter). As expected, the measured blood velocities are much lower in the vein (PSV) and the parenchyma (10 mm•s 21 and 1-5 mm•s 21 , respectively), as compared to ASA.

Imaging in the sagittal view using both fUS imaging and ULM confirmed over several lumbar levels the previous observation of strong signals in arterial elements but weaker in the dorsal horn. We also observed new sets of "y" structures formed by the CSA in its anteroposterior projections (Fig. 1F). These structures were observed at the level of the midline (lateral 5 0 mm). This was confirmed using longitudinal spinal sections stained in toto using DiI (red arrow, Fig. 1G). An analysis of the microbubble's direction pointed out an ascending flow from the ASA to the dorsal horn, and a descending flow to the ventral horn (Figs. 1D andF). Finally, the comparison of the DiI staining in the transverse and sagittal views allowed us to observe that the density of blood vessels is much denser in the gray matter than in the white matter (Figs. 1C andG).

Functional activations in the rat spinal cord in health and disease

The second aim of our study was to determine the physiological properties of the spinal hemodynamic response induced by noxious and innocuous peripheral stimulations. Indeed, as shown in the past by several authors, such stimulations induce the activation of spinal neurons in different laminae of the dorsal horn, with specific somatotopic distributions. In addition, through mechanisms of neurovascular coupling, [START_REF] Frostig | Cortical functional architecture and local coupling between neuronal activity and the microcirculation revealed by in vivo high-resolution optical imaging of intrinsic signals[END_REF][START_REF] Iadecola | The neurovascular unit coming of age: a journey through neurovascular coupling in health and disease[END_REF] this neuronal activation should theoretically lead to a local hemodynamic response.

Mechanical stimulation of the hind paw

To address this question, we first characterized the hemodynamic response in the dorsal horn of the spinal cord induced by natural nociceptive mechanical stimulations. Manual application of nociceptive mechanical stimulations with a calibrated mechanical clamp induced a strong statistically significant local hemodynamic response in the dorsal horn (Figs. 2A-C). Quantification across all animals and acquisitions showed that this effect is restricted to the dorsal horn (weakly statistically significant in the ventral horn, red square, Fig. 2C) and to the lumbar level (not statistically significant more rostrally or caudally, see the green square located 4 mm caudally or rostrally from the activated area on Figs. 2A-C).

Natural stimulation of Ab fibers using dynamic stimulation of the foot pad with a brush induced also a statistically significant increase of the local SBV in the dorsal horn (Fig. 2D), but this effect was of a lesser amplitude. This effect was also restricted to the dorsal horn (not statistically significant in the ventral horn, see the red square on Fig. 2D) and to the lumbar level (not statistically significant more caudally, see the green square on Fig. 2D).

Electrical stimulation of the hind paw

To study the HRF and the fiber selectivity of the spinal cord stimulation, we chose to apply noxious or innocuous stimulations using electrical stimulations of either Ab, Ad, or C fibers. In this particular case, although the stimulations applied are obviously less natural, they have the great advantage of being triggered, and provide high temporal accuracy in the evoked response.

Because to date the HRF function of the rat spinal cord was never established, we determined it experimentally to further use it for the deconvolution of the stimulus patterns and therefore obtain a more accurate calculation of the hemodynamic response induced by electrical stimulations.

Determination of the hemodynamic response function

To accurately compute the hemodynamic responses to short stimuli induced by electrical stimulations, we determined hemodynamic impulse response in the rat spinal cord experimentally, to use it as a transfer function.

Figure 2. Functional ultrasound imaging reveals that natural mechanical noxious or innocuous stimulations induce a localized hemodynamic response in the spinal dorsal horn in anesthetized rats. (A-C) Hemodynamic response induced by noxious mechanical stimulation, applied using calibrated forceps (,4.5 N). A shows a typical example of the significant spinal hemodynamic response induced by noxious mechanical stimulation, imaged in the sagittal plane using fUS imaging. Are superimposed in (A): the vascularization of the cord imaged in fUS (black and white background) and the z-score of the pixels in the image (red/green pixels) and finally 3 regions of interest: the dorsal horn (DH, blue), ventral horn (VH, red), and the caudal dorsal horn (CA, green). (B) Averaged temporal signals (over all stimulations: N 5 52, from 5 animals) over time for the ROIs in the dorsal, showing a strong and robust increased spinal blood volume (SBV) in the lumbar dorsal horn (blue line), but not caudal horn during the stimulation (green line). Black squares indicate the stimulations. The red ROI is not shown for better clarity. The y-axis unit is the relative increase in spinal blood volume compared to baseline (portion of the signal without stimulation). (C and D) Box plots and overimposed individual values of changes in spinal blood volume (SBV) in all individual activations induced by either noxious (C) or innocuous (D) peripheral mechanical stimulations. In both cases, a statistically significant hemodynamic response is observed locally in the dorsal horn (not significant in the ventral horn [red], nor 4 mm more caudally [green]). Note that the amplitude of response evoked by noxious stimulations is much greater than for innocuous stimulations. Colored dots represent individual DSBV obtained for all individual stimulations. (C) Baseline: N 5 65 from 5 rats-dorsal horn: N 5 52 from 5 rats-ventral horn: N 5 44 from 5 rats-caudal area: N 5 48 from 5 rats. (D) Baseline: N 5 30 from 4 rats-dorsal horn: N 5 24 from 4 rats-ventral horn: N 5 24 from 4 rats-caudal area: N 5 24 from 4 rats. Statistical analysis: The statistical analysis was performed using a LMM. In this set of experiments, the fixed effect factor is the condition (eg, location of hemodynamic response [dorsal horn, ventral horn] or the intensity of fiber stimulation), whereas the animal and acquisition are random effect factors. Inclusive results of this analysis are shown in details in Supplementary table 1, and results of the post hoc tests are presented in the Supplementary table 2 (available at http://links.lww. com/PAIN/B176). They show a lack of effect of the random factors but statistical effect of various conditions. ***P , 0.001. NS, nonstatistically significant. To do so, we applied the shortest supramaximal peripheral stimulation of C fibers that elicited a reproducible response (2 seconds at 5 mA). Five acquisitions per rat on 2 rats were made to determine the HRF for ultrafast ultrasound imaging in the spinal cord of the rat. All acquisitions were averaged (Figs. 3B andC), and an inverse gamma distribution was individually fitted to each acquisition. Finally, we averaged the parameters obtained across all the experiments to have the analytical HRF parameters (Fig. 3C). We found the following HRF expression (Equation 1):

HRF : t→A: b a GðaÞ :t 2 ða 1 1Þ exp 2 b t ; (1)
with A 5 6.1 6 1.8, a 5 2.5 6 1.3, and b 5 12.7 6 7.0. The onset time (threshold corresponding to an SBV increase reaching twice the amplitude of the SD of the baseline) was of 1.4 6 0.2 seconds after the stimulation, at 14% 6 4% of SBV increase, with a maximal hemodynamic response at 3.0 6 0.9 seconds after the onset. It is interesting to note that this transfer function is very close to the transfer function between neuronal and vascular responses recently described by Aydin et al 2 in the olfactory bulb. Also, despite the difference of location in the central nervous system, the physiological properties of this HRF (onset time and time to peak) were very similar with previous observations in the brain of anesthetized rats: (2.0 seconds 6 0.3/5 6 0.7, [START_REF] Errico | Transcranial functional ultrasound imaging of the brain using microbubble-enhanced ultrasensitive Doppler[END_REF] ) and awake nonhuman primates (1.4 seconds 6 0.6/2.5 [START_REF] Dizeux | Functional ultrasound imaging of the brain reveals propagation of task-related brain activity in behaving primates[END_REF] ). These results are also consistent with measurements described in the literature using functional magnetic resonance imaging (fMRI). [START_REF] Giulietti | Characterization of the functional response in the human spinal cord: impulse-response function and linearity[END_REF][START_REF] Silva | Early temporal characteristics of cerebral blood flow and deoxyhemoglobin changes during somatosensory stimulation[END_REF][START_REF] Uluda G | Linking brain vascular physiology to hemodynamic response in ultra-high field MRI[END_REF] Some authors reported very different time response (10 seconds) in the human spinal cord. The fact that this was for a motor task (and not sensory) might explain the difference. However, as we measured blood volume and not deoxyhemoglobin, neither initial dip nor poststimulus undershoot was present.

Functional activation in the spinal cord using electrical stimulation

To evaluate the ability of the different afferent modalities to induce an evoked hemodynamic response, we stimulated electrically different groups of peripheral afferents, using previously described parameters. [START_REF] Chang | A fMRI study of brain activations during non-noxious and noxious electrical stimulation of the sciatic nerve of rats[END_REF][START_REF] Bars | Diffuse noxious inhibitory controls (DNIC). I. Effects on dorsal horn convergent neurones in the rat[END_REF][START_REF] Silva | Early temporal characteristics of cerebral blood flow and deoxyhemoglobin changes during somatosensory stimulation[END_REF] Note that when the intensity of stimulation is increasing, fibers of lower threshold are also stimulated. Therefore, Ab fibers are also stimulated when Ad and/or C fibers are stimulated.

Electrical stimulations of Ab fibers (0.5 mΑ, [START_REF] Chang | A fMRI study of brain activations during non-noxious and noxious electrical stimulation of the sciatic nerve of rats[END_REF] ) did not evoke any hemodynamic response (20.4 6 6.0% SBV) compared to baseline (0.5 6 0.7% SBV, P . 0.05). Electrical stimulation of Ad (1.5 mΑ, [START_REF] Chang | A fMRI study of brain activations during non-noxious and noxious electrical stimulation of the sciatic nerve of rats[END_REF][START_REF] Silva | Early temporal characteristics of cerebral blood flow and deoxyhemoglobin changes during somatosensory stimulation[END_REF] ) and C fibers (3 mA, [START_REF] Chang | A fMRI study of brain activations during non-noxious and noxious electrical stimulation of the sciatic nerve of rats[END_REF][START_REF] Bars | Diffuse noxious inhibitory controls (DNIC). I. Effects on dorsal horn convergent neurones in the rat[END_REF] ) fibers, however, induced a statistically significant hemodynamic activation (Figs. 4A-C, 35.4 6 5.1% SBV, P , 0.001/73.4 6 7.5% SBV, P , 0.001, respectively). The amplitude of the response induced by C fibers was greater than the one induced by Ad fibers, suggesting that C fibers have a major role in this effect.

Finally, spinal topic application of ketamine (an NMDA receptor antagonist) completely blocked the hemodynamic response induced by C-fiber stimulation (20.4 6 1.5% SBV, P . 0.05 compared to baseline). These results suggest that (1) the evoked spinal neurovascular response is fiber-specific and (2) that glutamate (through the involvement of the NMDA receptor) plays a major role in the hemodynamic activation induced by C fibers. As further discussed in the discussion section, this hemodynamic response is likely to be the dual consequence of specific neuronal activation and neurovascular coupling.

Spatial localization of the spinal hemodynamic response

All neurovascular responses observed in the dorsal horn imaged in the sagittal plane were of heterogenous oval shape (Fig. 4A), typically (after C-fiber stimulation) 2 mm long.

The reason for this shape and size could be (1) the spinal areas where responsive neurons are located and (2) the spread of increased local blood supply, 2 hypotheses that are not mutually exclusive.

At each spinal level, activated primary sensory afferents enter the dorsal horn of the spinal cord ipsilaterally, terminating in different laminae of the dorsal horn, locally activating various subpopulation of dorsal horn neurons. Therefore, imaging the hemodynamic response induced by C fibers in the transverse plane could provide crucial clue on the specificity of the hemodynamic response measured in fUS towards neuronal activation. As shown in Figure 4D, despite the thin network of blood vessels and weak Doppler signal in laminae I-III (Fig. 1), a strong increased blood volume was observed in laminae I-V (Figs. 4D andE) of the dorsal horn ipsilateral to noxious stimulation (44.0 6 3.1% SBV vs 0.2 6 0.7% SBV for baseline, P , 0.001). This effect was not observed in the contralateral dorsal horn (1.4 6 1.8% SBV, P . 0.05 compared to baseline), suggesting a specific hemodynamic response where activated neurons are known to be located. Therefore, this restricted area of activation of approximately 1 mm suggests that fUS imaging faithfully measures the local increased blood volume induced by neuronal activation.

Spinal hemodynamic response respects the somatotopy of afferent entries in the dorsal horn

Compelling evidence brought forward from different but complementary techniques ranging from electrophysiological recordings [START_REF] Men Étrey | An analysis of response properties of spinal cord dorsal horn neurones to nonnoxious and noxious stimuli in the spinal rat[END_REF][START_REF] Woolf | Somatotopic organization of cutaneous afferent terminals and dorsal horn neuronal receptive fields in the superficial and deep laminae of the rat lumbar spinal cord[END_REF] to tracing experiments [START_REF] Swett | The somatotopic organization of primary afferent terminals in the superficial laminae of the dorsal horn of the rat spinal cord[END_REF][START_REF] Woolf | Somatotopic organization of cutaneous afferent terminals and dorsal horn neuronal receptive fields in the superficial and deep laminae of the rat lumbar spinal cord[END_REF] and Evans blue extravasations [START_REF] Takahashi | Dermatomes and the central organization of dermatomes and body surface regions in the spinal cord dorsal horn in rats[END_REF][START_REF] Takahashi | Dermatomes in the rat limbs as determined by antidromic stimulation of sensory C-fibers in spinal nerves[END_REF] have shown that entries of afferent fibers are somatotopically organized in the spinal cord, activating nociceptive specific or nonspecific neurons in the corresponding spinal cord levels. However, due to the difficulty to image in the spinal cord, no neuroimaging technique so far was able to verify such mechanisms in vivo. Combining the large field of view achievable in fUS imaging of the spinal cord, the high spatial resolution, and sensitivity of the technique, using electrically stimulated C fibers in different dermatomes consecutively (L3, L4, or L5) in the same animal, we analyzed the spatial shift of the activated areas. As shown in Figure 5A (example in one animal), the areas of activation shift along the neuroaxis caudally as the stimulated dermatomes shift from L3, L4, to L5. This phenomenon was quantified over several animals by calculating the weighted centroid of the activated area (white stars in Fig. 5B), using the location of L5 in each animal as the spatial reference. Note that surroundings of the activated areas are shown as correlation level curves in Figures 5B andC, using the following color code: red (L3), blue (L4), and green (L5). This analysis confirmed that the area functionally activated by C-fiber stimulation shifted caudally of approximately 1.2 6 0.4 mm between L5 and L4 and 2.5 6 0.3 mm between L5 and L3.

In a transverse plane of imaging, this analysis showed that the functionally activated areas were restricted to the whole laminae I, II, III, and the medial aspect of laminae IV and V, with a weighted centroid located in the medial part of lamina III (Fig. 5E). These are the spinal locations of nociceptive specific and nonspecific neurons.

Imaging of Ab fibers induced mechanical allodynia

Allodynia is an abnormal pain perception induced by the application of an innocuous stimulation. It is probably the most frequent and debilitating clinical symptom observed in chronic pain diseases (both in clinical condition and in animal models [START_REF] Jensen | Allodynia and hyperalgesia in neuropathic pain: clinical manifestations and mechanisms[END_REF] ). Because pathways mediating pain and touch are segregated in physiological conditions, there is a large consensus that allodynia might be due to the functional reorganization of spinal circuits. [START_REF] Kuner | Central mechanisms of pathological pain[END_REF][START_REF] Latremoliere | Central sensitization: a generator of pain hypersensitivity by central neural plasticity[END_REF] During the installation of peripheral inflammation, due to the tonic activation of sensitized primary afferent fibers, the spinal mechanism of central sensitization takes place, in which spinal neurons are in a state of facilitation, potentiation, augmentation, or amplification. [START_REF] Latremoliere | Central sensitization: a generator of pain hypersensitivity by central neural plasticity[END_REF] In this inflammatory condition, wide-dynamic-range neurons have an increased spontaneous firing rate, and the application of a lowthreshold mechanical stimulation induces both a pain sensation (allodynia) and an increased firing rate of these neurons, but not of nociceptive-specific neurons. [START_REF] Mcgaraughty | Characterization and pharmacological modulation of noci-responsive deep dorsal horn neurons across diverse rat models of pathological pain[END_REF] In this study, we reevaluated the neurobiological mechanisms underlying allodynia in inflammatory pain, imaging the neurovascular response to mechanical tactile allodynia in a model of central sensitization induced by short-term inflammatory pain (24 hours after unilateral subcutaneous injection of Complete Freund's Adjuvant [25 uL, Sigma Aldrich, France] in the hind paw). Using the large field of view and the temporal resolution of fUS imaging, we imaged the hemodynamic response evoked by stimulation of low-threshold mechanoreceptors (induced by brush stimulation on the hind paw) and observed an increased spinal hemodynamic response in inflamed animals compared to the controls (Fig. 6). This increased responsiveness was characterized by an increased evoked local SBV in the dorsal horn (Fig. 6D).

Discussion

Functional ultrasound imaging is a relatively new modality of neuroimaging that allows, through ultrafast imaging of blood flow, the study of functional activations [START_REF] Bimbard | Multi-scale mapping along the auditory hierarchy using high-resolution functional UltraSound in the awake ferret[END_REF][START_REF] Gesnik | 3D functional ultrasound imaging of the cerebral visual system in rodents[END_REF][START_REF] Osmanski | Functional ultrasound imaging reveals different odor-evoked patterns of vascular activity in the main olfactory bulb and the anterior piriform cortex[END_REF] and functional connectivity. [START_REF] Osmanski | Functional ultrasound imaging reveals different odor-evoked patterns of vascular activity in the main olfactory bulb and the anterior piriform cortex[END_REF][START_REF] Rahal | Ultrafast ultrasound imaging pattern analysis reveals distinctive dynamic brain states and potent sub-network alterations in arthritic animals[END_REF] Compared to other modalities of neuroimaging, it has a very good spatial resolution, competing with the best fMRI scanners. [START_REF] Cohen-Adad | Functional magnetic resonance imaging of the spinal cord: current status and future developments[END_REF][START_REF] Stroman | The current state-of-the-art of spinal cord imaging: Methods[END_REF][START_REF] Tsurugizawa | Awake functional MRI detects neural circuit dysfunction in a mouse model of autism[END_REF] The temporal resolution is remarkably high. However, because the physiological process measured by fUS imaging is happening at slower temporal and larger spatial scales than the actual imaging modality due to the neurovascular coupling, both the effective spatial and temporal resolutions are thereby limited. It has to be noted, nevertheless, that this fast temporal resolution enabled us in the past to measure directional propagation of hemodynamic signals within layers of the cortex in the primate brain, thanks to lag correlation at high sampling rates (,10 ms) faster than the cardiac cycle or the hemodynamic response. [START_REF] Dizeux | Functional ultrasound imaging of the brain reveals propagation of task-related brain activity in behaving primates[END_REF] The effective spatial resolution provided in fUS imaging was also shown to be comparable to the resolution of the CBV images in another recent work. [START_REF] Bimbard | Multi-scale mapping along the auditory hierarchy using high-resolution functional UltraSound in the awake ferret[END_REF] During auditory stimuli, Bimbard et al. [START_REF] Bimbard | Multi-scale mapping along the auditory hierarchy using high-resolution functional UltraSound in the awake ferret[END_REF] demonstrated in awake ferrets that fUS imaging was able to provide voxel-to-voxel independent information both in the auditory cortex and deepseated regions with an effective resolution of 100 mm for voxel responsiveness and 300 mm for voxel frequency tuning. Probably the 2 main advantages of fUS imaging are its high sensitivity (comparable to calcium imaging at the pixel level, as demonstrated recently [START_REF] Boido | Mesoscopic and microscopic imaging of sensory responses in the same animal[END_REF] ) and its portability compared to other techniques. Apart for very recent developments in 3D imaging, [START_REF] Rabut | 4D functional ultrasound imaging of whole-brain activity in rodents[END_REF] the main disadvantage of the technique (compared to others) is the imaging in 2 dimensions. Also, as ultrasounds are attenuated by bones, apart in conditions where the skull is thin (rat pups, mice [START_REF] Tiran | Transcranial functional ultrasound imaging in freely moving awake mice and anesthetized young rats without contrast agent[END_REF] ) or naturally (through the fontanel in infants [START_REF] Demene | Functional ultrasound imaging of brain activity in human newborns[END_REF] ) or surgically opened (perioperative interventions in adults [START_REF] Imbault | Intraoperative functional ultrasound imaging of human brain activity[END_REF] ), fUS imaging requires more invasive preparations (such as a laminectomy in this study).

Functional activations in the rat's spinal cord

It has long been established that peripheral application of a noxious stimulus induces in the spinal cord a concomitant increased neuronal activity, [START_REF] Coghill | Spatial distribution of nociceptive processing in the rat spinal cord[END_REF] energy consumption, [START_REF] Kadekaro | Differential effects of electrical stimulation of sciatic nerve on metabolic activity in spinal cord and dorsal root ganglion in the rat[END_REF] metabolic response, [START_REF] Marcus | Regulation of total and regional spinal cord blood flow[END_REF] and BOLD signal. [START_REF] Lawrence | Correlation of functional activation in the rat spinal cord with neuronal activation detected by immunohistochemistry[END_REF] Interestingly, the physiological properties of the neurovascular coupling in the brain and the spinal cord are very similar. [START_REF] Cohen-Adad | BOLD signal responses to controlled hypercapnia in human spinal cord[END_REF][START_REF] Nix | Comparison of vascular reactivity in spinal cord and brain[END_REF] One important difference, however, is the tight neurovascular coupling that exists in the spinal cord between neuronal activation and increased blood flow, independently of changes in mean arterial blood pressure, a feature not observed in the brain after C-fiber stimulations. [START_REF] Pich É | Tight neurovascular coupling in the spinal cord during nociceptive stimulation in intact and spinal rats[END_REF] Our study showed, using fUS imaging, a C-fiber-evoked hemodynamic response (elicited by activation of either natural noxious mechanical stimulations or electrically activated C fibers) in the ipsilateral dorsal horn to the stimulus and was restricted spatially to approximately 1 to 2 mm, after the somatotopy of the dermatome stimulated. The localization in the ipsilateral dorsal horn is consistent with previous observations by other groups using fMRI, [START_REF] Lilja | Blood oxygenation level-dependent visualization of synaptic relay stations of sensory pathways along the neuroaxis in response to graded sensory stimulation of a limb[END_REF][START_REF] Majcher | Simultaneous functional magnetic resonance imaging in the rat spinal cord and brain[END_REF][START_REF] Malisza | Functional MRI of the rat lumbar spinal cord involving painful stimulation and the effect of peripheral joint mobilization[END_REF][START_REF] Zhao | BOLD and blood volumeweighted fMRI of rat lumbar spinal cord during non-noxious and noxious electrical hindpaw stimulation[END_REF] intrinsic optical imaging, [START_REF] Sasaki | Optical imaging of intrinsic signals induced by peripheral nerve stimulation in the in vivo rat spinal cord[END_REF] and in light scattering imaging, [START_REF] He | Hemodynamic and light-scattering changes of rat spinal cord and primary somatosensory cortex in response to innocuous and noxious stimuli[END_REF] and suggests, as discussed by these authors, that the spinal neurovascular response observed is extremely locally restricted to ipsilateral site and to laminae in the dorsal horn where activated neurons are located. A large number of studies in human subjects also reported a contralateral activation after noxious heat stimulations [START_REF] Cahill | Mapping of neural activity produced by thermal pain in the healthy human spinal cord and brain stem: a functional magnetic resonance imaging study[END_REF][START_REF] Rempe | Spinal and supraspinal processing of thermal stimuli: an fMRI study: processing of Thermal Stimuli[END_REF][START_REF] Stroman | Spinal fMRI investigation of human spinal cord function over a range of innocuous thermal sensory stimuli and study-related emotional influences[END_REF][START_REF] Weber | Functional magnetic resonance imaging of the cervical spinal cord during thermal stimulation across consecutive runs[END_REF] (and sometimes innocuous ones [START_REF] Rempe | Spinal and supraspinal processing of thermal stimuli: an fMRI study: processing of Thermal Stimuli[END_REF] ). But this feature was not always observed in human subjects [START_REF] Brooks | Stimulus site and modality dependence of functional activity within the human spinal cord[END_REF][START_REF] Eippert | Direct evidence for spinal cord involvement in placebo analgesia[END_REF][START_REF] Nash | Functional magnetic resonance imaging identifies somatotopic organization of nociception in the human spinal cord[END_REF] and usually not observed in anesthetized rodents. The nature and intensity of the stimulations, as well as the presence of anesthesia, might explain the difference of results. The observations both in human [START_REF] Summers | A quantitative comparison of BOLD fMRI responses to noxious and innocuous stimuli in the human spinal cord[END_REF] and macaque, [START_REF] Yang | Differential fMRI activation patterns to noxious heat and tactile stimuli in the primate spinal cord[END_REF] in the same subjects of a unique ipsilateral activation induced by innocuous stimulation but both ipsilateral and contralateral activations by noxious heat, support this hypothesis.

Stimulation of the interdigit territories are known to activate diverse dermatomes, inducing the transmission of sensory signal in the tibial and (to a smaller extent) to the saphenous nerves to the spinal segments L3 to L5 where they activate postsynaptic neurons. [START_REF] Takahashi | Dermatomes and the central organization of dermatomes and body surface regions in the spinal cord dorsal horn in rats[END_REF][START_REF] Takahashi | Dermatomes in the rat limbs as determined by antidromic stimulation of sensory C-fibers in spinal nerves[END_REF] Study of the somatotopy of hemodynamic responses in the spinal cord using neuroimaging was performed previously in human subjects [START_REF] Brooks | Stimulus site and modality dependence of functional activity within the human spinal cord[END_REF][START_REF] Cahill | Mapping of neural activity produced by thermal pain in the healthy human spinal cord and brain stem: a functional magnetic resonance imaging study[END_REF][START_REF] Nash | Functional magnetic resonance imaging identifies somatotopic organization of nociception in the human spinal cord[END_REF][START_REF] Stroman | Mapping of neuronal function in the healthy and injured human spinal cord with spinal fMRI[END_REF][START_REF] Stroman | Somatotopic arrangement of thermal sensory regions in the healthy human spinal cord determined by means of spinal cord functional MRI[END_REF] but, to the best of our knowledge, it had never been performed in rodents. Despite the spatial heterogeneity of the activated area, we identified the centroid of this activation area and demonstrated its spatial shift along the lumbar segment when stimulations were successively applied to the L3, L4, or L5 dermatomes. The spatial translation of these activation areas along the lumbar segment (L3-L5: 2.5 6 0.3 mm, L4-L5: 1.2 6 0.4 mm) was found to be consistent with the known neuroanatomy of the rat spinal cord [START_REF] Waibl | the Albino Rat (Rattus norvegicus)[END_REF] and location of the activated areas in these specific segments. [START_REF] Swett | The somatotopic organization of primary afferent terminals in the superficial laminae of the dorsal horn of the rat spinal cord[END_REF][START_REF] Takahashi | Dermatomes and the central organization of dermatomes and body surface regions in the spinal cord dorsal horn in rats[END_REF][START_REF] Woolf | Somatotopic organization of cutaneous afferent terminals and dorsal horn neuronal receptive fields in the superficial and deep laminae of the rat lumbar spinal cord[END_REF] It is interesting to note that although the vasculature of the parenchyma could not be detected in fUS imaging (Fig. 1, and was only detected using ULM), functional activations (ie, measurements of Doppler signal in correlation with the stimulus) could be easily detected using fUS imaging, suggesting the high sensitivity of the technique to detect subtle changes of task-evoked SBV. This agrees with previous observations that fUS imaging is sensitive enough to measure neurovascular changes at the pixel level, within single trials. [START_REF] Boido | Mesoscopic and microscopic imaging of sensory responses in the same animal[END_REF] Surprisingly, the shape of the spinal area activated, observed in the sagittal axis, was spatially heterogeneous. Independently of the stimulus intensity, it was consistently formed by 2 to 3 discontinued ovals of 700-mm diameter in the anteroposterior axis. It can be hypothesized that this heterogenous form could be due to (1) the possible heterogeneity in either the afferent fibers and/or the spatial distribution of spinal neurons and/or (2) a heterogeneity in the diameter of the blood vessels irrigating the spinal segment. It has long been described that the somatotopic distribution of tibial, sural, and saphenous nerves terminates in a widespread manner in laminae I-V of the lumbar segments L3 to L5. [START_REF] Shortland | Morphology and somatotopic organization of the central terminals of hindlimb hair follicle merents in the rat lumbar spinal cord[END_REF][START_REF] Swett | The somatotopic organization of primary afferent terminals in the superficial laminae of the dorsal horn of the rat spinal cord[END_REF] Concomitantly, postsynaptic neurons are located homogeneously either in superficial laminae (laminae I, where nociceptive specific neurons are located) or in deep laminae (where nonspecific nociceptive neurons are located). [START_REF] Bullitt | Somatotopy of spinal nociceptive processing[END_REF][START_REF] Wilson | Plasticity of cutaneous primary afferent projections to the spinal dorsal horn[END_REF][START_REF] Woolf | Somatotopic organization of cutaneous afferent terminals and dorsal horn neuronal receptive fields in the superficial and deep laminae of the rat lumbar spinal cord[END_REF] Therefore, because both afferent and postsynaptic neurons are quite homogeneously distributed throughout each spinal segment, it is likely that the heterogenous shape of the activated areas has a vascular origin.

Intensity-specific hemodynamic response

The comparison of the efficacy of the various types of afferent fiber stimulation to induce an evoked spinal hemodynamic response showed that C fibers and Ad stimulation induce a robust hemodynamic response due to local neuronal activation through an NMDA receptor-dependent mechanism. Zhao et al. [START_REF] Zhao | fMRI investigation of the effect of local and systemic lidocaine on noxious electrical stimulation-induced activation in spinal cord[END_REF] showed previously that imaging of the spinal cord allows pharmacological studies in the rat. Stimulation of large myelinated fibers through innocuous stimulations gave rise to increased SBV, but to a smaller extent than other types of fibers. Such increased neurovascular coupling induced by innocuous stimulation has been detected in fMRI in humans, [START_REF] Brooks | Stimulus site and modality dependence of functional activity within the human spinal cord[END_REF][START_REF] Nash | Functional magnetic resonance imaging identifies somatotopic organization of nociception in the human spinal cord[END_REF][START_REF] Rempe | Spinal and supraspinal processing of thermal stimuli: an fMRI study: processing of Thermal Stimuli[END_REF][START_REF] Sprenger | Spinal cord-midbrain functional connectivity is related to perceived pain intensity: a combined spinocortical fMRI study[END_REF][START_REF] Stroman | Mapping of neuronal function in the healthy and injured human spinal cord with spinal fMRI[END_REF][START_REF] Stroman | Somatotopic arrangement of thermal sensory regions in the healthy human spinal cord determined by means of spinal cord functional MRI[END_REF] but not in animal models, even using intravenous injection of ultrasmall superparamagnetic iron oxide particles using a gradient echo. [START_REF] Zhao | BOLD and blood volumeweighted fMRI of rat lumbar spinal cord during non-noxious and noxious electrical hindpaw stimulation[END_REF] Only one study managed to detect it using light scattering, [START_REF] He | Hemodynamic and light-scattering changes of rat spinal cord and primary somatosensory cortex in response to innocuous and noxious stimuli[END_REF] but this approach lacks spatial discrimination.

Imaging spinal cord's central sensitization

Persistent pain associated with injury or diseases can result from both peripheral and central sensitization. Central sensitization is characterized by a state of hyperexcitability of the central nervous system, in which application of nonnociceptive stimuli (such as the gentle brush of the skin) is perceived as noxious, a clinical symptom termed "allodynia." [START_REF] Jensen | Allodynia and hyperalgesia in neuropathic pain: clinical manifestations and mechanisms[END_REF] Among the approaches usually used to study these mechanisms, electrophysiological recordings measure the neuronal activity with excellent temporal resolution, but in an extremely small area. Other approaches (immunohistochemistry) do not allow dynamic reading over time. Imaging of central sensitization using fMRI has previously reported an increased spinal responsiveness in healthy volunteers [START_REF] Weber | Functional magnetic resonance imaging of the cervical spinal cord during thermal stimulation across consecutive runs[END_REF] and diabetic rats. [START_REF] Malisza | Functional magnetic resonance imaging of the spinal cord during sensory stimulation in diabetic rats[END_REF] Our study measured the spinal responsiveness to allodynic stimulations in the contact of inflammatory pain. Consistent with observations of Rempe et al., [START_REF] Rempe | Spinal and supraspinal processing of thermal stimuli: an fMRI study: processing of Thermal Stimuli[END_REF] we observed an exacerbated response in the ipsilateral dorsal horn. These results, that suggest an increased neuronal activity in this model, are consistent with previous demonstrations that in inflammatory conditions, wide-dynamic-range neurons have both an enhanced spontaneous activity and exacerbated response to innocuous mechanical stimulations. [START_REF] Mcgaraughty | Characterization and pharmacological modulation of noci-responsive deep dorsal horn neurons across diverse rat models of pathological pain[END_REF] As evidenced in other parts of the CNS, it is very likely that the neuronal hyperexcitability, through a glutamate-mediated signaling mechanisms (such as nitric oxide [from neurons] and arachidonic acid derivatives [from astrocytes]), leads to this enhanced hemodynamic response. [START_REF] Attwell | Glial and neuronal control of brain blood flow[END_REF][START_REF] Lecrux | Neuronal networks and mediators of cortical neurovascular coupling responses in normal and altered brain states[END_REF] In conclusion, this study revealed the ability of fUS imaging to provide insightful information on the spinal physiology and its potential in the future for the development of analgesic drugs.

Figure 1 .

 1 Figure 1. Vascular anatomy of the rat spinal cord. (A) Diagram showing the major veins and arteries in the spinal cord and their position with respect to the sagittal and transverse imaging planes. (B-D) Imaging of the rat spinal blood volume in transverse view, using fUS imaging (B) or using ultrasound localization microscopy imaging (D), showing a strong spinal blood volume in the posterior spinal vein (PSV), a vein that runs dorsally between both right and left dorsal horns; the anterospinal artery (ASA) that, on the contrary, runs ventrally between both ventral horn; and the dorsal extension of the ASA from ventral to the dorsal horn. Although classic fUS imaging does not reveal any Doppler signal in the white matter and the most part of the dorsal horn (laminae I and II), ULM imaging (D) and in toto staining of the vascularization (C) reveals a thin vascular network of capillaries in these areas. (E-F) Functional ultrasound (E) and ULM imaging (F) of the spinal cord in sagittal view, at the level of the midline, showing highly contrasted "Y" structures, that originate from the ASA, whose vertical part is the CSA, while splitting in 2 branches at the level of the central canal, either dorsally in the dorsal horn (positive flux, red arterioles, E), or ventrally in the ventral horn (negative flux, blue arterioles, F). (D and F) Ultrasound localization microscopy imaging and analysis of the directionality of the flow. Ultrasound localization microscopy imaging was achieved with the intravenous injection of biocompatible microbubbles. The super-resolution image was reconstructed from microbubble tracking and false colorized. The colour indicates the directionality of the bubbles (blue: going ventrally, red: going dorsally) and the colour intensity indicates the velocity of the bubbles in these directions. (H) Quantitative measurements of both the diameter of some spinal blood vessels and the velocity of bubbles in these blood vessels. (C and G) In toto staining using vascular injection of DiI reveals a dense network of blood vessels in the spinal dorsal horn, and a sparser network in the ventral horn. The central canal (cc) that delineates the dorsal from ventral horns is clearly visible as a lack of staining. The arrow in G points at one of the "Y" structures, which is clearly a vascular structure of large diameter. ASA, anterior spinal artery; CSA, central sulcal artery; DH, dorsal horn; PSV, posterior spinal vein; VH, ventral horn; wm VH, white matter of the ventral horn; ULM, ultrasound localization microscopy.

  1052 J. Claron et al. • 162 (2021) 1047-1059 PAIN ®

Figure 3 .

 3 Figure 3. Determination of the spinal hemodynamic response function using functional ultrafast ultrasound imaging. (A) Image illustrating the implantation of the stimulation acupuncture electrodes in the rat hind paws. (B) Average changes (blue line) 6 standard error deviation of the mean of the spinal blood volume evoked by peripheral noxious electrical stimulations in all animals and stimulations (2 rats, 5 acquisitions per rat, N 5 20 individual stimulations). Black line illustrates the stimulation (2 seconds, 40 Hz, 5 mA). (C) Analytical determination of the fUS HRF using inverse gamma distribution. These acquisitions were fitted individually, and the final values are the mean of each fit. The figure shows the individual values (blue dots) and the fitted function (red line). The equation and parameters are indicated in the panel (A 5 6.1 6 1.8, a 5 2.5 6 1.3, b 5 12.7 6 7.0). (D-E) Example of activation maps in one representative animal, showing the Z-score of all the pixels in the image with a P-value , (0.05/image size in pixel, Bonferroni correction). (E) Without the use of the HRF in the correlation model, whereas (D) is with the use of our hemodynamic response function convoluted with the stimulation pattern. Orientation: caudal (C)-rostral (R) is indicated. fUS, functional ultrasound; HRF, hemodynamic response function.

Figure 4 .

 4 Figure 4. The spinal hemodynamic response induced by peripheral electrical stimulation is fiber-specific and localized ipsilaterally. (A) Typical example of functional activation induced in the dorsal horn in a naive rat by the electrical stimulation of C fibers (imaged in sagittal view). The functional activation was calculated using a global linear model, using the analytical HRF obtained previously. The color bar indicates the Z-score. (B) Averaged temporal dynamic of Doppler signal in relative augmentation of the spinal blood volume induced by electrical stimulation of C fibers in all stimulations (3 mA, N 5 84, from 9 rats). Black dots represent electrical stimulations (5 seconds). The blue line represents the average experimental Doppler signal over all stimulations (N 5 84), whereas the orange line represents the theoretical response to the stimulation pattern describing the spinal cord HRF. (C) Statistical analysis of hemodynamic responses induced by either 0.5 mA (Ab fibers), 1.5 mA (Ad fibers), or 3 mA (C fibers) stimulation showing that only Ad-and C-fiber stimulations induced a significant hemodynamic response. The C-fiber response was specifically blocked by NMDA-receptor antagonist, ketamine, suggesting that glutamatergic neurotransmission through the NMDA receptors is a necessary mechanism in this neurovascular response. (D-E) Functional ultrasound imaging of the hemodynamic response induced by C-fiber stimulations viewed in a transverse plane. (D) Example of a functional activation map in one rat showing that the area activated is in the ipsilateral dorsal horn. The scale applies to both panels of (D). (E) Statistical analysis of these functional activations over all individual stimulations (N 5 48, from 5 animals), showing that the response is statistically significant only in the ipsilateral dorsal horn, suggesting a very localized and specific neurovascular response to the C-fiber stimulations. Results were expressed as in (C). In (C and E), colored dots represent individual DSBV obtained for all individual stimulations (C: Baseline: N 5 105 from 9 rats-3 mA: N 5 84 from 9 rats-1.5 mA: N 5 48 from 5 rats-0.5 mA: N 5 36 from 9 rats-3 mA 1 ketamine: N 5 16 from 4 rats). (E) Baseline: N 5 60 from 5 rats-ipsilateral and contralateral: N 5 48 from 5 rats). *P , 0.05, ***P , 0.001. Orientation: caudal (C)-rostral (R), dorsal horn (DH), and ventral horn (VH) are indicated. HRF, hemodynamic response function; NS, nonstatistically significant.

5 .

 5 The spinal hemodynamic response respects the somatotopy of afferent entries in the dorsal horn. (A) Typical example in one animal of spatial shift the spinal area functionally activated (z-score map) by C-fiber stimulation in the dermatomes L3, L4, and L5. The scale bar (2 mm) applies to all panels in (A). (B-E) Quantification of the spatial areas functionally activated over several animals (and their spatial shift) by calculation of the weighted centroid of these areas, in the sagittal plane (B-D) and in the transverse plane (E). (B and E) The white stars are the weighted centroid of the correlations averaged over all the acquisitions (N 5 6 individual acquisitions in B-D; N 5 6 in E). (B-D) Color code of the different dermatomes: red: L3, blue: L4, and green: L5. The blue, red, and green curves represent contour lines corresponding to the first, second, and third quartiles of the correlations, respectively. (C-D) The centroid of L5 (star) was used as an invariant point to allow relative calculation for L4 (diamonds) and L3 (circle) positions on different animals. Only the horizontal position was taken into account; the vertical position was ignored. (N 5 4 rats; 4 for L5, 3 for L4, and 3 for L3). (D) The yellow star represents the centroid of the calculated relative positions for the dermatome L4 (diamonds, blue) and L3 (circle, red). The scale bar in (D) applies for (C and D). (E) Same as (B) in a transverse plane at lumbar level L4 (obtained from one animal and 6 activation maps). The bottom panel in E allows the localization of the laminae where the response is observed. Orientation: caudal (C)-rostral (R), dorsal horn (DH), ventral horn (VH), medial (M), and lateral (L). 1056 J. Claron et al. • 162 (2021) 1047-1059 PAIN ®

April 2021 • Volume 162 •Number 4 www.painjournalonline.com

Acknowledgements

The authors thank the CNRS, INSERM, and ESPCI Paris for their financial support, and Iconeus company for their collaboration and support on the software development. Supported by the Foundation Bettencourt-Schueller, the LABEX WIFI (Laboratory of Excellence within the French Program "Investments for the Future") under reference ANR-10-IDEX-0001-02 PSL, the National Agency for Research under the program "future investments" with the reference ANR-10-EQPX-15, and the European Research Council (ERC) Advanced Grant FUSIMAGINE. Line Rahal was supported by a fellowship from the "Fondation NRJ"-Institut de France award granted to Mickael Tanter. Finally, this work was supported by the Inserm ART (Technology Research Accelerator) in "Biomedical Ultrasound." The authors thank Miss Alice Fermigier for her early involvement in the project, Mr Alexandre Dizeux for his help in graphics design, and for his drawing of the cover, Mr J ér émy Ferrier for the critical reading of the manuscript, and Mrs No émie Dherbe and Nathalie Ialy-Radio for animal husbandry.

Conflict of interest statement

M.T. and T.D. are co-founders and shareholders of Iconeus company and are co-inventors of several patents in the field of neurofunctional ultrasound licenced to Iconeus company.The remaining authors have no conflicts of interest to declare.

Appendix A. Supplemental digital content

Supplemental digital content associated with this article can be found online at http://links.lww.com/PAIN/B176. Figure 6. After short-term peripheral inflammation, a peripheral mechanical stimulation (brush, ie, allodynic stimulation) induces an enhanced spinal hemodynamic response. (A) Drawing illustrating the experimental design. At T 5 0, animals received a unilateral intraplantar injection of either saline (control group) or inflamed group (Complete Freund's Adjuvant, CFA). The functional spinal response to a brush stimulation was studied 24 hours later. (B) Typical example of a functional activation in the dorsal horn during mechanical stimulation of the hind paw in one inflamed animal, using a brush. Quantification and presentation of the functional activation was determined as described above. (C) Changes in averaged Doppler signal over time (induced by brush stimulation) in control (blue, 7 rats, N 5 76 individual hemodynamic responses) and CFA inflamed rats (red, 7 rats, N 5 96 individual hemodynamic responses) in the ipsilateral dorsal horn. The black squares indicate the stimulations of the injected hind paw. The traces were averaged values over all responses in each group. (D) The statistical analysis of these functional activations reveals an enhanced responsiveness in the dorsal horn of inflamed animals (D). Colored dots represent individual DSBV values obtained for all individual stimulations (N 5 76 in control animals; N 5 96 in CFA animals). Statistical analysis (mixed linear model) shows a statistically significant effect of the treatment (Saline/CFA), but also of the animal, suggesting that this effect was not observed in all animals (see Supplementary