Hyper-Ballistic Superdiffusion of Competing Microswimmers

Abstract

Hyper-ballistic diffusion is shown to arise from a simple model of microswimmers moving through a porous media while competing for resources. By using a mean-field model where swimmers interact through the local concentration, we show that a non-linear Fokker–Planck equation arises. The solution exhibits hyper-ballistic superdiffusive motion, with a diffusion exponent of four. A microscopic simulation strategy is proposed, which shows excellent agreement with theoretical analysis.

1. Introduction

The transport of self-propelled particles through an obstacle-laden or porous environments is of importance to a broad range of scientific disciplines [1]. With applications ranging from the dispersion of contaminants in soils to transport of cells inside the body, and even medical applications in the context of micro- and nano-robotics, such processes provide a fruitful and important avenue of research from theoretical, experimental and industrial perspectives.

Transport of passive particles, such as Brownian colloids or tracers, through disordered porous media have been studied intensively in recent decades, and is, to a large extent, understood through theories such as the framework presented by Brenner [2]. Biological swimmers, on the other hand, display a plethora of phenomena not seen in the passive counterparts, most of which arise from their active nature, i.e., their ability to absorb energy from the environment [3]. In the majority of cases, it is assumed that this energy is used to produce directed persistent motion, after which it is dissipated back into the environment. This break both the classical fluctuation-dissipation theorem and detailed balance, making such systems non-equilibrium and fundamentally different from passive systems [4]. Self-propelled motion in porous media is an active area of research, displaying many intriguing phenomena such as enhanced motion and optimal swimming strategies [5,6,7,8,9,10,11], directional locking [12,13], and hydrodynamic trapping at obstacles [14,15].

Persistent motion of self-propelled particles is known to be the origin of many interesting phenomena in active systems. One example is the motility-induced accumulation of particles near solid obstacles or surfaces. Such accumulation can be purely dynamical in origin, where the finite-time correlations in the particle direction of motion give rise to extended durations, whereby a particle collides with the solid [1,16,17,18]. In porous media, this may give rise to long durations of trapping in dead-ends of the porous matrix. Combined with other effects of geometric confinement present in media with a high filling fraction, this can give rise to anomalous diffusion of the sub-diffusive type, whereby the particles mean square displacement scales in time as $\langle r^2\left(t\right)\rangle \sim t^{2\tau }$, with $0\le \tau \le 1/2$ [19]. This is typically attributed to power-law distributed trapping times, and is common for both passive and active systems in strongly heterogeneous environments [20,21,22,23,24,25,26].

Superdiffusion with $1/2\le \tau \le 1$ can also be observed, with typical examples including Lévy flights, animal migration patterns and tracers in turbulent flows [27,28,29,30,31]. Hyper-ballistic diffusion, on the other hand, where $\tau \ge 1$, is much more rare. Hyper-ballistic diffusion seems almost like a contradiction, since it suggests transport that is faster than motion without any change in the direction of motion. There are still a few systems that exhibit such behavior. The random acceleration process is a model where $\langle r^2\left(t\right)\rangle \sim t^3$, which has seen many applications [32]. For example, it can be mapped onto the motion of a semi-flexible polymer inside a tube [33]. Quantum interference effects give Gaussian distribution with $\langle r^2\left(t\right)\rangle \sim t^3$ [34,35]. Such effects have also been observed in optical experiments studying wave packets moving through random potentials [36], an effect which is related to Anderson localization [37]. Both in classical and quantum systems, hyper-ballistic transport may occur for particles that receive unbounded amounts of energy from some random potential. For example, Golubovic et al. finds superdiffusion with $\langle r^2\left(t\right)\rangle \sim t^{18/8}$ in a time-dependent random potential [38]. A simple random walk model with memory effects, called the “elephant random walk” because elephants have long memories, also gives rise to superdiffusion, but with sub-ballistic behavior [39].

One way in which hyper-ballistic diffusion can appear is when the particle velocity is allowed to grow without bounds [40]. For normal Langevin systems, such behavior is not possible, since the fluctuation-dissipation theorem ensures a balance between the fluctuations driving the system and the dissipation. For active systems, such as biological microswimmers, the fluctuation-dissipation theorem is famously broken, making these systems potential candidates for hyper-ballistic diffusion. In this paper, we propose a simple model where self-propelled swimmers move in a porous media while competing for resources. Hyper-ballistic motion is shown to arise as a direct consequence of a nutrient landscape that dynamically changes with the swimmer density.

2. Model of Competing Swimmers

The N swimmers move in an isotropic porous medium that defines a mean free path length $\lambda$ (see Figure 1). The swimmers are assumed to move in straight lines at a speed $\mathit{v}(\mathit{x},t)$ until they hit the walls of the porous media. As mentioned above, some cells generate more persistent straight-lined motion when starved or in the absence of external signals [41,42]. When the swimmers collide with the walls of the medium, after a characteristic distance $\lambda$, they randomly change direction. The porous medium is only represented in this minimal way, so that the effect of potential trapping in dead ends or near obstacle surfaces is ignored [19].

The swimmers also move through a nutrient concentration $C_N(\mathit{x},t)$ which can change dynamically in time due to the swimmers presence. The study of swimmers motion in various chemical gradients, chemotaxis, has been studied for decades in the biological and biophysical communities [43,44]. The motion of cells may also change drastically under starved conditions, when the nutrient concentration is low. While some cells respond to such conditions by deactivating flagellar motors [45], other observations include more persistent straight-lined motion [41,42]. More recent studies have modeled nutrient or activity landscapes implicitly as time- or space-dependent self-propulsion speeds in models of active matter [46,47,48,49]. Recently, a dynamical resource landscape was considered in a lattice model [50].

In the present context, we consider in the simplest case an inverse relationship $C_N(\mathit{x},t)$$\sim 1/C(\mathit{x},t)$ between the nutrient concentration $C_N(\mathit{x},t)$ and swimmer concentration $C(\mathit{x},t)$, implying that the presence of many swimmers deplete the nutrients in that region. For example, the nutrients can be supplied at a given rate from the walls in the medium, for example due to prior accumulation. We shall take the pore size to be much smaller than the nutrient diffusion length, so that the nutrient supply is limited by the release rate from the walls. Each swimmer eats at a rate $\sim C_N(\mathit{x},t)$ and we shall assume that the nutrient supply is almost depleted, so that $C_N(\mathit{x},t)$ is much smaller than the saturation limit, hence the competition. In a given pore volume, the food consumption will thus be $\sim C_N(\mathit{x},t)C(\mathit{x},t)$. This consumption is balanced by the constant supply from the walls so that the nutrient concentration scales as $C_N(\mathit{x},t)\sim 1/C(\mathit{x},t)$, as above.

Since the swimmers move against a Stokes drag force $\propto v$, the power they dissipate by swimming is proportional to $v^2$. Assuming that they convert a constant fraction of the energy they consume into this power, we arrive at a swimming velocity scaling as $v^2(\mathit{x},t)\sim C_N(\mathit{x},t)$. A similar scaling is also found in microscopic models for active particles where energy consumption is explicitly modeled [51,52,53]. Combining the above, we arrive at

$$\left|\mathbf{v}\left(C\right)\right|\equiv v\left(C\right)=v_0{\left(\frac{C}{C_0}\right)}^{-1/2}$$

(1)

where $v_0$ is the velocity at some reference concentration $C_0$. In other words, the swimmers will accelerate once their concentration is down and the food competition is reduced. Their motion is diffusive in the sense that swimmers perform a random walk with a concentration dependent step length $\lambda$. This concentration dependence represents an interaction with the neighboring swimmers that define the local concentration.

3. Simulation Model

In order to simulate the collective motion of the swimmers, we represent them as point particles with positions ${\mathbf{x}}_i$, $i=1,\dots N$ that are updated according to the algorithm

$$\begin{array}{cc}\hfill {\mathbf{x}}_i(t+dt)& ={\mathbf{x}}_i\left(t\right)+\mathbf{v}\left(C\left({\mathbf{x}}_i\right)\right)dt\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill {\mathbf{v}}_i(t+dt)& =\left\{\begin{array}{cc}\mathcal{R}\mathbf{v}\left(C\left({\mathbf{x}}_i\right)\right)\hfill & if\Delta x_i>\lambda ,\hfill \\ \mathbf{v}\left(C\left({\mathbf{x}}_i\right)\right)\hfill & otherwise,\hfill \end{array}\right.\hfill \end{array}$$

(3)

where $dt$ is fixed small timestep, $\mathcal{R}$ is a rotation operator representing a random collision that uniformly re-orients the swimmers direction of motion, and $\Delta x_i$ is the displacement, since the last application of the random rotation operator $\mathcal{R}$. We emphasize that, in the above algorithm, the velocity vector $\mathbf{v}\left(C\left({\mathit{x}}_i\right)\right)$ only depends on the local concentration through its magnitude through Equation (1). The random re-orientations due to collisions only affect the swimmers direction of motion. To sample the local concentration around particle i, we use a scheme similar to that in Ref. [40]. Here, a number of interaction partners $N_r$ is introduced, along an associated volume $V_r\left(\mathbf{x}\right)$ of a sphere that contains the $N_r$ nearest neighbors. This is illustrated in Figure 2. We emphasize that this type of interaction is non-local, in the sense that the range of interaction has no bounds. The local concentration used in Equations (2) and (3) is then sampled as

$$C\left({\mathbf{x}}_i\left(t\right)\right)=\frac{N_r}{V_r\left({\mathbf{x}}_i\left(t\right)\right)}$$

(4)

The operator $\mathcal{R}$ in Equation (3) rotates ${\mathbf{v}}_i$ into a new arbitrary direction without changing the speed, and the distance $\Delta x_i$ is calculated from this point, so that the next application of $\mathcal{R}$ may be determined. So, the swimmers are random walkers with a step length $\lambda$, and potentially, they spend several time steps to move $\lambda$ along straight lines.

3.1. The Fokker–Planck Equation

Following the classical methods [54,55], we derive the Fokker–Planck equation that governs the time evolution of the particle concentration. The general form of the master equation takes the standard form

$$\frac{\partial C(\mathbf{x},t)}{\partial t}=\int d^{3}r\left[C(\mathbf{x}-\mathbf{r},t)W(\mathbf{x}-\mathbf{r},\mathbf{r})-C(\mathbf{x},t)W(\mathbf{x},-\mathbf{r})\right],$$

(5)

where $W(\mathbf{x},r)$ is the jump rate associated with a step $\mathbf{x}\to \mathbf{x}+\mathbf{r}$. In the present case, the size of a step a swimmer takes depends on the nutrient concentration and, hence, also on the particle concentration. Assuming a slowly varying spatial dependence in the jump rates $W(\mathbf{x}-\mathbf{r},\mathbf{r})$, we may Taylor expand around $\mathbf{x}$ in its first argument. This results in the Fokker–Planck equation

$$\frac{\partial C(\mathbf{x},t)}{\partial t}=\frac{1}{2}\int d^{3}r{r}_i{r}_j\frac{{\partial }^2}{\partial x_i\partial x_j}\left[C(\mathbf{x},t)W(\mathbf{x},\mathbf{r})\right],$$

(6)

Rearranging the derivatives, we have

$$\frac{\partial C(\mathbf{x},t)}{\partial t}=\frac{1}{2}{\nabla }^2\left[a_2\left(\mathbf{x}\right)C(\mathbf{x},t)\right],$$

(7)

where $a_2\left(\mathbf{x}\right)$ is the second moment of the jump length per time,

$$a_2\left(\mathbf{x}\right)=\int d^{3}r\frac{{\mathbf{r}}^2}{3}W(\mathbf{x},\mathbf{r})=\frac{1}{3}\frac{{\lambda }^2}{\tau }=\frac{1}{3}\lambda v,$$

(8)

where $\tau$ is the local mean free time, and we have used that the local swimmer velocity is $v=\lambda /\tau$. The factor of $1/3$ comes from the identity ${\mathbf{r}}^2=3r_i^2$. Using $v=v_0{(C/C_0)}^{-1/2}$, we obtain

$$\frac{\partial C}{\partial t}=\frac{v_0\lambda C_0^{1/2}}{2}{\nabla }^2{C}^{1/2},$$

(9)

or

$$\frac{\partial C}{\partial t}=D_0\nabla ·\left({\left(\frac{C}{C_0}\right)}^{-1/2}\nabla C\right),$$

(10)

where $D_0=v_0\lambda /4$. It may be shown [56] that this equation has the normalizable solution

$$C(r,t)\propto {\left[y{(r,t)}^2+\frac{8{\pi }^4}{N^2}\right]}^{-2}{t}^{-2}$$

(11)

where $y(r,t)=r/\left(C_0{\left(D_{0}t\right)}^2\right)$. At large values of r, the solution has a tail $C(r,t)\sim r^{-4}$. For other properties of the solution to concentration dependent diffusivity, see Ref. [56].

3.2. Mean Square Displacement and Simulation Results

From the above solution, we can calculate the mean squared displacement $\langle r^2\left(t\right)\rangle$ or, equivalently, the root-mean-squared displacement (rms) $r_{rms}\equiv \sqrt{\langle r^2\left(t\right)\rangle }$, by the integrals

$$r_{rms}^2=\langle r^2\left(t\right)\rangle =\frac{\int dV{r}^{2}C(\mathbf{r},t)}{\int dVC(\mathbf{r},t)}=\frac{{\int }_0^{\infty }dr{r}^{4}C(r,t)}{{\int }_0^{\infty }dr{r}^{2}C(r,t)},$$

(12)

where $dV=dxdydz$ is the volume element in three dimensions, and where we changes to radial coordinates in the second equality, with r the radius from where the swimmers are released. To extract the leading temporal scaling at late times, we observe that from Equation (11) we have for any exponent $\alpha$ the scaling

$${\int }_0^{\infty }dr{r}^{\alpha }C(r,t)\sim t^{2+2\alpha }{\int }_0^{\infty }dy{\left[y^2+\frac{8{\pi }^4}{N^2}\right]}^{-2},$$

(13)

where we changed integration variable from r to $y(r,t)$, as given above. With this change in the variable, the time-dependence has been made explicit. Using this result together with Equation (12) immediately gives

$$r_{rms}^2\sim t^4.$$

(14)

This is a remarkable result, as it predicts a spreading rate which is beyond that of ballistic growth. Ballistic growth, which would result if the swimmers would move at constant speed and never change direction of motion, would give $\langle r^2\left(t\right)\rangle \sim t^2$.

Note that, since $r^{4}C(r,t)\sim r^0$ for large r, the integral in the numerator diverges. However, as was pointed out in Ref. [40], there will always be a maximal swimmer position $r_{max}$ in any numerical realization, which effectively acts as a cut-off. This means that the proportionality with $t^4$ survives and that the prefactor varies linearly with $r_{max}$. Since $r_{max}$ fluctuates in each simulation, the prefactor in Equation (14) is noisy, even with large ensemble numbers. Figure 3a. shows results from numerical simulations where particle concentration was initialized as a delta-peak. The results show that the swimmers indeed conform to the prediction of Equation (14). They are indeed noisy, as expected, since the largest swimmer positions contribute significantly to the $r_{rms}\left(t\right)$ evaluation.

Generalizing the model to $v\sim \lambda C^{-\gamma }$ where $\gamma \ne 1/2$, the solution becomes $C(r,t)\sim {(y^2+k)}^{-1/\gamma }$. Now, the integrals in Equation (12) converges when $\gamma <0.4$ and, in this case, the $\langle r^2\rangle \left(t\right)$ becomes significantly less noisy. The generalized solution leads to $r_{rms}\left(t\right)\sim t^{\tau }$, where [40]

$$\tau =\frac{1}{2-3\gamma }.$$

(15)

In order to validate both the analytical and computational model, we have carried out simulations using this generalized concentration dependence in the velocity. The results, which are shown in Figure 3b, demonstrate close agreement between simulations and theory.

4. Conclusions

We have studied the diffusive behavior of microswimmers competing for resources. Through simple arguments, we have derived a density-dependent Fokker–Planck equation that effectively models this scenario, and have proposed a microscopic numerical model that agrees with the theory. We find that the swimmers exhibit a hyper-ballistic superdiffusive behavior, where the growth of the mean squared displacement in time is quartic.