Helical and oscillatory microswimmer motility statistics from differential dynamic microscopy

TL;DR

This paper develops a theoretical framework to analyze the three-dimensional motion of micro-swimmers using Differential Dynamic Microscopy data, enabling extraction of motility parameters without advanced 3D imaging.

Contribution

The authors derive ISF expressions for helical and oscillatory swimmers and demonstrate how to extract 3D motility parameters from 2D DDM data, simplifying experimental analysis.

Findings

Derived ISF models for helical and oscillatory swimmers

Validated the framework with simulated swimmer data

Showed how to obtain 3D motility parameters from 2D DDM data

Abstract

The experimental characterisation of the swimming statistics of populations of microorganisms or artificially propelled particles is essential for understanding the physics of active systems and their exploitation. Here, we construct a theoretical framework to extract information on the three-dimensional motion of micro-swimmers from the Intermediate Scattering Function (ISF) obtained from Differential Dynamic Microscopy (DDM). We derive theoretical expressions for the ISF of helical and oscillatory breaststroke swimmers, and test the theoretical framework by applying it to video sequences generated from simulated swimmers with precisely-controlled dynamics. We then discuss how our theory can be applied to the experimental study of helical swimmers, such as active Janus colloids or suspensions of motile microalgae. In particular, we show how fitting DDM data to a simple, non-helical ISF model can be used to derive three-dimensional helical motility parameters, which can therefore be obtained without specialised 3D microscopy equipment. Finally, we discus how our results aid the study of active matter and describe applications of biological and ecological importance.

Equations

$$f(q,\tau )=⟨e^{i\mathbf{q}\cdot \Delta {\mathbf{r}}_j(t+\tau )}⟩$$

$$\mathbf{r}(\tau )={\mathbf{r}}_c+\delta {\mathbf{r}}_b,$$

$${\mathbf{r}}_c=v_p\tau {\mathbf{e}}_z+R{\mathbf{e}}_r,$$

$$\mathbf{p}={\mathbf{v}}_p/v_h=(v_p{\mathbf{e}}_z+{\omega }_{h}R{\mathbf{e}}_{\psi })/v_h,\text{mbox}where{v}_h=\sqrt{v_p^2+({\omega }_{h}R{)}^2}.$$

$$\delta {\mathbf{r}}_b=\delta r_b\left[\frac{v_p}{v_h}{\mathbf{e}}_z+\frac{{\omega }_{h}R}{v_h}{\mathbf{e}}_{\psi }\right],$$

$$\delta r_b(\tau )=A_b\sin ({\omega }_b\tau +{\varphi }_b).$$

$$\eta \equiv \mathbf{q}\cdot \Delta \mathbf{r}(\tau )=q{v}_p\tau \cos \theta$$

$$+q{A}_b\frac{v_p}{v_h}\left[\sin \left({\omega }_b\tau +{\varphi }_b\right)-\sin {\varphi }_b\right]\cos \theta$$

$$+qR\left[\cos \left({\omega }_h\tau +{\varphi }_h\right)-\cos {\varphi }_h\right]\sin \theta$$

$$-q{A}_b\frac{{\omega }_{h}R}{v_h}\left[\sin \left({\omega }_h\tau +{\varphi }_h\right)\sin \left({\omega }_b\tau +{\varphi }_b\right)-\sin {\varphi }_b\sin {\varphi }_h\right]\sin \theta ,$$

$$f(q,\tau )=⟨e^{i\eta }⟩\equiv \frac{\int \int d{\varphi }_{h}d{\varphi }_b{\int }_{{\mathbf{\Gamma }}_{\mathbf{s}}}P(v_p)e^{i\eta }{d}^3{\mathbf{v}}_p}{\int \int d{\varphi }_{h}d{\varphi }_b{\int }_{{\mathbf{\Gamma }}_{\mathbf{s}}}P(v_p)d^3{\mathbf{v}}_p},$$

$$f(q,\tau )=\frac{1}{8{\pi }^2}{\int }_0^{2\pi }d{\varphi }_b{\int }_0^{2\pi }d{\varphi }_h{\int }_0^{\pi }\sin \theta d\theta {\int }_0^{\infty }{P}_s(v_p)e^{i\eta }d{v}_p,$$

$$\overline{v_h}={\int }_0^{\infty }[v_p^2+({\omega }_{h}R{)}^2{]}^{1/2}{P}_s(v_p)d{v}_p,$$

$${\sigma }_h^2=\overline{v_h^2}-{\overline{v_h}}^2={\sigma }_p^2+{\overline{v_p}}^2-{\overline{v_h}}^2+({\omega }_{h}R{)}^2.$$

$$P_s(v_p)=\frac{1}{Z!}{\left(\frac{Z+1}{{\overline{v}}_p}\right)}^{Z+1}{e}^{-\left(\frac{Z+1}{{\overline{v}}_p}\right)v_p}{v}_p^Z,$$

$$\eta \approx q{v}_p\tau \cos \theta +X_b\cos \left(\frac{{\omega }_b\tau }{2}+{\varphi }_b\right)-X_h\sin \left(\frac{{\omega }_h\tau }{2}+{\varphi }_h\right),$$

$$X_b(\theta ,\tau )\equiv 2q{A}_b\sin \left(\frac{{\omega }_b\tau }{2}\right)\cos \theta \left(1-ϵ\sin {\varphi }_h\tan \theta \right)$$

$$X_h(\theta ,\tau )\equiv 2qR\sin \left(\frac{{\omega }_h\tau }{2}\right)\sin \theta ,$$

$$f(q,\tau )\approx \frac{1}{2}{\int }_0^{\pi }\sin \theta d\theta {\int }_0^{\infty }{P}_s(v_p)e^{iq{v}_p\tau \cos \theta }d{v}_p$$

$$\times {\int }_0^{2\pi }{e}^{i{X}_h\sin \left(\frac{{\omega }_h\tau }{2}+{\varphi }_h\right)}{J}_0(X_b)d{\varphi }_h$$

$$f(q,\tau )\approx {\int }_0^{1}W(\chi )J_0\left[2q{A}_b\sin \left(\frac{{\omega }_b\tau }{2}\right)\chi \right]$$

$$\times J_0\left[2qR\sin \left(\frac{{\omega }_h\tau }{2}\right)\sqrt{1-{\chi }^2}\right]d\chi .$$

$$W(\chi )\equiv \frac{\cos \left[(Z+1){\tan }^{-1}\left(\Lambda \chi \right)\right]}{{\left[1+{\left(\Lambda \chi \right)}^2\right]}^{(Z+1)/2}}.$$

$$f(q,\tau )\approx {\int }_0^{1}W(\chi )J_0\left[2qR\sin \left(\frac{{\omega }_h\tau }{2}\right)\sqrt{1-{\chi }^2}\right]d\chi ,$$

$$W_1(\chi )\equiv \cos (q{v}_p^1\tau \chi ).$$

$$f(q,\tau )\approx {\int }_0^{1}W(\chi )J_0\left[2q{A}_b\sin \left(\frac{{\omega }_b\tau }{2}\right)\chi \right]d\chi ,$$

$$f(q,\tau )\approx {\int }_0^{1}W(\chi )d\chi =\frac{1}{Z\Lambda }\frac{\sin (Z{\tan }^{-1}\Lambda )}{(1+{\Lambda }^2{)}^{Z/2}}.$$

$$f(q,\tau )\approx \mathrm{sinc}(q{v}^1\tau ).$$

$${\dot{\mathbf{r}}}_j=v_{h,j}(t){\mathbf{p}}_j,\text{mbox}{\dot{\mathbf{p}}}_j={\omega }_h{\mathbf{n}}_{0j}\times {\mathbf{p}}_j,$$

$$v_{h,j}(t)=v_{h,j}+{\omega }_b{A}_b\sin {\omega }_{b}t,$$

$$P(v_{p,j})=\frac{1}{\Gamma (\alpha )}\frac{1}{{\beta }^{\alpha }}{v}_{p,j}^{\alpha -1}{e}^{-v_{p,j}/\beta },$$

$$g(q,\tau )=A(q)\left[1-f(q,\tau )\right]+B(q),$$

$${\omega }_{h}R=\sqrt{({\sigma }_h^2-{\sigma }_p^2)+({\overline{v_h}}^2-{\overline{v_p}}^2)}.$$

Full text

• Manuscript draft date:

1 Introduction

The behaviour of swimming microorganisms and artificially propelled microscopic particles, collectively termed “microswimmers,” is of both fundamental and practical interest. On the one hand, suspensions of microswimmers reveal qualitatively distinct statistical mechanics [1, 2, 3] and fluid mechanics [4, 5, 6, 2, 3] from those of passive colloids. On the other hand, a quantitative understanding of microswimmer dynamics opens up exciting new possibilities for active material engineering [3, 7] and microbial biotechnologies [8], as well as for microbe-dependent environmental and climate science, such as in marine plankton population dynamics [9, 5].

Microorganisms and artificial swimmers typically are observed to swim along trajectories of a helical nature [10, 11, 12]. The helical motion arises from body rotation about an axis that differs from the swimming direction [13, 11] and is the expected outcome of systems that lack perfect symmetry. For biological swimmers such as microalgae [10, 14] or spermatozoa [12] the rotation originates from torques caused by non-planar flagellar motion [13, 15]. In the case of microalgae, the rotation is observed to have a biological role in that it allows cells to sample the light environment and move towards regions that are photosynthetically optimal (phototaxis) [16, 17]. For artificial microswimmers such as catalytic Janus particles the rotation is likely due to a combination of body and coating imperfections [11].

Measurements of the physical characteristics of swimmers permit the parameterization and improvement of theoretical models of active matter, which, if successfully predictive, can be used for innovative (bio)engineering design [8, 3]. Significantly, even in the absence of a theoretical framework, measurements of motility statistics also allow one to make direct inferences on the biological, ecological and biotechnological behaviour of microswimmers. For example, motility patterns, including helical swimming, and associated motility statistics change when heterotrophic microalgae prey on smaller phytoplankton [10]. Therefore, it is crucial to be able to characterise three-dimensional swimming motions for statistically significant population sizes.

Single-particle tracking [18] and ensemble-averaged techniques, such as dynamic light scattering (DLS) [19], have been until recently the main techniques used to probe the spatio-temporal dynamics of particle suspensions. Particle tracking in video microscopy, developed for the characterisation of passive colloidal dynamics [18], has been widely applied to study active systems, biological [20, 15] and synthetic [11, 21]. However, the characterisation of three-dimensional motion, such as helical swimming, with standard imaging microscopy is limited by the tracking depth of the microscope [18]. Specialised microscopy apparatus and image processing algorithms are required to extract three-dimensional motility information, such as multiple cameras [22, 23, 14, 24], exploitation of optical phase information in phase-contrast microscopy [25], digital holographic microscopy [10, 26] or ‘Lagrangian microscopes’ [27, 28, 29]. However, these techniques can be limited either in statistical accuracy and/or to low particle concentrations. For example, the use of multiple cameras is limited to relatively dilute samples because cross-correlation of camera outputs becomes challenging at high concentrations (due to particle trajectory overlaps) [14].

In DLS, fluctuations in the light scattered from a sample are collected in the far field at a given scattering vector $\mathbf{q}$, and analysed to infer the microscopic dynamics. DLS delivers statistical information for the dynamics of colloidal samples in three dimensions and plays a crucial role in the study of passive soft matter [19]. Importantly, DLS allows one to measure the Intermediate Scattering Function (ISF), also known as the Dynamic Structure Factor. The ISF is the Fourier component of the probability density function of particle displacements at a given time. It thus encodes full statistical information about the particle dynamics at a given length scale $l=2\pi/q$, with $q=|\mathbf{q}|$, and delay time $\tau$. While the potential for the application of DLS to suspensions of motile microorganisms was recognised early (e.g. for bacteria [30], microalgae [31, 32] and sperm [33]), standard DLS is restricted to large scattering angles, corresponding to small length-scales or large $q=|\mathbf{q}|$ values. At these small length-scales, many processes such as ballistic, rotational, and oscillatory motions all contribute to the ISF. Thus extracting dynamical information becomes impractical. Therefore, DLS is not suitable for the study of microswimmers. In particular, because of its limitation to motions on small scales, DLS studies [34, 35, 33] did not succeed in obtaining information about helical swimming trajectories, as this requires probing larger scale dynamics.

The discovery and development of differential dynamic microscopy (DDM; see [36] for a review on recent developments) has made the dynamics of active systems amenable to being probed by standard imaging microscopy [37, 38]. DDM yields the ISF and is particularly suited to low optical resolution imaging microscopy with a large field of view, thus giving access to the particle dynamics over large length-scales, i.e. more than one order of magnitude larger than DLS. Additionally, DDM is not restricted to low particle concentrations [39] and thus can more easily provide statistically significant information for dense microswimmer suspensions.

DDM has been applied to a range of microorganims to extract key motility parameters (including bacteria [37, 38, 40], algae [38] and spermatozoa [41]). For bacteria, DDM has also been used to clarify the interaction between motile and non-motile cells [42], characteristics of swimming in a polymer solution [43], and dynamics of concentrated suspensions [44]. Also, it has been employed to study biological active matter (e.g. [39, 45]) and artificial swimmers in quasi-two-dimensional geometries [21]. However, despite the fact that many of the above artificial and biological swimmers swim helically, to the best of our knowledge, no theoretical expression for the ISFs of helical swimmers has been derived to allow the use of DDM to study their full motion.

In this work we derive the ISF for swimmers with helical trajectories combined with progressive back-and-forth body motion. The latter is included so that the ISF can be used to describe biflagellate algae, such as the model species Chlamydomonas reinhardtii, which propel themselves by beating flagella with a breaststroke motion [46]. We then derive approximations to the ISF that facilitate extraction of helical and breaststroke swimming statistics from DDM data. The accuracy of these approximations is assessed with video sequences generated from simulated microswimmers. Finally, we discuss how our analysis suggests a new method to extract helical swimming parameters using DDM with standard microscopy setups and simple ISF models. Our method should allow the experimental study of more concentrated suspensions of active swimmers than afforded by current 3D methods.

2 Theory: intermediate scattering function and approximations

The ISF for independent (non-interacting) swimmers is given by:

[TABLE]

where $\Delta{\bf r}_{j}(t+\tau)={\bf r}_{j}(t+\tau)-{\bf r}_{j}(t)$ is the displacement of swimmer $j$, $\tau$ the delay time, q is the wavevector, with magnitude $q=|\mathbf{q}|$, probing the dynamics at a length-scale $l=2\pi/q$, and angled brackets denote averages over time $t$, and all swimmers and wavevector direction $\mathbf{q}/q$. The position vector of a helical swimmer with an oscillating back-and-forth component is given by

[TABLE]

where ${\bf r}_{c}$ is the position of the centre around which the particle centre oscillates back-and-forth, with a diplacement $\delta{\bf r}_{b}(\tau)$. Here and henceforth, we omit the swimmer index $j$ (e.g. ${\bf r}_{c,j}\to{\bf r}_{c}$) for clarity and with a view to later replacing sums over swimmers with integrals. The tip of the vector ${\bf r}_{c}$ traces a helical path, and $\delta{\bf r}_{b}$ models a back-and-forth oscillatory motion, such as the body rocking that results from the back-and-forth swimming of biflagellate algae like Chlamydomonas spp. [46], see figure 1. For each swimmer, we consider a reference frame with the $z$ axis coinciding with that of the helix, so that, adopting the cylindrical polar coordinate system $(r,\psi,z)$ shown in figure 1(b), we can decompose the helical motion as a superposition of translation along the helix axis and rotation around it

[TABLE]

where $v_{p}$ is the progressive speed, the projection of the cell velocity on the helical axis, and $R$ is the helical radius. The unit vector ${\bf e}_{r}$ rotates around the helical axis with angular speed $\omega_{h}=2\pi f_{h}$ (with $f_{h}$ the helical frequency), described by the azimuthal coordinate $\psi=\omega_{h}\tau+\phi_{h}$, where $\phi_{h}$ is a random phase, uniformly distributed in the interval $[0,2\pi]$, added to ensure the helical rotations of different cells are not synchronised. The back-and-forth motion is along the instantaneous swimming direction

[TABLE]

This is the speed along the helical trajectory. The back-and-forth displacement can then be written as

[TABLE]

where

[TABLE]

Here, $A_{b}$ and $\omega_{b}=2\pi f_{b}$ are the back-and-forth oscillatory amplitude and angular speed (with $f_{b}$ the back-and-forth oscillatory frequency), respectively, and $\phi_{b}$ a random phase added to avoid synchronisation, as for the helical motion.

With these approximations, the phase contribution to the ISF due to swimmers can be shown to be given by (see A)

[TABLE]

where $\theta$ is the angle between ${\bf q}$ and the helical axis along $z$. The first term in equation (7) is the contribution of progressive swimming, the second term corresponds to breastroke swimming oscillations, the third to helical rotation, and the forth is a cross-term reflecting the coupling between breastroke and helical swimming.

To evaluate the ISF given by (1) we need to average over all swimmers, which means averaging over all distributions of swimming parameters. To make analytical progress, henceforth we assume isotropic swimming, and consider a distribution $P(v_{p})$ only for swimming speeds: all other swimming parameters are approximated as having single values. The ISF for the algae is then a multiple integral over the volume of the swimming velocity space ${\bf\Gamma_{s}}$, as well as over the random phases $\phi_{h}$ and $\phi_{b}$

[TABLE]

where normalisation ensures $f=1$ for $q=0$ or $\tau=0$. Upon substituting (7) into (1) and choosing spherical polar coordinates for the velocity space integration we obtain

[TABLE]

where the isotropic speed distribution is defined as $P_{s}=4\pi P(v_{p})v_{p}^{2}$.

2.1 Speed distribution transformations

The reader will have noticed that the averaging over speed in the ISF integral (9) is evaluated using $P_{s}(v_{p})$, the distribution of progressive speeds $v_{p}$ (projected onto the helical axis), and not $P_{s}(v_{h})$, the distribution of along-helix speeds. In fact, either distribution can be used, as a simple change of variables in the integral demonstrates: $\int_{\omega_{h}R}^{\infty}h(v_{h})P_{s}(v_{h})dv_{h}=\int_{0}^{\infty}h(v_{h}(v_{p}))P_{s}(v_{p})dv_{p}$, where $h$ is a general function. This same result implies we can derive the along-helix motility statistics if we know the distribution of progressive speeds. Indeed, from equation (4) it follows that the mean along-helix speed average is given by

[TABLE]

where here and henceforth overbars denoted averages over the speed distribution. The speed average in equation (10) provides a value of the along-helix speed given values of $\omega_{h}$, $R$, and the distribution $P_{s}(v_{p})$. Similarly, recalling $P_{s}$ is normalised, it is easily shown that the second moment of the along-helix speed is given by $\overline{v_{h}^{2}}=\overline{v_{p}^{2}}+(\omega_{h}R)^{2}$. Using this relation and $\sigma_{p}^{2}=\overline{v_{p}^{2}}-\overline{v_{p}}^{2}$, the variance of the progressive speed distribution, the variance of the along-helix speed distribution can then be written as

[TABLE]

We will use transformations (10) and (11) in the analysis of our simulated results, and in particular to deduce three-dimensional swimming parameters from two-dimensional data. These relations apply generally for swimmers with swimming distributions with finite first and second moments. In what follows, as in [37, 38], we consider swimmers with a Schulz distribution [47] of progressive swimming speeds, such that

[TABLE]

where $\overline{v}_{p}$ is the progressive mean speed and $Z=(\overline{v}_{p}/\sigma_{p})^{2}-1$ is a parameter related to this mean and the variance of the progressive speed distribution. The Schulz distribution (12) is identical to the gamma distribution, modulo a simple transformation of parameters (see section 3.1). Physically it displays the correct general features (going through the origin and peaked), and mathematically it allows us to make analytical progress.

2.2 Approximate ISF expressions

We summarise here limiting expressions of the ISF (9), which will be employed in section 4 to analyse simulations of helical swimmers and aid interpretation of experiments with real helical swimmers. Deriving these from the full ISF is straightforward, as demonstrated in the next section for the case when back-and-forth and helical swimming motions are uncoupled.

2.2.1 Uncoupled back-and-forth and helical swimming motions.

The swimming of realistic swimmers, such as the biflagellate algae Chlamydomonas reinhardtii and Dunaliella salina, displays both back-and-forth and helical swimming motions. In this case, the full ISF (9) is required, which is hard to simplify analytically much further (integration is straightforward over either, but not both, random phases). For helical swimmers like C. reinhardtii that move along tight helices with well-separated helical ($1/\omega_{h}$) and back-and-forth ($1/\omega_{b}$) motion timescales [48, 31, 32], the ISF (9) can be written in terms of the non-dimensional small parameters $\epsilon:=\omega_{h}R/v_{p}$ and $\nu:=\omega_{h}/\omega_{b}$, so that (see B)

[TABLE]

where we have defined

[TABLE]

and

[TABLE]

and we have used the trigonometric identity $\cos(A)-\cos(B)=-2\sin[(A+B)/2]\sin[(A-B)/2]$ and similarly for $\sin(A)-\sin(B)$. Substituting equation (13) into the ISF (9), we can integrate over the random phase $\phi_{b}$ to obtain:

[TABLE]

where $J_{0}(X_{b})=(1/2\pi)\int_{0}^{2\pi}e^{iX_{b}\cos\left(\frac{\omega_{h}\tau}{2}+\phi_{b}\right)}d\phi_{b}$ is the zeroth order Bessel function of the first kind [49]. As shown in B, if $qA_{b}\epsilon$ is small, a condition always met at large scales where $qA_{b}\ll 1$, it is possible to further integrate the ISF over $\phi_{h}$. For swimmers with a Schulz distribution (12), we can also integrate over swimming speeds $v_{p}$ to obtain the ISF

[TABLE]

where we have changed variables to $\chi=\cos\theta$, and defined a ‘ballistic kernel’ function

[TABLE]

We have given it this name because integrating over $\chi$ provides the ballistic model [38], see Equation (22). We have also used definitions (14, 15) to make swimming parameter dependencies evident and defined the constant $\Lambda=q\overline{v}_{p}\tau/(Z+1)$, where we recall $Z=(\overline{v}_{p}/\sigma_{p})^{2}-1$ is a parameter from which the standard deviation $\sigma_{p}$ of the progressive speed distribution can be obtained.

2.2.2 Pure helical swimming (negligible breastroke motion).

Some microswimmers, such as dinoflagellates [10], are known to have negligible or nonexistent oscillatory back-and-forth motion ($A_{b}\approx 0$ for all $q$ values). In this case, since $J_{0}(0)=1$, the ISF (2.2.1) simplifies to

[TABLE]

where all quantities have been defined above. In what follows, to understand the effect of helical swimming on the ISF, we also consider swimmers with speed distribution $P_{s}(v_{p})=\delta(v_{p}-v^{1}_{p})$, where $v^{1}_{p}$ is the single progressive speed possessed by all swimmers. In this case, the ISF is still given by (19), but with the substitution $W\to W_{1}$, where

[TABLE]

is the single speed ballistic kernel.

2.2.3 Pure back-and-forth swimming (negligible helical motion).

It is also of interest to consider the case of negligible helical motion ($R\approx 0$ for all $q$). In this case the ISF (2.2.1) reduces to

[TABLE]

where $W$ is given by equation (18). However, we note that, in the absence of helical motion, $\Lambda=q\overline{v}\tau/(Z+1)$, where $\overline{v}$ is the average swimming speed and $Z=(\overline{v}/\sigma)^{2}-1$ provides the standard deviation $\sigma$. The distinction between progressive and along-helix speeds is not relevant for this model: there is only one speed, which we denote by $v$.

2.2.4 Ballistic swimming (negligible back-and-forth and helical motion).

Finally, we consider the limit of ballistic motion, where both back-and-forth and helical motion are neglected. It is then possible to integrate the ballistic kernel to obtain

[TABLE]

This ISF, just as the one for pure breastroke motion (21), provides the mean speed $\overline{v}$ and standard deviation $\sigma$. When applied to particle dynamics that are not purely ballistic, Eq. (22) yields effective quantities, that are scale-dependent, as we shall see later, and has been previously used to model the swimming of bacteria [37] and algae [38].

For swimmers with a single speed, as for the pure helical swimming case, we can make the substitution $W\to W_{1}$, where $W_{1}$ is given by equation (20) but with a speed $v^{1}$ replacing the progressive one. Then equation (22) integrates to

[TABLE]

3 Numerical methods and DDM analysis

We performed simulations of swimmers undergoing helical and/or back-and-forth motions (in 3D). These were carried out to demonstrate that it is possible to recover swimming parameters from the simulated ISFs using DDM analysis and the ISF models just presented. ISFs were obtained by applying such analysis to ‘microscopy-like’ pseudo-image sequences derived from the coordinates of simulated swimmers, as detailed below.

3.1 Simulations and pseudo-image sequence generation

We use an individual based model, adapted from [13], to simulate the dynamics of swimmers with both helical and back-and-forth motions. Each swimmer $j$ obeys the dynamical system

[TABLE]

where $\mathbf{r}_{j}$ is the position of the centre of oscillation of the $j$-th swimmer and $\mathbf{p}_{j}=\mathbf{v}_{h,j}/v_{h,j}$ its swimming direction, which rotates around the direction $\mathbf{n}_{0,j}$ (the helical axis direction for swimmer $j$, see figure 1). Dots denote time derivatives. The speed along the helix is given by

[TABLE]

a superposition of a linear translation with constant speed $v_{h,j}$ and a sinusoidal oscillation that models back-and-forth motion, e.g. for swimming algae. This is analogous to what we did for the derivation of the ISF model. Using Cartesian coordinates for $\mathbf{r}_{j}$ and sphericals for $\mathbf{p}_{j}$, equation (24) is expanded in the system of component ordinary differential equations, and these were solved numerically with MATLAB using swimming parameters realistic for algae such as C. reinhardtii (see C for details). As in the ISF model, all of these parameters were assumed to have a single value, with the exception of the swimming speeds. These were obtained from the distribution of speeds projected along the helix, $v_{p,j}$, provided by the Gamma distribution

[TABLE]

where $\alpha=Z+1$ and $\beta=\overline{v}_{p,j}/(Z+1)$ are the distribution parameters (substituting these into (26) yields the Schulz distribution (12)). Using relation (4) we then obtain the speed along the helix $v_{h,j}=\sqrt{v_{p,j}^{2}+(\omega_{h}R)^{2}}$, which can be substituted into (24) for each swimmer. The swimming direction, $\mathbf{p}_{j}$, and beat phase are initialised randomly (again, see C for details), while the initial values for the helix angle $\gamma_{j}$ between $\mathbf{p}_{j}$ and $\mathbf{n}_{0,j}$ can be obtained from the initial along-helix, $v_{h,j}$, and progressive, $v_{p,j}$, speeds. This is given by $\mathbf{p}_{j}\cdot\mathbf{n}_{0,j}=v_{p,j}/v_{h,j}=\cos\gamma_{j}$, and needs to be prescribed for each simulated swimmer, as it is a function of speed. Image sequences of the simulated helical swimmers were constructed from the simulations by generating a Gaussian pseudo-diffraction spot for each point position $(x_{i},y_{i},z_{i})$, as detailed in C.

3.2 Differential dynamic microscopy analysis

For each simulated movie, we calculate the differential image correlation functions (DICF), $g(\vec{q},\tau)$, i.e. the power-spectrum of the difference between pairs of images delayed by time $\tau$. Under suitable conditions, such as isotropic motion, $g(q,\tau)=\left<g(\vec{q},\tau)\right>_{\vec{q}}$ is related to the ISF $f(q,\tau)$, the $q^{\mathrm{th}}$ mode of the number density autocorrelation function, also called the normalised dynamic structure factor:

[TABLE]

with $B(q)$ the instrumental noise and $A(q)$ the signal amplitude. Details can be found elsewhere [38]. For a given simulated movie, the $g(q,\tau)$s were fitted as a function of $\tau$ either for each $q$ independently or simultaneously over a given $q$ range. A range of models were used. For simplified models that ignore the helical path, we found that applying a fitting weight at longer times allowed better recovery of the input parameters. Details are mentioned in each figure caption when relevant. For practicality, we present directly ISFs extracted from equation (27) using the fitted parameters $A(q)$ and $B(q)$.

4 Results

4.1 Pure helical swimming

We start by analysing simulations of swimmers undergoing pure helical swimming (no back-and-forth motion). To avoid the complications of swimming speed heterogeneity, we consider first a population of swimmers swimming with a single speed. Typical trajectories are shown in Figure 2. Figure 3a shows the ISFs as a function of delay time $\tau$ at several values of $q$ for a population of swimmers with a single along-the-helix speed $v^{1}_{h}=120$ $\mu$m/s, helical radius $R=2$ $\mu$m, and helical frequency $f_{h}=2$ Hz. The decay is qualitatively different depending on the value of $q$, which defined the length scale of interest. Indeed, equations (19-20) suggest that the characteristic times of ballistic and helical contribution scale as $\tau_{b}\sim 2\pi/qv_{p}$ and $\tau_{h}\sim 2/f_{h}$, respectively, so that there exists a crossover $q_{c}\sim\pi f_{h}/v_{p}\approx 0.1$ $\mu$m*-1* for which $\tau_{b}=\tau_{h}$. Thus, helical motion should provide the fastest contribution to the dynamics for $q\ll q_{c}$, and the slowest at $q\gg q_{c}$.

For small $q$, i.e. $q\ll q_{c}$ (large length-scales), the ISF displays two main decays, corresponding to ballistic and helical contributions to the swimming motion. In particular, a kink in the ISF can be seen at time $\tau\simeq 0.5$ s, independent of $q$, for $q\lesssim 0.1$ (see vertical dashed line in figure 3a), which corresponds to the characteristic period $1/f_{h}=0.5$ s of helical swimmers rotating at $f_{h}=2$ Hz (the simulation input frequency). At large values of $q$, i.e. $q\gg q_{c}$, the ISF shows a single decay, which fully decorrelates on time-scales $\tau\ll\tau_{h}$. This is characteristic of ballistic motion: large $q$ values correspond to small length scales, where helical trajectories are ballistic to a good approximation, so helical features in the ISF are not pronounced (or easily distinguished). Additionally, the amplitude of the helical process scales with $q$ and thus decreases with $q$ as expected from the $J_{0}$ term in equation (19).

The effect of helical rotation on the ISF, via the observation of a kink, is strongly sensitive to the helical swimming parameters. For example, the characteristic kink evident for single speed swimmers (figure 3(a)) is less apparent in the following instances when the radius $R$ of the helical path is decreased (figure 3b), or when swimmers obey a speed distribution with non-zero standard deviation (figure 3 (c) and (d)).

The signature of helical motion is nevertheless encoded in the swimmer ISF. We demonstrate in the next sections that fitting the ISFs with appropriate models derived in section 2 allows to measure the 3D motility parameters, with good accuracy, from 2D movies of the swimmers.

4.1.1 Motility parameters from the simulated swimmer ISF: single speed swimmers.

We first analyse the ISF of a population of helical swimmers with one swimming speed. The ISF is fitted using equation (19) with $W$ given by equation (20). Fitting for each $q$ independently, we can infer the swimming parameters as a function of $q$. Figure 4 displays the $q$ dependence of the progressive single swimming speed $v_{p}^{1}$, helix radius $R$ and frequency $f_{h}$ (panels a-c respectively). We see that parameters are recovered well with better results at large $R$. For small $R$, features of helical motion in the ISF are not pronounced (e.g. see figure 3(b) and (c) for $R=2$ $\mu$m), resulting in less accurate parameter recovery. However, global fits over given ranges of $q$ permit accurate recovery of input values, see dotted lines in figure 4 and figure caption.

It is also interesting to analyse the ISF of helical swimming with the simpler ballistic model (23), which does not account for helical swimming motions. The ballistic model yields a speed that varies with $q$ for non-zero $R$ (figure 4(d)). The ballistic model does not account for helical motion and thus provides an effective speed, with a meaning that depends on length scale. For small length scales (high $q$), it represents the along-helix speed $v^{1}_{h}$, matching the input values to $\leq 3\%$, where the effects of helical motion are not appreciable. At large length scales ($q\ll q_{c}$), the speed tends towards the progressive speed which, rearranging equation (4), is given by $v^{1}_{p}=\sqrt{(v^{1}_{h})^{2}-(\omega_{h}R)^{2}}$. We recall this is the projection of the along-helix speed onto the helix axis: at large length scales, helical motion is averaged out over several helical pitches.

4.1.2 Motility parameters from the simulated swimmer ISF: Schulz speed distribution.

We show in figure 3 typical ISFs for helical swimmers with a Schulz swimming speed distribution. As mentioned, a distribution of speeds makes helical features in the ISF more subtle and they qualitatively look similar to swimmers with pure ballistic motion ($R=0$) (figure 3(d)). Here, we fit the ISFs with a model for helical swimmers with a Schulz swimming speed distribution given by equation (12), using (19) for the ISF with $W$ given by the ballistic kernel (18). For smaller values of the helix radii, we found that global fitting the ISFs provides more reliable results than fitting each $q$ independently (data not shown). Figure 5(a) shows the ratio of the swimming parameters extracted from the global fit to simulation input values, as a function of the helix radius $R$. The parameters are denoted collectively as $X=\overline{v}_{p},\sigma,R,f_{h}$, the swimming speed, variance, helix radius and frequency, respectively (corresponding input values are denoted by $X_{input}$). It is clear from figure 5(a) that most parameters are recovered with good accuracy (within $\approx 10\%$) for all $R$. This demonstrates the applicability of the helical ISF model for swimmers with a Schulz distribution and shows how this model is sensitive to helical signatures not evident from inspection of the ISF (see figure 3(b)). At lower values of $R$, while $\overline{v}_{p}$ and $f_{h}$ are recovered to within a few percent across the range, the accuracy of recovery for $R$ is less accurate. This is because, at such small $R$ values, the helical signature in the ISF is too weak for even the global fit to pick up.

It is useful to also analyse the Schulz swimmer ISF using the ballistic model, as we did for single speed swimmers. In this case we can fit the ISF with equation (22) independently for each value of $q$. The results for the effective mean swimming speed $\overline{v}$ and distribution standard deviation $\sigma$ are shown in figure 5(b) and (c), respectively, as a function of $q$ for different values of $R$. As for the single speed case, fitting the ballistic model to the helical swimmer ISF results in a $q$-dependent speed $\overline{v}(q)$, whose value increases from its progressive value $\overline{v}_{p}$ at low $q$ to its along-helix value $\overline{v}_{h}$ at high $q$. A similar trend is observed for $\sigma$ although the data can be noisier.

As is clear from Figure 5b, the mean along-helix speed value grows with the helix radius $R$. The functional form for this variation, provided by equation (10), is explicitly plotted in Figure 5d and compared with the ratio between the high $q$ value of the mean speed from fits and the simulation input. Similarly, the theoretical and fitted ratios of the corresponding standard deviations, using equation (11) are plotted on the same graph. The speed and standard deviation data follow the expected variation with $R$, which provides a useful check that the dependencies predicted by our analytical theory are borne out by the simulated experiments.

4.2 Pure back-and-forth swimming

Next, we analyse simulations of swimmers with only back-and-forth motion (no helical swimming) and a Schulz distribution of swimming speeds using typical swimming parameters for biflagellate microalgae, like Chlamydomonas spp. or Dunaliella spp.. The ISF obtained from these simulations is shown for different values of $q$ in figure 6(a). As for helical swimming, a $q$-independent signature of back-and-forth motion is evident. This is oscillatory in time and its amplitude decreases with decreasing back-and-forth amplitude $A_{b}$, vanishing for $A_{b}=0$, see figure 6(b), as expected from the $J_{0}$ term in equation (21). The short-time decay of the ISF ($\tau\lesssim 0.02$ s) corresponds to the back-and-forth swimming oscillation, while longer time decay corresponds to ballistic motion. Subsequent oscillations at longer times in the decay, before the noise floor is reached, reflect oscillations in the Bessel function $J_{0}$ (see equation (21)). In contrast to the helical case, the back-and-forth signature in the ISF is more evident at high $q$ values, because of its small length-scale amplitude.

4.2.1 Motility parameters from the simulated swimmer ISF.

Analogously to the helical case, fitting the simulated back-and-forth swimmer ISF with equation (21) provides values for the swimming parameters. Figure 7(a)-(c) shows the variation with $q$ of the mean swimming speed $\overline{v}$ and standard deviation $\sigma$ (inset), back-and-forth amplitude $A_{b}$ and frequency $f_{b}$, respectively. For small values of $A_{b}$, the values of $\overline{v}$ and $\sigma$ deviate from the simulation inputs at low $q$: here, the back-and-forth signal in the ISF for small $A_{b}$ is weak, and thus harder to fit (similarly to the small $R$ case for helical swimmers). Unsurprisingly, on the other hand, since $f_{b}$ and $A_{b}$ are $q$-independent, their values can be accurately determined across the $q$ range. The effect of the magnitude of $A_{b}$ on parameter accuracy is shown in figure 7(d), where we plot the ratio of the high $q$ value of the swimming parameters to the input value. Note how the effect of $A_{b}$ on accuracy is analogous to that of $R$ for helical swimmers.

4.3 Back-and-forth and helical swimming

Finally, we consider simulations of swimmers with both helical and back-and-forth motions, and a Schulz distribution of speeds. Figure 8 shows the corresponding ISF decay for different $q$ values. The ISFs encode both pure helical and back-and-forth swimming although only back-and-forth features are clearly evident.

4.3.1 Motility parameters from the simulated swimmer ISF.

We fit the back-and-forth and helical swimming ISF with equation (2.2.1) to obtain swimming parameters. As for the pure helical case, a global fit over a selected $q$ range provides the best accuracy, with parameter value recovery $\leq 5\%$ (figure 9).

We have also fitted the ISF with the pure back-and-forth model (ignoring the helical path) given by equation (21), independently for each value of $q$. Figure 9(a) and (b) shows the mean effective swimming speed $\overline{v}$, distribution standard deviation $\sigma$, swimming amplitude $A_{b}$ and frequency $f_{b}$, respectively, as a function of $q$. We found that fitting with a stronger weight at longer times, thus ignoring short-time contributions, yielded a similar $q$ dependence we observed for the pure helical case. As in this case, this is an artefact of using a model that is unable to account for helical trajectories. As a result, $\overline{v}$ transitions from its progressive value, $\overline{v}_{p}$, at low $q$, to its along-helix value, $\overline{v}_{h}$, at high $q$, and similarly for the distribution standard deviation $\sigma$. As for the helical case, we can use equation (10) and (11) to relate along-helix and progressive values for a given helical angular frequency $\omega_{h}$ and radius $R$. More interestingly, when the helical parameters are not known (as is the case in experiment) we can estimate the product $\omega_{h}R$, a measure of the rotational speed around the helical axis. Indeed, simply rearranging equation (11) provides:

[TABLE]

Thus, along-helix (high $q$) and progressive (low $q$) values of the speed and distribution standard deviation can be used to estimate $\omega_{h}R$ based on a simple ISF model that does not account for helical motions. From figure 9, we measure the following values: $v_{p}\approx 122\mu m/s$, $\sigma_{p}\approx 24\mu m/s$, $v_{h}\approx 154\mu m/s$, and $\sigma_{h}\approx 21\mu m/s$. Thus, using equation (28), we obtain $(\omega_{h}R)\approx 95$ $\mu$m rad/s, only $6\%$ from the simulation input value. Finally, while oscillating components are only recovered towards higher $q$ values, fitting with a stronger weight at short-time allows determination of the oscillating component within $\leq 0.01\%$, see caption of figure 9.

5 Discussion and Conclusions

The experimental characterisation of the helical and oscillatory motions of microswimmers is essential for understanding their statistical mechanics [2, 3], fluid mechanics [5, 2, 3], biology [46] and ecology [10]. In this work, we developed a theoretical framework enabling the quantification of such motions by differential dynamic microscopy (DDM), a technique allowing to infer motility statistics from videomicroscopy without the need for specialised apparatus. We modelled the intermediate scattering function (ISF) corresponding to helical, back-and-forth, and combined helical and back-and-forth motions. Individual based model simulations of swimmers with theses same motions were used to generate artificial movie sequences, from which the ISFs were extracted, as in DDM [38] and fitted with our models. This study provides the basis for applying ensemble averaging techniques such as DDM to measure the motility statistics of microswimmers following a helical path or with a combination of helical and other motions.

The ISFs from our simulations allow to explore characteristic ISF features emerging from helical and back-and-forth swimmer motions. In particular, for swimmers with pure back-and-forth motion, a characteristic dip in the ISF is evident at the timescale corresponding to the oscillation (figure 6). A similar feature is evident for pure helical swimmers, corresponding to helical rotation. This is only apparent for swimmers with a single speed, and at large enough value of the helix radius. For swimmers with a speed distribution, the helical swimming feature is not conspicuous (see figure 3), because it is spread across different length-scales for swimmers with different swimming speeds.

Beyond predicting features to be expected experimentally, our simulated ISFs can be fitted using the models we have developed to test how accurately the simulation input parameters can be extracted using DDM. For all the models considered in this work, we find that it is possible to extract helical swimming statistics from simulated ISFs with good accuracy, even when qualitative features are not evident in the ISFs. We note that some of the assumptions we have made in our models (e.g. that all swimming parameters but speed are delta function distributions, different swimming motions can be decoupled) may not apply in real experimental systems. More detailed models of swimmer ISFs than considered here may be formulated, e.g. when several contributions to the ISFs from different swimming motions cannot be decoupled. In this case, however, the added details (e.g. multiple integrals, additional dynamical parameters, parameter distributions) may render fitting impractical.

In alternative to detailed models, our analysis also shows that it is also possible to extract information about complex swimming motions from simple ISF models. For example, when fitting the ISF of helical and back-and-forth swimmers with a simpler model ignoring helical motion, we find that a scale-dependent effective mean speed, $\overline{v}(q)$, and a speed distribution standard deviation, $\sigma(q)$, emerge. Their low $q$ (large lengthscale) values correspond to the mean progressive speed $\overline{v}_{p}$ and standard deviation $\sigma_{p}$, while their high $q$ (small lengthscale) values represent the mean along-helix speed $\overline{v_{h}}$ and standard deviation $\sigma_{h}$. These progressive and along-helix values can be combined to estimate the product $\omega_{h}R$ between the helical angular frequency $\omega_{h}$ and radius $R$, see equation (28), where $\omega_{h}R$ is a measure of the speed around the helical axis. Further, independent knowledge of $\omega_{h}$, obtainable e.g. from particle tracking, can provide an estimate of the (mean) helix radius $R$. For example, assuming a measured helical frequency $f_{h}=2$ Hz, we can use the value $\omega_{h}R\approx 95$ $\mu$m rad/s (estimated in the previous section from the data in figure 9) to obtain a helix radius $R\approx 7.6$ $\mu$m, which is within $5\%$ of the input value used for the simulation.

The helix radius is an intrinsically 3D property of helical swimming. Thus far, to acquire motility statistics for this parameter, specialised techniques have been used, e.g. microscopy involving multiple cameras [22, 23, 14, 50, 24] or holographic microscopy [10]. Such techniques are limited in statistical accuracy and/or to relatively dilute microswimmer systems, as explained in the introduction. Gurarie al [51] recently obtained 3D motility statistics of dilute suspensions of dinoflagellate microalgae from 2D projections of helical trajectories obtained using standard microscopy using a continuous stochastic model (CSM) of helical swimming [51]. In this study, the dinoflagellates swam with pronounced (large radius) helical trajectories, and little other body oscillations. It would be challenging, however, to apply the CSM method to swimmers with small helical radius and back-and-forth oscillations, such as C. reinhardtii. Instead, the hybrid DDM-tracking approach proposed here has great potential to enable measurements of the helix radii for moderately concentrated suspensions of helical swimmers.

The efficacy of our methods and analyses needs to be tested by applying differential dynamic microscopy to suspensions of real artificial and biological swimmers. As discussed, helical swimming motions are commonplace for many types of microswimmer [3]. Using the methods we have here developed, high quality 3D motility statistics can be collected with standard 2D-imaging microscopy. This provides economical experimental set-ups for the analysis of 3D microswimmer motility, making this important area of research accessible to researchers with a limited equipment budget, such as in schools and developing countries. The ability of our new methods to extract helical swimming parameters without the use of specialised microscopy setups could also be exploited in field work (e.g. to characterise unculturable helical swimmers, such as dinoflagellates, in situ). Future work could develop the models used in our analysis, adapting them to the large variety of interesting synthetic and biological swimmers, both currently known and to be discovered. This will allow theoretical and experimental progress in our statistical understanding of microswimmers, facilitating the development of microswimmer based biotechnologies, and aiding the characterisation of biological microswimmers in the environment.

6 Data availability.

The raw data (i.e. simulated movies) presented in this publication is available on the Edinburgh DataShare [52].

7 Acknowledgements

We thank J. Arlt for technical discussions about DDM processing, A. Brown for discussions and comments on the manuscript and acknowledge financial support from ERC-AdG-340877 PHYSAPS (V.A.M and W.C.K.P), EPSRC EP/N509620/1 (T.J.) and the Winton Programme for the Physics of Sustainability (T.J., O.A.C.).

Appendix A Phase of the full ISF

From equation (1) the phase of the ISF is defined as

[TABLE]

If we substitute equations (3), (5) and (6) from the main text into (2) then, writing ${\bf q}=q\sin\theta\,{\bf e}_{x}+q\cos\theta\,{\bf e}_{z}$, and noting ${\bf e}_{x}\cdot{\bf e}_{r}=\cos\psi(\tau)$ and ${\bf e}_{x}\cdot{\bf e}_{\psi}=-\sin\psi(\tau)$ (see figure 1), we obtain

[TABLE]

We recall that $\psi(\tau)=\omega_{h}\tau+\phi_{h}$ is the azimuthal coordinate of the swimmer. Substituting into equation (30) and setting $\tau=0$, gives ${\bf q}\cdot{\bf r}(0)=qR\sin\theta\cos\phi_{h}+qA_{b}\sin\phi_{b}\left[\frac{v_{p}}{v_{h}}\cos\theta-\frac{\omega_{h}R}{v_{h}}\sin\phi_{h}\sin\theta\right]$. Subtracting this from equation (30) we arrive at expression (7) for the ISF phase angle.

Appendix B Uncoupled helical and back-and-forth motions approximation

We define the small parameters $\epsilon:=\omega_{h}R/v_{p}$ and $\nu:=\omega_{h}/\omega_{b}$. Then, for the second and fourth terms in phase of the ISF (7) in the main text, we have: $v_{p}/v_{h}=[1+(\omega_{h}R/v_{p})^{2}]^{-1/2}=1-\epsilon^{2}/2+O(\epsilon^{4})$ and, similarly, $\omega_{h}R/v_{h}=(\omega_{h}R/v_{p})[1+(\omega_{h}R/v_{p})^{2}]^{-1/2}=\epsilon-\epsilon^{3}/2+O(\epsilon^{5})$. Further, rescaling time such that $\tau^{\prime}=\omega_{b}\tau$, $\sin\left(\nu\tau^{\prime}+\phi_{h}\right)=\nu\tau^{\prime}\cos\phi_{h}+\sin\phi_{h}+O(\nu^{2})$, so that equation (7) becomes

[TABLE]

If $\nu\epsilon\ll 1$, the above expression for $\eta$ is identical to that in equation (13) in the main text, once dimensional time is restored, and $X_{b}$ and $X_{h}$ defined in equations (14) and (15).

Next, we consider the approximation necessary to integrate (2.2.1) over $\phi_{h}$. Consider the Bessel function of the first kind $J_{0}(X_{b})$, with argument written as $X_{b}=X_{b}^{0}+\delta X_{b}(\phi_{h})$, where $X_{b}^{0}:=2\sin\left(\frac{\omega_{b}\tau}{2}\right)\cos\theta$ and $\delta X_{b}:=-\epsilon qA_{b}X_{b}^{0}\sin\phi_{h}\tan\theta$. Taylor expanding this function,

[TABLE]

where we have used $J_{0}^{\prime}=-J_{1}$. Thus, the integral over $\phi_{h}$ in (2.2.1) can be written as

[TABLE]

where we have used $J_{1}(X_{h})=(1/2\pi)\int_{0}^{2\pi}e^{iX_{h}\sin\left(\frac{\omega_{h}\tau}{2}+\phi_{h}\right)}\sin\phi_{h}d\phi_{h}$. Since $J_{1}(X_{h})|X_{b}^{0}|\lesssim 1$, contributions from the second term to the ISF integral are negligible when $\epsilon qA_{b}$ is small. If $\epsilon=O(1)$ (but still $\epsilon\nu\ll 1$, as above), the latter condition is still met at large lengthscales, where $qA_{b}\ll 1$. For the parameters used in our simulations, which correspond to C.reinhardtii-like swimmers, $\epsilon\approx 2\pi f_{h}R/\bar{v}_{p}=0.2$-$0.8$ and $\nu=f_{h}/f_{b}=0.04$, $\epsilon\nu$ is always small, making equation (B) a good approximation to the phase of the ISF (7). The condition $\epsilon qA_{b}\ll 1$ required for the approximation to the Bessel integral just discussed, however, requires small values of the helix radius $R$ to be strictly satisfied for all $q$ values. Successful recovery of helical parameters from simulations analysed using ISF expressions for swimmers with uncoupled helical and back-and-forth motions, based on the approximations just described, further justifies the appropriateness of these approximations.

Appendix C Simulation details

Expanded in component form, equations (24) become

[TABLE]

where $v_{j}(t)=v_{h,j}+\omega_{b}A_{b}\sin{\omega_{b}t}$ is the net along-helix speed including a sinusoidal contribution due to back-and-forth motions. The angles $\theta_{n_{0}}$ and $\phi_{n_{0}}$ determining the helix orientation are given by the constraint $\mathbf{p}_{j}\cdot\mathbf{n}_{0,j}=v_{j}/v_{h,j}=\cos\gamma_{j}$. This system of ODEs is solved numerically using MATLAB (Mathworks, Natick, MA, USA) with a Runge-Kutta-Fehlberg (RK45) method parallelised for $N=1000$ swimmers in a periodic box of size $L_{x}\times L_{y}\times L_{z}$, where $L_{x}=L_{y}=L_{z}=1000$ $\mu$m in figures 3, 6, 7 and 8. The box size was increased to $L_{x}=L_{y}=L_{z}=2000$ $\mu$m in figures 4, 5 and 9, where the analysis of swimming parameters requires reaching lower $q$ values. In the first case, a $500$ $\mu$m thick slice of the box was taken to generate microscopy videos, whereas the full box was used in the latter case. Microscopy-like videos were generated by assigning a Gaussian pseudo-diffraction spot to the position of each simulated swimmer in the slice volume. This is achieved by subtracting the intensity $I_{0}\exp(-((x-x_{i})^{2}+(y-y_{i})^{2})/(2\sigma^{2}))[1-(z^{o}_{i}/\delta)^{2}]$ from the background intensity $I_{B}$ for of each pixel $(x,y)$ within a cut-off distance (in $x$ and $y$ only) of each alga, where $z^{o}_{i}$ is the offset of the algae from the centre of the simulation box in z-direction. Here $I_{0}$ is the Gaussian peak intensity, $\sigma$ its standard deviation in the plane and $\delta$ its extent in depth. The following image parameters were used in all simulations: $I_{B}=255$, $I_{0}=50$, $\delta=L_{z}/2$ and $\sigma=1$. The videos simulated capture at $500$ Hz, resulting in sequences of duration $15.32$ s ($L_{z}=1000$ $\mu$m) or $32$ s ($L_{z}=2000$ $\mu$m). As in the ISF model, the simulations neglect orientational noise and other realistic effects discussed in [15].