Mathematical model
We developed a basic mechanistic model to capture the movement pattern of the self-propelled diatom cells at an individual level in a two-dimensional space. Adapted from the motion behavior of self-propelled non-living micro-rods (34), our discrete time model uses the abovementioned 5 movement parameters, assuming white noises with intensity \(D_{r}\) and \(D_{\theta}\) in translational diffusivity and rotational diffusivity respectively. The stochastic difference equations for the movement of a single diatom cell are given by
\(x\left(t+\Delta t\right)-x\left(t\right)=\kappa(t)V_{0}\cos{\theta(t)}\Delta t+\sqrt{2D_{r}\text{Δt}}\xi_{1},\)(1a)
\(y\left(t+\Delta t\right)-y(t)={\kappa(t)V}_{0}\sin{\theta(t)}\Delta t+\sqrt{2D_{r}\text{Δt}}\xi_{2}\),(1b)
\(\theta\left(t+\Delta t\right)-\theta(t)=\kappa(t)\omega\Delta t+\sqrt{2D_{\theta}\text{Δt}}\xi_{3},\)(1c)
\(\kappa\left(t+\Delta t\right)-\kappa\left(t\right)=-2\kappa\left(t\right)B\left(\nu\right)\),(1d)
where θ is the direction of the movement, x and y indicate the spatial coordinate and \(t\) is the time. \(\kappa\) is the rotational direction (\(\kappa=1\ \)for CCW and -1 for CW), and \(\omega\) is the angular speed. Reversal events are represented by the telegraph process, characterized by \(B\left(\nu\right)\) which is a Bernoulli random variable with success probability νΔt (\(\nu\) is reversal rate). The noise terms \(\xi_{1}\), \(\xi_{2}\) and \(\xi_{3}\) follow a standard normal distribution. The default values of the parameters were derived from the experimental data (Table 1).
The mathematical derivation of the effective diffusion coefficientD (as a measure of foraging efficiency) can be obtained by calculating the probability distribution functions\(\Psi_{\pm}(\mathbf{r},\theta,t)\) (34,35),
\(\frac{\partial\Psi_{\pm}}{\partial t}+\nabla\bullet(\dot{\mathbf{r}}\Psi_{\pm})+\frac{\partial}{\partial\theta}\left(\dot{\theta}\Psi_{\pm}\right)=\nu(\Psi_{\mp}-\Psi_{\pm})\), (2)
with\(\dot{\mathbf{r}}={\pm V}_{0}\mathbf{n(}\theta\mathbf{)}-\ D_{r}\nabla\log\Psi_{\pm}\),\(\dot{\theta}=\pm\omega-D_{\theta}\frac{\partial}{\partial\theta}\log\Psi_{\pm}\). Here, the probability density function of cells’ spatial position\(\mathbf{r}=(x,y)\) changes following the Fokker-Planck equation associated with the Langevin equations (Eqs. 1).
We can obtain the time-dependent expected change in orientation of the diatom cells by multiplying Eq. (2) by \(\cos{\theta}\) and separately by \(\sin{\theta}\) respectively, and then integrating both equations over \(\theta\) and \(\mathbf{r}\). By solving a linear system of ordinary differential equations for\(\left\langle\operatorname{\kappa\ cos}{\theta}\right\rangle(t)\) and\(\left\langle\kappa\sin{\theta}\right\rangle(t)\) (see Text 1 for details), the analytical prediction of temporal correlation of orientation \(\left\langle\cos{\theta}\right\rangle\) is given by
\(\left\langle\cos{\theta}\right\rangle\left(t\right)=e^{-\left(D_{\theta}+\nu\right)t}(\cos{\sqrt{\lambda}\text{νt}}-\frac{1}{\sqrt{\lambda}}\sin{\sqrt{\lambda}\text{νt}})\), (3)
where\(\left\langle\bullet\right\rangle=\int{d\mathbf{r}}\int_{0}^{2\pi}{\text{dθ}\left(\Psi_{+}+\Psi_{-}\right)}\)and \(\lambda=\left(\frac{\omega}{\nu}\right)^{2}-1\).
Theoretically, we can further obtain the analytical predictions of the time-dependent mean-squared displacements (MSD) and effective diffusion coefficient from the Fokker-Planck equation (Eq. 2). Using mathematical derivation, we can obtain the analytical expression of MSD\(\left\langle\mathbf{r}^{2}\right\rangle\left(t\right)\ \)as
\(\left\langle\mathbf{r}^{2}\right\rangle\left(t\right)=4Dt+\ \frac{2V_{0}^{2}}{\nu^{2}{(\lambda+\alpha^{2})}^{2}}\left[\left(\lambda+2\alpha-\alpha^{2}\right)\left(1-e^{-\alpha\nu t}\cos{\sqrt{\lambda}\text{νt}}\right)+\left(\lambda-2\alpha\lambda-\alpha^{2}\right)\frac{1}{\sqrt{\lambda}}e^{-\alpha\nu t}\sin{\sqrt{\lambda}\text{νt}}\right],\)(4)
where\({D=D}_{r}+\frac{V_{0}^{2}\left(\alpha-1\right)}{2\nu\left(\lambda+\alpha^{2}\right)}\), (5)
is the effective diffusion coefficient (or diffusivity for simplicity hereafter) with \(\lambda=\frac{\omega^{2}}{\nu^{2}}-1\) and\(\ \alpha=\frac{D_{\theta}}{\nu}+1\). Note that for noncircular motion, i.e. \(\omega=0,\nu=0\), our model (Eq. 5) defines a system of persistent random walks characterized by a diffusivity\(D=D_{r}+\frac{V_{0}^{2}}{(2D_{\theta})}\) (28, 36).