{'fraction_non_alphanumeric': 0.037767358211831206, 'fraction_numerical': 0.01872310360331428, 'mean_word_length': 4.667819439388424, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Multifocal imaging for precise, label-free tracking of fast biological processes in 3D

Many biological processes happen on a nano- to millimeter scale and within milliseconds. Established methods such as confocal microscopy are suitable for precise 3D recordings but lack the temporal or spatial resolution to resolve fast 3D processes and require labeled samples. Multifocal imaging (MFI) allows high-speed 3D imaging but is limited by the compromise between high spatial resolution and large field-of-view (FOV), and the requirement for bright fluorescent labels. Here, we provide an open-source 3D reconstruction algorithm for multi-focal images that allows using MFI for fast, precise, label-free tracking spherical and filamentous structures in a large FOV and across a high depth. We characterize fluid flow and flagellar beating of human and sea urchin sperm with a z-precision of 0.15 µm, in a volume of 240 × 260 × 21 µm, and at high speed (500 Hz). The sampling volume allowed to follow sperm trajectories while simultaneously recording their flagellar beat. Our MFI concept is cost-effective, can be easily implemented, and does not rely on object labeling, which renders it broadly applicable.


Introduction
Life happens in three dimensions (3D). Organisms, cells, and sub-cellular compartments continuously undergo 3D movements. Many biological processes happen on the micrometer to millimeter scale within milliseconds: in one second, insects flap their wings 100 to 400 times 4 optical path lengths 14,16,23 , and (5) lenses of different focal power 15,17,19 . Several approaches suffer from some disadvantages: optical gratings are wavelength-dependent and have low optical efficiency. Deformable mirrors are bulky, expensive, and highly wavelength-dependent, which prohibits their use for low-cost compact applications. Finally, dual-objective microscopes are expensive and complex -e.g., they require two tube lenses and cameras -and are limited to two focal planes only and require fluorescent labels. Therefore, we aim for an MFI system that combines approaches (4) and (5) and thus, can be accommodated in different experimental set-ups. We use a multifocal adapter that splits the light coming from the microscope into four light paths that are projected to different areas of the camera chip ( Supplementary Fig. 1a). By inserting four lenses of different focal power into the individual light paths, four different focal planes of the specimen are obtained (Supplementary Fig. 1b and   1d). Using only four focal planes minimizes the loss of SNR due to light splitting and keeps the system low in complexity. We established an alignment procedure based on recording images of a calibration grid (see Materials & Methods) to precisely overlay the four focal plane images ( Supplementary Fig. 2) and to correct for subtle magnification differences between the four focal planes (Supplementary Table 1, 2). Of note, the images of the calibration grid did not show any apparent comatic aberrations, spherical aberrations, field curvature, astigmatism, or image distortions ( Supplementary Fig. 1c, 2a). Using the thin-lens approximation, we can predict the defocusing of a pattern in the set-up based on the set of lenses and the magnification (Supplementary Fig. 1c and 1e). We assembled the multifocal adapter to study objects of submicrometer to millimeter size by changing the magnification (Supplementary Table 3). This flexibility enables studying fast-moving objects, ranging from subcellular organelles to whole animals.

Extended depth-of-field imaging of fast-moving objects
Long-term imaging of a fast-moving object with high spatial and temporal resolution is a common, yet challenging goal in microscopy. Imaging at high resolution with objectives of 5 high magnification and numerical aperture constrains the depth-of-field and the FOV, which increases the odds that the specimen exits the observation volume, thereby limiting the duration of image acquisition. This limitation can be overcome by combining images acquired at different depths using extended-depth-of-field (EDOF) algorithms 24 . However, this technique is not suited for fast-moving objects if different focal planes are acquired sequentially. MFI allows employing EDOF algorithms for fast-moving objects (Fig. 1). We combine MFI (  [1][2][3][4]. Of note, the EDOF method also reveals a depth map ( Fig.   1c) that, in combination with MFI, allows coarse 3D tracking at high speed.

High-precision 3D tracking with MFI
MFI has been applied for 3D tracking 13,[15][16][17][18][19] , but only with poor precision or a small sampling volume. The image of an object, i.e. the sharpness, size, and intensity of the image, depends on the object's z-position relative to the focal plane. This dependance is used in "depth-fromdefocus" methods to infer the z-position of an object based on the defocus of the object. Such methods have been applied to estimate the third dimension in physical optics [25][26][27] and, for biological specimens, from 2D microscopy images [28][29][30][31] . To determine the z-position, a precise calibration of the relationship between the image of the object and the z-position is required.
We combined the depth-from-defocus approach with MFI to achieve high-precision 3D tracking. We illustrate the power of this concept by reconstructing the 3D movement of 6 spherical objects (latex beads) and filamentous structures (sperm flagella) using label-free imaging by dark-field microscopy and a multifocal adapter.
To calibrate the relationship between the image and the z-position of a spherical object, we acquired multifocal z-stacks of non-moving 500-nm latex beads by moving the objective with a piezo across a depth of 21 µm; the step size was 0.1 µm. The bead image underwent two characteristic changes as the objective was moved: the radius of the bead image increased (Fig.   2a), and, as a result, the image intensity decreased (Fig. 2b). Both features are z-dependent and can be used to determine the bead z-position. The z-position of a particle cannot be inferred unequivocally from a single plane. For any given bead radius, two possible z-positions exist ( Fig. 2c). However, combining the information from multiple planes allows unequivocally localizing the particle in z (Fig. 2c) and, in addition, allows combining multiple measurements of the z-position in different planes by applying Bayesian inference probability 32 or other statistical methods to improve the z-precision. We estimated the z-precision for different zpositions and for different focal planes based on the relationship between bead radius and bead z-position (see Materials & Methods) (Fig. 2d). The z-position near the focus is imprecise because the radius of the bead image hardly changes when the bead's z-position is changed. By contrast, if the bead is slightly defocused, the relationship between the radius of the bead image and the z-position is steep, providing a better precision. When multiple focal planes are available to infer the z-position of a bead, the plane yielding the best precision can be selected to infer the z-position (Fig. 2e). Thereby, a better z-precision across the entire sampled volume is achieved. Increasing the number of planes allows increasing the range and z-precision (Fig.   2e). For the setup with four focal planes, we predicted a mean z-precision of 162 µm across a depth of 20 µm.
We next experimentally scrutinized the predicted z-precision. We recorded a new data set of multifocal z-stacks from immobilized latex beads using a piezoelectric-driven objective, measured the radius of the bead images in the different focal planes, and inferred the z-positions 7 of the beads using the relationship between radius and z-position (Fig. 2 a) (see also Materials & Methods). Comparison of the inferred and piezo z-positions revealed a linear relationship (mean ± standard deviation of the slope: 0.98 ± 0.02, n = 5 beads). We determined the zprecision based on linear fits with a slope of unity (Fig. 2f); the predicted and experimentally determined z-precision closely matched across the entire depth (Fig. 2g).
To gauge the validity of our 3D-localization approach, we studied the stochastic Brownian motion of beads (n = 81) that were tracked for a mean duration of 3.3 s (Fig. 2h) and we determined the x, y, and z distributions of bead displacement (Fig. 2i). The distributions indicated that stochastic bead motion was isotropic and normal distributed over all spatial coordinates. From the distributions of bead displacement (Fig. 2i) we inferred the x, y, and z precision of our 3D-localization approach based on the assumption that the distributions represent the positional variations resulting from diffusion superimposed on the positional variations resulting from the measurement precision (see Materials & Methods, Fig. 2j). We obtained a z-precision of 45 nm, 37 nm, and 154 nm in x, y, and z, respectively. Of note, the zprecision agrees with the z-precision predicted and measured during the calibration procedure ( Fig. 2g). We conclude that MFI is highly suitable for precise 3D tracking in a large volume.

High-precision 3D reconstruction of the flagella beat with MFI
Eukaryotic cells such as sperm or green algae deploy motile, lash-like appendages protruding from the cell surface -called flagella. The 3D flagellar beat propels the cell on a 3D swimming path 29,33 . Because the beat frequency is high (up to 100 Hz) 6,7 , the precise reconstruction of flagellar beating at high temporal and spatial resolution is challenging. The fastest calibrated 3D reconstructions of human sperm achieved temporal resolutions of 100 frames per second (fps) with DHM 34 and 180 fps using a piezo-driven objective 35 . This temporal resolution allows determining the beat frequency of human sperm that ranges, considering cell-to-cell variability, from 10 to 30 Hz 36 . However, this resolution is not sufficient to characterize higher harmonics of the beat 6,37 or to track and characterize the first harmonic of fast-beating flagella, 8 e.g. from Anguilla sperm (about 100 Hz) 7 . Detecting the second harmonic frequency (60 Hz) for a sperm beating at a fundamental frequency of 30 Hz will require, according to the Nyquist Shannon theorem, an acquisition rate of at least 120 fps. To determine the beat frequency of Anguilla sperm (about 100 Hz), an acquisition rate of at least 200 fps is required. Of note, the beat frequency of sperm is variable, rendering the determination of the flagellar beat frequency noisy. Thus, a temporal resolution of 120 fps would suffice only for regularly beating cells.
Recently, using holography, the flagellar beat of human sperm has been reconstructed with a high temporal resolution (2,000 fps) 38 . However, the small spectral power compared to the first harmonic and the low number of beat cycles sampled during the short recording time (0.5 s) did not allow for a clear detection of the second harmonic. Moreover, the tracking precision was not evaluated and the sampled volume was very small (ca. 30 x 45 x 16 µm) 38 compared to the flagellum length (50 µm), which compromises deriving the relation between flagellar beat and trajectory.
We aimed to establish MFI for reconstructing the 3D flagellar beat of human sperm at high speed (500 Hz), with high precision, and in a large FOV that allows tracking both the flagellar beat and trajectory over a long time.
We used a magnification of only 32x, which decreased the optical resolution compared to previous methods, employing 100x magnification 17,38 . However, this yielded a 58-times larger sampled volume of ca. 240 x 260 x 20 µm compared to a DHM-based method (30 x 45 x 16 µm) 38 and an 80-times larger sampled volume compared to an MFI-based method (ca. 80 x 35 x 5.6 µm) 17 . We compensated for the lower optical resolution by establishing a precise reconstruction method similar to that for bead tracking.
To characterize the relationship between image and z-position of flagella, we acquired multifocal z-stacks of immotile sperm using the piezo-driven objective. Similar to the latex beads, the flagellum image displayed a characteristic widening (Fig. 3a) and an intensity decrease ( Fig. 3b) when imaged out of focus.
We refined SpermQ, a software for analyzing motile cilia in 2D 39 Supplementary Fig. 3). To determine the z-precision, we inferred the z-position of an immotile sperm using SpermQ-MF in multifocal z-stacks recorded with the piezo-driven objective (Fig. 3c). Inferred and piezo z-positions showed a linear relationship with a slope of unity (Fig. 3c). From linear fits to the relationships of inferred and piezo z-position at different arc length positions (n = 3 sperm) (Fig. 3c), we determined a zprecision of ~1 µm at the head and midpiece and 0.6 µm at the principal piece (Fig. 3d).
We next tested whether our approach is applicable to flagella of other species by applying the same procedure to sperm from the sea urchin Arbacia punctulata ( Supplementary Fig. 4). Here, fits to determine the flagellar width at the neck region were compromised by the highly flexible neck and the high-intensity contrast between head and flagellum in sea urchin sperm. Thus, the z-precision at the flagellar neck region (arc length < 10 µm) was twofold less for sea urchin sperm (~2 µm, Supplementary Fig. 4c-d) than for human sperm (~1 µm, Fig. 3d). In the more distal part of the flagellum (arc lengths > 10 µm), sea urchin sperm were reconstructed with a better z-precision (~0.4 µm, Supplementary Fig. 4d) compared to human sperm (~0.6 µm, Fig.   3d).
In conclusion, our approach provides high temporal resolution (500 Hz), a large FOV of ca. 240 x 260 x 20 µm, and allows flagellar 3D reconstruction with a z-precision of down to ~0.4 µm, whereas for other methods the z-precision has not been reported 17,38 .

3D reconstruction of the flagellar beat with MFI
We next studied and compared the flagellar beat of free-swimming sperm from human and sea urchin. Using the calibrated MFI setup and SpermQ-MF, we recorded the 3D flagellar beat of n = 6 human and n = 10 sea urchin sperm swimming in a 150-µm deep chamber. For human sperm, the tracking duration was limited by the FOV to about 2-3 s, i.e. the time it takes a sperm cell crossing the FOV (Fig. 4a, Supplementary Movie 5). Due to the precise localization of the flagellum, we resolved the 3D flagellar beat in detail (Fig. 4b, c). The kymographs revealed a flagellar beat wave that propagates in 3D along the flagellum (Fig. 4b).
Making use of the large FOV of our MFI system, we recorded the flagellar beat and swimming path over a time span that allows relating beat pattern and trajectory ( Sea urchin sperm can spontaneously depart from circular swimming to produce a stereotypical motor response known as "turn-and-run" 40 . Making use of the large FOV of our method, we were able to record a sperm during such a motor response and correlate the swimming trajectory and the flagellar beat ( Supplementary Fig. 9). The flagellar beat transitioned from an asymmetric beat during the "turn" to a symmetric beat during the "run" (Supplementary Fig.   9).
We further quantified the 3D trajectory and the 3D flagellar beat (Fig. 5c-l). Human sperm swam at lower speed than sea urchin sperm ( Fig. 5c-d). Additionally, the flagellar shape of human sperm was characterized by a large out-of-plane component (Fig. 5a) while the flagellar shape of sea urchin sperm was rather planar (Fig 5b). We quantified the out-of-plane component using the non-planarity ratio 29 (see also Materials and Methods). This ratio ranges from zero to one. A ratio of zero indicates a perfectly planar beat, whereas a ratio of unity indicates a flagellar beat that has no preferential plane. The ratios largely fluctuated over time, especially for human sperm (Fig. 5e): The mean standard deviation over time of the non-planarity ratio was 0.17 ± 0.02 (mean ± standard deviation of n = 6 sperm) in human sperm and 0.08 ± 0.03 (mean ± standard deviation of n = 10 sperm) in sea urchin sperm. To exclude any bias by large fluctuations, we determined the median of the ratio over time as a reference for the nonplanarity per sperm cell (Fig. 5e, f). The median non-planarity ratio was 0.26 ± 0.04 (mean ± standard deviation of n = 6 sperm) in human sperm and 0.17 ± 0.04 (mean ± standard deviation of n = 10 sperm) in sea urchin sperm ( Fig. 5f), confirming a more planar beat of sea urchin compared to human sperm. In addition, the beat plane in human sperm rolled around the longitudinal sperm axis (Supplementary Movie 6), whereas no rolling was detected for sea urchin sperm (Fig. 5g, h). Our measurements of the non-planarity ratio (0.26) and the rolling velocity (3.3 Hz, representing 20.7 rad/s) agree with a previous analysis using phase-contrast microscopy 29 .
Sperm steering has been related to asymmetric flagellar beat patterns that can be described by  12,41,42 . We observed higher harmonics in the spectrum of human and sea urchin sperm (Fig. 5i). The curvature components C0 and C1 were much more prominent in sea urchin compared to human sperm, whereas the second harmonic component C2 was only marginally different between sea urchin and human ( Supplementary Fig. 10). We determined a second harmonic component of C2 = 0.006 ± 0.002 µm -1 (mean ± standard deviation) for human sperm. Notably, the widths of the frequency peaks were narrower for sea urchin than for human sperm (Fig. 5i, k), indicating that the beat frequency of human sperm is less constant than the beat frequency of sea urchin sperm. To confirm this observation, we analyzed the beat frequency of human and sea urchin sperm at different recording times, which revealed large variations in the beat frequency of human sperm over time but only small variations for sea urchin sperm (Fig. 5l). Such variations have not been reported in previous analysis, where a constant beat pattern with sharp frequency peaks was measured 6,37 . However, these measurements relied on tethering sperm at their head to establish sufficiently long observations for resolving the power-spectrum. When tethered, human sperm are prevented from rolling -they are forced into a stable condition and their beat plane is aligned with the glass surface by hydrodynamic interactions 43 . In contrast, our technique, due to precise tracking in a large FOV, allows recording freely-swimming human sperm under unconstraint conditions and over a sufficiently long time to detect variations in the flagellar beat.

3D particle imaging velocimetry to reconstruct fluid flow around a human sperm cell
The behaviour of ciliated cells and the underlying physics is of great interest for cell biology and biotechnological applications. We are only beginning to understand how sperm manage to navigate in a highly complex environment like the female genital tract 44

13
To visualize the 3D flow around human sperm with our MFI technique, we used latex beads whose relationship between radius and z-position had been previously calibrated. To follow the direction of the bead movement relative to the sperm flagellum, we tethered the sperm head to the glass surface.
Bead tracking revealed long-range hydrodynamic interactions between sperm and beads ( Fig. 6a, Fig. 6b, Supplementary Movie 7). At large distances from the sperm cell, beads displayed a random Brownian motion (Fig. 6c). Near the cell, beads displayed a Brownian motion with a drift towards the sperm head (Fig. 6d) and away from the flagellar tip (Fig. 6e).
Near the flagellum, beads described a fast 3D spiral-like movement resulting from the 3D flagellar beat (Fig. 6d-e).
We computed the 3D flow field around a human sperm based on a time-lapse recording of 40 s at 500 Hz, revealing 2,920,731 bead positions with a predicted median z-precision of 0.34 µm (prediction based on the approach shown in Fig. 2d; median used to exclude outliers in this calculation). The flow field shows that the fluid moves from a large distance towards the sperm cell and converges at the flagellum, where the flow speed steeply increases (Fig. 6f,

Discussion
State-of-the-art MFI systems that employ four focal planes have been limited to a depth < 6 µm and use high-magnification objectives (≥ 60x) to resolve the specimen, thereby constraining the FOV (≤ 80 x 35 µm) 15,17 . Although such a small FOV allows reconstructing the beat period of a freely-swimming human sperm (flagellum length: 50 µm), it is impractical to relate the beat pattern to the swimming trajectory because cells quickly leave the FOV.
Here, we present a new, depth-from-defocus-based 3D reconstruction algorithm that allows to computationally achieve an unprecedented z-precision in 4-plane multifocal images as high as 0.15 µm. Thereby, we are able to use a lower magnification (32x) and a wider plane spacing (planes spanned 8.8 µm). Using this algorithm, we precisely reconstruct beads and flagella across a 21-µm depth, a FOV of up to 240 x 260 µm, and at a high volume-sampling rate (500 Hz). This represents a 3.5-fold increase of depth, a 22-fold increase of FOV, and a 2.5-fold increase of sampling speed compared to other MFI systems 15,17 . The sufficiently large FOV allows relating swimming trajectories and beat pattern of sperm without compromising precision and speed.
For sperm, depth-from-defocus methods have been applied to roughly estimate the z-position in single-plane microscopy, without considering the z-precision 29, 39 . We demonstrate that the z-precision inferred from a single plane is not equal for different z-positions and poor when the object is located near focus. We have surmounted these shortcomings by using multiple focal planes to derive the z-position and show that this allows unambiguous determination of the zposition, extensive improvement of the z-precision, approximation of the z-precision over a large depth, and significant increase in depth compared to using only one focal plane. Of note, these improvements depend on a well-designed MFI setup. To establish a good tracking precision over a high depth range, the interplane distance needs to be chosen based on the relationship between object width and z-position for an individual plane image.
Our MFI system represents a cost-effective 3D imaging device that does not rely on any special equipment and uses a commercially available adapter that is compatible with most microscopes.
The system is flexible: it allows customizing the focal distance between image planes using different adapter lenses and it can be adapted to object sizes ranging from nano-to millimeter using different objectives. Thus, our MFI system is affordable and accessible to a large scientific community. To allow broad accessibility and applicability of our computational methods, we have implemented the whole image analysis pipelines presented in this study as Our method does not rely on a particular microscope configuration, which ensures broad applicability. We envisage three additional applications.
First, although we focused here on label-free samples using dark-field microscopy, our method is directly adaptable to fluorescently-labelled samples. Dark-field microscopy with a standard condenser, as applied here, limits the application to moderate NA objectives and consequentially also limits the optical resolution. However, using oil-immersion dark-field condensers or fluorescence-based MFI, high NA objectives could be applied, improving optical resolution.
Second, our MFI-based filament reconstruction method could be applied to larger objects, for instance to study rodent whisking during active vibrissal sensing 69 .
Third, we show that the concept of combining depth-from-defocus with multiple focal planes significantly improves depth estimates. Depth estimates are not only relevant for high-speed microscopy but also for research in 3D object detection by computer vision, including autonomous driving.

Characteristics of the multifocal adapter
The multi-channel imaging device QV-2 (Photometrics, USA) was used to split the light coming from the microscope into four light paths ( Supplementary Fig. 1a). The light within each optical path was projected to a different location on the camera chip. The focal length of each optical path was varied by an additional lens located at the position designed for holding filters in the QV-2. For measurements, a set of lenses with the following focal lengths was used (in mm): f1 = ∞ (no lens), f2 = 1,000, f3 = 750, and f4 = 500 ( Supplementary Fig. 1a).

Intensity normalization across planes
Local differences of image intensity in the four focal planes (Supplementary Fig. 12a) were measured by recording an image without a specimen and generating an intensity heat-map ( Supplementary Fig. 12b). For normalization, pixel intensity was divided by the respective fractional intensity value in the heat-map ( Supplementary Fig. 12c).

Alignment of plane images
The four focal planes were aligned at each time step using the ImageJ plugin Multi-Stack Reg. 71 . An alignment matrix was determined from the first time point and applied to all subsequent frames. Because the alignment algorithm uses the first plane as a template and corrects the other planes based on the first plane, this technique also allowed to correct for minute magnification differences appearing in planes 2-4 while not in plane 1 (Supplementary Table 1,   Supplementary Table 2).
For the images of freely-swimming sperm, an image of a calibration grid was recorded ( Supplementary Fig. 2a) and processed in ImageJ by image segmentation (threshold algorithm Supplementary Fig. 2b), skeletonization 73 ( Supplementary Fig. 2c), and alignment (Multi-Stack Reg 71 ; Supplementary Fig. 2d). The output alignment matrix was used for plane alignment.

Determination of the inter-plane distances
To determine the inter-plane distances, a series of multifocal images of a calibration grid at different positions set by a piezo was acquired. For each of the four planes, it was determined the piezo position at which a reference grid was in focus. The focus position was defined as the position at which the grid image featured the highest standard deviation of pixel intensity 74 .
From the differences in focus position between neighboring planes, the inter-plane distances were determined.
For low magnification (stereomicroscope, 4x), where the piezo range was insufficient for spanning all planes, the inter-plane distance was estimated manually by measuring the required focus displacement between two adjacent planes for sharp imaging of the specimen in each plane.

Simulating images produced by the multifocal adapter
Based on the inter-plane distances and the focal lengths of the lenses in the QV-2, we performed numerical calculations of images produced by the setup. The different imaging planes produced by the MFI setup were obtained with ray-optics calculations based on the thin-lens equation.
Defocusing was simulated by convolving the point spread function (PSF) with the sharpest image that we acquired of the calibration grid. We assumed a Gaussian PSF with transverse intensity distribution: In the previous equation, ( ) = ' *1 + -

Extended depth-of-field (EDOF)
Multifocal images were processed into extended depth-of-field (EDOF) images using the complex wavelet mode 24 of the ImageJ plugin "extended depth-of-field" by "Biomedical Imaging Group", EPFL, Lausanne. We created a customized version of this ImageJ plugin to enable automated processing of time series of image stacks.  Each bead's z-position was analyzed using a java-based custom-written ImageJ plugin using the "TrackMate"-derived xy and time coordinates as a template.

Tracking beads in multifocal images
For each of the focal planes, all pixels within a radius of 10 px (about 3.4 µm) were subjected to an unconstraint 2D radial fit to retrieve the radius. The fit was developed as follows. , and is the weight If the quality of the fit r 2 > 0.8, the radius was used to infer the z-position. If the fit to the radial profile did not match the criterion in more than two planes, the bead was excluded from analysis and no z-position was determined. The z-position was determined for each remaining bead as follows. In each focal plane, a difference function of the determined bead radius and the plane's calibrated relationship between bead radius and z-position was determined. Each difference function shows two minima, which represent the two most probable z-positions (one located above and another below the focal plane). Finally, the position of the bead was unequivocally ) , where sR is the standard deviation of the calibration for different beads.
Individual bead coordinates on each frame were connected to form tracks using custom-made software written in MATLAB. The software searched for the nearest particle across frames within a distance of less than 1.5 µm.
Determining the x-,y-, and z-precision based on Brownian motion of particles where * denotes the convolution operation. Assuming isotropic diffusion and normally distributed measurement errors, it can be readily found (use the convolution theorem for this) that the variance of the resulting distribution sM 2 is characterized by the sum of the variances: where the standard deviation sp characterizes the precision of the method, and sD 2 the variance of the bead displacement due to diffusion. For each dimension of space, the corresponding diffusive variance can be calculated provided that the diffusion coefficient D is known: For a spherical particle at high Reynolds numbers, the diffusion coefficient is well described by the Stokes-Einstein relation: where KB is the Boltzmann constant, T the temperature (295 K), R the particle radius (250 nm), and h(T) the dynamic viscosity of the embedding fluid (ca 0.95·10 -3 kg·m -1 ·s -1 for water at 295 K). Combining Eqs. 2-4, we find that the precision of the method to be: Using Eq. 5 with a frame period Dt of 2 ms, we obtain a z precision of 150 nm (N = 82 beads recorded for a mean duration of 3.3 s), which is in agreement with the precision found during the calibration procedure (Fig. 2g). Additionally, the precision in the XY plane can be estimated by this procedure to be 45 and 37 nm for the x and y axis, respectively.

Reconstruction of the flagellar beat
To reconstruct the flagellar beat, we developed a fully automated ImageJ plugin written in Java

Data availability
Exemplary high-speed video recordings to run the analysis routines, extracted data, and analysis supporting this study are freely accessible via the GitHub repository https://github.com/hansenjn/MultifocalImaging-AnalysisToolbox. Additional data that support the findings of this study are available from the corresponding author upon request.

Code and software availability
The developed ImageJ plugins along with the source code, manual, and example data are freely accessible on GitHub (https://github.com/hansenjn/MultifocalImaging-AnalysisToolbox). Like ImageJ, all the developed ImageJ plugins are open-source.