Adama Creppy
email: acreppy@imft.fr

Olivier Praud
email: opraud@imft.fr

Xavier Druart
email: xavier.druart@tours.inra.fr

Philippa L Kohnke
email: philippa.kohnke@tours.inra.frfplourab@imft.fr

Franck Plouraboué

Turbulence of swarming sperm

Collective motion of self-sustained swarming flows has recently provided examples of small-scale turbulence arising where viscous effects are dominant. We report the first observation of universal enstrophy cascade in concentrated swarming sperm consistent with a body of evidence built from various independent measurements. We found a well-defined k -3 power-law decay of a velocity field power spectrum and relative dispersion of small beads consistent with theoretical predictions in 2D turbulence. Concentrated living sperm displays long-range, correlated whirlpool structures of a size that provides an integral scale of turbulence. We propose a consistent explanation for this quasi-2D turbulence based on self-structured laminated flow forced by steric interactions and alignment, a state of active matter that we call "swarming liquid crystal." We develop scaling arguments consistent with this interpretation.

I. INTRODUCTION

Large-scale coherent dynamics in sessile drops of bacterial flows [1] or sperm [2,3] were reported many years ago. For spermatozoa confined between two plates [2], coherent vortex rings composed of as few as 10 individuals were observed to form. In the case of sperm, collective motility analysis based upon the observed turbidity in the vicinity of the sessile drop contact line has been used as a reliable index for male fertility scoring [3] in livestock. At a more fundamental level, there has been recent interest in the long-range spatiotemporal dynamics in concentrated suspensions of swimmers [2,[4][5][6][START_REF] Zhang | Proc. Natl. Acad. Sci. USA[END_REF][START_REF] Wensink | Proc. Natl. Acad. Sci. USA[END_REF][START_REF] Sokolov | [END_REF][10][11].

Collective motion of swarming flows is an emerging topic that involves challenging issues related to its continuum mesoscopic description, especially at high densities [12][13][14][15][16][17][18][19][20][21].

Recent observations in bacterial suspensions have shown that highly concentrated swarming flows display "turbulentlike" [4,[START_REF] Wensink | Proc. Natl. Acad. Sci. USA[END_REF]22] behavior without inertia. While the use of the term "turbulent" might seem inappropriate in the context of inertialess flow ("chaotic" or "intermittent" are arguably more appropriate), it has been used in the present article due to its use in the majority of previous studies in the field [4,[START_REF] Wensink | Proc. Natl. Acad. Sci. USA[END_REF]22].

Various turbulent phases are expected [START_REF] Wensink | Proc. Natl. Acad. Sci. USA[END_REF] depending on the confinement, the concentration, and the aspect ratio of the active swimmers, and different turbulent behaviors have indeed been observed [START_REF] Wensink | Proc. Natl. Acad. Sci. USA[END_REF]22]. However, many details remain to be considered to provide a clear and reliable picture of this complex state of active matter. This contribution suggests new evidence that turbulent swarming flows can indeed present clear self-similar behavior which exhibits many similarities with a direct enstrophy cascade.

Ram semen is highly concentrated with a volume fraction φ ≈ 50%. This volume fraction is calculated based on the concentration and the ram spermatozoon dimensions. The ram spermatozoon length is approximatively 60 µm with 50 µm of it being the flagellum. Its head has an ellipsoidal shape, flat at the bottom and round at the top with a thickness of 3 µm, a 10 µm length, and 5 µm width.

B. Camera recording and microscope setup

For gray-level images, we use a high-speed camera PCO Dimax S4 (Photon Lines, France) mounted on a phasecontrast microscope (Nikon, Japan). For fluorescent images, we use a pco.edge camera (Photon Lines, France) mounted on a fluorescent microscope (Olympus BH2, Japan). Both microscopes are equipped with temperature control systems and with ×4 magnification objectives. Recorded images are 2000 × 2000 pixels and have 12 bits of gray level. The frame rate is 50-100 Hz. Phase contrast is performed by a transparent phase plate in a phase-shift ring so as to enhance sample contrast. For the particle image velocity (PIV) measurement the pixel size is 1.08 µm × 1.08 µm, which gives recorded windows of 2 mm × 2 mm. The pixel size is 1.5 µm × 1.5 µm for particle tracking velocity (PTV) acquisitions giving 1.628 mm × 1.628 mm recorded windows. The ram sperm samples were confined in Leja (IMV) chambers of 100, 150, and 200 µm depth.

C. Particle image velocity (PIV) reconstruction

The standard image correlation PIV algorithm Davis TM has been used on the gray-level phase-contrast images. The projection of the Eulerian velocity field in the focus plane is reconstructed on a two-dimensional (2D) grid using a cross-correlation PIV algorithm. The algorithm estimates the most probable displacement of grid-centered small windows (identification window) between two successive frames. Iterative calculations are done to take into account the deformation of the identification window in order to increase the accuracy and reduce peak locking.

The spatial resolution is dictated by the size of this window, which in our case is 16 pixels × 16 pixels, which corresponds to 17 µm × 17 µm. This size is small enough to resolve spatial structures on the order of a few spermatozoa. The velocity field is obtained on a grid of size 126 × 126.

D. Particle detection, tracking, and PTV

Fluorescent microparticles (Spherotech, FH-2056-2) were mixed with semen at a final concentration of 3.9 × 10 6 beads/ml. The fluorescent microparticles have high-intensity Nile red fluorescence (excitation/emission 530/605 nm), size 1.7-2.2 µm and density 1.05 g/cm 3 . They are spherical and uniform in shape.

Particle detection algorithms have been coded in house using MATLAB TM . Individual beads are tracked along time from PTV of the fluorescent images as illustrated in Fig. 1.

Detection

Individual particles are first detected using a particle mask convolution analysis [23]. A particle template which consists of a 2D Gaussian distribution G having an isotropic standard deviation of σ = 5 pixels is scanned over the entire image in order to detect the peaks of image intensity which correspond to the particles. We compute the normalized convolution

C( i , j ) = i,j I (i,j )G(i + i )G(j + j ) i,j I 2 (i,j ) i,j (GiGj ) 2
, where G(i) = exp(-i 2 /2σ 2 )/ √ 2πσ is the normalized Gaussian, and I (i,j ) is the original image gray-level intensity. A 0.45 threshold is then applied to C( i , j ) in order to provide an isolated binary island for each event detection. The center of the detected particle is then computed by performing a 2D Gaussian regression [24] into the convolution field C( i , j ) to provide a subpixel accuracy of the position x n p of particle p of image n.

Tracking

Since the mean displacement of the particle images is small compared to the mean inter-particle distance, the matching particles pairs between two successive frames can be found by pairing each particle in the first frame with its closest neighbor in the second frame. For each detected position x n p of image n we find any particle p ′ of image n + 1, whose position x n+1 p ′ lies within a radius of R = 10 pixels, i.e., |x n+1 p ′ -x n p | < R. We then select the closest event for the pairing.

Particle tracking velocity (PTV)

For a particle p whose position, x n p , in image n, moves to position x n+1 p in the subsequent image n + 1, its velocity is v n p = (x n+1 p -x n p )/ t where t is the time between frames. By considering in each image n the relative distances

D 2 = |x n p -x n p ′ | 2 and relative velocities V 2 = |v n p -v n p ′ | 2 of
each (p,p ′ ) particle pair, we can build a very robust statistics of the relationship between D and V averaged over all pairs and all images. These data are shown in Fig. 10 below.

Comparison of PIV and PTV

Finally, we investigate whether phase-contrast gray-level texture PIV gives a sensible measurement of the real flow field from a direct comparison with PTV. Since small particles beads are perfect inertialess Lagrangian tracers of the flow, they offer the most reliable estimate of the local flow velocity.

In addition to the fact that PTV and PIV fluctuations have similar PDFs as shown in Fig. 2, we directly investigate for each event, the relative error between PIV and PTV. On the same sequence of images of pure semen seeded with Lagrangian PMMA-Rhodamine B-particles of 1-35µm diameter, a comparison was performed by interpolating the PIV to the beads particle locations. We found the relative error between PIV and PTV to systematically drop from 50% for particle diameters of 10-15 µm to 25% for particle diameters of 30-35 µm as shown in Fig. 3. This is consistent both with the optical averaging of the velocity field within the microscope : cubic interpolation; : spline interpolation. This is consistent with the optical field depth of the microscope δ = 25 ± 10 µm.

field depth δ = 25 ± 10 µm, as well as the PIV window size and the resulting interpolation error. This observation provides strong evidence to support that the velocity field measured from the image gray-level texture correlation is indeed very close to the depth-averaged true velocity of the flow over the microscope field depth.

E. Sample collection and scoring

Ejaculates from adult Lacaune rams were collected with an artificial vagina [START_REF] Baril | Manuel de formation pour linsémination artificielle chez les ovins et les caprins[END_REF]. All ejaculates used in this study scored over three for Massal mobility (MM), had a thick and creamy appearance, and were free of blood and urine contamination. Wave motion of semen was observed in a 10 µl drop of raw semen on the warming plate of a BH-2 Olympus microscope with a 4× phase-contrast objective. The MM score was performed based on the Herman definition (13). The MM score was refined from its original definition to subjectively assess the rotation speed of dark waves (5 being the maximal speed). The sperm concentration was also assessed by measurement of absorbance using spectrophotometry. Such measurements provide a precalibrated estimation of the cell concentration inside a cylindrical tube (the calibration has been validated against other techniques such as manual counting). From the concentration, and knowing the volume of each individual, we infer the volume fraction.

F. Rheological behavior of pure and dilute semen

We perform a rheological analysis of pure and dilute fresh samples. Between 1 and 4 h after collection, samples maintained at 15 • C were placed into the lower plate of a coneplate rheometer (Mars III Fischer Scientific). All dilutions were performed with isotonic NaCl. We mainly use a 1 • angle, 35 mm diameter cone, since the quantity of fluid involved is obviously limited for semen collections (typically 1.5 ml). Nevertheless, some results have been confirmed and validated with a 1 • angle, 60 mm cone diameter. Each measurement was performed during a steady state of the torque. The time lag associated with its stabilization is 180 s, and for each γ for pure semen and seminal plasma and their dilutions. Full symbols are for semen and hollow ones for seminal plasma. The color codes pertain to ( , ), pure-( , ), ×10 dilution, ( , ), ×100 dilution. The viscosity value decreases as the dilution rate increases for both semen and plasma until reaching the water viscosity for high shear rate. For each dilution rate the semen is more viscous than the plasma but still displays the same shear-thinning behavior with approximatively the same slope. measurement, the integration time for the torque evaluation is 30 s. Each experiment has been performed at 37 • while preventing the evaporation of the fluid during the measurement process (about 30 min) by confining the cone/plate setup into a controlled saturated water vapor enclosure. The shear rate γ has been varied from 10 -3 to 2 × 10 2 s -1 .

Figure 4 illustrates the semen viscosity, which displays a typical shear-thinning power-law behavior. Since a very similar curve is observed for seminal plasma, it can be deduced that this shear-thinning behavior is typical of the plasma protein suspensions and weakly dependent on the presence of spermatozoa. This result is consistent with previous rheological analysis for other species [START_REF] Severa | [END_REF]27]. This is also confirmed from the curve obtained for diluted plasma where we obtain the viscosity of water at large shear rate and the same power-law behavior at lower shear rates. Spermatozoon is a pusher-like swimmer, the presence of which decreases the semen viscosity [28][29][30], unlike to puller [31]. Due to the very reduced lifetime of spermatozoon activity, we were unable to perform viscosity measurements on active semen. All the presented rheology measurements were performed on weakly active or inactive spermatozoa suspensions. As found for suspensions of B. subtilis (pusher-like swimmer) [29], our results show that for nonactive cells the viscosity increases with concentration.

III. EXPERIMENTAL RESULTS

We experimentally investigate the statistical turbulence of semen with various complementary indicators.

After collection, a pure, fresh ram semen sample was introduced into a confined transparent glass chamber of 1 cm in length in each horizontal direction and a few hundred microns (H = 100,150,200 µm) in thickness. H is much larger than cell size since the ram sperm head size is a ≃ 10 µm. Direct

(x) = ∇ × u ′ = ∂ x v ′ (x) -∂ y u ′ (x)
normalized by the enstrophy ′ = |ω ′2 |/2 where • ≡ S •dS/S and S is the total horizontal surface. We see that the largest flow structure is typically smaller than the recorded windows but much larger than individual cells (see movie S3 in Ref. [32]).

In order to further analyze this integral length scale the spatiotemporal correlation functions were evaluated. This provides insight into the structure cascade and its dynamics [START_REF] Zhang | Proc. Natl. Acad. Sci. USA[END_REF][START_REF] Wensink | Proc. Natl. Acad. Sci. USA[END_REF][START_REF] Sokolov | [END_REF][10][11]. The energy-to-enstrophy ratio defines an integral length scale [START_REF] Wensink | Proc. Natl. Acad. Sci. USA[END_REF]

L ⋆ = E ′ / ′ , (1)
which is shown in Fig. 7(b). The relationship between ′ and E ′ displays a linear trend, leading to a well-defined integral length scale L ⋆ that is independent of the initial energy of the active flow. As expected, L ⋆ is close to the observed maximum whirlpool sizes. weak dependence is found, the significance of which will be discussed in Sec. IV.

We wish to compare this integral length scale L ⋆ with the whirlpool structures present in Figs. 5(a) and 5(b). For this we need another quantification independent of the flow field, but related with image texture, so as to detect and quantify a typical length scale for the whirlpools. For this, gray-level phase-contrast images are thresholded at 70% of the gray-level mode. The resulting binary images are then dilated (using dilation morphological operator) with 5 × 5 square structural elements. The resulting images are filtered with a median filter whose size is between 31 and 41 pixels and again binarized with the same threshold. A typical example is illustrated in Fig. 6(b), which can be compared with the original gray-level display in Fig. 6(a). Connected components associated with the whirlpools are then detected. An ellipse is extracted, the properties of which preserve the second moments of the connected component pixel region [cf. Fig. 6(b)]. The center of each ellipse corresponds to the barycenter of the connected component. We then use the major axis of the ellipse as a typical length scale l ⋆ of the whirlpools and compare it with the integral scale L ⋆ defined in Sec. III A. L ⋆ and l ⋆ are compared in Table II where one sees they have similar values. An average 6% difference is found between both estimation of the integral length scale. This is an interesting observation since one is obtained from kinematics properties of the flow, and the other one results from the analysis and segmentation of gray-level texture.

B. Two-point correlation of the flow field

The PIV and PTV statistics of the velocity fluctuations, u ′ (x) = u(x) -u , is found to be isotropic and close to Gaussian (see Fig. 2) as in Refs. [START_REF] Wensink | Proc. Natl. Acad. Sci. USA[END_REF]22]. The two-point correlation between point x and point x + ξ of u ′ (x) is analyzed along the direction parallel to ξ , u ′ = u ′ • ξ/|ξ |. The correlation function was found to be isotropic. Hence we consider only its dependence with distance |ξ |, as illustrated in Fig. 7(a). Considering the dimensionless distance |ξ |/L ⋆ , provides a nice collapse of the correlation functions in Fig. 7(c). We see that the integral length scale L ⋆ also defines the typical correlation length that is independent of the initial energy of the active fluid [see Fig. 7(b)]. It is interesting to mention that an integral time scale τ ⋆ can be also extracted from our measurements, as illustrated in Fig. 7(e). A nice collapse of the time-correlation profile is also obtained. This integral time scale displays a linear trend with the typical mesoscopic time scale built upon the enstrophy ′ = |ω 2 |/2 , τ ∼ 1/ √ ′ as shown in Fig. 7(d).

We now analyze the "internal" length scale of the velocity field from the swimmer size a to the integral scale L ⋆ . For this purpose, we consider the energy spectrum, which is the signature of the cascade in the kinetic energy. From the observation of an isotropic velocity Fourier transform û′ (k) in the wavelength k plane, we evaluate the averaged azimuthal kinetic energy spectrum in Fig. 8.

This spectrum displays three distinct regimes: at small k, it saturates to a white noise plateau, whereas, above a definite critical wavelength k c , it displays a clear power law decay. individual experiment, we estimate the error bars on the power-law decay to be -3 ± 0.4648. Such k -3 power-law decay is indeed found in quasi-2D turbulent flows. Then the last white noise plateau regime for k > k a = 2π/a is consistent with size of the order of individual swimmers. These suggest that the velocity correlations build up by coarse graining from the individual swimmer scale to reach broader collective scales within the range [a,L ⋆ ]. For dilute samples, the structures persist at large scale but turn cloudy at small scale (see movie S4 [32]). This qualitative observation is confirmed by the computation of the integral scale L ⋆ obtained from relation (1), while a linear trend is indeed observed between enstropy and energy (not shown) during the time decay of the sample motion. This integral scale L ⋆ does not change with the dilution as we found L ⋆ = 118 µm for pure, L ⋆ = 90 µm for ×2 dilution, and L ⋆ = 114 µm for ×4 dilution. Hence, volume fraction hardly affects the integral scale. This is confirmed by the computation of the power spectrum provided in Fig. 9, where it can be seen that the normalized spectrum amplitude is increased at small wavelength and decreased at larger wavelengths. This is the signature of a decrease of small-scale velocity correlations and lower energy at small scales. The effect of the dilution on the small-scale structures is also visible in the movies (see movie S4 [32]).

Finally, it is important to report that semen dilution drastically decreases the optical phase contrast of the suspension, so that texture correlation computed by PIV is not reliable or difficult to handle. In many cases, we have not been able to get reliable PIV measurements from many diluted samples. For dilution larger than ×8 the number of false vectors in the resulting PIV field becomes larger than 10% so that velocity field can no longer reliably be estimated. Hence, for diluted samples, the statistics of the presented results are less substantial than for pure semen ones.

C. Quasi-2D properties of the flow

The superposition of normalized instantaneous vorticity and streamlines, presented in Fig. 5(d), are in good agreement, qualitatively suggesting that the flow field is quasi-2D. This fundamental feature of the flow is more deeply investigated with two complementary indicators: (1) the ratio of the in-plane 2D kinetic energy to the out-of-plane kinetic energy and ( 2) the ratio between Lagrangian tracers' curvilinear path length in the observation plane to the microscope field depth.

(1) We first evaluate the mean 2D kinetic energy E = |u| 2 /2 = (u 2 + v 2 )/2 and examine the time over which the semen energy decays. A slow decay is observed over a few minutes. During recording sequences of 40 s, the kinetic energy typically decreases by about 40%. Computing the 2D divergence of the flow, we evaluate the out-of-plane fluctuating kinetic energy. The divergence of the flow field can indeed provide an estimate of the vertical velocity at the edges of the microscope depth field. Since locally fluid incompressibility holds, the divergence of the 2D flow components is

div u ′ = -∂ z w ′ , (2)
where w ′ is the vertical component of the fluctuating velocity field. Integrating (2) in z from z = -δ/2 to z = +δ/2, where δ is the microscope depth field, and assuming that the actual measured in-plane velocity is the depth-average velocity leads to div

δ/2 -δ/2 u ′ dz = -δw ′ , (3)
where δw ′ = w ′ (+δ/2) -w ′ (-δ/2) is the out-of-plane fluctuation velocity component differences at the edges of the microscope field depth δ. Thus, by computing the divergence of the 2D PIV in-plane flow, we can estimate the fluctuating component of the vertical flow and the vertical fluctuating kinetic energy,

E ′ w = |δw ′ | 2 /2 . Comparing E ′
w to E ′ averaged for all time over 3000 images and six different experiments, we found that the ratio E ′ w /E ′ = 8 × 10 -3 . This is a clear indication that the out-of-plane to the inplane velocity ratio |w ′ |/|u ′ | ∼ 0.0845 ∼ 1/12 is typical of quasi-2D flows. (1) We seed the semen with highly diluted fluorescent particles (volume fraction equals 1.2 × 10 -5 ). By tracking the particle positions within the flow (see Fig. 1 and Sec. II D) we evaluate their curvilinear path length in the horizontal plane ℓ 2D and then compute its average value ℓ 2D for all particles denoted. We find ℓ 2D = 356 ± 80 µm, which is much larger than the estimated microscope field depth δ = 25 ± 10 µm, as δ/ℓ 2D ∼ 0.0844 ∼ 1/12. This is consistent with the previous estimate |w ′ |/|u ′ | ∼ 1/12. One measurement is illustrated with data gathered from 3000 frames, with, in average, 575 detected particles per image pair. Symbols represent experiment on different males ( : male no. 1; : male no. 2; ⋆: male no. 3), and solid black lines indicate a power-law fit whose slope is equal 1.05 consistent with V ∼ D scaling.

D. Separation law

Here we analyze the specific scaling of the particle tracer trajectory statistics. One interesting prediction in turbulent flows is related to the so-called separation law, also called relative dispersion [33], which is the relation between the average relative distance and the average relative velocity for each particle pair. This law for three-dimensional (3D) Kolmogorov turbulence differs from that from the enstrophy cascade [33]. For 3D Kolmogorov turbulence, the relative velocity V [each particle pair (p,p ′ ), V 2 = |v p -v p ′ | 2 ; see Sec. II D] for two particles in the flow scales as V ∼ D 1/3 where D is the relative distance [34,35]. In 2D turbulence associated with an enstropy cascade scaling, the prediction becomes V ∼ D [33,36]. From PTV (Sec. II D), we evaluate the relative velocity V of a very large number (≈3 × 10 7 ) of particle pairs in order to build the histogram of the relative velocity versus relative distance D. By computing the average velocity V over the interval D,D + δD, we can test the scaling prediction for relative particle dispersion. Figure 10 shows that the linear prediction V ∼ D is consistent with our observations within more than one decade in distances. The quality of the experimental curve relies on the very large number of events: an average of 575 particles found correlated for each image pair (image pairs are needed to compute each particle velocity) and statistics are computed over 3000 images to give an excellent average. The scaling obtained in the log-log plot of Fig. 10 are highly reproducible from one sample to another (only three of them are illustrated). This very robust behavior is a reliable test for the enstrophy-cascade-like velocity spectrum, which is consistent with the spectrum obtained in Fig. 8.

IV. DISCUSSION

A. Scaling derivation

The development of continuum equations for concentrated suspensions is still an active field of research [12][13][14][15][16][17]22,[37][38][39][40]. Here we provide simple insights for the leading asymptotic structure of a confined active suspension to examine whether the quasi-2D structure of the flow is imposed by the Hele-Shaw chamber or by some internal self-structured layering.

In the following, we show that self-structured layering gives a description consistent with our experimental findings, as a vertical length scale a ≪ H is shown to be a relevant length-scale. Since this length scale a also equals the size of a spermatozoa head we show that layered sheets of horizontally aligned spermatozoa are slipping one over another and forming quasi-2D turbulent layers.

Obviously the standard 2D asymptotic reduction of Stokes flow leads to the Darcy-Brinkman equations [41]. But in the case of active fluid suspensions, the driving mechanism comes from a supplementary stress associated with the influence of active swimmers which drive the flow through dipolar forcing generates.

Fluid incompressibility associated with velocity field (u,v,w) ≡ (u,w) reads

∇ • (u,v,w) = ∇ • u + ∂ z w = 0, (4)
where ∇ ≡ (∂ x ,∂ y ). When considering an integral length scale L ⋆ for the typical in-plane variations along (x,y), a typical swimming velocity of spermatozoa u ∼ v ∼ U s , and a typical vertical length a, then incompressibility provides

w ∼ a U s /L ⋆ ∼ ǫU s (5)
and the definition of a small parameter ǫ = a/L ⋆ . The momentum balance equations which govern the motion of an active suspension on sufficiently large length and time scales at zero Reynolds number [37,38] are

µ (u,v,w) -∇p = -∇ • σ E = Cσ 0 ∇ • Q, Q = pp - I 3 , (6)
σ 0 = 3φµU s ℓ 2 4π a 3 ,
where the extra-stress tensor σ E = -Cσ 0 Q is composed of C, a positive constant for spermatozoa which are pushers, Q, the stresslet-nematic tensor averaged over the orientation of the swimmers, and σ 0 , a dimensional active stress constant. σ 0 is composed of two lengths: a, associated with the spermatozoa head, and ℓ, which is the flagellum length, and φ, the volume fraction. The extra-stress tensor thus scale as σ E ∼ µU s ℓ 2 /a 3 . Using dimensionless length scale L ⋆ ( x, ỹ,ǫ z), gradient ∇ ≡ 1/a(ǫ∂ x ,ǫ∂ ỹ ,∂ z), and velocity scales U s ( ũ, ṽ,ǫ w), the leading order dimensionless form of ( 6) is

µU s a 2 ∂ 2 z ũ -∂ x p = ǫ µU s ℓ 2 a 4 φ ∂ x σ x x + ∂ ỹ σ x ỹ + ∂ z σ x z , µU s a 2 ∂ 2 z ṽ -∂ ỹ p = ǫ µU s ℓ 2 a 4 φ(∂ x σ ỹ x + ∂ ỹ σ ỹ ỹ + ∂ z σ ỹ z), ǫ µU s a 2 ∂ 2 z w - 1 ǫ 2 ∂ z p = ǫ 2 µU s ℓ 2 a 4 φ(∂ x σz x + ∂ ỹ σz ỹ + ∂ z σzz ). (7)
The dimensionless pressure is chosen to balance the in-plane viscous forces, so at leading order p ≡ (µU s L ⋆ /a 2 ) p, which is invariant along the gap span from the third momentum balance component, ∂ z p = 0. Balancing the viscous and extra-stress leading order contributions in each (x,y) direction provides the same consistent scaling ǫ ∼ a 2 φℓ 2 = a L ⋆ so that one finds an estimate for the longitudinal length scale L ⋆ ∼ φℓ 2 a . One should also note that the dimensionless right-hand side of [START_REF] Zhang | Proc. Natl. Acad. Sci. USA[END_REF] implies that the horizontal orientation gradients ∂ p are of same order of vertical one ∂ z p z . This indicates that the averaged orientation is mainly horizontal in the (x,y) plane, i.e., p ∼ p x ∼ p y ∼ O(1) and weakly vertical p z ∼ ǫ. One can also evaluate the resulting scaling for the in-plane kinetic energy E ∼ U 2 s and the enstrophy ∼ (U s /L ⋆ ) 2 , so that L ⋆ is also given by L ⋆ ∼ E ′ / ′ , as evaluated in Fig. 7(b). This analysis is consistent with the following: (1) The ratio between the vertical and horizontal velocity |w|/|u| ∼ |w|/|U s | ∼ a/L ⋆ ∼ ǫ, which also equals the ratio of vertical to horizontal bead displacements. Furthermore, since we measure |w|/|u| ≃ 1/12, it gives ǫ ≃ 1/12, and L ⋆ ≃ 12a. Choosing a ≃ 10 µm for ram sperm leads to L ⋆ ≃ 120 µm, which is consistent with the observed range L ⋆ ≃ 120-145 µm. Furthermore, this leading order scaling predicts that L ⋆ is independent of the width H consistent with the observation of a weak experimental variation as H varies (Table II).

(2) Integral scale estimate L ⋆ ∼ φℓ 2 /a with a ≃ 10 µm and ℓ ≃ 50 µm for ram sperm, φ ≃ 1/2, predicts L ⋆ ≃ 125 µm, which is also consistent with the observed range L ⋆ ≃ 100-150 µm.

(3) The ratio of p z / p ∼ ǫ ∼ a/L ⋆ indicates a flow structured into preferentially horizontal layers the thickness of which is a, reminiscent of swarming liquid crystals [19][20][21].

Finally, it is important to recall that the above derivation also includes the hypothesis of strongly confined layering of the suspension, which is consistent with most concentrated semen studied here, but not diluted ones. This observation provides strong evidence to support that the velocity field measured from the image gray-level texture correlation is indeed very close to the depth-averaged true velocity of the flow over the microscope field depth.

B. Discussion and perspectives

Let us first summarize the points raised by this study. We analyze the velocity of concentrated sperm suspensions with PIV reconstruction applied to gray-level textures obtained from phase-contrast microscope high-speed image recordings (as in Refs. [1,4,[START_REF] Wensink | Proc. Natl. Acad. Sci. USA[END_REF][START_REF] Sokolov | [END_REF]). In Sec. II D, we obtain a consistent comparison between PIV measurements on the optical fielddepth averaged image and passive fluorescent beads, real local flow tracers velocity. This is strong and direct experimental evidence that PIV provides a measurement of the real velocity field. Since the precise physical origin of gray-level texture associated with phase-contrast images is unclear, such direct evidence was a necessary step. As a side remark, we should also mention that we have indirect evidence for a close link between the phase-contrast signal and local orientational order of nonaxisymmetric sperm cells. Furthermore, we put forward additional coherent evidence:

(1) The instantaneous streamline and isovorticity lines are nearly identical.

(2) Direct and indirect (computing the flow divergence; see Sec. III C) PIV measurements provided a ratio between vertical to horizontal velocity close to 1/12.

(3) This ratio has been found to be consistent with the ratio between the horizontal Lagrangian tracers meandering and the microscopic field depth.

(4) The integral scale associated with the largest structures of the flow is found to be much larger than the individual swimmmers (see Sec. III A) and weakly sensitive to the gap span of the experimental chambers.

(5) We found consistent values of the kinematic estimation of the integral scale (provided by a clear experimental linear relationship between the energy and the enstrophy) and that image provided by based extraction from whirlpool characteristics, we found consistent values (see Sec. III A). This is experimental evidence that the typical maximal size of gray-level whirlpools depict the decorrelated units of the flow structure.

(6) The horizontal energy spectrum is found to scale like an enstrophy cascade with k -3±0.4648 power law (see Sec. III B).

(7) When decreasing the volume fraction, the large structures of this turbulent flow are mainly unaffected, while the small scales display less energy. [START_REF] Wensink | Proc. Natl. Acad. Sci. USA[END_REF] The physical evidence for the flow field velocity to display such a specific power-law spectrum has been challenged with an independent experimental measurement based on the separation law of particle beads (see Sec. III D). By estimating the relative velocity V of each pair of beads with their relative distance D we found clear evidence for a V ∼ D scaling. Since the theoretical prediction for a V ∼ D scaling is based upon a 2D k -3 power law of the energy spectrum, this is consistent with the k -3±0.4648 behavior that we measure for the horizontal energy spectrum.

(9) Finally, using a very simple model (see Sec. IV A) for active swimmers in a Newtonian flow where active swimming provides an extra stress, we found that a typical vertical scale close to the sperm size leads to a consistent scaling for the integral scale.

All these points provide a solid body of evidence that a layered 2D stratified turbulent flow is taking place in the described experiments. We believe that this feature is generic for concentrated suspensions of active and elongated swimmers. This layering can be explained by the combined effect of steric-interaction alignment and boundary influence. Highly concentrated elongated self-propelled particles, close to packing limit have no other choice than local alignment. This steric alignment mechanism applies in the horizontal direction because this is the only possible direction near the upper and lower glass plate. We believe that horizontal alignment propagates from upper and lower boundaries, down into the fluid bulk (the thickness of the chosen chambers are between 10 to 20 cell layers). A simple scaling analysis indicates that if the thickness of the layers are close to the micro-swimmers size a, as suggested by a steric alignment mechanism of this layering, a horizontal integral scale can be predicted, the estimate of which is consistent with our measurements.

Since the proposed mechanism to explain the observed lamination is a steric repulsion alignment, we believe there might be some analogy between the system under study and "swarming nematic" systems (e.g., see Fig. 3 of [19], [START_REF] Zhou | Proc. Natl. Acad. Sci[END_REF]). While we do not wish to push further this analogy in more detail, we mention it for completeness. As mentioned earlier, previous studies had already documented the statistical properties of swarming flows of active suspensions [4,5,11]. Most of studies were conducted at volume fraction much lower than ours with the exception of Ref. [START_REF] Wensink | Proc. Natl. Acad. Sci. USA[END_REF], mainly using bacterial suspensions. Furthermore, the energy spectrum found in these studies rarely presents a simple behavior. Reference [10] report a velocity spectrum which increases at low wave vectors and decreases at large ones. Moreover, the observed decay fitted with a k -8/3 at large k over a fraction of decade is indeed quite similar with the k -3 reported here. This similarity supports the theory that the reported turbulence is quite generic to stratified or layered anisotropic concentrated active matter.

Finally, it is important to mention that we report in Sec. II F that the rheology of semen displays a power-law shear-thinning behavior. Such non-Newtonian rheology adds an important amount of complexity to the reported measurements. However, we believe that most qualitative features of the reported observations weakly depend on the fluid rheology. As a matter of fact, short-range (possibly lubricated) steric interactions between self-propelled elongated objects will always produce alignment. This is the alignment mechanism which produced at very high concentration the observed layering effect. Hence, the steric induced horizontal layering scenario evoked earlier which is responsible for the resulting 2D behavior of the flow would remain a salient property of anisotropic concentrated active matter flows.

Nevertheless, it is difficult to assess how much fluid rheology influences the observed power-law decay of the energy spectrum. It would certainly be of interest to better understand how the plasma rheology influences the observed dynamics. However, one should bear in mind that separating active cells from plasma is not easy, while, at the same time, preserving their mobility property. It will be experimentally challenging to access the impact of the separation methodology on the apparent changes of active suspension motion.

C. Conclusion

We have performed a systematic experimental analysis of the collective motion of fresh semen confined horizontally between two glass plates. PIV flow fields show "turbulent-like" behavior between the scale associated with the swimmer (spermatozoa) to a well-defined integral scale L ⋆ . A k -3±0.4648 behavior of the energy spectrum is reported which is consistent with the V ∼ D scaling of the separation law. That is typical of k -3 2D energy spectrum. Other indirect quantitative indicators have been put forward, all consistent with an apparent 2D turbulent flow in this concentrated active suspension. The reported results are interpreted to be associated with a laminated turbulent layer flow.

FIG. 1 .

 1 FIG. 1. (Color online) Trajectories of fluorescent spheric shape particles. Each color represents the curvilinear path of a particle. x and y units are in microns.

1 FIG. 2 .

 12 FIG. 2. (Color online) The normalized PDFs of the PIV ( ) and PTV ( ) fluctuations velocity have a Gaussian similar fit. ( ) Distribution.

5 FIG. 3 .

 53 FIG. 3. (Color online) The decrease of the relative error against beads diameter d. Comparison between PIV and PTV for different interpolants of the PIV on the particle positions : linear interpolation;: cubic interpolation; : spline interpolation. This is consistent with the optical field depth of the microscope δ = 25 ± 10 µm.

1 FIG. 4 .

 14 FIG. 4. (Color online)Viscosity µ dependence on shear rate γ for pure semen and seminal plasma and their dilutions. Full symbols are for semen and hollow ones for seminal plasma. The color codes pertain to ( , ), pure-( , ), ×10 dilution, ( , ), ×100 dilution. The viscosity value decreases as the dilution rate increases for both semen and plasma until reaching the water viscosity for high shear rate. For each dilution rate the semen is more viscous than the plasma but still displays the same shear-thinning behavior with approximatively the same slope.

FIG. 5 .

 5 FIG. 5. (Color online) (a, b) Phase-contrast microscope images of fresh semen placed between two glass plates separated by 100 and 150 µm, respectively (movies S1 and S2 in Ref. [32]). (c) PIV field obtained from (a). (d) Normalized color-coded vorticity field ω ′ (x,y)/ √ ′

FIG. 6 .

 6 FIG. 6. (Color online) (a) Phase contrast gray-level images (b) Whirlpool detection, segmentation, and ellipse extraction.

FIG. 7 .

 7 FIG. 7. (Color online) (a) Longitudinal correlation function of the fluctuation velocity field, ρ(ξ ) = u ′ (x)u ′ (x + ξ ) / u ′2 obtained from PIV measurements. The different set of curves corresponds to different semen samples recorded within chambers whose depths are : H = 100 µm; : H = 150 µm; : H = 200 µm. (b) The 2D fluctuating kinetic energy E ′ = |u ′ | 2 /2 is plotted versus the fluctuating enstrophy ′ = |ω ′ | 2 /2 for six different experiments each of more than 2000 images. Each circlular symbol represents an experiment in a 100 µm chamber where error bars represent the enstrophy standard deviation on each hole sequence. The dotted line shows the linear trend associated with a constant ratio E ′ / ′ , the square root of which provides the integral length L ⋆ = E ′ / ′ . (c) The collapse of the curves when considering dimensionless distance ξ/L ⋆ along the horizontal axis. (d) The linear trend between integral time scale τ ⋆ obtained from the first zero of the temporal correlation function ρ(τ ) = u ′ (t)u ′ (t + τ ) / u ′2 and the mesoscopic time scale τ obtained from enstrophy. (e) Temporal correlation function ρ(τ ) plotted versus the dimensionless time τ/τ ⋆ , where τ ⋆ defines the characteristic time obtained from the first zero of the temporal correlation function.

FIG. 8 .

 8 FIG. 8. Fluctuating energy spectrum S(k) = 2π 0 | û′ (k)| 2 k dk θ averaged over the azimuthal angle dk θ in the wave-vector plane k with polar representation k ≡ (k,k θ ) plotted versus the modulus of k, k = |k|, using bilogarithmic scales. The wave-vector modulus k is rescaled by a critical wavelength k c = 2π/L ⋆ based upon the integral scale L ⋆ . Data gathered from six different experiments in the 100 µm chamber are averaged over 2000 velocity fields reconstructed from PIV. The spectrum displays three regimes: k within [k c , k a ] where the motion is well correlated and presents the k -3 power law (gray dotted line); k > k a which defines the uncorrelated noise (at the individual scale); and k < k c where we have at large scale coherent but uncorrelated structures shown in Fig. 5.

2 FIG. 9 .

 29 FIG. 9. (Color online) Fluctuating energy spectrum S(k) versus wave-vector modulus k for a single sample of concentration 4.4 × 10 9 cells/ml and its dilutions. , pure, , ×2 dilution, , ×4 dilution.

3 FIG. 10

 310 FIG.10. (Color online) Relative velocity V of two particles versus their relative distance D is plotted in bilogarithmic scales in order to challenge turbulent particle dispersion scaling predictions. One measurement is illustrated with data gathered from 3000 frames, with, in average, 575 detected particles per image pair. Symbols represent experiment on different males ( : male no. 1; : male no. 2; ⋆: male no. 3), and solid black lines indicate a power-law fit whose slope is equal 1.05 consistent with V ∼ D scaling.

  Table I shows L ⋆ for different values of H . A

TABLE I .

 I The integral length scale L ⋆ for different gap spans of the glass chambers H . N is the sample number and φ the evaluated volume fraction.

	N	H (µm)	L ⋆ (µm)	C (10 9 cells/ml)	φ(%)
	6 1	100 150	126 ± 16 133	5.2 ± 1.1 5.4	50 ± 11 52
	1	200	146	5.4

TABLE II .

 II The integral length scale L ⋆ and the whirlpools characteristic length l ⋆ for different values of H .

	H (µm)	l ⋆ (µm)	L ⋆ (µm)
	100	119	109
	100	125	126
	100	112	108
	100	129	128
	100	150	137
	100	130	150
	150	175	133
	200	183	146

ACKNOWLEDGMENTS

This work has been supported by the French Agence Nationale pour la Recherche (ANR) in the frame of contract MOTIMO (ANR-11-MONU-009-01). We thank Pierre Degond, Eric Climent, Laurent Lacaze, and Frédéric Moulin for interesting discussions. We thank Sébastien Cazin for the imaging technical support. We thank Alexis Viel and Christophe Korbuly for their help in preliminary experiments.


Dominique Lavenier
email: dominique.lavenier@inria.fr

Charles Deltel

David Furodet

Jean-François François Roy

BLAST on UPMEM

Keywords: BLAST, genomic banks, Processing--in--Memory BLAST, banques génomique, Processing--in--Memory

This paper presents the implementation of the BLAST software on the UPMEM architecture. BLAST is a well--known molecular biology software to rapidly scan DNA or protein genomic banks. It is daily used by thousands of biologists. Given a query sequence, BLAST reports all the sequences of the bank where similarities are found. The result is a list of alignments that detail the similarities. UPMEM's Processing--In--Memory (PIM) solution consist of adding processing units into the DRAM, to minimize data access time and maximize bandwidth, in order to drastically accelerate data--consuming algorithms. A 16 GBytes UPMEM--DIMM module comes then with 256 UPMEM DRAM Processing Units (named DPU). To find similarities, BLAST proceeds in 3 steps: (1) search of common words between the query sequence and the bank sequences; (2) evaluation of local similarity on the neighbourhood of these words; and (3) computation of the final alignment. As the 2 first steps are limited by memory bandwidth, and represent the majority of time of the overall computation, they have been massively parallelized on UPMEM DPUs. The 3 rd step is performed on the host processor and is overlapped with the UPMEM processing. Experimentation on real datasets shows a speed--up of 25 when using UPMEM configuration, versus a standard server running 20 Intel cores.

RR N° 8878

Blast Software

BLAST is a popular bioinformatics software aiming at finding similarities between DNA or protein sequences [START_REF] Altschul | Basic local alignment search tool[END_REF] [START_REF] Altschul | Gapped BLAST and PSI-BLAST: a new generation of protein database search programs[END_REF]. BLAST is used every day by thousands of biologists or bioinformaticians to scan the genomic databases or to perform intensive sequence comparison between large genomic data sets. It is a fundamental processing tool, used in many bioinformatics pipelines. BLAST software consists of a family of software, each one dealing with different types of sequences.

With the fast evolution of sequencing techniques, BLAST (as other software) must face with the explosion of data. Comparing large sets of data coming from sequencing projects with regular standard banks is becoming extremely time-consuming. And, in more and more cases, becomes a serious bottleneck. There is a strong need to speed up this task.

The primary goal of BLAST is to search genomic banks to report similarities between sequences. More precisely, from a set of query sequences, BLAST reports all alignments between theses sequences and all the sequences of a reference bank. The output of BLAST is a list of local alignments, as represented below:

A local alignment details the way two portions of DNA or protein sequences are matching. Attached to an alignment is an expected value (probabilistic information) allowing BLAST to select only relevant ones. When large genomic banks are processed, an alternative output format summarizing the features of the alignment is preferred. This compressed format is better suited for post processing tasks.

The BLAST algorithm is based on a powerful heuristic, assuming that relevant alignments contain similar sequences with common words on N characters. These words are called seeds and are primarily searched. The structure of the BLAST algorithm can thus be decomposed into 3 main steps:

1. Seed search 2. Seed extension 3. Alignment computation To search seeds, the query sequence is first indexed on N-word characters basis. Then, the reference bank is scanned, words by words, to find similar N-words. When there is a match between a query and a RR N° 8878 reference word, the neighborhood of the seed is explored, as step 2. This step simply computes a local similarity. If it exceeds a threshold value, then a complete alignment is computed.

Profiling the BLAST code indicates that most of the computation time is performed in step 1 & 2. As a matter of fact, in many cases, the match of two seeds arises by chance and does not lead to any alignments.

Step 2 works as a filter that rejects most of the match proposals.

Comparing one query sequence with a full reference bank generates millions of seed matches. But all these matches can be processed independently. In other words, there is a great opportunity for massive parallelism, as offered by the UPMEM technology.

UPMEM Architecture Overview

UPMEM technology is based on the concept of Processing-in-Memory (PIM). The basic idea is to add processing elements next to the data, i. e. in the DRAM, to maximize the bandwidth and minimize the latency. Host processor is acting as an orchestrator: It performs read/write operations to the memories and controls the co-processors embedded in the DRAM. This data-centric model of distributed processing is optimal for data-consuming algorithms.

A 16 GBytes UPMEM DIMM module comes with 256 processors: One processor every 64 MBytes of DRAM. Each processor can run its own independent program. In addition, to hide memory latency, these processors are highly multithreaded (up to 24 threads can be run simultaneously) in such a way that the context is switched at every clock cycle between threads.

The UPMEM DPU is a triadic RISC processor with 24 32-bits registers per thread. In addition to memory instructions, it comes with built-in atomic instructions and conditional branching bundled with arithmetic and logic operations.

From a programing point of view, two different programs must be specified: (1) the host program that will dispatch the data to co-processors memory, sends commands, input data, and retrieve the results; (2) the program that will specify the treatment on the data stored in the DPU memory. This is often a short program performing basic operations on the data. It is called a tasklet. Note however, that the architecture of the UPMEM DPU allows different tasklets to be specified and run concurrently on different blocks of data.

Depending on the server configuration (i. e. the number of 16 GBytes UPMEM modules), a large number of DPU can process data in parallel. Each DPU only accesses 64 MBytes and cannot directly communicate with its neighbors. Data exchanges, if needed, must go through the host processor. A DPU has a fast working memory (64 Kbytes) acting as cache/scratchpad memory and shared by all tasklets (threads) running on the same DPU. This memory working space can be used to transfer blocks of data from the DRAM, and can be explicitly managed by the programmer.

To sum up, programing an application consists in writing a main program (run on the host processor) and one or several tasklets that will be executed on the DPUs. The main program has to synchronize the data transfer to/from the DPUs, as well as the tasklet execution. Note that the tasklet execution can be run asynchronously with the host program, allowing host tasks to be overlapped with DPU tasks.

BLAST implementation on UPMEM

In this section, we present the implementation of BLASTp/BLASTx that compares a set of query proteins or DNA sequences with a protein databank. The difference between BLASTp and BLASTx is that BLASTx first translates the DNA sequences into the six possible protein reading frames. For each query, an index is computed. Figure 1 details the index structure. A first table (Index1) memorizes for each seed its number of occurrences in the query (N part) and an entry number (P part) in a second table (index2). This second table stores the coordinates where the seed occurs in the query. Inria reported. To speed--up the process, this task is optimized with SSE instructions and run in parallel using all available cores of the processor. HOST/DPU task overlapping When a large set of query sequences are compared to a reference bank, host and DPU tasks can be advantageously overlapped as shown below. The host has the charge of computing the query index and the final alignment step, while the DPUs perform a parallel scanning of the reference bank. In that case, the overall execution time will be constrained by the longest task. An optimal implementation will have to provide a good balance between these tasks that can be adjusted either by the number of DPUs or the number of processor cores.

Performance evaluation

Speed up evaluation

Performances of the BLAST implementation have been evaluated on a DELL server (Xeon Processor E5-2670, 40 cores 2.5 GHz, 64 GBytes RAM) running Linux Fedora 20.

The two following dataset have been considered:

For dataset 1, the ncbi BLASTp program (version 2.2.29+) has been run with 20 threads and with an evalue set to 10 -6 :

blastp -query query_bank -db ref_bank -evalue 1e-6 -comp_based_stats F -outfmt 6 -max_target_seqs 500 -num_threads 20 -out query_ref.m8

For dataset 2, the ncbi BLASTx program (version 2.2.29+) has been run with 20 threads and with an evalue set to 10 -6 :

blastx -query query_bank -db ref_bank -evalue 1e-6 -comp_based_stats F -outfmt 6 -max_target_seqs 500 -num_threads 20 -out query_ref.m8

Dataset 1

The DPU frequency is set to 750 MHz The speed-up is given by: Blast time / max (DPU time, HOST time)

The speed up cannot exceed 25 because the host contribution is the limiting factor for this dataset.

On this dataset, in average, a DPU requires 280 cycles to process one amino acid (AA) on the reference bank, that is 373 ns. Knowing that the reference bank is equitably split over all the DPUs, the overall DPU execution time is given by: DPU execution time = NB_QUERY x NB_AA_REF / NB_DPU x 373 ns NB_QUERY is the number of blocks of queries sent to each DPU (4658) NB_AA_REF is the number of amino acids of the reference bank (194 x 10 6 )

We also need to add the communication time between the host and the DPUs. Intel mlc tool has been used to evaluate the memory bandwidth of the server (see detail on Annex 2). We choose the 3:1 Reads/Writes configuration providing a bandwidth of 46.8 GBytes/sec. Host/DPU communication can be decomposed into 4 steps as described below (annex 3 details each steps).

1 -Load of the database (only once at the beginning of the process). The reference bank load time (BT) is given by: BT = SIZE_REF_BANK / MEMORY_BANDWIDTH = 194 x 10 6 / 46.8 x 10 9 = 4.14 ms

• SIZE_REF_BANK is the size of the reference bank. One amino acid is stored on 1 Byte.

2 -Load of the queries. Several queries are packed together in order to optimize the query index. The size of the index is equal to 6.4 Kbytes. The same index is broadcasted to all DPUs. The query load time (QT) is given by :

Dataset 2

The DPU frequency is set to 750 MHz The speed-up is given by: Blast time / max (DPU time, Host time)

The speed up maximum is 28.5. The number of DPU is the limiting factor for this dataset.

On this dataset, in average, a DPU requires 220 cycles to process one amino acid (AA) on the reference bank, that is 293 ns. Knowing that the reference bank is equitably split over all the DPU, the overall DPU execution time is given by: • The seed search step considers words of only 2 amino acids only. The seed size of BLAST is 3, but is leveraged by the substitution matrix. The UPMEM implementation is thus more sensitive at this step. • The seed extension step is restricted to a window of 16 amino acids right next to the seed. It counts the number of identical pseudo amino acids. The amino acid alphabet has been reduced to 14 according to the work of [START_REF] Li | Reduction of protein sequence complexity by residue grouping[END_REF]. When a specific number of matches is detected, we consider that a complete alignment computation can be fired. Overall, the sensitivity of UPMEM BLAST is comparable to the NCBI BLAST. Both of them are based on the same general heuristic. However, step 1 & 2 of UPMEM BLAST have better sensitivity, but it can be slightly damaged by the SSE constraint of step 3. Practically, both software find 95 to 97 % of identical alignments. The rest belong to the grey zone depending of many optimization threshold acting all along the production of an alignment, including statistical consideration allowing an alignment to be discarded or not.

Comparison with CUDA BLASTP

CUDA BLASTP is one of the best implementation of BLAST on GPU. In [START_REF] Liu | CUDA-BLASTP: Accelerating BLASTP on CUDA-Enabled Graphics Hardware[END_REF], the authors test their implementation with the NCBI BLAST software by comparing protein sequences of different length with the GenBank Non-redundant Protein Database (3,163,461,953 Amino Acids in 9,230,955 Sequences).

The GPU board has the following features: NVIDIA GeForce GTX 295: dual GPU,

Estimation on a 512-DPU UPMEM system

We can estimate the performances of a 512 DPU architecture as follows:

DPU execution time = NB_AA_REF / NB_DPU x T_AA

The number of cycles for a DPU to process one amino acid of the bank depends of the query size. There is a fixe cost to scan the reference (about 150 cycles) then a cost that linearly depends of the query size, up to 250 cycles for the largest sequence (2026 AA). Thus T_AA = NB_CYCLE / 0.75 x 10 9 (frequency set to 750 MHz).

For the host execution time, we apply again a ratio of 1/25 on the NCBI BLASTP execution time.

The following table gives an estimation of the UPMEM execution time (and speed up compared to NCBI BLAST) on the dataset used by the authors of CUDA-BLASTP.

1 .

 1 Blast Software ......................................................................................................................................................................... 2. UPMEM Architecture Overview ............................................................................................................................................. 3. BLAST implementation on UPMEM ...................................................................................................................................... 4. Performance evaluation .......................................................................................................................................................... 5. Comparison with CUDA BLASTP ......................................................................................................................................... Bibliography ................................................................................................................................................................................ Annex 1: Tasklet Program .......................................................................................................................................................... Annex 2: Intel MLC report ......................................................................................................................................................... Annex 3: Host/DPU data transfer ............................................................................................................................................... Inria

Inria

  As introduced in section 1, the BLAST algorithm can be decomposed into 3 steps: (1) seed search; (2) seed extension; (3) alignment computation. Only the two first steps, that are the most time consuming, are offloaded to the UPMEM DPUs.The protein bank is split across all DPUs. Hence, each DPU stores only a part of the protein bank. Every query is broadcasted to all DPUs to perform both the search seed step and the search extension step. Inside a DPU, the task is again split into similar subtasks that will only have to look after a restricted part of the protein bank devoted to that DPU. The general implementation can be described as follows:Dispatching the reference bank into the DPUs is done in a round-robin fashion to avoid burst of similar sequences into the same DPU. This is done once, at the beginning of the treatment. Data are read from the disk and the reference bank is stored both on the host DRAM and on the UPMEM DRAM. Actually, in the current implementation, the host processor needs to completely access the bank to perform the final step. Storing the bank only on the UPMEM side would require extra time dedicated to data exchange, an overhead estimated to 10% of the overall computation. The advantage of the current implementation is that the time devoted to exchange the data between the host processor and all the DPUs is negligible compared to the effective computation time (see host/DPU exchange estimation time in section 4).

Figure 1 :

 1 Figure 1: structure of the query index

  Dispatch the reference bank inside DPU For each query sequence: Index the query Broadcast the query and its index to the DPUs Perform parallel scanning on UPMEM Get results from UPMEM Compute alignments

  QT = (NB_BLOCK x SIZE_INDEX x NB_DPU)/MEMORY_BANDWIDTH QT = (4658 x 6400 x 1024)/46.8 x 10 9 = 652 ms 3 -Load of parameters needed by tasklets. Two integers are sent for each tasklet. The parameter load time (PT) is given by: PT = (NB_BLOCK x NB_DPU x NB_TASKLET x 8)/MEMORY_BANDWIDTH PT = (4658 x 1024 x 10 x 8)/46.8 x 10 9 = 8.of the results. Each DPU send back to the host the coordinates (2 integers) of the sequences when a match has been detected. In average, for this dataset, and for a block of queries a DPU generate 200 couples of integers. The results transfer time (RT) is given by: RT = (NB_BLOCK x NB_DPU x 200 x 8)/ 46.8 x 10 9 = 163 ms The complete host/DPU communication time is equal to: BT + QT + PT + RT = 828 ms. It is included in the host execution time in the above table.

  DPU execution time = NB_QUERY x NB_AA_REF / NB_DPU x 293 ns NB_QUERY is the number of blocks of queries sent to each DPU (163) NB_AA_REF is the number of amino acids of the reference bank (18 x 10 9 ) The Host/DPU computation time is evaluated using the same bases as Dataset 1 and provide the following times: BT = SIZE_REF_BANK / MEMORY_BANDWIDTH = 18 x 10 9 / 46.8 x 10 9 = 385 ms QT = NB_BLOCK x SIZE_INDEX x NB_DPU)/MEMORY_BANDWIDTH = 24.2 ms PT = (NB_BLOCK x NB_DPU x NB_TASKLET x 8)/MEMORY_BANDWIDTH = 0.3 ms RT = (NB_BLOCK x NB_DPU x 300 x 8) = 9 ms Quality Evaluation Alignments computed by the UPMEM version of BLAST have been compared with the BLAST NCBI software. Results slightly differ for the following reasons:

  The table below gives the elapsed execution time reported with the Linux time command.

			Dataset 1 -BLASTp	Dataset 2 -BLASTx
	BLAST Exec. time	2h05		1h40
	The execution time of the tasklets (step 1 & 2) is estimated from the Cycle Accurate Simulator (CAS)
	developed by the UMEM Company. The execution time of the banded Smith & Waterman algorithm (step
	3) is directly measured from an optimized code developed by the GenScale IRISA/INRIA research group.
	HOST	Index Query	Index Quer		Compute step 3	Index Quer	Compute step 3	Index Quer
		1	y 2		Query 1	y 3	Query 2	y 4
			Compute step 1&2	Compute step 1&2	Compute step 1&2
	DPU		Query 1		Query 2	Query 3
								time
			Query bank		Reference bank
	Dataset 1	Arabidopsis Thaliana proteom 31477 proteins	SwissProt, version xxx 545 536 proteins (194 x10 6 AA)
	Dataset 2	Marine meta-transcriptomic 1000 RNA sequences	RefSeq, version 74 50 496 795 proteins (18 x 10 9 AA)
	RR N° 8878

•

  The alignment computation is done on restrictive banded Smith & Waterman algorithm using the SSE technology available on today processors. The size of the band depends of the technology: SSE4, AVX and AVX2 provide respectively a band size of 16, 32 or 64. BLAST has no such restriction. For fragmented alignments (i.e. alignment with many long gaps) the impact is that several alignments are generated instead of a single one.

	# DPU	512	1024	2048	4096
	DPU t execution time (sec)	1680	840	420	210
	HOST execution time (sec)	180	180	180	180
	Speed up	3.5	7.1	14.2	28.5
	RR N° 8878

  1,242 MHz processor clock speed, 60 SMs (30 per GPU), 1,792 MB GDDR3 device memory (896 MB per GPU). Tests have been conducted with these cards installed in a PC with an Intel Quad-Core i7-920 2.66 GHz CPU, 12 GB RAM running Linux Fedora 10. The following table reports the performance speed up of GPU for different query sizes:

	Query size (AA)	127	254	517	1054	2026
	NCBI BLAST Exec. Time (sec)	34.5	36.1	63.8	131.7	155.3
	CUDA-BLASTP Exec. Time (sec)	5.8	5.9	11.9	37.7	29.2
	Speed up	5.9	6,1	5.3	3.5	5,3

Inria

RR N° 8878

For performance purpose, the index must fit in the UPMEM working memory and has a limited size. The query sequences must not exceed 3000 characters. Fortunately, the average size of a protein sequence range from 300 to 400 amino acids. We then optimized the process by packing several sequences in the index when possible. The query sequence (or the group of queries) and its index are then broadcasted to all DPUs. Each DPU starts its tasklets that perform independent searches. As in BLAST, the reference is scanned word by word. Each word queries the index that indicates (or not) a potential match (step 1 of the algorithm). If so, the second step is fired and a score corresponding to a sequence similarity is calculated (see tasklet program below). When the score exceeds a predefined threshold value, coordinates of both the query and reference sequences are memorized in a list. At the end of the scan, each tasklets send back this list that represents the coordinates of all successful seed-extensions. When all DPUs have finished, the host computes the final alignments based on this list of coordinates. Note that this execution scheme is sequential: step 1 & 2 need to be performed before step 3 starts. To maximize the overall throughput an overlap strategy is implemented: As there are many queries to process, overlapping step 1 & 2 with step 3 is straightforward. The total execution time is thus given by the maximum of step (1 + 2) and step 3.

Tasklet program

The tasklet program performs step 1 & 2. It scans, sequence by sequence, a small part of the protein bank. A sequence is systematically decomposed into words of N characters (seeds) that are compared to the query index. For each seed match, a local similarity is computed, by counting the number of identical characters in a restricted neighborhood (here, a window of 16 characters). If the number of identical characters exceeds a threshold value, then the coordinates of both the query and the bank sequences are sent as results. Below is a simplified version of the algorithm. The complete code is given in annex 1.

Alignment program

The alignment program performs step 3 and is executed on the host processor. It receives from all the DPUs a list of coordinates where significant seed extensions have been found between the query sequence and some sequences of the bank. From these coordinates, a banded Smith & Waterman algorithm [3][4]