Fully automated digital holographic processing for monitoring the dynamics of a vesicle suspension under shear flow

: We investigate the dynamics of a vesicle suspension under shear flow between plates using DHM with a spatially reduced coherent source. Holograms are grabbed at a frequency of 24 frames/sec. The distribution of the vesicle suspension is obtained after numerical processing of the digital holograms sequence resulting in a 4D distribution. Obtaining this distribution is not straightforward and requires special processing to automate the analysis. We present an original method that fully automates the analysis and provides distributions that are further analyzed to extract physical properties of the fluid. Details of the numerical implementation, as well as sample experimental results are presented.


Introduction
Optical microscopy is limited by the small depths of focus due to high numerical apertures and the high magnification ratios. The extension of the depth of focus has become an important goal in optical microscopy. In this way, Digital Holographic Microscopy (DHM) provides a natural and elegant solution to this problem by giving the capability to reconstruct numerically different slices inside an experimental volume. On the basis of the knowledge of the complex amplitude in one plane, one can easily propagate this complex amplitude to different planes along the optical axis by implementing numerically the propagation equations of Kirchhoff-Fresnel. Furthermore, DHM provides a quantitative optical phase which is particularly useful for pure phase objects (like living cells). DHM has been successfully used for different applications in fluid physics [1,2] but also in biological applications such as living cell culture monitoring, observation of biological samples [3][4][5][6][7], classification of microscopic living organisms [8,9], living cells tracking [10. For the last decade, we have developed DHM working with reduced spatial coherent sources [11]. It has been shown that the reduction of the spatial coherence strongly improves the quality of the holograms [12], especially when the experimental cell under investigation cannot be optimized with respect to the optical quality. DHM can be a powerful measuring instrument but is often problem dependent and requires the implementation of specific holographic information processing techniques.
Channel flows of complex fluids containing deformable particles such as Red Blood Cells (RBCs) are influenced by many phenomena such as hydrodynamic interactions between particles and between channel walls and particles, phase separation at bifurcations, which all depend on the mechanical properties of the particles. While blood flows have been observed for nearly two centuries [13], their understanding is still incomplete. In particular, understanding the concentration inhomogeneities between vessels or across one given vessel, or the interplay between RBCs, White Blood Cells an platelets still requires huge efforts from the community, in terms of experiments or numerical simulations.
Two phenomena are of particular interest in confined shear-flow in the understanding of blood rheology: when sheared close to a wall, deformable objects such as vesicles or RBCs undergo a lift force of viscous origin which pushes them away from walls [14,18,39] while on the other hand hydrodynamic interactions and collisions within the suspension lead to shearinduced diffusion that tends to homogenize the concentration. When these two effects balance, a cell-free layer is formed near walls and plays a decisive role in biological functions (cell adhesion and aggregation, immunitary response, rheology of blood in confined situations [40]) and may be of interest in technological applications such as cell sorting devices, etc… In this paper, we aim at describing an automated protocol to study the 4D behavior of a vesicle suspension in flow. We illustrate the procedure by shearing a suspension between two plates. The bottom plate is fixed while the top one is rotating creating a linear profile of velocity of the fluid between the plates. The thickness of the channel is about 200 µm. The distribution of the vesicle suspension over the thickness of the channel is obtained by monitoring the suspension with a spatially reduced coherence DHM. Depending on the concentration of the samples, the duration of the experiments (between 20 and 100 sec at 24 frames/sec), the number of experiments needed to get reliable statistical data, the different experimental parameters to be investigated etc … the number of holograms to be processed becomes rapidly huge and there is a real need to automate the processing. This automation is not straightforward as it requires detecting the tri-dimensional position of each vesicle present in the field of view but also specific features such as its size and shape. The proposed technique combines several algorithms in an original way. Because vesicles are pure phase objects, they become invisible when they are in the focus plane. To solve the problem, we take benefit of the full interferometric information by working on the intensity image and especially on the phase map provided quantitatively by DHM.
The optical system is briefly described in section 2. The experimental protocol is detailed in section 3. Section 4 reports the processing sequence of the holograms. Illustrative examples are given in section 5 and concluding remarks in section 6.

Optical setup description
This section describes the DHM setup used for the experiments. The optical setup is described in Fig. 1. It is a Mach-Zehnder configuration with microscope lenses (x20) to provide the magnification. The interferometer is using a reduced spatial coherent source provided by a laser source (monomode laser diode, λ = 635 nm) focused by lens L1 on a rotating ground glass (RGG). By adjusting the position of the rotating ground glass on the optical axis, one can easily increase the spatial extension of the source and consequently decrease the spatial coherence. The use of spatial coherence reduces the coherent noise and provides a temporal-like effect that eliminates the coherent multiple reflection issue [12]. Holograms are acquired with a charge-coupled device (CCD) camera. We used a JAI CV-M4 camera with an array of 1280 x 1024 pixels corresponding to an optical field of view of 524 µm x 420 µm. The angle between the reference and the object beam (passing through the sample) is adjusted to work in an off-axis configuration allowing a fast acquisition (24 frames/sec) with a very short exposure time (200 µm) and avoid blurring when using high flow rates.

Objective and experimental protocol
The main objective of the experiment is to get dynamic properties of a vesicle suspension under shear flow. As outlined in the introduction, two physical effects govern the behavior of the suspension: a lift force due to the presence of walls and shear induced diffusion due to interactions between vesicles. At shear rates relevant to practical situations, sedimentation cannot be neglected and significantly affects the results. To remove the screening effects of gravity (the sedimentation speed of vesicles is of the same order as the lift velocity we investigate), experiments were conducted in microgravity (on parabolic flight and in sounding rockets) with a Couette shear flow chamber designed and manufactured by the Swedish Space Corporation (SSC). As shown in Fig. 2, it is made of two parallel glass discs, with a diameter of 10 cm and a gap of 170 µm. Vesicles were prepared by the electroformation method [41] from lipids (1,2 -Dioleoylsn-glycero-3-Phosphocholine (DOPC)) deposited on ITO glass plates and a solution to be encapsulated made of 80% water and 20% glycerol w:w, where 300 mM of sucrose are dissolved. Each electroformation provides a polydisperse sample with vesicles diameter ranging from 5 to 100 µm (median size between 15 and 25 µm), a volume fraction of a few percent and a total volume of 2 ml. Samples were then diluted into a slightly hyperosmotic solution made of 80% water and 20% glycerol w:w, and 350 mM of glucose creating a refractive index change Δn of 0.0061.

Digital processing of numerical holograms
Digital holograms are grabbed and stored on a hard disk for further processing. The huge amount of data stored requires a fully automated processing minimizing any intervention of an operator and getting the maximum reliability of the processing. Automated threedimensional detection is generally problem dependent and there is no miracle solution that solves the automatic detection problem. Different methods have been proposed based on the accumulation of interferometric information during the numerical propagation along the optical axis [8,42,43].
The processing sequence can be summarized as follow: for each hologram, find the position, the size and the shape of each vesicle present in the considered hologram. In the following, we will consider the processing of one hologram, all the holograms of the sequence are processed in the same way and do not depend on their temporal position in the sequence. The very first step is to extract the complex amplitude U(x,y) = A(x,y)exp(iφ(x,y)) (amplitude A(x,y) and phase φ(x,y), i being the imaginary number) from the hologram using the Fourier method [43][44][45]. On the basis of the complex amplitude U(x,y), the task can be divided in two parts: first find the position xyz of each vesicle; secondly extract the features of each vesicle in its corresponding z plane. Those two operations are described in details in the following subsections.

Automatic XYZ detection
Vesicles are pure phase objects becoming transparent in intensity when they are in the focus plane. Intuitively, one can easily understand that in the vicinity of a vesicle, the xy gradient is minimal at the focus plane. This property could be exploited to scan the experimental volume along the z direction and, for each plane, compute the xy gradient of the reconstructed intensity. The xyz detection of the vesicles would then be performed by thresholding the different gradient of intensity in the different reconstructed planes. This approach is intuitively attractive but suffers from two drawbacks. First, it requires propagating the whole hologram along the z axis (covering the whole experimental volume) with a sufficiently small increment Δz to detect the smallest vesicles. This becomes rapidly heavy in computation time. Secondly, the method is very sensitive to noise and consequently provides poor results. The number of false detections is important and the method is not sufficiently reliable. For those two reasons, it has been rapidly discarded.
We propose a novel approach that consists in decoupling the xy and the z detection. The xy detection is performed on the phase image φ(x,y) rather than the intensity one. The in-depth z detection is performed locally within a region of interest (ROI) surrounding each detected vesicle. Indeed, although the image of a vesicle undergoes strong changes in intensity while propagated along the z direction, it keeps a very slowly changing shape on the phase image.
As an example, Fig. 3 compares the intensity and the phase image for one vesicle (diameter of about 20 µm) in three different planes (respectively −70 µm, 46µm -focus plane -and + 100 µm). Excepted for very small vesicles, a simple threshold on the phase image (corresponding to the focus plane of the microscope) detects almost all the vesicles. Nevertheless, due to phase jumps induced by the periodicity of the phase, one cannot work directly on the wrapped phase image. Usually, the phase map is unwrapped to get a continuous phase distribution. Unwrapping might be a tough task and is generally time consuming [46]. In our case, the phase change induced by vesicles is, excepted for very big vesicles (which are not considered in our study), below 2π. The removal of phase jumps can thus be easily made by subtracting (modulo 2π) the background phase [47]. Because of the symmetry between the reference and the object beam of the interferometer, we used a compensation method which consists in fitting the background phase with a bi-dimensional quadratic phase map. Due to phase jumps, the fitting cannot be performed directly on the phase image. For this reason, we perform the fitting on the derivative of the phase map [15,48]. Figure 4 illustrates the compensation method with a phase map (a), the fitted background map (b) and the subtracted phase map (or compensated phase map) (c). Unlike the original phase map, the compensated phase presents an almost uniform background on which vesicles appear as bright patches. This phase map compensation is indispensable as will be highlighted later for the extraction of geometrical features of vesicles. Nevertheless, the compensation method used here is based on a fitting and does not perfectly compensate the background phase. As a consequence, the compensated phase does not present a perfectly flat background allowing a simple threshold. Putting a threshold a little bit higher than the background value does not solve completely the problem as it might swallow up small vesicles. We remove those background fluctuations (which can be of several grey levels) by subtracting the mean compensated phase image (computed on all the compensated phase maps of the sequence under processing -it corresponds to a minimum of 600 images). A simple filtering of the thresholded image removes the small local fluctuations that can be considered as false detections (we remove all the connected area smaller than 200 pixels which corresponds to vesicles having a diameter smaller than 6 µm). This detection method is very simple and presents several advantages. First, the detection is very fast and very efficient (more than 95% of the vesicles are correctly detected). Secondly, the proposed detection method performs a first segmentation of the vesicles which gives a rough determination of their size and shape. This point is very important for the z detection. Figure 5 shows an example of xy detection and pre-segmentation of the vesicles. Note that the smallest vesicles, correctly detected, present very low contrasts demonstrating the efficiency of the proposed method. Once the xy position of each vesicle is known, the z detection is performed locally on a ROI centered on the vesicle. The complex amplitude and the corresponding compensated phase are cropped locally and the z scanning is performed using the numerical propagation inside this ROI. The dimensions of the ROI are determined by the pre-segmentation performed during the xy detection. We used ROIs whose dimensions equal twice the width and the height of the segmented vesicle. In this way, ROIs are selected differently for each vesicle, minimizing the dimensions of the Fourier Transforms needed for the numerical propagations. To avoid border diffractions due to the ROI borders, we extend the ROI size with values that minimize the borders diffraction [49]. In this way, the ROI is extended to the smallest value that fits the requirement of Fast Fourier Transforms (the size must be a power of 2). For example, for a vesicle whose width and height are respectively 80 and 105 pixels, we will crop the complex amplitude with an ROI of 160 x 210 pixels that will then be enlarged to 256 x 256 pixels. Considering enlarged ROI of 256 x 256 pixels, the numerical propagation inside ROI is faster than the whole hologram propagation as long as the number of vesicles present in the hologram is smaller than 88. For a larger number of vesicles, the whole hologram propagation becomes faster. Nevertheless, the highest concentrations studied didn't exceed 50 vesicles per hologram.
Digital Holography allows to numerically investigate a volume by refocusing the different slices but does not provide any information about the focus status. In [50], we have proposed a focus metrics based on the integral of the modulus of the complex amplitude. It has been shown that this metrics is maximal for phase objects when the reconstructed distance reaches the best focus plane. Nevertheless, this criterion is based on an integral and gives a value that is influenced by all the objects distributed in the volume. If, following the discussion of the previous paragraph, another vesicle is present in the ROI, the z value will be influenced by both vesicles and will give a mean value which does not correspond to the z position of any of the two vesicles (unless they are in the same plane). To cope with all those requirements, we use a mask that determines the region where the integral must be calculated. This mask is computed on the basis of the pre-segmented shape provided by the xy detection. We enlarge the pre-segmented shape by a border of 10 pixels (this value is chosen after successive trials). This is made simply by applying 10 times a dilatation operation. Indeed, it has been demonstrated [51] that the flux density of the complex amplitude distribution is proportional to the gradient of the complex amplitude, and consequently, its evolution direction is perpendicular to the object border during the propagation. Figure 6 shows an example of a vesicle with the different ROI considered and the mask defining the domain of integration. To illustrate the efficiency of the processing, Fig. 7 shows an example of a region containing 4 vesicles with close xy position but different z positions. Table 1    Of course, when the distance between two (or more) vesicles is too small, it becomes impossible to consider them separately and the discussion above fails into default. We will see in the next subsection how to remove those occurrences from the output data.
The processing time is directly linked to the increment Δz used when propagating the cropped amplitude along the z axis but determines also the accuracy on the z axis. Decreasing Δz will increase the accuracy of the measured z position but will also increase the processing time. Instead of using a dichotomous approach, we take benefit of a previous theoretical analysis of this metrics [52]. Indeed, it has been demonstrated that around the focus point, the curve of the integrated modulus of the complex amplitude was of parabolic shape. The accuracy is then obtained by fitting the curve with a parabola around the maximum. The best focus plane is then derived analytically using the fitting parameters. In this way, we can use larger increment Δz while keeping a high accuracy in z. In the experiments, we propagate numerically along z over 180 µm (covering the thickness of the experimental cell) with an increment Δz of 5 µm. The parabolic fitting is performed on 9 points equally distributed around the curve maximum. In this way, we obtain an accuracy along z of 1 µm (which is much better than the classical depth of focus of 6 µm). This accuracy is determined by measuring, with the same setup, the standard deviation of the linear-fitted curve z(t) of the sedimentation of a vesicle (which is linear in time).

Feature extraction
The xyz position of each vesicle obtained, size and shape must be determined. For this purpose, vesicles must be segmented on the phase map in order to delimit their border. The pre-segmentation obtained during the xy detection is rough and not sufficiently precise to get reliable parameters of the shape. This is due to the fact that the pre-segmentation is performed in the focus plane of the microscope which is, most of the time, not the one of the vesicles (see for example Fig. 7). The compensated phase is thus propagated to the corresponding z plane before to be segmented.
Boundaries of vesicles are smooth and require more sophisticated segmentation techniques like active contours for feature extraction. Active contour techniques are iterative processes that attempt to minimize a defined function and are much more suitable for applications with smooth boundaries. Starting from an initial contour centered on the object, the successive iterations deform the contour until the function is minimized and the deformed contour fits the boundaries of the object. We used an active contour algorithm based on the Mumford-Shah functional minimization [53]. Our motivation for choosing the active contour technique is guided by the following: we know the position (x,y) of the object and its approximate shape (from the pre-segmentation), around which we have cropped the ROI. Thus, as an initial contour, we can reasonably choose a circle with a diameter extracted from the presegmentation. In this way, the processing time of the segmentation is drastically reduced and much more reliable as the Mumford-Shah algorithm is very sensitive to the initial contour. The segmentation process delimits the boundaries of the vesicle. Feature extraction is then performed by fitting the segmentation contour with an ellipse. Size of vesicles is given by the big and small axis (respectively a 1 and a 2 ) of the fitted ellipse. Shape is given by the aspect ratio AR = a 1 /a 2 related to the deflation of the vesicle [14]. Figure 8 illustrates the segmentation and feature extraction process of two vesicles of almost the same equivalent diameter (around 20 µm) but different deflations.  Table 2 describes the different parameters of Fig. 8. As we can see, other parameters are extracted from the fitted ellipse such as the orientation of the big axis with respect to the flow direction (from top to bottom of the images) and the mean square error of the fitting. Orientation should be considered only for elliptical vesicles with an aspect ratio different of 1. Those additional parameters are useful to detect and remove false detections (like aggregates) or segmentation errors (mainly on the smallest and very low contrasted vesicles). Figure 9 shows example of false detections.  In the results presented in the following section, we have considered as false detections, MSE values above 1.1 µm. This value has been determined by selecting manually 100 different false detections and among those, measuring the smallest MSE. Further filtering of the data has been made on the orientation by removing orientations bigger than 15° when the AR exceeds 1.3. Indeed, the measurement of the orientation is relevant only for significantly deflated vesicles, whose long axis should be aligned with the flow direction. If this not the case, the detected object is an aggregate as depicted in the first row of Fig. 9. By this procedure, and depending on the concentration of the sample, we discard up to 5% of the data but strongly increase the reliability. Note that specific methods exist to process and separate aggregates [51,54]. Nevertheless, we did not implement them since for our study only the dynamics of single vesicles is considered and the reliability criterion consists in ensuring that no false detection is made. Tracking of individual vesicles is not of interest in our case as we are essentially interested in the description of the temporal evolution of the suspension. Obtaining the trajectories of individual vesicles would become meaningful for a deeper study of the diffusion process. However, this would require a much higher frame rate and a different setup, since in the range of shear rates explored, the velocity of vesicles is of order a few mm/s and the field of view only 0.37 mm long, which does not allow observing a significant portion of individual vesicle trajectories or even interaction processes.
Applying the processing presented in the two previous sub-sections on all the holograms of an experiment results in a 4D distribution describing the evolution of the tri-dimensional distribution. From this data, physical properties of the suspension can be extracted such as the lift velocity, shear-induced diffusion, population … Note that this algorithm is not restricted to vesicles and has been successfully applied for RBC dynamics under shear flow [55].

Illustrative examples
The hydrodynamic interactions between deformable particles and walls or other particles strongly depend not only on their size but also on their ability to deform. For objects with incompressible membranes, such as lipid vesicles or red blood cells, this ability strongly depends on their initial volume over surface ratio, which stays constant in time under flow. At rest, the excess area may however not be visible: the bending modulus is typically of order a few kT, which favors local folding of the membrane. The characterization of deflation within a given population therefore need to be done under flow. In simple shear flow, vesicles adopt an ellipsoidal shape whose main axis is in the shear plane (xy). Our analysis provides the small axis a 2 in the vorticity direction (z) as well as the apparent long axis a 1 in the flow direction (x), which is the long axis of the projection of the vesicle shape onto the (xz) plane. The ratio a 1 /a 2 characterizes the deflation of the vesicle. On Figs. 10 and 11, we show the size and apparent aspect ratio distribution within a vesicle suspension. Electroformation method provides a vesicle suspension with a strong size disparity. This suspension was sorted out thanks to a pinched flow fractionation method [56,57] and two subpopulations were mixed as to get a bi-disperse suspension. The distribution was extracted from a sequence of 1476 holograms in which 24362 vesicles were correctly detected. Depending on the size of the vesicles, the processing times may slightly vary. In average, on a recent processor, the algorithm processes between 700 and 800 vesicles per hour. Note that the process can be easily spread on multiprocessor platforms.  An evolution of the z position of all vesicles in the field of view along a shearing experiment is shown on Figs. 12 and 13. The data was extracted from two sequences of 600 holograms in which 2257 (Fig. 12) and 7575 (Fig. 13) vesicles were correctly detected. The vesicles initially lay on the bottom plate and are sheared during a microgravity phase, so that only hydrodynamic interactions between themselves and with the walls are present. The initial lift-off stage where the top plate is too far to play a role was described in a preceding paper [14] and is in agreement with several numerical and theoretical papers [27,31,36]. Lift velocity decreases with distance to the wall with a power law with exponent -2. At longer time, vesicles start to feel the influence of the top plate and eventually reach, on average, the middle plane of the cell. As shown on the graphs, this mean behavior is well described by assuming that the effects of both walls are additive, so that the mean velocity along the (z) axis scales as 1/z 2 -1/(h-z) 2 . As the vesicles migrate towards the center, their distribution broadens. This diffusion cannot be a Brownian diffusion due to the large size of the vesicles, but it is due to the multiple interactions between the vesicles. For not too large concentrations, pair interactions are predominant, so that the diffusion constant depends linearly on concentration [55,58]. This implies in particular that the width should increase with time following a power law with exponent 1/3 [55,59]. This is illustrated on Fig. 14, where standard deviations in the vesicle distribution are plotted as a function of time.   These examples demonstrate the accuracy of our procedure and will be embedded in a more complete set of data and theoretical modeling to be presented elsewhere, where the diffusion coefficients will be determined, as well as the cross-diffusion coefficients in the case of a bi-disperse suspension.

Conclusions
In this paper, we have proposed a fully automated processing of digital holograms for monitoring the dynamics of a vesicle suspension under shear flow. The understanding of complex fluids in confined flows is of importance as they have direct consequences on the rheology of biological fluids like blood and especially on the transport of oxygen by red blood cells in the smallest capillaries. The novel approach presented combines in elegant way different holographic techniques while keeping reasonable processing times with a good reliability. Automated detection is most of the time problem dependant, especially when we add the in-depth dimension. We separated the tri-dimensional detection by an in-plane followed by an in-depth detection. This approach presents several advantages notably a presegmentation of the objects indispensable for a good in-depth detection and feature extraction. Experimental results show a very good agreement with theoretical models demonstrating that DHM coupled with the proposed method is an adequate and powerful tool for studying complex fluids.