A geometric criterion for the optimal spreading of active polymers in porous media

Abstract

Efficient navigation through disordered, porous environments poses a major challenge for swimming microorganisms and future synthetic cargo-carriers. We perform Brownian dynamics simulations of active stiff polymers undergoing run-reverse dynamics, and so mimic bacterial swimming, in porous media. In accord with experiments of Escherichia coli, the polymer dynamics are characterized by trapping phases interrupted by directed hopping motion through the pores. Our findings show that the spreading of active agents in porous media can be optimized by tuning their run lengths, which we rationalize using a coarse-grained model. More significantly, we discover a geometric criterion for the optimal spreading, which emerges when their run lengths are comparable to the longest straight path available in the porous medium. Our criterion unifies results for porous media with disparate pore sizes and shapes and for run-and-tumble polymers. It thus provides a fundamental principle for optimal transport of active agents in densely-packed biological and environmental settings.

Introduction

Microorganisms display agile motility features to optimize their survival strategies and efficiently navigate through their natural disordered and porous habitats^1,2,3,4,5. While locomotion by swimming represents the most prominent, inevitable transport feature of many microorganisms, sudden changes of their swimming direction are also an essential tool for their efficient search for nutrients⁶ or escape from harmful environments⁷. These reorientation events are generated by intrinsic biophysical mechanisms and generate different swimming modes, such as the run-and-tumble motion of Escherichia coli³ or Bacillus subtilis⁸, run-reverse(-flick) patterns of diverse bacteria^9,10, sharp turns in swimming algae¹¹, and run-reverse behavior of different species of archaea¹². In unconfined media these transport features lead to trajectories reminiscent of a random walk, yet their consequences for the navigation through real, porous environments, characteristic of a wide variety of biological, biomedical, and environmental contexts, such as biological gels and tissues or environmental soils and sediments, remain largely unexplored. Understanding the underlying physical mechanisms is thus paramount for revealing fundamental microbiological processes, such as biofilm formation and community ecology^13,14, and has significant potential to enable novel nanotechnological applications^15,16.

Engineering the propulsive mechanisms of microorganisms has proven to be a promising route towards the design and development of smart, self-propelled cargo-carriers^17,18 that overcome several limitations of their passive counterparts (e.g., ordinary colloids). Yet, their ability to self-propel might not suffice to make them generally suitable for performing complex tasks in biomedical and environmental settings, where, for example, they may be expected to deliver drugs to a specific target¹⁹, penetrate the porous structure of tumors^20,21, or find and induce degradation of contaminants^16,22. In fact, randomizing the swimming direction of these autonomous agents could be an efficient strategy for reaching a target. To date, however, experimental realizations of controlled reorientation of self-propelled synthetic agents are sparse^23,24.

Experiments of biological microswimmers^25,26,27,28 and synthetic, active agents^29,30 in confined, disordered environments are often concerned with near-surface motility. These studies display a range of unusual phenomena, ranging from the circular swimming motion of bacteria near walls^25,26, hydrodynamic trapping^29,30, and enhancement of bacterial transport near surfaces due to the presence of obstacles²⁷. Similarly, theoretical studies on active transport in crowded environments mainly focus on 2D models for active, point-like particles moving in periodic structures³¹, on a lattice with obstacles³², or disordered environments³³. Accounting for elongated shapes, Brownian dynamics simulations of self-propelled flexible polymers have revealed subdiffusive motion in 2D porous media³⁴. Quantitative studies of active transport in 3D porous media, however, are sparse and it was only recently that the hop and trap mechanism of individual E. coli cells moving in a 3D porous structure was identified^4,5.

While such studies shed light on how pore-scale confinement influences bacterial motility, it is still unclear what motility patterns are optimal for spreading in porous media. A clue comes from the seminal work of Wolfe and Berg¹, who studied the spreading of engineered bacteria, which lacked the ability to sense chemical gradients and whose tumbling rate could be controlled chemically. Their experiments indicated that smooth swimming strains of E. coli get stuck in the porous structure of the semi-solid agar, similarly to incessantly tumbling cells, while at an intermediate tumbling rate bacterial transport appeared more efficient. A non-monotonic transport behavior as a function of the tumbling rate has been found also theoretically in 2D systems^32,35,36, where the effects of pore shape³⁶ and mobile obstacles³² have been addressed. However, the underlying optimal transport mechanism, dictated by the interplay of swimming characteristics and geometric features of the 3D porous medium, remains an open question.

In this work, we elucidate the spreading of self-propelled stiff polymers, as model systems for elongated microorganisms, in porous media by performing Brownian dynamics simulations and developing a coarse-grained theory. We demonstrate that reorientation mechanisms are indispensable for efficient dispersion through porous media, as intrinsic reversals enhance the overall diffusivities by up to two orders in magnitude. We identify a competition between the pore length, a direct measurement of straight pathways available in the pore space, and the run length of self-propelled polymers. In particular, the hopping lengths of rarely and frequently reversing polymers are constrained by the pore length and their intrinsic run length, respectively. Most importantly, maximal spreading occurs when the intrinsic run length of the polymers is comparable to the longest pore length of the porous medium, which allows us to introduce a simple but robust geometric criterion for optimal transport. Subsequently, we rationalize the non-monotonic transport behavior in terms of a renewal theory.

While such a non-monotonic behavior is predicted in refs. ^32,35,36, our study unravels the underlying mechanism and demonstrates that this large-scale non-monotonic behavior persists irrespective of the pore shapes and reorientation mechanism of the polymers and is dictated by the maximal pore length only. These findings together with our geometric criterion provide fundamental, physical insight into earlier experimental observations¹ and thereby should guide the future design of synthetic cargo-carriers, applicable in biomedical and environmental settings.

Results

Model: run-reverse polymer in a porous environment

We model the elongated shape of a bacterial cell by a stiff polymer with aspect ratio L/σ = 5, where L denotes the polymer length and σ the diameter of the individual monomers. We employ Brownian dynamics simulations of the discretized polymer, where the monomers are connected via stiff springs according to the well-established bead-spring model (see Fig. 1a and Methods). The stiffness of the polymer is characterized by a persistence length ℓ_p, which exceeds its contour length L, ℓ_p/L ≫ 1. The self-propulsion of the polymer is modeled by active forces of magnitude F_p acting on each monomer in the direction tangential to its backbone³⁷ (Fig. 1a). The velocity along the contour of the polymer is then determined by the friction coefficient ζ of a single monomer, v_c = F_p/ζ, and each monomer is subject to translational diffusion with diffusivity D₀ = k_BT/ζ, where k_B is the Boltzmann constant and T the temperature. The diffusive time scale of a single monomer is τ₀ = σ²/D₀.

Fig. 1: Model set-up of a run-reverse polymer in a porous environment.

a Sketch of a self-propelled stiff polymer composed of monomers (gray spheres) with active forces acting on each monomer in the direction tangential to the backbone of the polymer. (Zoom) The polymer chain is modeled using a bead-spring model, where r_i is the position of bead i with diameter σ, and t_i,i+1 is the tangent vector between bead i and i + 1. The beads are connected by elastic springs. Furthermore, ${{{{{{{{\bf{F}}}}}}}}}_{p}^{(i)}={F}_{p}({{{{{{{{\bf{t}}}}}}}}}_{i-1,i}+{{{{{{{{\bf{t}}}}}}}}}_{i,i+1})$ denotes the active force acting on bead i. b Schematic of the run-reverse mechanism of the active polymer. At the run-reverse event the polymer randomly reverses its swimming direction, in particular, the active forces randomly change sign, ${{{{{{{{\bf{F}}}}}}}}}_{p}^{(i)}\to -{{{{{{{{\bf{F}}}}}}}}}_{p}^{(i)}$. c Simulation snap-shot of several active stiff polymers immersed in a porous environment composed of overlapping spheres. d Slice of the 3D porous environment. Red circles indicate 2D slices of the porous structure and white areas correspond to the open pore space. e Chord-length distributions ${\varphi }_{{L}_{c}}(\ell )$ for porous environments composed of a different number of N overlapping spheres. Solid lines correspond to an exponential fit of the data. Source data are provided as a Source Data file.

In addition, the active polymer randomly reverses its swimming direction at exponentially distributed times with reversal rate λ (Fig. 1b). At these ‘pseudo’ ‘reversal events’ the polymer either instantaneously changes its swimming direction and moves along the opposite direction or continues swimming along the same direction. The run length of the polymer is then defined as the length the polymer moves before it can reverse: ℓ_run ≡ v_c/λ. A run-reverse mechanism is employed by several microorganisms^9,12, yet for studying the large-scale spreading of active agents in a porous environment we anticipate that it can be used to model run-and-tumble bacteria^3,8 since often the only way to escape narrow pores is by reversing the swimming direction^4,5. We demonstrate this by complementing our findings for run-reverse polymers by run-and-tumble polymers, which change their swimming direction randomly^3,8 (see Methods).

The swimming characteristics of run-reverse polymers can be described by two dimensionless parameters: (1) the Péclet number Pe ≡ v_cL/D₀, which measures the self-propulsion strength relative to diffusion, and (2) the reversal rate λ with respect to the characteristic diffusive time scale of a single monomer τ₀, i.e., λτ₀.

To model the porous medium, we generate a disordered, monodisperse, porous structure composed of N overlapping spheres of size 4σ within a cubic box of length 30σ [Figs. 1c and d for N = 10³]. The micro-architecture of the porous medium can be characterized by the distribution of the straight paths, referred to as chord lengths, available in the medium, ${\varphi }_{{L}_{c}}(\ell )$³⁸. It is shown for different N in Fig. 1e. We further introduce the maximal chord length ${L}_{c,\max }$ via ${\varphi }_{{L}_{c}}({L}_{c,\max })=1{0}^{-5}/\sigma$. Throughout the manuscript we mainly focus on very densely packed environments and hence N = 10³, unless stated otherwise. For this case, the maximal chord length is ${L}_{c,\max }\simeq 20\sigma$ and the medium is characterized by narrow channels of average pore diameter σ_p comparable in size to the individual monomer σ_p/σ ≃ 3.5 (or the polymer σ_p/L ≃ 0.7), as is typical of many natural environments³⁹ (see Methods).

Dynamics: mean-square displacement

To investigate the dynamics of a self-propelled polymer in a porous environment, we measure the mean-square displacement (MSD) of the center monomer $\langle {|{{\Delta }}{{{{{{{{\bf{r}}}}}}}}}_{c}(t)|}^{2}\rangle$ with Δr_c(t) = r_c(t) − r_c(0) and ${{\Delta }}{{{{{{{{\bf{r}}}}}}}}}_{c}(t)={[{{\Delta }}{x}_{c}(t),{{\Delta }}{y}_{c}(t),{{\Delta }}{z}_{c}(t)]}^{T}$. We keep the Péclet number fixed, Pe = 50, and tune the rate of reversal events, λ (Fig. 2). At short times $t\lesssim {\tau }_{{{{{\mbox{diff}}}}}}\equiv {D}_{0}/{v}_{c}^{2}$, the MSDs of active agents with different reversal rates λ collapse and display a linear increase reflecting their diffusive motion at short times, which remains independent of the porous structure of the environment. At intermediate times, the directed motion of the polymers dominates and the MSDs exhibit a superdiffusive increase (for λτ₀ ≳ 5 · 10⁻⁴), which varies for different λ. The MSDs eventually cross over to a linear regime, which characterizes the effective diffusive behavior of the run-reverse polymers in a porous medium.

Fig. 2: Mean-square displacements, $\left\langle {\left|{{\Delta }}{{{{{{{{\bf{r}}}}}}}}}_{c}(t)\right|}^{2}\right\rangle$, of the center monomer of a polymer with different reversal rates λ.

Here, the polymer has 5 monomers and the Péclet number is Pe = 50. Further, τ₀ denotes the diffusion time, ${\tau }_{\,{{{{\mbox{rot}}}}}\,}^{0}$ the rotational relaxation time of a non-tumbling polymer, τ_diff the cross-over time between short-time diffusion and directed swimming motion, and L the polymer length. The colored triangles indicate the characteristic time of tumbling, 1/λ divided by τ₀, for different reversal rates. Source data are provided as a Source Data file.

In contrast, non-reversing polymers can only reorient due to the interplay of thermal fluctuations, activity, and conformational chain dynamics³⁷. The corresponding rotational relaxation time of the polymer is denoted by ${\tau }_{\,{{{{\mbox{rot}}}}}\,}^{0}$ and can be extracted by measuring the fluctuation of the end-to-end direction of the polymer, e(t) = [r₅(t) − r₁(t)]/∣r₅(t) − r₁(t)∣, which evolves as $\langle | {{\Delta }}{{{{{{{\bf{e}}}}}}}}(t){| }^{2}\rangle =2-2\exp (-t/{\tau }_{\,{{{{\mbox{rot}}}}}\,}^{0})$⁴⁰. The MSD for rarely reversing polymers, $\lambda\tau_0=5\cdot10^{-4}$ (and polymers with λτ₀ = 5 ·10⁻³), displays a subdiffusive behavior at intermediate times, $t \sim {\tau }_{\,{{{{\mbox{rot}}}}}\,}^{0}$, where the effect of the porous environment becomes apparent as the polymer slows down. At long times, the polymer has merely moved a distance comparable to its own length over the whole simulation time, which is a fingerprint of confined transport and indicates that motion in tight narrow spaces becomes hindered if active agents cannot reverse efficiently.

Moreover, we find that the cross-over time to long-time diffusion is determined by the faster of the two times: the reorientation time due to reversing 1/λ and the orientational relaxation time ${\tau }_{\,{{{{\mbox{rot}}}}}\,}^{0}$ (indicated in Fig. 2). To quantify the long-time behavior, we extract the effective diffusivities via

$${D}_{\,{{{{\mbox{eff}}}}}}^{{{{{\mbox{sim}}}}}\,}\equiv \mathop{\lim }\limits_{t\to \infty }\frac{\langle | {{\Delta }}{{{{{{{{\bf{r}}}}}}}}}_{c}(t){| }^{2}\rangle }{6t}.$$

(1)

Our findings demonstrate that the long-time effective diffusivities can be enhanced by more than two orders of magnitude (over the whole simulation time of 10³τ₀) upon introducing a reversal mechanism, ${D}_{\,{{{{\mbox{eff}}}}}}^{{{{{\mbox{sim}}}}}}(\lambda {\tau }_{0}=0.5)/{D}_{{{{{\mbox{eff}}}}}}^{{{{{\mbox{sim}}}}}\,}(\lambda {\tau }_{0}=0.0005)\simeq 4\cdot 1{0}^{2}$. This is in stark contrast to the motion of run-reverse polymers in a dilute environment, whose long time transport becomes suppressed upon increasing the reversal rate, ${D}_{\,{{{{\mbox{eff}}}}}}^{{{{{\mbox{sim}}}}}}(\lambda {\tau }_{0}=0.5)/{D}_{{{{{\mbox{eff}}}}}}^{{{{{\mbox{sim}}}}}}(\lambda {\tau }_{0}=0.0005)\simeq 3\cdot 1{0}^{-1}$.

Non-monotonic transport behavior: effective diffusion

Most prominently, the MSDs indicate that the long-time diffusivities ${D}_{\,{{{{\mbox{eff}}}}}}^{{{{{\mbox{sim}}}}}\,}$ of run-reverse polymers display a non-monotonic behavior with respect to the reversal rate λ. We further introduce the scaled path length${{\Lambda }}={L}_{c,\max }/{\ell }_{{{{{\mbox{run}}}}}}$ = ${L}_{c,\max }\lambda /{v}_{c}$, which characterizes the trade-off between the maximal pore length of the medium and the run length of the polymer. We find that the non-monotonic behavior persists for a broad range of Péclet numbers Pe = 25, …75 (Fig. 3a filled symbols and solid lines). Most significantly, the effective diffusivities for all Péclet numbers display a prominent maximum, where the scaled path length is

$${{\Lambda }}=\frac{{L}_{c,\max }}{{\ell }_{{{{{\mbox{run}}}}}}}={{{{{{{\mathcal{O}}}}}}}}(1).$$

(2)

Hence, we propose that equation (2) serves as geometric criterion for optimal transport of active agents in porous media, which occurs when the run length ℓ_run is comparable to the maximal pore length ${L}_{c,\max }$ characteristic of the porous the environment.

Fig. 3: Effective diffusivities of active polymers spreading in porous media.

a Effective diffusivities, ${D}_{\,{{{{\mbox{eff}}}}}}^{{{{{\mbox{sim}}}}}\,}$, extracted from the simulations (filled symbols, solid lines) as a function of the scaled path length Λ determined by the trade-off between the maximal chord length ${L}_{c,\max }$ and the run length of the polymer ℓ_run = v_c/λ, for different Péclet numbers. Dashed lines correspond to the theoretical predictions [equation (5)] of the hop-and-trap model with parameters extracted from individual trajectories. b Effective diffusivities, ${D}_{\,{{{{\mbox{eff}}}}}}^{{{{{\mbox{sim}}}}}\,}$, as a function of the reversal rate λτ₀ for porous environments with different number of obstacles N. c Rescaled effective diffusivities, ${D}_{\,{{{{\mbox{eff}}}}}}^{{{{{\mbox{sim}}}}}\,}$, of run-reverse (RR) and run-and-tumble (RT) polymers for Pe = 50 and different N. The data are rescaled by the maximal pore length ${L}_{c,\max }$ extracted from the chord-length distributions in Fig. 1e. Source data are provided as a Source Data file.

Further, we found that the non-monotonic transport behavior persists for porous media with fewer obstacles, N = 500, 750, but becomes rather weak for dilute environments, N ≲ 500, where it approaches the monotonic behavior of active polymers in an unconfined environment, see Fig. 3b. The chord-length distributions, ${\varphi }_{{L}_{c}}(\ell )$ (Fig. 1e), show a strong decay for densely packed environments.

Furthermore, we have addressed the effect of run-and-tumble motion on the overall spreading of active polymers, where, instead of reversing, the swimming direction after the tumbling event is random. Our results (Fig. 3c) demonstrate that our geometric criterion remains valid. We observe that the overall diffusivities increase for N = 1000 with respect to those of run-reverse polymers, as 3D tumbling allows polymers to spread further. Most importantly, in dense porous environments (N = 1200) the effective diffusivities of run-and-tumble polymers collapse with those of run-reverse polymers. This indicates that the overall spreading in dense porous environments is characterized by hop-and-trap dynamics, which are fully determined by the geometry, irrespective of the re-orientation mechanism. It further emphasizes that indeed the longest pore length of the environment is the characteristic length scale dictating the transport behavior of active polymers and furthermore confirms our geometric criterion.

In contrast to previous work³⁶, we have also found non-monotonic spreading of active agents in porous environments with concave pore shapes (Fig. 5 in Methods). This highlights that solely the interplay of length scales dictates this intricate behavior and not the pore shape.

Our findings offer an interpretation for the observations of the seminal experimental work by Wolfe and Berg¹, who used engineered strains of non-chemotactic E. coli. The tumbling rate and thus the run length of these mutants can be tuned by varying the concentration of an external inducer (isopropyl β-D-thiogalactoside (IPTG)) in the medium. Monitoring the flagella bundles of cells tethered to a surface suggested an increase of the tumbling rate upon increasing the concentration of IPTG. Qualitatively, the experiments of Wolfe and Berg revealed that bacterial swarms that tumble at an intermediate rate spread more than either incessantly tumbling cells or smooth swimming, i.e., non-tumbling, strains. However, the dynamical behavior of individual cells has not been quantified. Therefore, we anticipate that our simulations, which qualitatively confirm these experimental findings, enable us to elucidate the physics of the non-monotonic behavior and provide predictions for the optimal spreading.

Individual trajectories

To elucidate this non-monotonic behavior of the diffusivities, we investigate individual trajectories for Pe = 50 (Fig. 4a). For rarely reversing agents, λτ₀ = 5 · 10⁻⁴, the trajectories in the xy plane indicate highly confined motion, whereas the trajectories of a polymer reversing at higher rates, λτ₀ = 0.5 and 15, cover significantly larger areas including many pores [Fig. 4a (inset)]. By inspecting the associated time evolution of the displacements, we observe that the motion of rarely reversing polymers λτ₀ = 5 · 10⁻⁴, i.e., the black curve in Fig. 4, is characterized by fast hopping events, at which the polymer moves through the pore space, and long, extended phases of trapping, where the polymer is trapped inside a pore. This motility pattern agrees with experimental findings of wild-type E. coli moving in a porous structure^4,5. We further observe that upon increasing the reversal rate the trapping events become shorter, see blue curve for λτ₀ = 0.5 in Fig. 4a. In fact, they disappear for large reversal rates, λτ₀ = 15, where the trajectory is dominated by hopping events, i.e., the red curve in Fig. 4a.

Fig. 4: Polymer trajectories and distributions for the trapping time and hop length.

a Representative (1D) displacements, Δx_c(t), of polymers with different reversal rates, λ, as a function of time for Pe = 50. Horizontal solid lines indicate the trapping phases for different λ. (Inset) Particle trajectories of the center monomer of a rarely reversing polymer, λτ₀ = 5 · 10⁻⁴ (black ; magnified by a factor of 4), and moderately reversing polymers, λτ₀ = 0.5 and 15 (blue and red, respectively). The trajectories are shown in the xy plane, (Δx_c(t), Δy_c(t)). b Trapping time distribution φ_T(t) and c hop length distribution ${\varphi }_{{\ell }_{H}}(\ell )$ for Pe = 50 and different reversal rates, λ, extracted from individual trajectories. The inset in b shows the fraction of time spent hopping p = τ_H/(τ_H + τ_T) and the average trapping time τ_T/τ₀ as a function of λτ₀. The black solid line in panel c indicates the chord length distribution of the medium ${\varphi }_{{L}_{c}}(\ell )$. The inset in c shows the average hopping length 〈ℓ_H〉 normalized by the average chord length 〈L_c〉 of the porous medium. Source data are provided as a Source Data file.

Now, we explain the mechanism for the trapping and hopping, which identifies the main features necessary for an active agent to explore a porous medium. Suppose the active agent is exploring the porous environment with an intrinsic run time τ_run ~ 1/λ and after a while it enters a dead-end pore at a time t. Consequently, the agent will be trapped within the dead-end pore for the remaining time interval [t, τ_run]. This implies that the trapping time of a self-propelled agent can potentially be reduced by increasing the reversal rate λ. The exemplary trajectory for λτ₀ = 0.5, which corresponds to the optimal reversal rate, exhibits the largest displacements with few trapping events (Fig. 4a). Moreover, we observe that despite the sparse and short trapping events of frequently reversing polymers, $\lambda\tau_0=15$, their hops become significantly shorter.

We quantify this behavior by extracting the distributions (1) φ_T(t) for the trapping time, i.e., the time the polymer spends trapped inside a pore, and (2) ${\varphi }_{{\ell }_{H}}(\ell )$ for the hop length, i.e., the distance the polymer moves from one trapping event to another or from one trapping event to the next reversal, from the individual trajectories (see Methods).

Trapping time distribution

We find that the trapping time distributions for active polymers with Pe = 50 exhibit a power-law scaling at long times with an exponent that increases from 2 to 7/2 for increasing reversal rates λ (Fig. 4b). To rationalize this behavior, we extend the concept of an ‘entropic trap’ model introduced in the study of diffusing long polymers escaping small outlets⁴¹ and recently for E. coli bacteria moving in a porous medium^4,5. In this model, the escape of an active polymer from an obstruction in the porous medium is determined by the number of orientations (Ω_t) keeping it trapped inside the pore and the number of orientations (Ω_e) allowing the polymer to leave it. In short (see Methods), our model predicts a power-law trapping time distribution φ_T(t) ~ t^−(1+β) at long times. Here, the exponent β = X/C₀ is related to the active energy X, which has dimensions of F_pL, and the average depth of the entropic trap C₀, which can be quantified by the average free energy difference between the two states: ${C}_{0}={k}_{B}T\langle {{{{{{\mathrm{ln}}}}}}}\,({{{\Omega }}}_{t}/{{{\Omega }}}_{e})\rangle$, where the brackets correspond to the ensemble average over all pores.

In accord with this phenomenological model, we find that the exponent β increases for increasing reversal rate λ (Fig. 4b). Specifically, at a fixed Pe (corresponding to X = const.), a larger reversal rate increases the probability Ω_e of leaving the (e.g., dead-end) pores and leads to a relatively lower trap depth C₀, and thus a larger exponent β = X/C₀. Further, we note that for β≤1, corresponding to the case where the active energy becomes smaller than the trap depth X ≤ C₀, the entropic trap model predicts a divergence of the mean trapping duration τ_T ~ (β−1)⁻¹. This becomes evident in the monotonic increase of the mean trapping duration by orders of magnitude with decreasing reversal rates (Fig. 4b (inset) open circles).

Hopping length distribution

We further extract the distributions of the hopping length of the agents ${\varphi }_{{\ell }_{H}}(\ell )$ (Fig. 4c for Pe = 50) and rationalize that the hop lengths are determined by the interplay of the intrinsic run length of the active polymers ℓ_run ≡ v_c/λ and the pore geometry. In particular, we find that the chord length distribution (Fig. 4c black solid line) is approached by the hopping length distribution of rarely reversing polymers, λτ₀ = 5 · 10⁻⁴ (and 5 · 10⁻³). These agents have an intrinsic run length of ℓ_run ≳ 10³σ (and 10²σ), which is much longer than the longest available straight pathways, ${L}_{{{c}},\max }\simeq 20\sigma$, and therefore their hopping motion is fully determined by the chord lengths of the porous medium (Fig. 4c). The distributions also show longer hops, however, the probabilities are rather low. Further, the average hop length is comparable to the average chord length for small λ, $\langle\ell_H\rangle\sim\langle L_c\rangle$ [Fig. 4c (inset)]. For large reversal rates λτ₀ ≳ 5 the distribution decays faster than that of polymers with an intermediate reversal rate λτ₀ = 0.5 and the average hop length approaches the run length of the polymer, 〈ℓ_H〉 → ℓ_run.

Our findings further show that polymers with an intermediate reversal rate, e.g., λτ₀ = 0.5, can follow the straight path available in the porous structure and continue to explore another successive pore without getting trapped, which leads to hopping lengths longer than the longest chord length. Furthermore, we find that the probability for longer hops becomes larger at an intermediate reversal rate than the chord-length distribution, which indicates that the polymer explores these more often than the shorter pores of the medium.

Our results also demonstrate that the probability for long hops and also the average hop length 〈ℓ_H〉 [Fig. 4c (inset)] vary non-monotonically with λτ₀: they are lowest for λτ₀ = 5 · 10⁻⁴ and λτ₀ = 15 and highest for λτ₀ = 0.5. This strong non-monotonic variation of the average hop length could explain the optimal spreading and thus the maximal long-time effective diffusivity of run-reverse polymers in porous media (Fig. 3a). In particular, agents with Péclet number Pe = 50 and reversal rate λτ₀ = 0.5 have an intrinsic run length of ℓ_run = 20σ, comparable to the longest pore length ${L}_{{{c}},\max }$. Our findings reveal that the active polymer continuously explores the pore space without getting trapped too often (Fig. 4a), which suggests the qualitative picture: the agent moves through the pores until it reaches an obstruction, where subsequent reversals allow it to continue to explore the pore space.

Coarse-grained dynamics

To quantitatively characterize the long-time effective diffusivities, we develop a coarse-grained model for the 3D dynamics of the alternating hopping and trapping phases. During the hopping phase the agent, modeled as a point particle, moves straight at effective velocity v and reverses its swimming direction at rate λ with propagator ${{\mathbb{P}}}_{H}({{\Delta }}{{{{{{{\bf{r}}}}}}}},t)$, which denotes the probability that the particle displaces Δr during time t while hopping. During the trapping phase the particle cannot move and after the trapping event it hops along a new, random direction. The time the particle spends in a hopping or trapping phase is determined by the hopping and trapping time distributions extracted from the simulations (Fig. 4b, c). As the swim speed of the polymer during the hopping events remains roughly constant, the hopping time, i.e., the duration of a hopping phase of length ℓ, t ~ ℓ/v, also follows an exponential distribution ${\varphi }_{H}(t)=\exp (-t/{\tau }_{H})/{\tau }_{H}$ with mean duration τ_H. The trapping time distribution obeys a power-law behavior φ_T(t) = β(1+t/τ)^−1−β/τ with mean duration τ_T = τ(β−1)⁻¹, where τ is a characteristic time scale for trapping. This power-law behavior suggests that the dynamics are non-Markovian and therefore we use a renewal theory^42,43.

The probability density P(Δr, t) for the particle to have displaced Δr during lag time t is the sum of the probability densities for the particle to be in a hopping or a trapping phase: P(Δr, t) = P_H(Δr, t) + P_T(Δr, t). We further introduce the probability densities (per time) for the particle to start a hopping phase and a trapping phase at displacement Δr and at lag time t by H(Δr, t) and T(Δr, t), respectively. Then the probability P_T(Δr, t) that the particle is trapped at Δr for lag time t is obtained as the sum of the probability of never having hopped before ${P}_{T}^{(0)}({{\Delta }}{{{{{{{\bf{r}}}}}}}},t)$ and the probability of having hopped at least once before getting trapped at Δr and at an arbitrary earlier time $t-t^{\prime}$:

$${P}_{T}({{\Delta }}{{{{{{{\bf{r}}}}}}}},t)=\; {P}_{T}^{(0)}({{\Delta }}{{{{{{{\bf{r}}}}}}}},t) +\int\nolimits_{0}^{t}dt^{\prime} T({{\Delta }}{{{{{{{\bf{r}}}}}}}},t-t^{\prime} ){\varphi }_{T}^{(0)}(t^{\prime} ).$$

(3)

Here, ${\varphi }_{T}^{(0)}(t)=\int\nolimits_{t}^{\infty }dt^{\prime} \,{\varphi }_{T}(t^{\prime} )={[\tau /(\tau +t)]}^{\beta }$ denotes the probability that the trapping time exceeds t. The probability to become trapped at Δr and time t reads

$$T({{\Delta }}{{{{{{{\bf{r}}}}}}}},t)= \; {T}^{(1)}({{\Delta }}{{{{{{{\bf{r}}}}}}}},t)\\ +{\int}_{{\!\!\!\!{\mathbb{R}}}^{3}}\ d{{{{{{{\boldsymbol{\ell }}}}}}}} \int\nolimits_{0}^{t}dt^{\prime} H({{\Delta }}{{{{{{{\bf{r}}}}}}}}-{{{{{{{\boldsymbol{\ell }}}}}}}},t- t^{\prime} ){{\mathbb{P}}}_{H}({{{{{{{\boldsymbol{\ell }}}}}}}},t^{\prime} ){\varphi }_{H}(t^{\prime} ),$$

(4)

where T⁽¹⁾(Δr, t) denotes the probability for the first trapping event. The second term corresponds to the sum over the probabilities that the agent has escaped a trap and started to hop at Δr − ℓ at an earlier time, $t-t^{\prime}$, until it gets trapped again at Δr and t. Similar equations hold for P_H(Δr, t) and H(Δr, t) (see equations (11)–(12) in Methods). We can derive an analytic solution for the probability density in Fourier–Laplace space, $P(k,s)=\int_{0}^{\infty }\ dt\,{e}^{-st}{\int}_{{{{{{{\mathbb{R}}}}}}}^{3}}\ d{{\Delta }}{{{{{{{\bf{r}}}}}}}}\,{e}^{-i{{{{{{{\bf{k}}}}}}}}\cdot {{\Delta }}{{{{{{{\bf{r}}}}}}}}}\,P({{\Delta }}{{{{{{{\bf{r}}}}}}}},t)\equiv \int\nolimits_{0}^{\infty }\ \ dt\,{e}^{-st}\langle {e}^{-i{{{{{{{\bf{k}}}}}}}}\cdot {{\Delta }}{{{{{{{\bf{r}}}}}}}}}\rangle$⁴². By isotropy, it can be expanded for small wave numbers k up to ${{{{{{{\mathcal{O}}}}}}}}({k}^{4})$ via $P(k,s)\simeq \int\nolimits_{0}^{\infty } dt\,{e}^{-st} \left[1 -{k}^{2}\langle | {{\Delta }}{{{{{{{\bf{r}}}}}}}}(t){| }^{2}\rangle /3!\right]={s}^{-1}-{k}^{2}\langle | {{\Delta }}{{{{{{{\bf{r}}}}}}}}(s){| }^{2}\rangle /3!,$ which allows for an analytical derivation of the Laplace transform of the mean-square displacement $\langle | {{\Delta }}{{{{{{{\bf{r}}}}}}}}(s){| }^{2}\rangle.$

Finally, we obtain the long-time transport behavior by taking the limit of $\langle | {{\Delta }}{{{{{{{\bf{r}}}}}}}}(s){| }^{2}\rangle$ for s → 0. We find $\mathop{\lim }\limits_{s\to 0}\langle | {{\Delta }}{{{{{{{\bf{r}}}}}}}}(s){| }^{2}\rangle \simeq 6{D}_{\,{{{{\mbox{eff}}}}}}^{{{{{\mbox{theo}}}}}\,}{s}^{-2}$, which corresponds to a long-time diffusive behavior $\langle | {{\Delta }}{{{{{{{\bf{r}}}}}}}}(t){| }^{2}\rangle \simeq 6{D}_{\,{{{{\mbox{eff}}}}}}^{{{{{\mbox{theo}}}}}\,}t$ for t → ∞ with effective diffusivity:

$${D}_{\,{{{{\mbox{eff}}}}}}^{{{{{\mbox{theo}}}}}\,}=\frac{{v}^{2}{\tau }_{H}^{2}}{3({\tau }_{H}+{\tau }_{T})\left[1+(1-\cos {\vartheta }_{0})\lambda {\tau }_{H}\right]}.$$

(5)

It is determined by the fraction of time spent hopping between traps, τ_H/(τ_H + τ_T), the hopping time, τ_H, the effective velocity, v, and the tumbling rate λ. The tumbling angle for run-reverse motion, as employed here, is⁹ ϑ₀ = − π. The expression [equation (5)] depends on the exponent β of the trapping time distribution via the average trapping time τ_T. We note that for a power-law distribution truncated at rate γ, ${\varphi }_{T}(t) \sim \exp (-\gamma t){(1+t/\gamma )}^{-1-\beta }$, which may be more appropriate to describe our data (Fig. 4b), the effective diffusivity assumes the same form as in equation (5) with average trapping time τ_T ≡ τ_T(γ, τ, β) depending on the parameters γ, τ, and β (see Methods).

To compare the simulation data with our coarse-grained theory, we extract the average hopping and trapping times, τ_H and τ_T, from the individual trajectories and obtain the effective diffusivity from the hop-and-trap model [equation (5)]. We observe that our model captures the non-monotonic trend of the simulation results, ${D}_{\,{{{\mbox{eff}}}}}^{{{{\mbox{sim}}}}} {\sim} {D}_{{{{\mbox{eff}}}}}^{{{{\mbox{theo}}}}}$ (Fig. 3a dashed lines). It is interesting that the hopping and trapping times, used as inputs for the model, can semi-quantitatively predict the non-monotonic trend of the effective diffusivities of a complex system. The deviations at small reversal rates, corresponding to the scaled path length Λ ≲ 10⁻¹, are expected because at such small reversal rates the polymers remain trapped most of the time. Therefore, the extracted mean trapping times may be underestimated and even longer simulations would be required.

We rationalize the non-monotonic behavior by inspecting the fraction of time spent hopping, p = τ_H/(τ_H + τ_T) (Fig. 4b inset). As expected, for rare reversal events the fraction of time spent hopping is small p ≪ 1/2. This corresponds to τ_H ≪ τ_T and therefore equation (5) reduces to ${D}_{{{{{\mbox{eff}}}}}}^{{{\mbox{theo}}}}={v}^{2}{\tau }_{H}^{2}/(3{\tau }_{T})\to 0$, as large trapping times suppress transport. In contrast, the trapping times for frequently reversing polymers are negligible compared to the hopping times, τ_T ≪ τ_H as p ≫ 1/2. Therefore, the behavior of the effective diffusivity is dominated by the reversal rate and equation (5) simplifies to $D_{{{\mbox{eff}}}}^{{{\mbox{theo}}}}\sim v^2/\lambda\to 0$, which vanishes as the reversal rate increases. Hence, maximal transport occurs at an intermediate reversal rate λ, in agreement with our simulation results for the long-time effective diffusivities and the average hopping lengths.

Our results indicate that the non-monotonic behavior of the effective diffusivities is a clear signature of hop-and-trap dynamics. This behavior vanishes for dilute environments, where the dynamics are solely determined by the run-reverse motion of the polymers. The transition between both motility modes is shown in Fig. 3b.

Discussion

To rationalize the seminal experiments of Wolfe and Berg¹, which characterized spreading in a porous environment of bacteria with different tumbling rates, we have performed Brownian dynamics simulations of run-reverse stiff polymers in a porous environment and find a non-monotonic transport behavior as a function of the reversal rate. We introduce, for the first time, a geometric criterion for the optimal spreading of active agents in a porous environment, which occurs when the intrinsic run length is comparable to the longest straight path available in the environment. In particular, we demonstrate that this criterion remains valid in dense environments for different reorientation mechanisms of the active polymers. We further show that individual polymer trajectories exhibit a hop-and-trap mechanism, in accord with recent experiments of E. coli^4,5. Our results suggest that optimal transport at an intermediate reversal rate is characterized by a maximal average hop length. In contrast, the motion of rarely and frequently reversing polymers is set by the pore length or the run-length, respectively. We corroborate these findings using a renewal theory for the coarse-grained hop-and-trap dynamics.

Our results demonstrate that this non-monotonic transport behavior of active agents in a porous environment persists irrespective of the details of the re-orientation mechanism and the shape of the pores. These findings indicate that the ‘size of the pores, not their shape, matters’ for this large-scale spreading and we therefore anticipate that this behavior is universal for densely packed environments. Nevertheless, this non-monotonic behavior fades for low packing fraction, which corroborates earlier predictions for dilute 2D porous environments with concave pore shapes³⁶.

In the future, it will be interesting to study the effect of swimmer shape anisotropy on active transport in a porous environment. In particular, recent work has shown that the non-monotonic behavior persists even for point particles on a 2D lattice with obstacles³². However, the interplay between pore geometry and run length remains to be explored. Most importantly, our study and recent experiments^4,5 predict power-law distributions for the trapping times, which could be amplified by particle shape. Taking into account this memory dependence in a coarse-grained description goes beyond earlier theoretical predictions^32,35.

We anticipate that our theoretical findings can be tested in various biological systems, such as bacteria^3,8,9,10, algae¹¹, or archaea¹². This could shed light on fundamental microbiological processes, which include the adaption of the tumbling behavior under varying environmental conditions. Experimental observations of Bacillus subtilis⁴⁴ have shown that these cells can reverse their swimming direction once they encounter an obstacle and, similarly, Spirochetes⁴⁵ display an increase of reversal rate while invading a heterogeneous fibrous medium. It would be interesting to elucidate if these spatially varying reversal mechanisms allow them to optimize their motion in a 3D porous medium. Similarly, our results could guide the design of future synthetic swimmers with, e.g., specific magnetic properties⁴⁶. Their spreading could be enhanced by reorienting them via externally applied magnetic fields. On the macroscale, our findings could find application in instructing robots during search and rescue operations in disaster zones³⁶.

Beyond porous media, our results lay the foundation for studying transport of active semiflexible polymers in dynamically re-arranging environments composed of other polymers, such as the interior of cells^40,47. The strongly interacting, crowded environment could lead to a relaxation of the extended trapping events of the active agents and entail complex entanglement effects of the individual constituents.

Methods

Brownian dynamics simulations of a stiff, run-reverse polymer in a porous environment

We model the active polymer chain in terms of the well-known bead-spring model in 3D³⁷. In particular, the polymer is a chain of contour length L composed of N_p spherical monomers with diameter σ that are connected by springs (Fig. 1a). The monomers have positions r_i and the distance between two monomers is denoted by r_i,j = ∣r_j − r_i∣, where i, j = 1, …N_p. We further introduce the tangent vector between two monomers by ${{{{\bf{t}}}}}_{i,i+1}=\left({{{{\bf{r}}}}}_{i+1}-{{{{\bf{r}}}}}_{i}\right)/{r}_{i,i+1}$. The dynamics of the semiflexible polymer are described by the equations of motion for each monomer,

$$\zeta \frac{d{{{{{{{{\bf{r}}}}}}}}}_{i}}{dt}=-{\nabla }_{i}U+{{{{{{{{\bf{F}}}}}}}}}_{p}^{(i)}+{{{{{{{{\bf{F}}}}}}}}}_{r}^{(i)}\quad \,{{\mbox{for}}}\,i=1,\ldots {N}_{p},$$

(6)

with friction coefficient ζ, active forces ${{{{{{{{\bf{F}}}}}}}}}_{p}^{(i)}={F}_{p}({{{{{{{{\bf{t}}}}}}}}}_{i-1,i}+{{{{{{{{\bf{t}}}}}}}}}_{i,i+1})$, and stochastic forces ${{{{{{{{\bf{F}}}}}}}}}_{r}^{(i)}$ characterized by zero mean $\langle {{{{{{{{\bf{F}}}}}}}}}_{r}^{(i)}\rangle ={{{{{{{\bf{0}}}}}}}}$ and variance $\langle {F}_{r,j}^{(i)}(t){F}_{r,k}^{(i)}(t^{\prime} )\rangle =2{k}_{B}T\zeta {\delta }_{jk}\delta (t-t^{\prime} )$. The interaction energy, U, is characterized by the interaction between neighboring beads and the interaction between non-neighboring beads of the polymer chain itself and with the porous environment (see later section on details about the porous environment).

The elasticity of the chain is characterized by the well-established wormlike chain (WLC) model pioneered by Kratky and Porod⁴⁸. The discretized interaction energy is

$$\frac{{U}^{{{{{\mbox{WLC}}}}}}}{{k}_{B}T}=-\frac{{\ell }_{p}N}{L}\mathop{\sum }\limits_{i=1}^{{N}_{p}-2}\left(1-{{{{{{{{\bf{t}}}}}}}}}_{i,i+1}\cdot {{{{{{{{\bf{t}}}}}}}}}_{i+1,i+2}\right),$$

(7)

where we have introduced the persistence length ℓ_p of a semiflexible polymer. It is a measure of the decay length of its tangent-tangent correlations and allows distinguishing between flexible, ℓ_p/L ≲ 1, semiflexible, ℓ_p/L ≃ 1, and stiff polymers, ℓ_p/L ≳ 1.

Interactions between neighboring beads are modeled by the finitely extensible non-linear elastic (FENE) potential,

$$\frac{{U}^{{{{{\mbox{FENE}}}}}}}{{k}_{B}T}=-{\epsilon }_{{{{{\mbox{FENE}}}}}}\mathop{\sum }\limits_{i=1}^{{N}_{p}-1}{{{{{{\mathrm{ln}}}}}}}\,\left[1-\left(\frac{{r}_{i,i+1}}{\delta}-\frac{L}{{N}_{p}\delta }\right)^2\right],$$

(8)

which ensures a finite, maximal distance between two monomers, L/N_p + δ with δ = σ/4. For monomer–monomer distances larger than L/N_p + δ the interaction potential diverges, U^FENE = ∞. Consequently, the entire polymer chain has a maximal length of L + N_pδ. The interactions between non-neighboring polymer beads and the interactions between the polymer and the surrounding, porous environment are modeled using the Weeks–Chandler–Anderson (WCA) potential⁴⁹,

$$\frac{{U}_{ij}^{\,{{{{\mbox{WCA}}}}}}}{{k}_{B}T}={\epsilon }_{{{{{\mbox{WCA}}}}}}\left\{\begin{array}{ll}4\left[{\left(\frac{d}{{r}_{ij}}\right)}^{12}-{\left(\frac{d}{{r}_{ij}}\right)}^{6}\right]&r < d\\ 0,\hfill&\,{{\mbox{else}}}\,.\end{array}\right.$$

(9)

Here, d = σ for the interaction of two monomers and d = (σ + σ_s)/2 for the interaction between a monomer and an obstacle of the porous environment with diameter σ_s (see later section on the Porous environment). The full contribution is the sum over all monomers and beads of the porous environment ${U}^{{{{{\mbox{WCA}}}}}}={\sum }_{i,j;i\ne j}{U}_{ij}^{{{{{\mbox{WCA}}}}}\,}$.

The total interaction energy is then U = U^WLC + U^FENE + U^WCA, which includes (amongst others) three important dimensionless parameters: the persistence length of the polymer chain relative to the contour length ℓ_p/L, and ϵ_FENE and ϵ_WCA, which measure the relative importance of the interaction energies with respect to thermal energy. In the present study, we simulate stiff polymers with ℓ_p/L = 50. To assure minimal polymer extension and prevent overlap of monomers and the environment, we set ϵ_FENE = 10³ and ϵ_WCA = 5, respectively. Furthermore, we use a Brownian time step of τ_B = 10⁻⁶τ₀ to integrate equation (6) using a modified version of LAMMPS⁵⁰ with τ₀ the diffusive time scale of a monomer and wait for 10⁹τ_B before taking measurements.

In addition to self-propulsion, the polymer performs a pseudo run-reverse motion, as sketched in Fig. 1b. The reversal events occur at randomly drawn times, t, which follow an exponential distribution $\lambda \exp (-\lambda t)$ with reversal rate λ. At the run-reverse event, the active forces, acting along the polymer chain, randomly change sign, ${{{{{{{{\bf{F}}}}}}}}}_{p}^{(i)}\to \alpha {{{{{{{{\bf{F}}}}}}}}}_{p}^{(i)}$ with α randomly chosen from { − 1, 1}. In our framework the reversal events occur instantaneously, so that the polymer does not stop before moving into the new direction.

Porous environment

The porous structure of the environment is generated by randomly distributed obstacles of diameter σ_s, which are allowed to overlap. By varying the diameter of the obstacles, we tune the average pore diameter, σ_p, of the medium, inspired by the experimental set-up of refs. ^4,5. The porous environment is characterized by the average pore diameter and the chord length. For ensemble averaging, we use 20 statistically independent structures.

We extract the average pore diameter of the medium by monitoring the transport behavior of Brownian particles with radius σ_B in the porous environment. Therefore, we measure the mean-square displacement, $\langle {|{{\Delta }}{{{{{{{{\bf{r}}}}}}}}}_{B}(t)|}^{2}\rangle$ with Δr_B(t) = r_B(t) − r_B(0), and calculate the local exponent, $\alpha (t)=d{{{{{{\mathrm{ln}}}}}}}\,\left[\langle | {{\Delta }}{{{{{{{{\bf{r}}}}}}}}}_{B}(t){| }^{2}\rangle \right]/d{{{{{{\mathrm{ln}}}}}}}\,t$, as a function of time. The local exponent displays a minimum at an intermediate time ${\tau }_{\min }$, which allows introducing the pore diameter by ${\sigma }_{p}={\sigma }_{B}+{\langle | {{\Delta }}{{{{{{{{\bf{r}}}}}}}}}_{B}({\tau }_{\min }){| }^{2}\rangle }^{1/2}$. For simplicity, we choose σ_B = σ.

To measure the chord length distribution ${\varphi }_{{L}_{c}}(\ell )$, we image several two dimensional planes randomly passing through the simulated porous medium (Fig. 1d). These images provide a map of the pore space. We then binarize these images where two phases represent solid obstacles and open pores, respectively. We calculate the distribution of chords of length ℓ, which fit within each binarized pore space image. This protocol yields a direct measurement of straight pathways available in the pore space.

Effect of pore shape

We have addressed the effect of pore shape on the large-scale spreading of active agents. In particular, we have replaced the WCA potential [equation (9)], corresponding to pores with convex boundaries, and modeled the interaction between the polymer and concave pores by the interaction potential:

$$\frac{{U}_{i}^{C}}{{k}_{B}T}=-{\epsilon }_{C}{e}^{-{({R}_{i}/a)}^{12}},$$

(10)

where R_i is the distance between monomer i of the polymer and the nearest obstacle. We choose ϵ_C = 100 and a = 3.7, corresponding to concave voids of diameter ~ 7.4σ, as illustrated in Fig. 5c. We find that also for pores with concave walls the long-time behavior of the mean-square displacement is diffusive, see Fig. 5a. Extracting the long-time effective diffusivities shows that the non-monotonic behavior persists as a function of the reversal rate (Fig. 5b). The maximal diffusivity occurs at a reversal rate of λτ₀ = 0.5 and a run length of ℓ_run = 20σ. This optimal behavior agrees with the findings for a porous environments with convex pore shapes, as discussed in the main text. Overall, the effective diffusivities are smaller than for a convex environment, which may indicate that the polymer explores the individual pores for a longer time and should depend on the energy depth ϵ_C.

Fig. 5: Active polymers in a porous environment with concave pore shapes.

a Mean-square displacements $\langle |\Delta{{{\bf{r}}}}_c(t)|^2\rangle$ of active polymers as a function of time for different reversal rates, λτ₀, and Péclet number Pe = 50. b Effective diffusivities ${D}_{\,{{{{\mbox{eff}}}}}}^{{{{{\mbox{sim}}}}}}$ as a function of ${\Lambda }={L}_{c,\max }/{\ell }_{{{{\mbox{run}}}}}$. c Slice of the 3D porous environment. Red areas indicate 2D slices of the porous structure and white areas correspond to the open pore space. Source data are provided as a Source Data file.

Effect of the re-orientation mechanism

We further study the effect of the re-orientation mechanism of active polymers on their spreading in a porous environment. To complement our findings for run-reverse polymers, we consider a run-and-tumble polymer, which randomly changes its swimming direction at exponentially distributed tumbling events with tumbling rate λ (Fig. 6a). At the tumbling event the active forces may change sign F_p → αF_p with α randomly chosen from { − 1, 1} and the polymer is subject to tumbling torques. In particular, only at the tumbling event (hence during one time step) opposite forces ${{{{{{{{\bf{F}}}}}}}}}_{t}^{0,1}$ and ${{{{{{{{\bf{F}}}}}}}}}_{t}^{3,4}$ are applied to the zeroth  and first  and third and forth monomers, respectively (Fig. 6b), so that the polymer gets a random ‘kick’ leading to its re-orientation.

Fig. 6: Run-and-tumble polymers in a dilute environment.

Sketch of the a run-and-tumble motion and b tumbling mechanism of a stiff polymer. Here, the unit vectors ${{{{{{{{\bf{n}}}}}}}}}_{0,1}^{(1)}$ and ${{{{{{{{\bf{n}}}}}}}}}_{0,1}^{(2)}$ span the plane normal to the tangent t_0,1. The normal vector n_0,1 lies in this plane. c Trajectories of three run-and-tumble polymers (solid lines) and run-reverse polymers (dotted lines) with vanishing translational diffusivity, D₀ = 0, and run-length ℓ_run/σ = 20. Source data are provided as a Source Data file.

These ‘random’ tumbling forces are chosen (on average) perpendicular to the backbone of the polymer: ${{{{{{{{\bf{F}}}}}}}}}_{t}^{0,1}={F}_{t}({{{{{{{{\bf{n}}}}}}}}}_{0,1}+{{{{{{{{\bf{n}}}}}}}}}_{1,2})/2$ and ${{{{{{{{\bf{F}}}}}}}}}_{t}^{3,4}=-{F}_{t}({{{{{{{{\bf{n}}}}}}}}}_{2,3}+{{{{{{{{\bf{n}}}}}}}}}_{3,4})/2$, where the force magnitude is F_t = 2 ⋅ 10⁵k_BT/σ ≫ F_p. As the polymer is stiff, we note that the forces approximately balance ${{{{{{{{\bf{F}}}}}}}}}_{t}^{0,1}+{{{{{{{{\bf{F}}}}}}}}}_{t}^{3,4}\simeq {{{{{{{\bf{0}}}}}}}}$.

To obtain the component of the force direction n_0,1 (and similarly n_1,2, n_2,3, n_3,4), we first define the plane normal to the tangent vector t_0,1, which is spanned by the unit normal vectors ${{{{{{{{\bf{n}}}}}}}}}_{0,1}^{(1)}$ and ${{{{{{{{\bf{n}}}}}}}}}_{0,1}^{(2)}={{{{{{{{\bf{t}}}}}}}}}_{0,1}\times {{{{{{{{\bf{n}}}}}}}}}_{0,1}^{(1)}$ (Fig. 6b). Then we choose n_0,1 as a random direction in this plane via ${{{{{{{{\bf{n}}}}}}}}}_{0,1}=\cos (\beta ){{{{{{{{\bf{n}}}}}}}}}_{0,1}^{(1)}+\sin (\beta ){{{{{{{{\bf{n}}}}}}}}}_{0,1}^{(2)}$ with β drawn from a uniform distribution ${{{{{{{{\mathcal{U}}}}}}}}}_{[0,2\pi ]}$.

To illustrate our algorithm, we show typical trajectories for polymers with vanishing translational diffusivity, D₀ = 0, and run-length ℓ_run/σ = 20, in Fig. 6c. The trajectories (solid lines) demonstrate that the polymer changes its swimming direction randomly at tumbling events and tumbles in 3D. We note that in this case reversing polymers only move back and forth along a straight line (dotted lines).

Data analysis of the individual trajectories

We extract the distributions for the trapping time, hopping time, and the hopping length from the individual trajectories. We note that the hopping time is defined as the time between hopping from one trap to the next or from one trapping event to the next reversal and the hopping length is the length the polymer moved during this time. To extract these quantities, we follow the approach from refs. ^4,5 and first measure the average velocity of a non-tumbling polymer in a free environment, 〈v〉. Then we calculate the instantaneous velocities of the center monomer, v_i = ∣r_c(t_i+1) − r_c(t_i)∣/(t_i+1 − t_i), where t_i corresponds to the i-th time step, of individual particle trajectories. Thus, we can classify a hopping phase by v_i≥〈v〉/3 and a trapping phase by v_i < 〈v〉/3, which allows extracting the trapping and hopping time distributions, φ_T(t) and φ_H(t) with mean durations, τ_T and τ_H. In addition, we keep track of the reversal times, which are input to our simulations, and compare it with the hopping duration. This provides the hopping length distribution, ${\varphi }_{{\ell }_{H}}(t)$ with mean hopping length 〈ℓ_H〉. In particular, if the reversal time since the last trapping event t_λ is shorter than the hopping phase t_hop, the hopping length corresponds to the length displaced until the reversal event. For t_λ > t_hop the hopping length is the displacement from one trapping event to the next. For the comparison of the effective diffusivities extracted from simulations, ${D}_{\,{{{{\mbox{eff}}}}}}^{{{{{\mbox{sim}}}}}\,}$, to the predictions of the hop-and-trap model [equation (5)], we further use as effective velocity the cut-off velocity, v = 〈v〉/3. We can fully recover the non-monotonic behavior and explain the simulation data up to a constant pre-factor.

To test our approach, we have varied the cut-off velocity between 〈v〉/5 to 〈v〉/2 and found that it does not change our conclusions: the power-law behavior of the trapping time distributions, the behavior of the hopping length distributions, and the semi-quantitative agreement of the effective diffusivities of the coarse-grained model and the simulations remain preserved.

Entropic trap model for the trapping time distribution

Motivated by other disordered media^4,5,41, the probability for the entropic trap C is assumed to follow $P(C)={C}_{0}^{-1}\exp (-C/{C}_{0})$ with average trap depth C₀. By analogy to equilibrium physics, the probability of an active polymer to escape a trap of depth C is assumed to obey an Arrhenius-like relation. Then the trapping duration is given by $t=\tau (\exp (C/X)-1)$, where X characterizes the active energy due to its swimming motion and τ corresponds to a characteristic time scale for trapping. In particular, the trapping duration vanishes, t → 0, for small entropic traps, i.e., C/X → 0. For a passive polymer the active energy X is replaced by the thermal energy k_BT. The probability distribution of the trapping times can be obtained as φ_T(t) = P(C)(∂t/∂C)⁻¹ = β(1+t/τ)^−1−β/τ with β = X/C₀.

Renewal theory for the hop-and-trap dynamics

The probability density for a particle to be in a hopping phase follows:

$${P}_{H}({{\Delta }}{{{{{{{\bf{r}}}}}}}},t)={P}_{H}^{(0)}({{\Delta }}{{{{{{{\bf{r}}}}}}}},t)+ \int_{{{\mathbb{R}}}^{3}} d{{{{{{{\boldsymbol{\ell}}}}}}}}\ \int_{0}^{t}dt^{\prime} H({{\Delta }}{{{{{{{\bf{r}}}}}}}}-{{{{{{{\boldsymbol{\ell }}}}}}}},t-t^{\prime} ){{\mathbb{P}}}_{H}({{{{{{{\boldsymbol{\ell}}}}}}}},t^{\prime} ){\varphi }_{H}^{(0)}(t^{\prime} ).$$

(11)

Here, ${P}_{H}^{(0)}({{\Delta }}{{{{{{{\bf{r}}}}}}}},t)$ denotes the probability that the particle has never been trapped before and the second term corresponds to the sum over all hopping phases, which started after at least one trapping event. Further, ${\varphi }_{H}^{(0)}(t)=\int\nolimits_{t}^{\infty }dt^{\prime} \,{\varphi }_{H}(t^{\prime} )$ is the probability that the hopping time exceeds t. The probability density (per time) that a new hopping phase starts obeys the equation of

$$H({{\Delta }}{{{{{{{\bf{r}}}}}}}},t)={H}^{(1)}({{\Delta }}{{{{{{{\bf{r}}}}}}}},t)+\int\nolimits_{0}^{t}\ \ dt^{\prime} \,T({{\Delta }}{{{{{{{\bf{r}}}}}}}},t-t^{\prime} ){\varphi }_{T}(t^{\prime} ),$$

(12)

where H⁽¹⁾(Δr, t) is the probability that the particle starts the first hop. After a Fourier transform of the probability densities, Δr → k, and by the convolution theorem, the renewal equations [Eqs. (3), (4), (11), (12)] simplify to

$${P}_{T}(k,t)={P}_{T}^{(0)}(k,t)+\int\nolimits_{0}^{t}dt^{\prime} T(k,t-t^{\prime} ){\varphi }_{T}^{(0)}(t^{\prime} ),$$

(13a)

$${P}_{H}(k,t)={P}_{H}^{(0)}(k,t)+\int\nolimits_{0}^{t}dt^{\prime} H(k,t-t^{\prime} ){{\mathbb{P}}}_{H}(k,t^{\prime} ){\varphi }_{H}^{(0)}(t^{\prime} ),$$

(13b)

$$T(k,t)={T}^{(1)}(k,t)+\int\nolimits_{0}^{t} dt^{\prime} H(k,t-t^{\prime} ){{\mathbb{P}}}_{H}(k,t^{\prime} ){\varphi }_{H}(t^{\prime} ),$$

(13c)

$$H(k,t)={H}^{(1)}(k,t)+\int\nolimits_{0}^{t}dt^{\prime} T(k,t-t^{\prime} ){\varphi }_{T}(t^{\prime} ).$$

(13d)

We further need to specify the probability densities

$${P}_{T}^{(0)}(k,t)=(1-p)\int\nolimits_{t}^{\infty }dt^{\prime} \,{\varphi }_{T}(t^{\prime} )(t^{\prime} -t)/{\tau }_{T},$$

(14a)

$${P}_{H}^{(0)}(k,t)=p\ {{\mathbb{P}}}_{H}(k,t)\int\nolimits_{t}^{\infty }dt^{\prime} \,{\varphi }_{H}(t^{\prime} )(t^{\prime} -t)/{\tau }_{H},$$

(14b)

with p = τ_H/(τ_H + τ_T), which account for the fact that the system starts in a stationary state⁴². In particular, the probability to have never hopped before, ${P}_{T}^{(0)}(k,t)$, depends on the probability that the particle is in a trapped state, 1 − p, and on the time integral, which represents the probability that the trapping phase exceeds time t. It can be rationalized as follows: The probability density that the trapping phase is of length $t^{\prime}$ is given by $t^{\prime} {\varphi }_{T}(t^{\prime} )/{\tau }_{T}$. The probability that after lag time t the particle is still trapped is $(t^{\prime} -t){{\Theta }}(t^{\prime} -t)/t^{\prime}$, where Θ( ⋅ ) denotes the Heaviside function. Then the probability that the particle has not yet started a hopping phase at time t is obtained by integrating over all durations $t^{\prime}$: $\int\nolimits_{t}^{\infty }dt^{\prime} \,{\varphi }_{T}(t^{\prime} )(t^{\prime} -t)/{\tau }_{T}$. Similarly, we can derive ${P}_{H}^{(0)}(k,t)$.

Moreover, the probability densities (per time) for the first trapping and hopping event are

$${T}^{(1)}(k,t)=p\ {{\mathbb{P}}}_{H}(k,t)\int\nolimits_{t}^{\infty }dt^{\prime} \,{\varphi }_{H}(t^{\prime} )/{\tau }_{H},$$

(15a)

$${H}^{(1)}(k,t)=(1-p)\int\nolimits_{t}^{\infty }dt^{\prime} \,{\varphi }_{T}(t^{\prime} )/{\tau }_{T}.$$

(15b)

Here, the probability for the first trapping event T⁽¹⁾(k, t) depends on the probability that the particle is in a hopping state p and has hopped for a time t with propagator ${{\mathbb{P}}}_{H}(k,t)$. Further, the probability density for a hopping phase to be of length $t^{\prime}$ is given by $t^{\prime} {\varphi }_{H}(t^{\prime} )/{\tau }_{H}$. The probability for the lag time t = 0 to be in the same interval is uniformly distributed, $1/t^{\prime}$. Then the probability density that the trapping phase starts at t given a hopping phase of length $t^{\prime}$ is obtained by the integral over all possible $t^{\prime}$, $\int\nolimits_{t}^{\infty }dt^{\prime} \ {\varphi }_{H}(t^{\prime} )/{\tau }_{H}$. Similar considerations hold for H⁽¹⁾(k, t).

Finally, we specify the propagator in the hopping phase ${{\mathbb{P}}}_{H}(k,t)$. Since we are interested in terms up to ${{{{{{{\mathcal{O}}}}}}}}({k}^{2})$, we expand it in k: ${{\mathbb{P}}}_{H}(k,t)\simeq 1-{k}^{2}{\langle | {{\Delta }}{{{{{{{\bf{r}}}}}}}}(t){| }^{2}\rangle }_{{{{{\mbox{RR}}}}}}/3!$, where the mean-square displacement of a run-reverse particle is ${\langle | {{\Delta }}{{{{{{{\bf{r}}}}}}}}(t){| }^{2}\rangle }_{{{{{\mbox{RR}}}}}}=2{v}^{2}/{\tilde{\lambda }}^{2}(\exp (-\tilde{\lambda }t)+\tilde{\lambda }t-1)$. The effective rate depends on the turning angle ϑ₀ via $\tilde{\lambda }=\lambda (1-\cos {\vartheta }_{0})$⁹. Subsequently, we perform a Laplace transform, t → s, of equation (13a)–(13d) and use the convolution theorem to derive an analytical solution for the renewal equations in Fourier–Laplace space⁴². The formal solution has been presented elsewhere⁴². We insert the expansion of ${{\mathbb{P}}}_{H}$ and the hop and trapping time distributions, φ_H and φ_T, into the theoretical predictions in Fourier–Laplace space and keep only terms up to ${{{{{{{\mathcal{O}}}}}}}}({k}^{4})$. Using the expansion from the main text, we derive an analytical solution for the mean-square displacement in Laplace space, 〈∣Δr(s)∣²〉. An analytical backtransform to time space is not possible, however, we can extract the long-time (corresponding to s → 0) behavior, see main text.

We further note that the long-time effective diffusivity can also be derived analytically for a truncated power-law distribution, ${\varphi }_{T}(t)=\exp (-\gamma t){(1+t/\tau )}^{-1-\beta }f(\gamma \tau ,\beta )$ with $f(\gamma \tau ,\beta )={\left[{e}^{\gamma \tau }\beta {{{{{{{\rm{ExpIntE}}}}}}}}(1+\beta ,\gamma \tau )\right]}^{-1}$, where ExpIntE( ⋅ , ⋅ ) denotes the generalized exponential integral function⁵¹. It assumes the same form as equation (5) with average trapping time ${\tau }_{T}=\int\nolimits_{0}^{\infty }t{\varphi }_{T}(t)\ {{{{{{{\rm{d}}}}}}}}t=\tau \left[{{{{{{{\rm{ExpIntE}}}}}}}}(\beta ,\gamma \tau )/{{{{{{{\rm{ExpIntE}}}}}}}}(1+\beta ,\gamma \tau )-1\right]$. We note that the average trapping time of the truncated power-law distribution reduces for γ = 0 to that of the power-law distribution used in our manuscript. Details of the distribution become apparent in the short-time behavior of the mean-square displacement, but do not affect our data analysis for the long-time effective diffusivities.