Impact of bacteria motility in the encounter rates with bacteriophage in mucus

Bacteriophages—or phages—are viruses that infect bacteria and are present in large concentrations in the mucosa that cover the internal organs of animals. Immunoglobulin (Ig) domains on the phage surface interact with mucin molecules, and this has been attributed to an increase in the encounter rates of phage with bacteria in mucus. However, the physical mechanism behind this phenomenon remains unclear. A continuous time random walk (CTRW) model simulating the diffusion due to mucin-T4 phage interactions was developed and calibrated to empirical data. A Langevin stochastic method for Escherichia coli (E. coli) run-and-tumble motility was combined with the phage CTRW model to describe phage-bacteria encounter rates in mucus for different mucus concentrations. Contrary to previous theoretical analyses, the emergent subdiffusion of T4 in mucus did not enhance the encounter rate of T4 against bacteria. Instead, for static E. coli, the diffusive T4 mutant lacking Ig domains outperformed the subdiffusive T4 wild type. E. coli’s motility dominated the encounter rates with both phage types in mucus. It is proposed, that the local fluid-flow generated by E. coli’s motility combined with T4 interacting with mucins may be the mechanism for increasing the encounter rates between the T4 phage and E. coli bacteria.

by the Langevin equation: where ξ ξ ξ ( ( (t t t) ) ) = 0 and ξ i (t)ξ j (t ) = δ i j δ (t − t ). This first order stochastic equation is linear on v v v B , and τ relax = m B /γ B represents the relaxation time in the exponential decay of the velocity when perturbed out of equilibrium (homogeneous equation). To estimate this relaxation time, one can use Stokes relation for the coefficient of friction γ B = 6πR B ν. Taking the viscosity of water ν ≈ 10 −3 Pa sec, the mass of E. coli m B ≈ 1 pg, and approximating E. coli as a sphere of radius R B ≈ 1 µm, this leads to a relaxation time of τ relax ≈ 0.05 µsec. This is four orders of magnitude smaller than the typical run time for E.
coli, τ R ∼ 1 sec. This suggests that the inertial term is not relevant, as far as the perturbation in the velocity is not drastic.
To estimate perturbation in the velocity the typical velocity change, ∆v B , induced by the kinetic energy of the fluid on the bacteria, one can assume a fluid similar to water, which is common in biology. The density of water is ρ H 2 O ≈ 1 kg dm −3 . The mass of a water molecule is m H 2 O = 18 g mol −1 . Thus, the typical distance between fluid molecules is l ≈ 3 · 10 −10 m. In Then, the final position of a bacterium after each run is where the mean run rate is The term x B (0) is the position of the bacteria at the beginning of the run.
In the tumble regime, E. coli stops swimming and rotates the flagella in the opposite direction, disassembling the bundle and moving erratically. The mean tumble time is about τ T ∼ 0.1 sec, which is an order of magnitude shorter than the running time.
The longitudinal axis of E. coli reorients on average θ T ≈ 70 • . While tumbling, the bacteria can diffuse due to the thermal fluctuations of the fluid. Using Einstein's relation from the dissipation-fluctuation theorem one can estimate the diffusion of the non-swimming bacteria, D B = kT /γ ≈ 0.22 µm 2 m sec −1 . The mean-squared displacement of the bacteria during the tumbling period is proportional to the bacterial diffusion and the tumbling time, MSD = 6D B τ T ≈ 0.13 µm 2 sec −1 . The typical S.2/S.16 displacement of the bacteria during the tumbling period is ∆r ∼ √ MSD ≈ 0.36 µm. This is an order of magnitude smaller than the actual size of E. coli. Thus, the displacement during the tumbling regime is negligible, and it is not included in the model. For the sake of computation efficiently, instead of using a Langevin equation of motion for the rotation, the tumbling process is generated probabilistically. The mean tumbling time is an order of magnitude shorter than the run time. In a first approximation, the tumbling regime is evaluated instantaneously. The change of direction is approximated from a uniform distribution between 0 and 2π, Thus, in the computational model, the updated directional change of the bacterium after a tumble event is where δt T → 0 indicates that the step is virtually instantaneous.

Walks
Based on experimental data, it is plausible to assume that the mechanism underlying phage subdiffusion in mucus emerges from the interaction of phage with mucus fibers and the subsequent diffusive excursion until encountering a new fiber. This mechanistic model can be studied using probabilistic theory to obtain a theoretical expression of the phage subdiffusion that facilitates the interpretation of experiments and the calibration of computational models.
The time of one production step is τ = τ W + τ D , where τ W is the waiting time spent by the phage adhering to a mucus fiber, and τ D is the time spent diffusing between the release from the mucus fiber and binding to another fiber. In this context, the probability of finding a phage between x and x + dx at time t is which is the basic CTRW transport equation. Here, t = ∑ n i=1 τ i is the time after n production steps, P(t; n) is the probability of producing t in n steps, and P n (x) is the probability density function (pdf ) of being at x after n diffusive steps. To evaluate the probability of finding a phage at x at time t, one needs to consider all possible trajectories, that is, any possible number of steps n, which leads to the sum in Eq. (S.2). The goal is to estimate the mean square displacement, which will require to estimate the first and second moments of x, x(t) = ∞ −∞ xP(x,t)dt and x 2 (t) = ∞ −∞ x 2 P(x,t)dt. As shown below, it is possible to obtain an analytic expression for the MSD in the asymptotic limit, that is, for t large.
Renewal theory allows to express the pdf P n (x) and the probability P(t; n) in terms of underlying physical processes.
The total displacement, x, is the sum of n displacements, ∆x so that x = ∑ n i=1 ∆x i . If the displacements ∆x i are independent,

S.3/S.16
identically distributed processes (i.i.d.) generated by the pdf f (∆x), then P n (x) is given by Applying the convolution theorem for Fourier transforms givesP n (k) = f (k) n , wherẽ . This result shows in Fourier space the pdf of the total displacement is just the product of the underlying pdf s for individual displacements. Thus, the goal pdf P(x,t) in Eq. (S.2) can be expressed in Fourier space as Recall that P(t; n) is the probability of having n events within (0,t), that is, 0 < t 1 < t 2 < · · · < t n < t < t n+1 , where the waiting time is defined by consecutive events variables with pdf ψ(τ), which is called the renewal process. The probability P(t; n) can be expressed as an indicator function The stochastic variable t is larger than zero, so it is convenient to use the Laplace transform to achievê In the first step of the above result, the linear properties of the average ensemble operator, · is applied and the fact that the indicator function is one within t n and t n+1 but zero elsewhere. Noticing that t n = ∑ n i=1 τ i for the i.i.d. random variables τ i , the average terms on the right hand side can be expressed in terms of the Laplace transform of the renewal process: where the same applies for the term t n+1 so that e −st n+1 = (ψ(s)) n+1 . Applying this result in Eq. (S.3), giveŝ Now, it is possible to express the goal pdf, P(x,t), in Eq. (S.2), in terms of the underlying mechanistic processes. Applying the Fourier and Laplace transforms to P(x,t) leads tỗ which is a geometric series that converges under the reasonable assumption |ψ(s)f (k)| < 1. The initial pdf s, ψ(t) and f (∆x) are normalized, which implies |ψ(s)| ≤ 1 and |f (k)| ≤ 1, where the only cases have the equality would be for delta functions.
Thus, except for single delta distributions the geometric series above will converge, leading tỗ .
which can be used to calculate the moments for the total displacement, which is the main goal of this derivation.
The expansion of the exponential kernel of the Fourier transform allows to express a pdf in terms of the moments of the variable. The expression is given bŷ This expression also allows to express the moments in terms of the derivatives of the transformed pdf : The particular expressions for the first and second moments are, respectively, Before diving in the calculation of these moments, it is important to emphasize the following properties off (k), which emerge from being the Fourier transform of the pdf f (∆x)): Now, applying the expression ofP(k, s) from Eq. (S.5) in the equations for the first moment x(s) and using the properties stated above, one obtains The prefactor related to the transformed renewal process,ψ(s), can be related to the average number of steps taken, n(s) . The average number of renewal steps within time (0,t) is related to the probability P(t; n) as n(t) = ∑ ∞ n=0 nP(t; n) with Laplace transform: This is a known arithmetico-geometric sequence, ∑ ∞ n=0 nr n = r/ (1 − r) 2 , where r = |ψ(s)| and convergences when |ψ(s)| < 1, S.5/S.16 which, as discussed above, is a common condition for most renewal processes. Then, Going back to the expression for x(s) , this leads to x(s) = n(s) ∆x which can be inverted using the Laplace transform to obtain the first moment of the total displacement which is a very intuitive result. On average, the total displacement is the average number of steps produced within (0,t) times the average individual displacement. In the absence of drift, that is, ∆x = 0, the total average displacement is zero, This is a relevant result for the research presented here.
Let's address now the second moment, which is necessary to determine the MSD and obtain a mechanistic expression for subdiffusion in the context of CTRW. By applying the expression ofP(k, s) given by Eq. (S.5) in the equation for x(s) , The prefactor of the first term of in the right hand side is n(s) , as derived above. The prefactor in the second term can be also related to averages ofn(s). The following technique illustrates this and is also convenient when generating higher order moments ofx(s). The key is to focus on the following generating function h(r) = 1/1 − r = 1 + r + r 2 · · · = ∑ ∞ n=0 r n , which can be obtained by expanding h(r) in a Taylor series around r = 0. Taking the first derivative of h(r) and multiplying both sides by r gives ∑ ∞ n=0 nr n = r/(1 − r) 2 , which is the arithmetico-geometric series used to obtain the first moment ofn(s) above. An analogous process with the second derivative of h(r) leads to ∑ ∞ n=0 n(n − 1)r n = 2r 2 /(1 − r) 3 . Using these expressions in the equation for x 2 (s) , one obtains which can be expressed in terms of series with the identities derived above The Laplace transform of this expression leads to which was the main goal of this calculation. Notice that this general result has a neat physical interpretation. The second moment of the displacement is determined by two contributions. The first term is related to the fluctuations of the local displacement ∆x and the second term to the fluctuations of the number of steps produced within (0,t).
Now the only step missing is to assign the type of local displacement, that is, the pdf f (∆x), as well as the renewal process, that is, the pdf ψ(τ). In the research problem addressed here, the excursion of phages are assumed to be diffusive and have no drift. Then where D is the diffusion constant and τ D is a set diffusive time, which, as shown later, it needs to be adjusted to reproduce the empirical subdiffusion of phages. Thus, the moments of x(t) reduce to x(t) = 0 and x 2 (t) = 2Dτ D n(t) . In this case, the MSD is The final task of the calculation reduces to estimate n(t) . The goal is to obtain a mechanistic expression for the subdiffusion, that is, where one can relate the subdiffusion exponent α and the prefactor K α with the underlying microscopic properties of the system. From the expression derived above for MSD(t), the objective is to obtain an average number of steps within (0,t) that is proportional to time to the α, that is: n(t) ∝ t α . By reverse engineering, one can show that, in the asymptotic limit (large t), this result emerges from a renewal process following a power law distribution, that is,

S.7/S.16
Taking the above as an ansatz it is possible to derive the expressions. The renewal process ψ(τ) is a pdf which must be normalized. The subdiffusion exponent is positive by construction, α > 0. This forces to define a minimum time, τ 0 , to guarantee the normalization of ψ(τ) expressed as ∞ τ 0 ψ(τ)dτ = 1. Then, where A = ατ α 0 . Mathematically it is more neat to take τ 0 as the unit time. But physically, it is convenient to keep it explicit during the derivation to keep track of the units and interpret better the assumptions that will be involved. The average number of steps in the Laplace space, n(s) , depends only on the Laplace transform of the renewal process, Eq. (S.7). In this case, the Laplace transform of ψ(t) isψ(s) = ατ α 0 ∞ τ 0 τ −1−α e −sτ dτ. If the moments of the renewal process were defined, one could use an analogous expansion of the exponential kernel to expressψ(s) as done for the Fourier transform case. Unfortunately, for the cases of interest here 0 < α < 1, all the moments diverge. Alternatively, it is possible to find an asymptotic form that followŝ ψ(s) ∼ 1 +Cs α o(s) for small s which is key to obtain the subdiffusion form of the MSD. To achieve this, it is convenient to express this integral in terms of gamma functions by applying the change of variables u = sτ and integrate by parts to get .
Applying this in the renewal process transform, one obtainŝ In the asymptotic limit for sτ 0 small (large times), one can expand the exponential in a Taylor series and compare terms. The exponent is constrained, 0 < α < 1, so the smallest order of the expansion is (sτ 0 ) α . Then the asymptotic expression for the renewal process iŝ

S.8/S.16
Applying this result in Eq. (S.7), one obtains the asymptotic result n(s) ∼ 1 The inverse Laplace transform of n(s) leads to Using the the reflection relation The unit term on the right hand of the first equation is relevant only when the power exponent α is infinitely small: lim α→0 n(t) ∼ 0. Thus, when the renewal process has an exponent that tends to one, the effective number of average steps is null. Notice that ψ(t) = A/τ can't be normalized because its integral is a logarithm and its upper limit does not decay to zero, that is, it is infinite. Physically, this translates into infinite waiting times, that is, no production steps.
For situations when t τ 0 , it follows (τ 0 /t) α ≈ 0 so that sinc (α) (τ 0 /t) α . In this case, Now, in the case for when t = τ D , by construction, the particle diffuses one step ( n(t) = 1) giving the relation This provides all the elements to obtain the mean squared displacement in the absence of drift, Eq. (S.8),

S.9/S.16
This result has been calculated for a single dimension. If the same process applies in multiple dimensions then, one must add a similar term for each dimension, for example, in 3D r 2 (t) = x 2 (t) + y 2 (t) + z 2 (t) = 3 x 2 (t) = d x 2 (t) with d = 3. Then, This provides a mechanistic way to interpret the generalized diffusion coefficient, K α , in Eq. (S.9), Solving for τ 0 in the above via Eq. (S.11) gives the remarkable formula In this expression, in the limit α → 1, the generalized diffusion recovers the regular diffusion prefactor in the MSD, that is, lim α→1 D α = 2dD. In this limit, the time corresponds to the number of steps diffusing, that is, t → nτ D . Then, This provides a very neat interpretation of the subdiffusion coefficient: It is just the regular diffusion coefficient divided by the diffusion time to the alpha. This provides also a mechanistic way to write the MSD subdiffusion equation: (S.14) This offers a very insightful look at subdiffusion: It's just regular diffusion rescaled by a power function of α for time in units of the diffusion time. The dimension (d) and diffusion (D) can usually be obtained by independent experiments or first principles.
Then, by fitting the MSD data to the equation above, one can obtain an empirical value for the diffusive time, τ D . This, in turn, is related to the minimum waiting time, τ 0 by Eq. (S.11). Thus, Eq. (S.14) provides a direct way to extract the molecular times of a subdiffusive process assuming CTRW. These timescales can be combined with molecular models (binding, polymeric separation) to extract even more valuable information from (sub)diffusion data.
Let's now show that the results obtained above are self-consistent. In particular, for α → 1, it was claimed that t → nτ D , that is, the majority of steps tend to be diffusive. One way to confirm this, is to calculate the probability of generating a waiting S.10/S.16 time, τ w , smaller than the diffusion time, τ D : where the identify in Eq. (S.11) has been used in the second equality. The probability of generating a waiting time is smooth monotonous function that grows from zero (α → 0) to one (α → 1). Thus, as the subdiffusive exponent approaches to one, the diffusive time will dominate, showing that the identify derived above for both time scales is self-consistent.
Another straightforward way to confirm self consistency is to use the identity in Eq. (S.11) to solve for τ 0 in Eq. (S.10) to get Letting α → 1 gives t ∼ n(t) τ D which shows the identity derived above for both time scales is self-consistent.

Supplementary Diffusive Encounter Rate Approximation of Two Motile Species
Here we develop a diffusive encounter rate model for phage-bacteria interactions based on their motility patterns. We assume the bacteria undergoes a random walk w w w B (t) while the phage has a random walk w w w P (t). As illustrated in panel (b), adopting a rest frame centered on the bacteria, we can form an equivalent random walk w w w(t) = w w w B (t) − w w w P (t). The square of the distance of the phage from the bacteria is r 2 = (x B − x P ) 2 + (y B − y P ) 2 + (z B − z P ) 2 . Without loss of generality, the mean of r 2 can be evaluated by first expanding the x-term: (x B − x P ) 2 and taking the mean · to get Assuming independence for both processes, we have x B x P = x B x P = 0 and one gets where v and τ is the average run velocity and time. For the phage x 2 P = D α t α where D α and α are the phage's diffusion coefficient and exponents. Thus, in three dimensions, we have The time rate of change of the variance of w w w is the diffusivity D(t) which is calculated as, Under the rest frame described above, an encounter model for the phage bacteria system can be treated exactly analogously to the flux of diffusing material to a perfectly absorbing sphere. So, if P(r,t) is the concentration of phage as a function of distance r from the center of the bacteria's spherical capture zone with radius R C then (S.20) subject to the boundary conditions that P(R C ,t) = 0 (all phage are consumed at the surface of the bacteria's capture zone) and P(r → ∞,t) = P ∞ , that at distances far from the bacteria the phage concentration approaches a constant background value. The initial conditions are P(r, 0) = 0 for r < R C and P(r, 0) = P ∞ when r ≥ R C .
Solution. In Eq. (S.20), since D(t) does not depend on r, it can be moved outside the differentials. The problem now more formally reads as with specified boundary and initial conditions: To transform the 2-D problem into a more workable 1-D framework, we use the substitution v(r,t) = P(r,t)r. Now, taking first and second partial derivatives with respect to r gives Taking the derivative of v with respect to time in our substitution gives Expanding Eq. (S.21), and substituting Eq's (S.22) and (S.23) gives with specified boundary and initial conditions: A transform of variables; z = r − R C (dz = dr), gives v(r,t) → w(z,t) and z ≥ 0, if and only if r ≥ R C . For the initial conditions we have w(z,t) = P ∞ z + P ∞ R C for z ≥ 0. The problem is now transformed into solving the diffusion equation on a semi-infinite interval (z ≥ 0). A standard way for solving the above problem, which can be found in most texts on partial differential equations, is to first introduce an odd extension of w(z,t) as −w(−z,t) = w(z,t) in the region z < 0. Next, using the method of Fourier transforms, one eventually arrives at the fundamental solution where we specially note, While we have solved the problem, we need to transform the solution back into a usable form (into an expression involving the phage concentration). The physical solution of the phage concentration exists when r ≥ R C which corresponds to z ≥ 0. So, it makes sense to inspect our solution by separating the integrals in Eq. (S.25), to get
S.14/S.16 √ π/2. So, we are left with Finally, since v(r,t) = P(r,t)r we have P(r,t) = P ∞ 1 − R C r + P ∞ R C r erf R C − r 4λ (t) The phage flux (encounter rate) J (t) is then given by for which we can evaluate at r = R C to get the time dependent solution, This solution shows the phage encountered varies λ (t) 1/2 for small t and for t 1 asymptotically tends to the steady state encounter rate where D ≈ v 2 B τ R /3 for α < 1 and D ≈ v 2 B τ R /3 + D for α = 1. Upon dividing by P ∞ gives the encounter kernel (S.29) which is essentially the maximum clearance rate of the phage. S.16/S.16