4D tracking of clinical seminal samples for quantitative characterization of motility parameters

: In this paper we investigate the use of a digital holographic microscope, with partial spatial coherent illumination, for the automated detection and tracking of spermatozoa. This in vitro technique for the analysis of quantitative parameters is useful for assessment of semen quality. In fact, thanks to the capabilities of digital holography, the developed algorithm allows us to resolve in-focus amplitude and phase maps of the cells under study, independently of focal plane of the sample image. We have characterized cell motility on clinical samples of seminal fluid. In particular, anomalous sperm cells were characterized and the quantitative motility parameters were compared to those of normal sperm.


Introduction
The human spermatozoon is a polarized motile cell, which delivers the haploid male genome to the oocyte, introduces the centrosome and triggers the oocyte into activity. Abnormal sperm behavior is one of the most common important indicators for clinical male infertility [1]. Of special relevance is sperm motility, i.e. curvilinear and linear velocity, wobble etc. In fact, both the concentration of progressively motile spermatozoa and their characteristics movements are significant factors in determining the outcome of homologous tests of human sperm-cervical mucus interaction [2]. For this reason, there is growing interest in understanding the kinematics and dynamics of swimming spermatozoa [3][4][5]. A commercial semi-automated device is currently available to perform sperm motility analyses [6]. This widely used instrument for clinical diagnostics is limited since the analysis is performed only in two dimensions and provides a partial in-plane representation of the motility features. In fact, 3D spatial motion is difficult to track using traditional microscopy as samples quickly move out of focus. A recent procedure for sperm cell tracking based on lens-less stereo microscopy has been demonstrated but, because of the low spatial resolution, it is incapable of imaging single cell features [7]. We propose an approach for particle detection and threedimensional tracking based on digital holographic (DH) microscopy. Digital holography allows the recording and the numerical reconstruction of the phase and amplitude of the specimen's optical wavefront [8] and the refocusing by numerical propagation from a single recorded hologram without realigning of the optical imaging system. Consequently a 3D volumetric field can be reconstructed by means of a single image (the hologram). DH enables high-resolution lensless imaging [9], nano sized particle detection [10] and the extension of the depth of focus [11]. Furthermore, DH may provide detection of three-dimensional features of biological microstructures [12], and allows quantitatively retrieving, in far field region, the amplitude and phase of the wavefront interacting with the structures themselves [13][14][15]. DH has been successfully employed to perform morphological analysis on bovine [16-18] and human [19, 20] sperm cells and the unique potentialities of DH for 3D particle tracking have been already demonstrated [9,21]. In this paper, we investigate the use of DH working with a partial spatial coherent source to detect automatically the 4D tracking (the 3D spatial motion over time) of organisms recorded on holograms, and to deduce from the retrieved path the relative motility parameters. The developed tracking algorithm has been initially implemented on easily available specimens characterized by complex motions, i.e. Paramecia moving in a free 3D liquid and a Giardia lamblia flowing in a microfluidic channel. Once tested and optimized, the DH tracking routine was applied to a clinical seminal sample moving in a free 3D space. The developed algorithm allowed us to collect resolved in-focus amplitude and phase maps of the cells under study, independently of the focus condition of the sample in the acquired image.

Experimental setup
The optical scheme relative to the DH is shown in Fig. 1. A major issue of the full coherent illumination is the speckle noise that arises on the unavoidable defects of the optical system, which can result in some cases in poor optical quality. The speckle results from the coherent superposition of random contributions that originate from out of focus locations. A way to reduce the strength of the speckle noise is to reduce the spatial coherence of the beam illuminating the sample. It has been demonstrated that a coherent beam incident on a moving ground glass behaves like a partial spatial coherent source, provided that the ground-glass motion is statistically averaging the optical properties of the outgoing laser beam [22]. The coherent source, (a mono-mode laser diode, λ = 635nm) is then transformed into a partially spatial coherent source by focusing the beam, by lens L1, close to the plane, with a rotating ground glass (RGG). The spatial coherence of the light can be varied by changing the position of the focused spot with respect to the RGG plane [22]. The Lens L2 (focal length f 2 = 50mm) collimates the beam that is divided by a beam-splitter BS1. As in the previous case, the object beam illuminates the sample S by transmission. When the RGG is stopped the speckle size can be measured in the object plane to control the spatial coherence. A plane of the sample is imaged by the couple of lenses L3-L5 on the CCD camera sensor, interfering with the reference beam. L3 is a microscope objective. An optical flat, not indicated in Fig. 1, is placed in front of the lens L4 to equalize the optical path of the reference and object beam. The camera is a JAI CV-M4 with a CCD array of 1280 × 1024 pixels (subsequently cropped to 1024 × 1024 to match the Fast Fourier Transform computation). The camera exposure time is 100 μs. The complex field of the object beam is reconstructed numerically from the frequency spectrum of the acquired hologram. The off-axis configuration of the employed set-up, i.e. a small angle is introduced between the two interfering beams, allows a spatial separation of the real and conjugate images due to the holographic reconstruction. Thus, the real image (or the conjugate one) can be separated from the whole spatial frequency spectrum with a bandwidth filter and shifted to the origin of the plane. In this way the complex spectrum of the field, except for a constant [23] is obtained and a quantitative phase measurement can be performed. Partial coherent illumination also removes the multiple reflections that can occur with coherent illumination when microscope slides on which biological specimen lying, are used. In fact, we assume that a reflection introduces an increase of the optical path d. If the distance d introduces a significant decorrelation of the speckle pattern, the contrast of the interference fringe pattern between the reflected beam and the direct beam is reduced.

Theoretical background
Once the whole field is known, it is possible to reconstruct the optical wavefront at different distances from the plane of acquisition, applying the Fourier formulation of the Fresnel-Kirchhoff diffraction formula [24]. An interesting approach based on the operator algebra, compatible with our propagation range, has been proposed by Shamir [25]. Fresnel diffraction is described by replacing the Fresnel-Kirchhoff integral, the lens transfer factor, and other operations by operators. The resulting operator algebra leads to the description of Fourier optics in a simple and compact way, bypassing the cumbersome integral calculus. The detail of the formalism can be found in [25]. By means of this approach the propagated field ( , ) prop O ξ η as a function of the initial field can be rewritten as: where N is the number of pixels in both directions, Δ is the sampling distance (i.e. the pixel dimension), m, n, U, V, h and j are integer numbers varying from 0 to N-1. The discretized Fourier transform is defined as Intensity and phase distributions can be retrieved starting from the propagated field according to the following relations: , t a n , . , The possibility offered by DH to manage the phase of the reconstructed image allows removing and/or compensating the unwanted wavefront variations (such as optical aberrations, slide deformations etc.) [27]. In this paper, a double exposure technique [28] has been used. A first exposure is relative to the object under investigation, whereas the second one refers to a flat reference region in proximity of the object. This second hologram contains information about all the aberrations introduced by the optical components, i.e. the defocusing due to the microscope objective. By numerically manipulating the two holograms, it is possible to compensate for these aberrations. As described by Eq. (2), DH provides a numerical investigation of the third dimension by performing a plane-by-plane refocusing. Indeed, if digital holographic reconstruction can refocus a sample slice-by-slice as the focusing stage of a classical imaging system, the refocusing degree of an object has to be determined by an external criterion. Several different solutions have been developed [29-31], the one we use in this study is based on the invariance of both energy and amplitude [32,33]. Those invariance properties allow building two focus criteria, respectively for pure amplitude and pure phase objects, based on the score of the integrated amplitude modulus. It has been demonstrated in [32] that the criterion is minimized for amplitude object while it is maximized for phase object. Biological cells are mostly transparent and can be considered as pure phase objects.

Tracking procedure
In order to detect automatically the 4D tracking (the 3D spatial motion over time) of the swimming samples a set of holograms has to be acquired at a constant sample -microscope objective distance. On the phase map, retrieved by means of the afore-mentioned numerical approach, the threshold is computed by applying the Otsu's method [34] and, after the automatic selection [12] of the region of interest, the area of the object under test is evaluated. Each following phase map is thus converted to a binary mask. The transversal position is evaluated by performing the cross correlation function between the selected object and the whole compensated phase map. In other words, the movement of the specimens from one image to following is tracked by finding the maximum value of a correlation coefficient between each pair of imaged cells. Therefore, the central position of the cell corresponds to the maximum value of the correlation function and its position is recorded as the X and Y coordinates of the cell. Finally the Z position (focus distance) is evaluated, by applying the self-focusing function (SFF) described above. In order to avoid a time consuming calculation of the SFF in the whole range of possible reconstruction distances and to overcome the limit of a chosen step for the evaluation of the SFF, an algorithm looking for the minimum of a function within a fixed interval has been used; The algorithm is based on golden section search and parabolic interpolation [35]. Reasonably the resolution of hundredth of micron for the evaluation of the focus plane distance (Z) has been chosen [36]. The criterion is repeated for the whole set of flux images, evaluating the X and Y coordinates and then the Z position.
The afore-described tracking algorithm has been initially tested on specimens with complex motions and more easily available than the seminal clinical samples. In particular, the movements of a Giardia lamblia flowing in a microfluidic channel (width of 120 μm and depth of 80 μm) have been initially tracked. These data have been captured using a 63 × microscope objective (NA = 0.70, Depth of Field = 1.2 μm) at 24 fps. 123 holograms of a Giardia lamblia have been acquired and the retrieved phase map is shown in Fig. 2(a). Once the phase irregularities due to the edges of the microchannel have been removed by means of a mask (Fig. 2(b)), the region of interest can be then automatically selected (Fig. 2(c)) and, by using Otsu's criterion, a binary image of the object is obtained (Fig. 2(d)). This image is used to find the same object in each phase map image. In fact, the central position of the cell corresponds to the maximum value of the correlation function ( Fig. 2(e)) between the selected object and the whole compensated phase map. The position is recorded as the X and Y coordinates of the cell (Fig. 2(f)).  The same 4D tracking algorithm has been tested on the movement of paramecia swimming in a small reservoir filled with stagnant water. For the following measurements the spatial coherence of the light, and consequently the speckle noise, has been increased to test the robustness of the algorithm. These unicellular specimens, with dimension of tenths of μm, are characterized by fast 3D movement and constitute the perfect sample for testing our tracking procedure on multiple particles. These data have been captured using a 10 × microscope objective (NA = 0.22) at 5 fps. The holograms acquired are noisier of the previous ones, because of the thickness of the liquid quantity (about 2 cm) and its composition: water, paramecia and suspended vegetable waste. So, a procedure for equalization has been performed in order to apply the 4D tracking procedure. In particular, the average of all the phase maps has been evaluated and it has been used as a background to remove by means of the relation: being φ the phase to be equalized ( Fig. 4(a)), φ mean the background obtained by meaning the reconstructed phase maps for different frames (Fig. 4(b)) and φ Eq the equalized phase ( Fig.   4(c)). Another solution proposed to remove the background noise consists in subtracting two following phase maps φ 1 and φ 2 (Fig. 4(d)-4(e)). The black and the white spots ( Fig. 4(f)) correspond to the position of the paramecium in the phase map φ 1 and φ 2 respectively. This latter approach has proved particularly useful for the tracking of multiple particles. In fact, when the tracking protocol described in previous section was applied to the whole phase map, it mixed up the different particles and a wrong tracking was estimated ( Fig. 5(a)). By means of a proximity criterion, the position in the (n + 1)-th frame has been searched in a reasonable neighborhood of the n-th frame position. The result of the proximity correction of the tracking procedure is reported in Fig. 5(b) and 5(c) (frame from Media 2). It is worth to note that the two objects in Fig. 5(b) are swimming on different planes even if, from a transversal perspective, they are apparently crossing each other. In order to increase the reliability of the tracking algorithm, the proximity criterion was applied to detect the Z position, too. In other words, the proximity criterion has been extended to a 3D box, allowing the tracking of objects whose phase map is not enough contrasted to apply successfully the Otsu's threshold method, as in the case, for example, of a blurred out of focus particle.

4D tracking of spermatozoa
The algorithm providing the extended proximity criterion has been finally used to track human sperm cells from a seminal clinical sample. The ejaculate has been collected by masturbation from a male scheduled to undergo in vitro fecundation at the Center for Assisted Fertilization -CFA, Naples (Italy). Donors provided informed consent under Local Ethical Approval. Sperm samples were collected by masturbation after abstinence for three to seven days and examined after liquefaction. All samples were washed using a silicon-based gradient of 40% overlaid over an 80% silicon solution (COOK Sperm Gradient, Ireland). The sample was centrifuged for 20 min at 1000 G, followed by a wash in Hams F-10. The final precipitate was resuspended to a final concentration of 1 × 106 sperm/ml and conserved in an atmosphere of 37°C and 6% CO2 until required. The analysis has been performed within 3 hours after the sample collection. These data have been captured using a 40 × microscope objective (NA = 0.65) at 5 fps.
Our algorithm has been applied to highlight spermatozoa anomalous behaviors. The cell in Fig. 6 is affected a morphological anomaly, known as "bent tail", which causes the nonlinear out of plane motion ( Fig. 6(a)-6(b), frame from Media 3). Using the holograms collected to track the spermatozoon motion, we can reconstruct the phase map of the cell (Fig.  6(c)) and highlight the morphological defect inducing the anomalous out of plane path. The unique feature of the Digital holography to retrieve quantitative phase information, allows using the tracking algorithm to evaluate many parameters useful for a reliable assessment of semen quality. In the following, the same parameters defined in [6] to provide quantitative motion values, have been estimated. The velocity values that are determined are the curvilinear velocity (VCL), straight-line velocity (VSL), and average path velocity (VAP). The VCL refers to the total distance that the sperm head covers in the observation period and it is always the highest of the 3 velocity values. The VSL is determined from the straight-line distance between the first and last points of the trajectory and gives the net space gain in the observation period. This is always the lowest of the 3 velocity values for any spermatozoon. The VAP is the distance the spermatozoon has travelled in the average direction of movement in the observation period: in cases in which the sperm head's trajectory is very regular and linear, with very little lateral movement, then the VAP is almost the same as the VSL. However, with irregular trajectories, such as those that are not linear, or where there is a high degree of lateral deviation of the head about the direction of movement, then the VAP will be much higher than the VSL. The average path is thus determined by adaptively smoothing the VSL. To describe the trajectory further, three velocity ratio values have been developed. These are linearity (LIN), a comparison of the straight-line and curvilinear paths, straightness (STR), a comparison of the straight-line and average paths, and wobble (WOB) a comparison of the average and curvilinear paths. They are thus defined as follows:  Linearity is an expression of the relationship between the three-dimensional path taken by the spermatozoon (i.e. curvilinear path) and its net space gain. A circling trajectory would have low linearity, because the curvilinear path (i.e. the circumference of the circle) would be much higher than the net space gain (i.e. the distance between the first and last points of the trajectory). A high linearity trajectory is one where the curvilinear path has a relatively low amplitude of lateral head displacement (see below), and the general direction of movement is the same as that of the straight-line path. Straightness gives an indication of the relationship between the net space gain and the general trajectory of the spermatozoon. A trajectory with evenly spaced track points and low amplitude would have a high STR value, since the average path would approximate the straight-line path. A circling track would be expected to have a low STR because the average path is the average of the curvilinear path, so the STR would be higher than the LIN, but still low. Wobble is the expression of the relationship between the average and curvilinear paths. WOB would be low for a track with a wide trajectory, but high for a circling track, since the curvilinear and average paths would be similar. For the cell in Fig. 6, VCL = 90.9 ± 0.5 μm/s, VSL = 1.2 ± 0.5 μm/s, VAP = 90.8 ± 0.5 μm/s, and consequently LIN = 1, STR = 1, WOB = 99.9.
The tracking of multiple spermatozoa, moving on different focal plane, is made possible by the full 3D approach described above. In Fig. 7(a)-7(b) five different spermatozoa are successfully tracked. It's worth noting that the sequence of holograms is acquired at a constant sample -microscope objective distance, allowing an in vitro analysis and subsequently a numerical approach allows the 3D volumetric field reconstruction. In the spermatozoa set shown in Fig. 7, an anomalous sperm cell is present, whose movement is plotted in green (frame from Media 4). In fact, while every other cell moves in parallel, swimming against a slight flow direction, the anomalous cell advances slower, along a broken track and on a tilted direction. Retrieving the motility parameters can highlight this anomalous behavior (Table 1). In particular the VSL measured for Cell1 is lower than for every other cell and describes effectively the inefficient cell movement. Accordingly, the wobble is pretty uniform for Cell2-5, varying between 97 and 99 and representing a motion with reduced oscillation around the average path. This value is instead sensibly lower for Cell1, providing a quantitative description of the wide fluctuation of the spermatozoon head. The reported data are consisted with the values reported in literature [6].

Conclusions
In this paper a label-free approach for spermatozoa tracking has been described. The algorithm, based on digital holography, allows acquiring a sequence of images at a constant specimen-microscope objective distance. The proposed algorithm has been developed and optimized in different experimental conditions, demonstrating a high robustness and versatility. Finally, in collaboration with CFA-ITALY, we studied the 4D tracking of human sperm cells, with the intent of providing a new and more complete approach for clinical investigations in the field of male fertility. Anomalous sperm cells were characterized and the retrieved quantitative motility parameters were compared to those of normal cells. The algorithm is potentially usable in every application where a three dimensional tracking may provide interesting information (i.e. migration analysis of cells in cancer research). Preliminary results of the experiments described in the present paper have already been reported in [37] and [38], and are the first reports of four-dimensional sperm cells tracking.