Cardiac pulsatility mapping and vessel type identification using laser speckle contrast imaging

: Systemic ﬂow variations caused by the cardiac cycle can play a role or be an important marker in both normal and pathological conditions. The shape, magnitude and propagation speed of the ﬂow pulse reﬂect mechanical properties of the vasculature and are known to vary signiﬁcantly with vascular diseases. Most conventional techniques are not capable of imaging cardiac activity in the microcirculation due to spatial and/or temporal resolution limitations and instead make inferences about propagation speed by making measurements at two points along an artery. Here, we apply laser speckle contrast imaging to images with high spatial resolution in the high frequency harmonics of cardiac activity in the cerebral cortex of a mouse. We reveal vessel dependent variation in the cardiac pulse activity and use this information to automatically identify arteries and veins.


Introduction
Various homeostatic processes in the brain and other organs depend on blood flow regulation.Although the mean flow and pressure are the most important parameters regulating blood flow, systematic variations caused by the cardiac cycle, vasomotion or other rhythmic activities can also play a significant role both in normal and pathological conditions [1,2].In the brain, the largest flow variation comes from the cardiac activity [3], which also governs dynamics of the intracranial pressure (ICP).It has been recognized that cardiac pulsatility in pressure and flow can change with disease, for example in dementia [4,5]; liver cirrhosis [6]; age-related macular degeneration, in which intraocular pulsatility increases with disease severity [7]; and in peripheral vascular disease [8].In the brain, arterial pulsation is considered to be the force in the glymphatic system driving cerebrospinal fluid flow through the interstitial space and playing an important role in waste product removal [9].Pulsatility characteristics, such as frequency, power and shape reflect vascular and overall intracranial compliance and can become a valuable tool in disease assessment and prognosis.
Unlike ICP, which propagates trough the brain at the speed of sound, flow pulsations can vary significantly depending on the location, thus requiring an imaging approach to permit the spatial analysis [3].Traditionally, both in clinical and preclinical studies, either transcranial Doppler ultrasound (TDU) or magnetic resonance imaging (MRI) is used to capture flow variations [10,11].Full-field imaging of microcirculation pulsatility has proven to be difficult due to the high requirement for both spatial and temporal resolution.This is particularly challenging in mouse animal models where cardiac activity can approach 10 beats per second.With recent developments in CMOS sensors, laser speckle contrast imaging (LSCI) provides a unique opportunity to monitor flow dynamics at frame-rates exceeding 100 Hz with high spatial resolution, making it possible to map high frequency pulsatility of the blood flow.LSCI provides a rapid wide field characterization of the motion of light scattering particles that allows monitoring of changes in blood flow [12].Over the past fifteen years it has become one of the major techniques applied for cerebral blood flow imaging [12] as well as for imaging blood flow in the skin [13], retina [14], kidney [1] and mesentery [2].Recent approval by the Food and Drug Administration of an LSCI based device "Laser speckle flowgraphy LSFG-NAVI" for human retinal imaging, has resulted in a significant rise in clinical applications [15,16].While most studies employ LSCI to monitor response-invoked changes in blood flow, it was also applied to analyze rhythmical flow activity caused by vasomotion, tubular-glomerular feedback and myogenic response [1,2,14,17].
In the present study we apply conventional LSCI to analyze and map cardiac pulsatility in the cerebral cortex of the mouse.We show that LSCI provides the temporal resolution and signal quality sufficient to capture local differences in the shape and timing of the pulse wave and that observed differences can be used to automatically separate veins and arteries.

Methods
Surgical Procedures.All animal procedures were approved by the Boston University Institutional Animal Care and Use Committee and were conducted following the Guide for the Care and Use of Laboratory Animals.
For this study, to test the performance of the vessel classification algorithm in a variety of conditions, mice of different strains, ages and cranial window preparations were used.All surgeries were performed when mice were approximately 12 weeks old.In one animal (ID: A) (C57BL6, C57BL/6J-Tg(Thy1-GCaMP6s)GP4.3Dkim/J,The Jackson Laboratory), a glass plug (made by fusing three 3-mm coverslips with one 5-mm coverslip) was inserted into the craniotomy, gently touching the brain.The plug was fixed with dental acrylic, along with an aluminium bar attached to the skull for head fixation during imaging (For details, see [18]).In two animals (ID: B and C) (B6C3, The Jackson Laboratory), the same surgical procedure was applied, but the animals were imaged approximately 6 months after the surgery.In the last animal (ID: D) (C57BL6, C57BL/6J-Tg(Thy1-GCaMP6s)GP4.3Dkim/J,The Jackson Laboratory), the skull over the right hemisphere was removed with craniotomy and the brain surface was then covered with a half-skull-shaped curved glass (Crystal Skull, LabMaker, Germany) [19].The plug was fixed with dental acrylic, along with an aluminium bar attached to the skull for head fixation during imaging.The animal was imaged approximately 6 weeks after surgery.In all animals imaging was done under isoflurane anesthesia (3% induction, 1-1.5% maintenance, in 1L/min oxygen) and fixing the head to a custom-made stage.Imaging lasted for approximately 20 minutes, then the animals were allowed to recover after the experiment and returned to their cages.No additional surgery was performed for the imaging experiment.Heart rate and oxygen saturation was non-invasively monitored (Mouse Stat Jr, Kent Scientific) and all noted measurements were within the expected physiological range.
Laser speckle imaging A CMOS camera (Basler acA2040-180kmNIR, 2048x2048 pixels, 4.4 µm pixel size) was used to record the backscattered light collected by a 5x objective with NA= 0.14 (Mitutoyo, Japan).A polarizer was placed in front of the objective.A subset of pixels was used -896x896 for animals ID:A-C and 1400x896 for animal ID:D.Recordings lasting 300 seconds were acquired with a 5ms exposure time and at a frame-rate of 114 frames per second.Coherent light was delivered to the object using a volume holographic grating (VHG) stabilized laser diode (785nm, LP785-SAV50, Tholrabs) with a power density of approximately 20mW/cm 2 , providing an optimal signal to noise ratio [20,21].The laser was operated with the recommended settings.

Visual identification of the cranial vasculature
The cranial windows were imaged under white light using a surgical microscope (Nikon SMZ-745T) and an iPhone 7 camera with a LabCam adapter (iDu Optics) immediately before the laser speckle imaging.Using these images as reference, a two researchers experienced in the cranial vascular anatomy of mice classified pial vessels as artery or vein, depending on visual cues like vascular orientation, diameter, color, and superficiality.

Data analysis
The recorded data were processed to get local spatial contrast [12].We preferred spatial over temporal contrast since it provides the high temporal resolution required to detect the cardiac pulsatility.To get the flow dynamics we have applied the commonly used model of the blood flow index BF I: where K is contrast, σ s and I s are standard deviation and mean of the pixels in a 7x7 neighbourhood.
The BFI was calculated for each pixel at each moment of time and then analyzed using a fast Fourier transform (FFT).The power spectrum for each pixel was calculated over subsets of 256 frames (2.24 seconds).The frequency of pulsatility harmonics were identified by averaging spectra over the whole field of view and identifying the peaks within the expected frequency bands.The phase at the first harmonic frequency was used to estimate pulse timing for each pixel.For each pixel and each subset of frames, we calculate absolute H and relative R peak-to-noise height as: where P( f i ) is power at frequency f of harmonic number i and P noise is the estimated power of the noise, which is calculated as the average power at non-cardiac frequencies.For all results, unless stated otherwise, instead of averaging power spectra over time, we averaged the H and R values.This allowed us to avoid issues with temporal variability in the cardiac frequency.

Automated identification of cranial vasculature
To automatically identify whether a pixel belongs to a vein or an artery, we have developed a simple algorithm that utilizes H and R data.The algorithm is based on the fact that the flow pulsation decreases from arteries to veins and is intermediate in the micro-vessels in between in the parenchyma.The main steps of the algorithm are described in Table 1.The algorithm was implemented in Matlab and is available for download [22].
Negative T corresponds to veins, positive to arteries.The closer |T| value to 0the higher uncertainty is.

Results and discussion
An example of BFI map with highlighted regions of interest belonging to the vein, artery and parenchyma is shown in Fig 1, A. From the corresponding power specters (Fig 1, B) one can see that up to 5 harmonics of cardiac activity can be clearly distinguished in arterial flow.The magnitude of slow flow oscillations in the artery and vein is much higher compared to parenchyma leading to higher values of absolute peak-to-noise height: H arter y = 20, H vein = 15 compared to H par enchyma = 2.5.This is, however, not true for higher order harmonics, which are almost absent in the vein, but are still evident in the parenchyma and the artery.Calculating relative and absolute peak-to-noise heights (see Eq. ( 2),( 3)) for each pixel allows us to plot activity maps for different harmonics.In • In arteries, pulsatility is governed only by the forward propagating pressure wave of the cardiac pulse.
• Oscillations in parenchyma and small vessels are affected by two factors: (i) the forward propagating pulse pressure wave changing the flow speed; and (ii) motion of the static tissue caused by pressure wave propagation through the tissue, particularly by that originating from the circle of Willis.
• The pulse in the vein is formed by: (i) the forward propagating pulse pressure wave, (ii) surrounding tissue motion and (iii) a backward propagating ("pulling-back") cardiac wave [23].
The fact that the parenchyma delay is independent of location and no gradient is observed from artery to vein (see Fig 2(E),(F)), indicates that bulk tissue motion is the dominating source of oscillations in parenchyma.Similar tissue motion has been observed in optical coherence tomography studies [24].Variation in the propagation of different harmonics through different tissues and vessels has been noticed previously and can potentially be used as a marker in pathologies related to changes in arterial stiffness or intracranial compliance [3].In the artery, H is higher compared to the vein and parenchyma.H in the vein is higher than in parenchyma for the 1st harmonic but not for the second.C,D -relative peak-to-noise height.The vein has reduced R compared to artery and parenchyma for the 1st harmonic, which becomes even more evident for the second harmonic.E -delay of the pulse occurrence compared to the artery.Reference region of the artery is marked as a red dot.
The delay is about 14ms for parenchyma and 20ms for the vein.F -relative change of BF I calculated by averaging multiple pulses.Solid red and blue lines correspond to artery and vein respectively.Dashed lines correspond to pixels in parenchyma: blue -along the large vein, red -along the artery, green -far from major vessels.While the difference between artery, vein and parenchyma is evident -there is no observable difference in parenchyma dynamics or delay in different regions.
While cardiac pulsatility mapping can be used as a marker of vascular stiffness in pathological studies -one the most obvious and straightforward applications is to use this information to separate veins and arteries.Incorrect identification of vessel type can lead to venous sacrifice which in turn can led to serious complications during surgeries [25][26][27], and thus, any additional information is valuable that helps a surgeon to make the right decision and possible early intervention.This becomes particularly important in pathologies which alter the vascular anatomy, making it difficult for a specialist to recognize vessel type [28].In some cases, techniques like fluorescent videoangiography, high-resolution 3D and 4D magnetic resonance angiography and contrast-enhanced ultrasound angiography have been used to identify the vasculature [28,29].Some of these techniques cannot provide real-time intraoperative imaging or require injection of agents and/or use of sterile imaging probes that need to be positioned on the tissue many times to get consecutive images.In other cases, particularly for processing large volumes of vessel specific response data, automated vessel type identification is required.Pulsatility characteristics can become a basis for such identification or can be combined with known white light characteristics [30] to reduce uncertainty in the identification.
Oscillations are attenuated while travelling from arterials through the vascular network, with oscillations in the veins mainly originating from the pulsation of the surrounding tissue and/or a backward propagating cardiac pulse wave.One can expect that this will result in significantly different characteristics of pulsatility in those types of vessels.In Fig. 2, the difference between venous and arterial pulsatility is evident.We have calculated the mean and standard deviation of H and R for ROIs belonging to artery, vein and parenchyma (see Fig. 3 ).It is evident that these vessel types can be distinguished based on the difference significance with less than 100 seconds of data acquisition.Based on the observations shown in Fig. 3, we have devised an algorithm for automated identification of the vessel type (see Table 1).This algorithm was applied to LSCI data recorded from mice of different strains, ages and cranial window preparations (n=4) and compared to the visual identification of the vessel type performed by a trained specialist (see Fig. 4).In animals ID:B,C,D, large veins were identified following step 3 and in animal ID:A, large arteries were identified at step 4 of the algorithm.Vessels identified at these steps do not require other vessels as a reference since they can be identified by comparison with the surrounding tissue.Vessels which type was not identified during steps 3 and 4 (arteries in ID:B,C,D, veins in ID:A as well as smaller vessels of both types) were identified with less confidence (steps 5 and 6) using the vessels identified at steps 3 and 4 as a reference.Concurrence between manual and automated segmentation was observed for all vessels.Fig. 4. Application of automated identification of the vessel type on a per pixel basis.Top row -algorithm output.Dark red or dark blue corresponds to high confidence in the vessel type identification (see Table 1 steps 3-4), lighter colors -to the pixels which were identified with less confidence (Table 1 steps 5-6).Bottom row -white light images with major vessels types identified by a specialist.

Conclusion
We have shown that LSCI can be used for per pixel mapping of cardiac pulsatility, revealing 5th order harmonics.Power spectrum characteristics are different in arteries, veins and parenchyma, with arteries generally having a strongly pronounced first harmonic peak and veins having the second harmonic peak significantly less pronounced compared to arteries or the parenchyma.Interestingly, the relative peak height R for the first harmonic is not significantly different between veins and parenchyma, and the second harmonic is not significantly different for arteries and parenchyma.This indicates that the shape of the cardiac signal, as well as the amplitude, varies significantly depending on vessel type.
Phase analysis shows that oscillations in the parenchyma and veins are delayed compared to arteries.This can be explained by the origin of the oscillations.The spatial uniformity of the delay in the parenchyma suggests that the oscillations are be dominated by tissue motion, rather than flow motion.
We have demonstrated that one can use differences in cardiac activity to identify arteries and veins.The cardiac activity mapping by itself provides information that can help in manual vessel type identification.We further explored the idea by suggesting an algorithm for automated or semi-automated vessel classification on a per pixel basis.In 4 tested animals, the algorithm provided 100% concurrence with vessel type identified by a specialist.An acquisition time of less than 100 seconds was sufficient to obtain a significant difference between veins and arteries.Unlike a previously published suggestion to identify vessel type based on the temporal minimum of intensity [31], our approach is not limiting the light source wavelength to 600-640 nm and can be applied in any setting that uses LSCI.For large veins, however, one can expect increased pulsatility from the backward propagating cardiac wave.To properly identify such vessels as veins it would be advised to incorporate use of other parameters, such as contrast, vessel diameter or pulse delay.
The overall approach can be directly applied to vascular networks in other organs with minimal changes in the vessel type identification algorithm.Particularly it might become useful as an intraoperational guide for surgeons or as an additional metric in emerging clinical studies that use laser speckle flowgraphy.Application of a deep learning approach or advanced masking algorithms [30] can further decrease uncertainty and improve segmentation results.

Fig. 1 .
Fig. 1.A -An example of an average BFI map, animal ID: A. Regions of interest (ROIs) over vein, artery and parenchyma are highlighted with black, red and green polygons respectively.Each ROI has 2000 pixels.B -Fourier power spectrum of BFI in the ROIs, averaged over 36 seconds.One can see up to 6 harmonics in the artery, but not in the vein which has a comparable signal magnitude.
Fig 2(A)-(D), example H and R maps are shown for the first and second harmonics.One can notice that not only absolute and relative peak heights are different depending on the vessel type, but their behaviour also depends on the harmonic order.The phase delay map shown in Fig 2(E), shows that the pulse arrives to the parenchyma and vein by 14 and 20 ms later compared to the artery.The observed differences are possibly driven by the following mechanisms:

Fig. 2 .
Fig.2.Mapping of cardiac activity for 1st and 2nd harmonics.A,B -absolute peak-to-noise height.In the artery, H is higher compared to the vein and parenchyma.H in the vein is higher than in parenchyma for the 1st harmonic but not for the second.C,D -relative peak-to-noise height.The vein has reduced R compared to artery and parenchyma for the 1st harmonic, which becomes even more evident for the second harmonic.E -delay of the pulse occurrence compared to the artery.Reference region of the artery is marked as a red dot.The delay is about 14ms for parenchyma and 20ms for the vein.F -relative change of BF I calculated by averaging multiple pulses.Solid red and blue lines correspond to artery and vein respectively.Dashed lines correspond to pixels in parenchyma: blue -along the large vein, red -along the artery, green -far from major vessels.While the difference between artery, vein and parenchyma is evident -there is no observable difference in parenchyma dynamics or delay in different regions.

Fig. 3 .
Fig. 3. Mean and standard deviation of H and R for ROIs shown in Fig. 1(A) depending on the signal acquisition time.A,C -first harmonic, B,D -second harmonic.All metrics show significant differences with p < 0.01 starting from the shortest acquisition time, thus it is not shown in the plot.Difference increases further for longer acquisition times.

Table 1 .
Algorithm description.* Step 4 is performed in the same loop as step 3. ** Steps 5-6 should only be performed when presence of both types of vessels is assumed (or identified as a result of steps 3-4).*** Steps 7-8 purpose is to increase robustness of the method.(x, y) notation in the third column represents a matrix of pixels a and R v parameters to 90th and 10th percentile of non-parenchyma R respectively Input: R(x, y), T(x, y), Output: R a , R v 6** Loop over non-parenchyma pixels that were not classified.Set T