2.1.

Single Plane Illumination Microscopy Setup

The SPIM setup^18,14 employs an argon−krypton laser (Melles Griot), expanded $1:4$ and focused into the back aperture of an illumination objective (Olympus, UMPLFLN 10XW, $10\times$, $\mathrm{NA}=0.3$, $\text{working distance}=3.5\mathrm{mm}$). The fluorescence emission is selected by a bandpass filter and is detected by an EMCCD detector (Cascade II, Photometrics) through an Olympus MPLFLN objective ($20\times$, $\mathrm{NA}=0.5$, $\text{working distance}=3.5\mathrm{mm}$; see Fig. 1 for details). The maximum acquisition rate was $4\mathrm{ms}/\text{frame}$ on reduced regions of $40\times 250\mu {\mathrm{m}}^2$.

Fig. 1

(a) Photo and (b) sketch of the SPIM setup. The source is an argon laser (514-nm emission line). A beam expander (1, focal lengths 25 and 50 mm) expands the incoming laser beam. A variable rectangular slit (2) creates a thin vertical aperture in the optical path. The cylindrical lens (3) is positioned at the focal length (50 mm) from the slit. The beam reducer (4, focal lengths 50 and 25 mm) reduces the incoming laser beam to fill the objective back aperture. The excitation objective (5) is an Olympus, water dipping, $10\times$, 0.3 NA. (6) The positioning system (for sample rotation and 3-D translation) and the sample holder. The detection objective (7) is an Olympus, water dipping, $20\times$, 0.5 NA. The tube lens (8) has focal length 150 mm. The detector (9) is an EMCCD detector (Cascade II—photometrics).

2.2.

Zebrafish Embryos

The measurements were done on embryos of zebrafish (3 days postfertilization, d.p.f.). We used the transgenic line ${\mathrm{mitfa}}^{\mathrm{w}2/\mathrm{w}2}$; ${\mathrm{roy}}^{\mathrm{a}9/\mathrm{a}9}$; $\mathrm{Tg}{(\mathrm{kdrl}:\mathrm{EGFP})}^{\mathrm{S}843}$; $\mathrm{Tg}\text{-}{(\mathrm{gata}1:\mathrm{dsRed})}^{\mathrm{sd}2}$, carrying green-labeled epithelium (EGFP) and red-labeled (dsRed) red blood cells.¹⁹ The zebrafish embryos were anesthetized²⁰ with tricaine [$40\mathrm{mg}/\mathrm{L}$ tricaine(ethyl 3-aminobenzene methansulfonate, Sigma–Aldrich Corporation, St. Louis, Missouri)], positioned in a 2-mm diameter fluorinated ethylene propylene tube (FT2X3, Adtech Polymer Engineering, Frampton Mansell, UK), and then immersed in a water cell.^14,21

2.3.

Spatiotemporal Image Correlation Spectroscopy Images Analysis

The STICS images were computed on a time stack of $M$ images $S={\{I_j\}}_{j=0,\dots ,M}$ taken at $M$ times spaced by $\delta t$ ranging typically from 13 to 19 ms (corresponding to 74 to $52\text{frames}/\mathrm{s}$ of the EMCCD, the fastest rates compatible with the chosen field of view). We take the SPIM plane as the $x-y$ plane. We assume that the flow velocity lies in this plane and is $\overrightarrow{v}=(v_x,v_y)$. Each image has a size of $X\times Y\text{pixels}$. The image autocorrelation function was computed on the fluctuation matrices obtained by subtracting pixel by pixel the average signal ${⟨I_k(x,y)⟩}_k=\sum _{k=1}^M{I}_k(x,y)/M$ to the content of each image of the stack, $\delta I_j(x,y)=I_j(x,y)-{⟨I_k(x,y)⟩}_k$, by the FFT algorithm as follows:

Each image of the time stack was padded with zeros up to double its size before performing the FFT.

On each correlation image, corresponding to the lag time $\tau$, we automatically selected the position $({\xi }_{\max },{\eta }_{\max })$ of the maximum by fitting the image to a 2-D Gaussian profile. The amplitude of the STICS function was derived from the amplitude of the Gaussian fit function, subtracted for the possible background.

Even though the morphological images are available, we do not want to presume any a priori knowledge of the flow direction. We have then simply computed the STICS image on rectangular regions of interest (ROIs) aligned with the axes of the image as it is acquired. To evaluate the components of the flow velocity in the SPIM plane, we follow the position of the STICS image with the lag time over at least 5 to 10 lag time points, up to a lag time ${\tau }_{\mathrm{fit}}$. We have checked that the condition to recover the STICS amplitude is that the ROI size is at least $S_{\mathrm{ROI}}\cong 2v_{xy}{\tau }_{\mathrm{fit}}$. A flow of the order $400\mu \mathrm{m}/\mathrm{s}$ and ${\tau }_{\mathrm{fit}}\cong 40\mathrm{ms}$ corresponds to a size of $16\mu \mathrm{m}$, or at least 20 pixels. This is our spatial limitation in the reconstruction of the field flow.

3.1.

Spatiotemporal Image Correlation Spectroscopy Allows to Follow the Flow Along Curved Vessels

We want to devise methods to map the blood flow in tiny branched capillaries. The flow should be characterized in terms of direction and amplitude of the velocity, and we will need to evaluate the flow branching ratio at the ramifications along the network. The rationale of the methods that we are describing is to autocorrelate and cross correlate in space and time stacks of images acquired through a wide field SPIM. In general, we could also map the velocity field by tracking, and then averaging, many trajectories of individual red blood cells. This is not always feasible, for example, when the vessel morphology cannot already be singled out from the images (initial angiogenesis) or the heart pulsation affects the flow, such as in arteries. STICS offers the advantage of directly providing an average flow map in which the heart pulsation effect is averaged out. With STICS, we obtain from a time stack of images a series of autocorrelation maps [Eq. (1)], $G(\xi ,\eta ;\tau )$, function of the lag time $\tau$. For a uniform drift, the map assumes the shape of a 2-D quasi-Gaussian^13,14 function whose maximum occurs at the position ${\overrightarrow{r}}_{\max }=({\xi }_{\max },{\eta }_{\max })$, measured with respect to the origin of the autocorrelation map and it is a function of $\tau$. Here, we want to exploit the possibility of following both the position and amplitude of the autocorrelation map maximum as a measure of the flow velocity and relative amount of cells taking either of the arms of a branching point. We will then extend this concept to cross correlation of different ROIs in the field of view.

We first assess the problem to reconstruct the flow along a curved vessel independent of possible ramifications. We assume that flow is occurring mainly in the SPIM plane and analyze the time stacks of images according to the STICS algorithm on an ROI (rectangular in shape) chosen in such a way to encompass a small ($40\mu \mathrm{m}$) tract of the vessel. The particular anatomy of the zebrafish embryo justifies our assumption. We took care to orientate the embryo in such a way to image in the same light sheet ($<12\mu \mathrm{m}$) both major vessels of the embryo. It is known (see, for example, the repository of image²²) from the anatomy of the embryo that most of the branching occurs in the same sagittal plane.

We first evaluate, at different positions along the putative vessel axis of the STICS image at increasing lag times and fix the lag time at the maximum value, ${\tau }_{\mathrm{fit}}$, for which the maximum ${\overrightarrow{r}}_{\max }=({\xi }_{\max },{\eta }_{\max })$ of the STICS image lies within the STICS image by at least 5 to 6 pixels. The average velocity of the cells is obtained as $\overrightarrow{v}=(-{\xi }_{\max }/{\tau }_{\mathrm{fit}},-{\eta }_{\max }/{\tau }_{\mathrm{fit}})$.^14,23,24 The speed is computed, by assuming nonconstrained motion, as $d_{\max }=\sqrt{{\xi }_{\max }^2+{\eta }_{\max }^2}$.

We need to find a compromise between the finesse of the spatial reconstruction of the flow, the computation time, and the spatial resolution of the STICS image. The latter depends on the brightness of the tracers. In our hands, for red blood cells in living embryos, the minimum displacement of the STICS maximum for which the corresponding velocity uncertainty is smaller than 15% is $|{\overrightarrow{r}}_{\max }|\cong 4\text{pixels}$ (that is $3.2\mu \mathrm{m}$ in the sample plane). The lower the signal/noise is, the larger the minimum value of $|{\overrightarrow{r}}_{\max }|$ is and the coarser the reconstructed flow map is.

Once the ROI size and the fitting lag time have been chosen, the algorithm computes sequentially a set of STICS maps, updating at each step the investigation of ROI with a rigid translation along the vector $\overrightarrow{T}=\alpha ({\xi }_{\max },{\eta }_{\max })$. The choice of small propagation factor $\alpha <1$ allows oversampling of the flow map. We typically adopt $\alpha =0.5$. A flow line can be built by assigning each flow velocity vector to the center of the corresponding ROI and by representing it with the Quiver python class.²⁵ As can be seen from Fig. 2(a), we can follow the blood flow on a wide field of view ($750\times 750\mu {\mathrm{m}}^2$) with a good resolution, $\cong 5$ to $10\mu \mathrm{m}$ [Figs. 2(c) and 2(d)].

Fig. 2

Evaluation of the red blood cells flow in zebrafish embryo (3 dpf) vessels based on STICS maps from ROIs of $(14\times 14)\text{pixels}$ ($\text{pixel size}=2.9\mu \mathrm{m}$). All the flow analyses have been done with propagation factor $\alpha =0.5$. (a) A wide field picture ($256\times 256\text{pixels}=742\times 742\mu {\mathrm{m}}^2$) of blood circulation in the abdomen of the embryo. The flow field is reported in four linear tracts of the vessel network as thin arrows whose length is proportional to the speed. At two branching points, the dominant and marginal flow components are reported as green and red filled arrows, respectively. (b) The STICS map for the vessel branching point marked by a dotted white box in (a), at lag time $\tau =311\mathrm{ms}$. (c and d) Zoomed images on one intersomitic vessel in which the trajectories computed from nine different starting points have been reported together with the average flow curve (see text). (e) The flow speed along the reconstructed average flow (green) line and the adjacent (blue) line reported in (c).

When approaching a branch in the vessel network, the STICS image assumes a double component character [see Fig. 2(b)]. At each step along the path (i.e., on each ROI), we evaluate the relative amplitude of the two peaks and displace the ROI along the direction that corresponds to the most relevant component (according to the maximum of each component detected in the STICS image). In this way, starting from two different points along the vessel network, we are able to evaluate the flow velocity even at the branching points [Fig. 2(a)] and follow it along the individual arms [Figs. 2(c) and 2(d)]. We observed a dependence of the resulting trajectory on the evolution parameter $\alpha$. We have found that for values $\alpha \cong 0.5$ of the evolution parameter, the flow trajectories closely follow the vessel axis, which can be measured from the morphological image [Fig. 2(a)]. There is still a dependence on the position of the starting point (the center of the first ROI of the sequence) along the putative cross section of the vessel, which is shown in Figs. 2(c) and 2(d) (different discontinuous colored lines). However, the trajectories coalesce on a common path that we have selected by averaging the points that correspond to the position of the individual flow lines on each cross section of the vessel (at steps of $20\mu \mathrm{m}$). The resulting line, reported as a thick semitransparent green line in Figs. 2(b) and 2(c), can be considered the common flow path in the selected vessel. There is no appreciable dependence of the amplitude of the flow speed on the starting point of the trajectory, which can be judged from the plot of the flow speed along two flow lines as shown in Fig. 2(e). This is to be expected since, here, we are tracking red blood cells that are approximately the size of the capillary diameter.

3.2.

Spatiotemporal Image Correlation Spectroscopy Allows Estimation of Flow Branching Ratios at the Vessel Ramifications

The application of the STICS algorithm on a ROI that encompasses a ramification of the vessel network, such as in the black rectangle of Fig. 3(a), results in correlation images, in which two components can be singled out [Fig. 3(b)]. These components can be ascribed to the cells that take either of the two branches. The rapidity of their displacement from the center of the correlation image, ${\overrightarrow{r}}_{\max }^{(1)}=[{\xi }_{\max }^{(1)},{\eta }_{\max }^{(1)}]$ and ${\overrightarrow{r}}_{\max }^{(2)}=[{\xi }_{\max }^{(2)},{\eta }_{\max }^{(2)}]$, as a function of the lag times is the measure of the flow direction and the flow speed for cells flowing in each arm. The direction of the velocity vector follows very well the morphology of the vessel ramification, which can be judged from the superposition of the correlation images with the morphological images [Fig. 3(b)]. To this purpose, the correlation images shown in Fig. 3(b) have been inverted by a central symmetry with respect to the center of the correlation space.

Fig. 3

Analysis of the flow ramification at branching points. (a) Morphological image taken on a zebrafish embryo (3 dpf) obtained as the standard deviation projection of the time stack of images used for the STICS analysis. (b) STICS maps computed on the ROI selected in (a) (black rectangle) for increasing lag time, from 164 to 348 ms from left to right and from top to bottom (flow speed coded by a “jet” LUT; morphological images superimposed in light blue). (c) Linear fit of position of the maximum of the STICS images as a function of the lag time. Open (blue) symbols refer to the prevalent component in each STICS image. Open (green) symbols refer to the marginal component (visible for lag times $\tau \ge 170\mathrm{ms}$). The linear fit to the data (dashed lines, same color code) gives velocities $v_1=55\pm 1.5\mu \mathrm{m}/\mathrm{s}$ and $v_2=69\pm 3\mu \mathrm{m}/\mathrm{s}$ for the prevalent and the marginal components, respectively. (d) The Gaussian fitting of the STICS map ($\text{lag time}=348\mathrm{ms}$) together with the profile plots along the $\xi$ and $\eta$ axes.

In the correlation maps [Fig. 3(b)], one peak was found to always be larger in magnitude irrespective of the lag time. Moreover, the standard deviation projection image (each pixel contains the standard deviation value found through the time stack) showed very similar signals on the two arms, which indicates that the emission of the individual red blood cells passing through the two arms is very similar. We assume, therefore, that the maximum amplitude of the two correlation components is proportional to the number of red blood cells flowing in one or the other ramification of the vessel branching. For example, as shown in Fig. 3(d), the brightest component of the STICS map corresponds to the red blood cells flowing upward, amounting to $66\%\pm 1\%$ of the overall amplitude. The marginal component represents cells flowing from left to right, amounting to $34\%\pm 1\%$ of the overall amplitude. The marginal component can be singled out with a reasonable signal/noise ratio (10%) only for lag times $\tau \ge 170\mathrm{ms}$. Two values of the flow speed can be computed from the linear fit of the displacements of the maxima of each component as a function of the lag time, as shown in Fig. 3(c): we measure $v_1=55\pm 1.5\mu \mathrm{m}/\mathrm{s}$ and $v_2=69\pm 3\mu \mathrm{m}/\mathrm{s}$ for the prevalent and the marginal components, respectively. We find then that the marginal component is slightly faster than the prevalent one. This corresponds quite well to the observation that the left–right vessel arm is visibly larger [see Fig. 3(a)] than the bottom-up arm and the fact that the vessels are approximately the size of the red blood cells diameter. It must be noticed that the algorithm for the computation of the flow velocity used here is customized to the case of right-angle crossing between vessels as in Fig. 3(a) (box). Since the flow is constrained to occur along the vessel axis, the displacement of the maximum of the autocorrelation map may, in general, not be along a line passing through the center of the map. In this case, we should compute the flow speed along the vessel axis through the incremental algorithm $|{\overrightarrow{r}}_{\max }(\tau +\delta \tau )-{\overrightarrow{r}}_{\max }(\tau )|=|\overrightarrow{v}|\delta \tau$.

3.3.

Spatiotemporal Cross-Correlation Analysis Enhance Mapping Features

The possibility of spatially mapping the flow relies on the choice of a small ROI on which we apply the spatiotemporal correlation algorithm. The space resolution of the reconstructed flow depends on the size of the ROI (half of its $\text{linear size}=0.5\times S_{\mathrm{ROI}}$). However, a second key parameter for obtaining high-precision flow speed estimations is the max tracer speed (K) inside the investigation ROI. If we assume that at least five lag time points are needed to linearly fit the relation between the displacement $|{\overrightarrow{r}}_{\max }|$ and the lag time, K can be estimated from the ROI size ($S_{\mathrm{ROI}}$) and the CCD frame rate (fps) as $K=S_{\mathrm{ROI}}\times \mathrm{fps}/5$. In the experiments shown in Figs. 1 and 2, $\mathrm{fps}=74\text{frames}/\mathrm{s}$ and $S_{\mathrm{ROI}}=16\mu \mathrm{m}$ and we compute $K\cong 240\mu \mathrm{m}/\mathrm{s}$. If we had larger speed values, we should increase the ROI size, losing then in spatial resolution of the flow mapping. Moreover, if the flow is not aligned with the ROI axes, we would need to increase the ROI size in both the $\xi$ and $\eta$ directions. This would increase substantially the computation time.

To solve these issues, we propose to adopt an approach consisting of two-stage analysis:

Fig. 4

Cross-correlation spatiotemporal analysis. (a) The time average of the stack of images. The light and dark gray rectangular boxes are the original and translated ROI, respectively. The white arrow indicates the translation vector. (b) Three autocorrelation images computed on the original ROI for lag times $\tau =18.2$, 36.4, and 54.6 ms. The white arrow indicates in each STICS image the displacement vector ${\overrightarrow{r}}_{\max }=({\xi }_{\max },{\eta }_{\max })$. (c) The spatiotemporal cross-correlation maps between the original and translated ROIs for absolute lag times $\tau =18.2$, 54.6, and 91 ms. The light green arrow indicates in each cross-correlation image the relative displacement vector ${\overrightarrow{r}}_{\max }^{\prime }=({\xi }_{\max }^{\prime },{\eta }_{\max }^{\prime })$. (d) The plot of the displacement $d_{\max }=\sqrt{{\xi }_{\max }^2+{\eta }_{\max }^2}$ as a function of the absolute lag time $\tau$ (blue color) and $d_{\max }^{\prime }=\sqrt{{({\xi }_{\max }^{\prime })}^2+{({\eta }_{\max }^{\prime })}^2}$ as a function of the reduced lag time, $\tau -{\tau }_0$ (green color, ${\tau }_0=18\mathrm{ms}$). The dashed line is the best linear fit to the data and corresponds to the speed $v=380\pm 9\mu \mathrm{m}/\mathrm{s}$.

In this way, we can increase the dynamic range of the speeds that can be detected. A flow speed almost 4 times larger than in Fig. 3 can be measured on the vessel reported in Fig. 4(a) with 3% accuracy at $20\text{-}\mu \mathrm{m}$ size resolution ($S_{\mathrm{ROI}}=20\mu \mathrm{m}$). At the same time, we are able to keep the computational costs at the same level as the original autocorrelation algorithm. The details of the algorithm can be gained by comparing Fig. 4(b) (STICS correlation) with Fig. 4(c) (STICCS, cross correlation). We evaluate the flow velocity on the cross-correlation images by translating the correlation frame of reference $(\xi ,\eta )$ at the position $({\xi }_0,{\eta }_0)$ at which the cross-correlation peak is first detected and the lag time $\tau$ axis to the corresponding lag time, ${\tau }_0$ [see Fig. 4(c)]. The two components of the flow velocity can then be evaluated as the inverse of slope of $\xi -{\xi }_0$ and $\eta -{\eta }_0$ [the corresponding vector is reported as a white arrow in Fig. 4(c)] as a function of $\tau -{\tau }_0$. This estimate is in very close agreement with the one obtained on the autocorrelation images, which can be judged from Fig. 4(d) (open and filled symbols).

We notice that the analysis presented here is limited to 2-D flows. We could extend it to 3-D flows by following, in addition to the position of the correlation map in the correlation space, the decay with the lag time of its amplitude. The main limitation in this case could be the computation time, more than the acquisition time. In fact, we would work out the information about the out-of-plane motion from the amplitude decay of the cross correlation computed on the time stack of images of a single light sheet. This possibility is currently under investigation in our laboratory.