Nutrient exposure of chemotactic organisms in small-scale turbulent flows

Micro-organisms living in a turbulent fluid environment often use directed motility to locate regions of higher than average nutrient concentrations. Here, we consider a simple continuum model for the distribution of such chemotactic particles when the particles and the chemoattractant are both advected by a turbulent flow. The influence of chemotactic sensitivity on the spatial distribution of the particles is characterized for different types of advected chemical fields. Using an effective diffusion approximation, we obtain an analytical expression for the nutrient exposure resulting from the chemotactic activity of the particles, generalizing previous results obtained for the case of phototaxis in flows. We show that the biological advantage of chemotaxis in such systems is determined by the spatial variability of the averaged chemoattractant field and the effective diffusivity of the turbulent flow.


Introduction
The motion of micro-organisms and cells in response to chemical signals plays an important role in a wide range of biological processes. For example, marine bacteria are capable of locating nutrients by detecting and moving towards higher concentrations following the nutrient gradients produced by phytoplankton photosynthetic products, decaying algal cells, etc [1,2]. However, complex oceanic environments, where the wind, convection and gravitational forces create turbulent flows, affect the bacterial search by altering both the concentration distribution surrounding the microscale nutrient sources and the motion of chemotactic bacteria. In recent years, some progress has been made in understanding the bacterial chemotactic behavior in this marine context. For example, Bowen et al [3] and Luchsinger et al [4] proposed theoretical models of bacterial clustering in a turbulent environment. In a recent work, Stocker et al [5] used microfluidic experiments to create a complex environment consistent with that expected in the ocean and found that bacterial chemotaxis can exploit patchily distributed nutrients and nutrient plumes. This confirms that microscale motility can have significant effects on net large-scale biogeochemical fluxes in the ocean.
In order to understand how the organisms effectively respond in different small-scale turbulent flows and the role of the spatial distribution of the nutrients in the organisms' exploitation of them, in this paper we consider a mathematical model that describes the distribution of a population of chemotactic particles that swim with a velocity proportional to the local gradient of a chemoattractant in a time-dependent chaotic flow. Here, we study the general case where the chemoattractant distribution is also transported by the moving fluid medium, extending previous results on phototactic swimming micro-organisms advected by a turbulent velocity field [6] where a static external field independent of the flow (i.e. an illumination field) was considered. In order to describe the temporal evolution of the nutrient concentrations and the distribution of micro-organisms at the population level, we consider a continuum model that is similar to those used in the study of chemotactic aggregation [7,8] (for a recent review, see [9]). In contrast to these works where the chemoattractant is assumed to be produced by the cells, creating an instability that induces the aggregation, here the evolution of the nutrients has external sources and is independent of the micro-organisms. A further significant difference with most continuum models of chemotaxis is that in this work we consider the effect of advection, which can dramatically alter the features of particle clustering [10].
We characterize the main statistical properties of the spatiotemporal distributions of the chemotactic particles advected by a velocity field generated using a stochastic model of a spatially smooth and temporally chaotic unsteady flow, representing microscale flows 3 corresponding to length scales below the Kolmogorov scale. Our analytical description, based on the effective turbulent diffusion approximation, allows us to reproduce the simulations and quantify the nutrient exposure experienced by the particles, describing how this depends on the particles' chemotactic sensitivity and the features of the chemoattractant source.

The model
We consider swimming particles that detect local instantaneous gradients of a chemoattractant (e.g. some nutrients) and move in the direction of increasing concentration with a swimming velocity proportional to the detected gradient. Although there can be a great variety of detailed microscopic mechanisms for the chemotactic response, this description can be considered as the simplest representation consistent with an effective chemotactic flux generated at the population level [9,11,12]. At the same time, both the particles and the chemoattractant are advected by a time-dependent incompressible flow and may diffuse. We describe the particle density and chemical concentration at any point r and time t by the continuous scalar fields ρ(r, t) and c(r, t), respectively, and assume that the chemoattractant is produced by a fixed localized source and degraded with a characteristic rate b. The simplest continuum model describing the temporal evolution of the system that includes these mechanisms is given by where V(r, t) is the carrier flow and D ρ and D c are constant diffusion coefficients. Typically the particles' diffusivity is smaller than the chemical diffusivity in real systems; however, as we will see below, the value of these coefficients will not significantly affect our principal results because the main effective mechanism for diffusion is the turbulent flow. The parameter χ is the chemotactic sensitivity and defines the strength of the particles' response to the chemical substance that is produced in the medium, at a rate following a stationary source distribution S(r). For our simulations we consider a Gaussian distribution, although these results do not depend on the exact functional form of S(r). Thus, the concentration field c(r, t) depends only on the properties of the carrier flow and the characteristics of the production and degradation. When the temporal scale of the degradation and production (b −1 ) is much faster than the characteristic timescale associated with the flow, the concentration field is not affected significantly by advection and can be approximated by the stationary solution of equation (2), that is, c(r) ≈ S(r). In this case, the system is reduced to the phototactic problem studied previously in [6]. However, when the production-degradation rate b is comparable or smaller than the inverse flow timescale, advection has a strong influence on the concentration field, creating a more irregular spatial structure and temporal variability. Thus, by varying the parameter b we consider different regimes where the chemical field is determined to differing degrees, by the source distribution or transport processes. The properties of the particle distribution depend on the interplay between the particle motility induced by the chemotactic term, and the advection. Assuming no-flux or periodic boundary conditions the total mass of the biological component is conserved and can be characterized by the integral of the spatial density, ρ 0 . After relaxation of transients, the chemical field also reaches a statistically stationary state with the total concentration, c 0 , equal to the spatial integral of S(r).

4
For the velocity field, we use a two-dimensional isotropic and homogeneous synthetic turbulence model introduced in [13]. The flow V(r, t) is defined through the stream function, ψ, as V = ∇ × (0, 0, ψ), where ψ is the solution of the Ornstein-Uhlenbeck process produced by the stochastic partial differential equation [6,14] ∂ψ Here ν is a damping parameter that determines the temporal correlations of the spatial structures in the flow, W represents a Wiener process and ξ is the amplitude of the noise that is chosen for each Fourier mode with wavenumber k according to a prescribed kinetic energy spectrum for the flow that verifies E(k) = |v k | 2 = ξ/2ν, wherev k is the velocity field amplitude for the corresponding mode of the streamfunction (for more details, see [14]). To represent small-scale turbulence below the Kolmogorov scale, we use the Kraichnan spectrum [15] which has a single peak at k 0 and an exponential tail, thus producing a spatially smooth time-dependent velocity field. We set the lengthscale of the highest energy mode to twothirds of the total system size defined on the unit square (k 0 = 3π/2). The temporal scale in equations (1) and (2) was also normalized to have a root-mean-square velocity of unity. We used periodic boundary conditions and a lattice grid (256 × 256 nodes) with spatial and temporal steps x = y = 4 × 10 −3 and t = 1 × 10 −4 , respectively, checking that the results do not differ significantly for smaller space and time steps. The normalized spatial integral of the particle density was set to ρ 0 = 1, which remains constant during the simulation. For the source distribution of the chemoattractant we used a Gaussian distribution with unit amplitude and width σ = 1/6 following S(r) = exp[−r 2 /(2σ 2 )], and normalized diffusion coefficients D ρ = D c = 0.002. Thus, after proper rescaling, the non-dimensional parameters b and χ can be both varied to study the different regimes for the solutions of equations (1) and (2). The numerical simulations were performed using a fourth-order Runge-Kutta method combined with a cubic interpolation semi-Lagrangian scheme [16] to describe the advection. For the spatial integration of the particle density, we used a novel discretization method proposed by Grima and Newman [17] for advection-diffusion equations that allows an accurate analysis of the evolution of the system without the dissipative effects of other schemes.

Simulations
In general, in the absence of a source and chemotaxis, the variance of a passive concentration field dispersed by an ergodic incompressible flow decays and asymptotically approaches a spatially uniform distribution; however, for equations (1) and (2), the spatial concentrations of c and ρ may remain highly non-uniform and their characteristic structure depends on the chemical degradation-production rate and the chemotactic sensitivity. An example of the temporal evolution of c and ρ for b = 5 and χ = 0.1 is shown in figure 1 (see also the supplementary movies available from stacks.iop.org/NJP/12/103043/mmedia). As we can see here the chemical field is dispersed by the turbulent flow in the region around the source. The instantaneous chemical field changes in an irregular fashion, but on average the concentration of c is larger at the center of the domain. The distribution of the chemotactic particles also has a higher density  (on average) near the center, but typically it is more concentrated into thin filamental structures that also extend into regions with lower concentration of the chemoattractant. The temporal evolution and spatial structure of the chemical field are independent of ρ and depend on the degradation-production rate, b [18]. In figure 2, typical snapshots of the concentration field, c, are shown after relaxation of the transients for different values of b. In general, the turbulent flow produces a chemoattractant distribution different from the source distribution, S(r). However, when the temporal scale associated with the degradation-production rate of c (given by b −1 ) is much smaller than the characteristic timescale of the turbulent flow (that is, of the order of unity), its shape is only slightly modified by the flow and remains quite similar to S(r) (see figures 2(a) and (b)). For values of b −1 of the order of or larger than unity, the distribution of chemoattractant becomes more homogeneous and is spread over a larger area compared to the source function (see figures 2(c) and (d)).  Figure 3 shows the particle density after relaxation of the transients for a set of different values of the degradation-production rate, b, and chemotactic sensitivity, χ. The steady state reflects a statistical equilibrium resulting from the competition between the chemotaxis, which tends to accumulate all particles in regions with the largest concentration of c, and the advection by the turbulent flow, which tends to homogenize the particle density. The temporal evolution of ρ depends on the time-dependent concentration field, c, that is dispersed over a larger area as b is decreased, having a similar effect on the particle density. Furthermore, the particle distribution concentrates into thin filaments that become more pronounced as χ is increased. Thus, for large values of χ, the particles tend to accumulate in thin structures separated by sparse regions following the nutrient concentration. Although one could expect that for the cases with almost homogeneous chemoattractant distributions (small values of b) the nutrient exposure of the particles cannot increase with the chemotactic sensitivity (χ ), as we will see below, even for small values of b, the simultaneous overlap (without spatial or temporal delay) between ρ and c fields increases as the chemotactic sensitivity χ is increased.
In order to characterize the statistical properties of the ρ and c fields in the steady state, we obtain the mean distributions,ρ(r) andc(r), by averaging the values of the concentration  and density at each point over a number of statistically equivalent realizations of the stochastic flow velocity field. Since in our case the source function S(r) is radially symmetric and the flow is isotropic, the averaged system presents an approximate radial symmetry that is slightly broken by the finite rectangular domain, but in general it is a good approximation for our simulations. In figure 4(a), we show the radial profile of the averaged chemical concentration in the steady state for different values of b. As expected, the concentration is more uniform and the maximum is smaller with decreasing b. The particle density also becomes more homogeneous with decreasing b (not shown) or with decreasing chemotactic coefficient χ , as shown in figure 4(b), where the averaged radial profile of the stationary particle density is plotted for different values of χ and a fixed value of b. We also note that the instantaneous particle density is more irregular for large values of χ with strong fluctuations and large standard deviations for the statistical data. The increased benefit of chemotactic motility for the micro-organisms can be characterized by the nutrient exposure of the population by evaluating which is a measure of the overlap between the density field and the high-concentration regions. This magnitude is shown in figure 5(a) for different values of b and χ. Here we see that the nutrient exposure increases with the chemotactic sensitivity of the particles, but the rate of increase strongly depends on the properties of the concentration field, controlled by the degradation-production rate b. In the limit of passive non-motile particles (χ → 0), their density distribution is not affected by the concentration field and is spatially uniform (being ρ 0 in all of the considered unit spatial domain). In this case, the nutrient exposure is proportional to the Since in the stationary state the total chemoattractant concentration c 0 = c(r) dr = S(r) dr is independent of b, for the considered case with ρ 0 = 1 we have 0 = c 0 = 0.17. On the other hand, for large χ all the particles are clustered into a small blob that tries to follow the maximum concentration. In this case the nutrient exposure reaches a plateau value, s , that is proportional (with a constant of proportionality ρ 0 ) to the maximum value of the averaged c field located at r = 0, which reads s = ρ 0c (r = 0). Thus, since on decreasing the value of b the maximum of c is also diminished (see figure 4(a)), the saturation value of the nutrient exposure is reduced for small values of b. In other words, the maximum benefit of chemotaxis becomes smaller when the chemical field is more strongly influenced by advection.

Mean field modeling
As seen above, although the instantaneous distributions change in an irregular fashion following the turbulent flow, averaging over time or different realizations of the random flow we can characterize the stationary distributions in the physical space. To obtain a transport equation for the averaged fields ρ and c, the effect of dispersion by the flow V can be described by using a turbulent diffusion approximation [6]. From equations (1) and (2) we have ∂ρ(r, t) ∂t = D e ρ ∇ 2ρ (r, t) − χ∇[ρ(r, t)∇c(r, t)], where the diffusion and advection terms are replaced by an effective diffusion for the particle and chemical densities D e ρ and D e c , respectively, and the overlined magnitudes represent their averaged values over different realizations of the stochastic flow. Certainly, for homogeneous flows and non-reactive scalar fields, one expects that such an effective diffusion will only depend on the statistical properties of the turbulent flow. Thus, for example, in [19] the effective diffusion is derived for a passive scalar advected by a stochastic velocity field. In the case of reactive fields, the effective diffusion coefficient also depends on the reaction time as shown by Plumb [20] and used, for example, by Pasquero to study the population dynamics in a plankton ecosystem model under turbulent transport [21]. To compute the effective diffusion and some of the statistical properties of the synthetic flow, we follow [13,19] where the spatial spread of a passive scalar field is evaluated using an approximate δ-function as the initial condition. Averaging over different realizations at different times we can relate the spatial variance Var(δ) with the finite-time effective diffusivity D e . For a two-dimensional system, this relation reads D e (t) = Var(δ)/(4t), where D e is the sum of a purely molecular diffusion and the contribution due to the turbulent flow. For the turbulent flow considered here, the effective diffusion obtained as described above is shown in figure 5(b). For short times the ballistic regime dominates and D e ∝ t, however for times larger than the characteristic time of the turbulent flow, T f , defined as the time required to reach the Brownian regime, a constant effective diffusivity, D e 0 , is attained. In the case of the flow used in the numerical simulations, we obtain D e 0 ≈ 0.034 and T f ≈ 0.5. For passive scalars or reactive systems in which the reaction timescales are larger than the characteristic flow time, the effective dispersion due to turbulent mixing is well described by D e 0 . In the case when the tracer reacts before the Brownian regime is reached, the effective spreading rate is smaller than D e 0 , and using the theoretical results of Plumb the effective diffusivity can be estimated by the value of the finite time diffusion coefficient D e (b −1 ) corresponding to the reaction time, as was done in [21]. We use this approximation to find the stationary solution of equation (8) for different values of b.
Since the profiles of the source function and the averaged chemical field are close to zero in the boundaries, we can extend the domain to infinity and use the Hankel transformation to obtain the stationary state ofc from equation (8) when a radially symmetric source is considered. The Hankel transform is just the two-dimensional Fourier transform for radially symmetric problems whose kernel is the Bessel function. The zeroth-order Hankel transform and its inverse transform of a radial function f are defined, respectively, by [22] H where J 0 is the zero-order first kind Bessel function. Equating to zero and applying the Hankel transform to the rhs of equation (8), we obtain Since the Hankel transform of the Gaussian source function is H 0 {S(r )} (s) = σ 2 exp(−s 2 σ 2 /2), it is straightforward to obtainc(r ) applying the inverse Hankel transform to equation (10). This reads  (11): (i) the maximum of the integral is always located at r = 0 since the maximum of J 0 (r ) is located at r = 0; (ii) the numerator in the integral in equation (11) is independent of the reaction time, but the denominator increases with decreasing b; thusc(0) decreases as b decreases; (iii) since the total amount of c does not depend on the value of b, the shape ofc(r ) must be spread over a wider region with decreasing b. This last point also means that the approximate analytical solutions are less accurate for small values of b sincec(r ) does not tend to 0 at the boundaries (see figure 4(a) for b = 1).
The chemotactic term tends to collect particles to the regions of high chemical concentrations but it is counteracted by the turbulent advection (represented by the first term on the rhs of equation (7)), which tends to smooth out the inhomogeneities in the particle density. Thus, even for relatively large values of χ, the consequence of turbulent advection is that the fluctuations in the density and the chemical fields are not correlated. This has been confirmed numerically by showing that the radial profiles of the magnitudes ρc(r ) andρ(r )c(r ) averaged over statistically independent realizations match ( figure 6(a)). Since the chemotaxis is driven by gradients of chemoattractant concentration, it is also interesting to compare the radial component of ρ∇ r c(r ) and its mean field approximationρ(r )∇ rc (r ). As presented in figure 6(b), the two profiles are very similar, showing that their simultaneous correlations are also weak and decrease with decreasing χ . This property can be used in equation (7) to obtain the stationary average density of the chemotactic particles from D e ρ ∇ 2ρ (r) − χ∇ [ρ(r)∇c(r)] = 0. (12) The solution of equation (12) is the Boltzmann distribution where Z is a normalization factor that is analogous to a partition function. Similarly to the equation forc(r), the effective diffusion in equation (7) may depend on the characteristic timescale corresponding to the chemotactic term. In this case, since we are using an extended nutrient source (with small gradients of the order of unity) and relatively small values of χ, this temporal scale is larger than the characteristic timescale of the flow and we have D e ρ = D e 0 . The analytical solution for the averaged particle distribution, from equation (13), is shown in figure 4(b) for different values of the chemotactic coefficient χ. This approximation works better for small values of χ when the particle and chemical fields are less correlated and the ρ field is more uniform.
As shown in figure 6(a), we can approximate ρc(r) byρ(r)c(r) and use equations (11) and (13) to obtain the following expression for the nutrient exposure: where we used the 'partition function' Z defined in (13). Using only the analytical expression for the averaged chemoattractant field given by equation (11) and the value of the effective diffusivities, we can obtain the nutrient exposure as plotted in figure 5(a) by dashed lines for different values of b.
For small values of χ, there is a linear relationship between and the chemotactic sensitivity. This is easy to see if we consider Z as the moment-generating function of the random variablec(r) with a uniform probability density function. In this case, its logarithm is the cumulant-generating function and satisfies log Z = ∞ n=1 χ/D e ρ n κ n /n!, where κ n are the cumulants ofc(r). Using this relation we can write equation (14) as Since the average chemical concentration is constant (independently of b), we have that κ 1 = c 0 = S(r) dr = 0.17. The second cumulant is the spatial variance of the averaged concentration field, κ 2 = Var(c). Thus, in this linear regime (i.e. weak chemotaxis), the biological advantage of directed motility (given by the increase in nutrient exposure) is proportional to the spatial variability of the averaged chemical field and inversely proportional to the turbulent diffusivity following the equation In order to estimate analytically the spatial variance of the averaged concentration field, we can use Parseval's theorem for the Hankel transform and equation (10) to obtain  (19).
where Ei(x) is the exponential-integral function for the real argument x. For large arguments (large values of σ 2 b/D e c ), we can expand the exponential-integral function as Ei(x) = e −x /x N −1 n=0 n!/(−x) n + O N !x N and the spatial variance is reduced to With this, we can approximate the nutrient exposure by the following linear function, which is strictly only valid for small values of χ in the limit of large b and σ . In figure 7, we represent this approximation for small values of χ and three different values of b. As we can see here, equation (19) represents accurately the exact solution, equation (14), for large values of b and shows explicitly that the nutrient exposure decreases with decreasing b in agreement with the numerical simulations. For small values of b, this approximation is less accurate and eventually, as occurs for b = 1, may become incorrect, resulting in negative slopes for equation (19).

Summary and discussion
We have studied the spatial distribution of a population of chemotactic micro-organisms when both the chemoattractant field and the micro-organisms are embedded in a moving microscale turbulent flow. The temporal evolution of the density and concentration fields has been described using a continuum model of a pair of coupled partial differential equations based on a variant of the standard Keller-Segel model. Numerical simulations and analytical results have allowed us to describe some important features of the considered system. Specifically, using a mean-field approach based on the properties of the averaged concentration field and an effective diffusion approximation for the dispersion by the flow we obtained an accurate description for the average distribution of the chemoattractant and the particles, and the nutrient exposure of the population. The obtained solution generalizes previous results for the case of a static external field and shows that the particle density distribution is crucially determined by the averaged chemoattractant field. As a consequence of the irregular character of the fluctuations due to the turbulent flow that inhibits the buildup of strong simultaneous correlations between both fields, on average the particles effectively follow the mean chemical field and cannot exploit the instantaneous concentration fluctuations. Thus, even for relatively large values of the chemotactic sensitivity, the turbulent flow decorrelates the particle and the concentration fields and the filamental structure of the underlying chemical field may be abstracted into a coarse-grained effective gradient. The obtained results show that the nutrient exposure of the population increases linearly with the chemotactic sensitivity for weakly chemotactic organisms and saturates for large chemotactic responses. The rate of increase and the saturation value both depend on the properties of the advected chemical field that is determined by the degradation-production rate, i.e. to what extent the chemical field is influenced by the carrier flow. Thus, for example, the benefit of chemotaxis becomes smaller when the chemical field is more strongly influenced by advection. Importantly, the increase in nutrient exposure also depends on the spatial variability of the averaged chemical field, i.e. a more spatially heterogeneous source field induces a larger nutrient exposure for organisms that effectively detect and respond to the resource.
Employing this continuous description, further work might take into account quorum sensing as a possible mechanism for a more effective search strategy [23,24] and possible three-dimensional effects, for example in the case of sinking marine snow particles. In addition, the presented model may also be applied to describe phoretic particles, which move due to the gradient of a solute. This phenomenon called diffusiophoresis is obtained by slaving the dynamics of large particles to a dilute salt [25]. It might be interesting to apply some of the obtained results to the study of the effect of small-scale turbulent flows on the migration of these particles towards regions of higher salt concentration.