The experimental system we are mimicking are bacterial cells confined in a box with the turbine at its center (see Figure 1B). The cells can move over the substrate via pili binding and unbinding and additionally interact with the turbine by pili attachment and excluded volume effects.

While Neisseria gonorrhoeae cells typically have a diplococcus shape [27], we simplify the in silico cell shape of bacteria as two-dimensional circular disks with a radius R. Please note, however, that a diplococcus shape can be also considered [19, 23–25]. Each cell possesses exactly N _p pili that are homogeneously distributed on the cell outline (see Figure 1A). We approximate pili as straight lines connecting two points: their start point (also called anchor point), located at the circular cell surface, and their end point. The distance between these two points is called the contour length l _c of the pilus.

The turbine is described by N _turb × 2 disks (with radius R) arranged in a double row array (see Figure 1B). Since the only way how the turbine can move is by rotation, the relative locations of the turbine disks towards each other are fixed. The orientation of the turbine is described by the angle α with respect to the y-axis and the turbine center is located at the position r ^(turb) (see Figure 1B).

Initially, a pilus protrudes with the velocity v _p in a direction perpendicular to the cell surface. When a pilus reaches a specific length, drawn from an exponential length distribution with mean length l _p [17, 18], it starts to retract with a velocity v r ( 0 ) . This retraction continues until the pilus has a contour length l _c = 0. In that case, the pilus is removed and immediately, a new pilus protrudes from the cell from the same position.

A pilus binds stochastically to the substrate or the turbine disks, independently whether it is protruding or retracting. In both cases, a pilus can only bind with its tip and it can bind only to either the substrate or the turbine disks. After binding, a pilus immediately starts to retract [28]. The binding is modeled as a Poisson process with the binding rate κ a t t ( s u b ) to the substrate and κ a t t ( t u r b ) to any of the turbine disks. Since the model is two-dimensional, we ignore that, for pili to bind to a substrate, they first need to be long enough to reach the substrate with their tips. In 3D, this leads to a delay in the initial binding event of a newly protruding pilus, but is less relevant for following attachments of pili. Hence, initially, the binding rate might be smaller to allow a pilus to protrude to reach the substrate. To account for this process, pili that bind to the substrate the first time do so with a rate κ a t t ( s u b , 0 ) = v p / l p . Again, the binding is modelled as a Poisson process. The rate κ a t t ( s u b , 0 ) corresponds to the binding of a pilus that protrude with the velocity v _p until it reaches a length determined by an exponential distribution with mean length l _p, in which case it binds. On average, this takes the time 1 / κ a t t ( s u b , 0 ) .

An attached pilus will generate a pulling force that acts on the cell and, if attached to a turbine disk, also on the turbine. Each pilus is modeled as a Hookean spring with the spring constant k _p. After attachment, a pilus is stretched due to its retraction and hence, mediates a pulling force. This force is proportional to the difference between the contour length l _c and the length of the pilus that it would have if it was not attached, here called the free length l _f. Here, we consider the case where a pilus can only generate a pulling force and no pushing force, thus the force is zero if l _c ≤ l _f.

Experimentally, it has been shown that the pilus force affects the retraction velocity by v r = m a x 0 , v r ( 0 ) 1 − F F s , (1)

with the stalling force F _s [22]. The force of a pilus also affects the unbinding from the substrate. Pilus detachment is modelled by a Poisson process with the rate κ d e t ( s u b ) F = 1 τ d e t ( s u b , 1 ) ⁡ exp − F F d e t ( s u b , 1 ) + τ d e t ( s u b , 2 ) ⁡ exp − F F d e t ( s u b , 2 ) , (2)

for detachment from the substrate, as motivated by [18]. Here, we introduce the detachment times τ d e t ( s u b , 1 ) and τ d e t ( s u b , 2 ) and the detachment forces F d e t ( s u b , 1 ) and F d e t ( s u b , 2 ) . For detachment from a turbine disk, we describe the rate by a simpler relation κ d e t ( t u r b ) = 1 τ d e t ( t u r b ) exp F F d e t ( t u r b ) , (3)

with the detachment time τ d e t ( t u r b ) and detachment force F d e t ( t u r b ) . We picked a different expression for the pili-turbine detachment (in comparison to the pili-substrate detachment) because in general, the turbine does not necessarily have to consist of the same material as the substrate, which determines pili binding properties [16, 17]. Below, we will explore turbine rotation in dependence on pili-turbine interactions, thus favouring a simpler single exponential form.

We model pili as Hookean springs with the spring constant k _p. For a cell i at location r i ( c e l l ) , an attached pilus j causes a force F i ( p i l u s j ) at the pilus anchor point r i ( a n c h o r j ) on the surface of the cell in the direction of the point r i ( a t t j ) where the pilus tip is attached (see Figure 1C).

Next to the active pili force, passive excluded volume forces are acting on the cells. Cells are located in a two-dimensional box with size L × L and a cell i that overlaps with the boundary wall is exposed to a repulsive force F i ( w a l l ) (see Figure 1C). This force is modelled as a harmonic interaction with spring constant k ^(wall) and acts on the cell center in the normal direction of the boundary wall if the overlap is smaller than the cell radius R. Additionally, intersections of 2 cells i and k lead to a repulsive force F i k ( c − c ) of the centers of both cells with the spring constant k ^(c−c) (see Figure 1C). A similar type of repulsive force F i ( t u r b m ) is acting between the cells and the turbine disks m with the spring constant k ^(turb) (see Figure 1C).

The total force of the cell i is given by F i ( t o t ) = ∑ attached pili j F i ( p i l u s j ) + F i ( w a l l ) + ∑ cells k F i k ( c − c ) + ∑ turbine disks m F i ( t u r b m ) . (4)

Additionally, the total torque acting on a cell is T i ( t o t ) = ∑ attached pili j r i ( a n c h o r j ) − r i ( c e l l ) × F i ( p i l u s j ) (5)

In the overdamped limit [29], a force mediates a translational motion of the cell with the velocity d r i ( c e l l ) d t = μ t F i ( t o t ) (6)

and the torque leads to a rotation with the angular velocity vector ω i ( c e l l ) = μ r T i ( t o t ) (7)

Here, we introduce the translational mobility μ _t and the rotational mobility μ _r. The same forces and torques cause an equivalent displacement of the pili anchor points r i ( a n c h o r j ) and attachment points r i ( a t t j ) in case of attachment to a turbine disk. The mobilities introduced here are a result of friction with the viscous solvent and the substrate on which the cells are moving.

The turbine is only able to undergo rotational motion. The total torque acting on the turbine is given by T ( t u r b ) = ∑ attached pili j r i ( a t t j ) − r ( t u r b ) × F i ( p i l u s j ) + ∑ cells i ∑ turbine disks m r i ( c e l l ) − r ( t u r b ) × − F i ( t u r b m ) . (8)

The rotation of the turbine is modelled in the overdamped limit with a mobility μ _turb. In that case, the turbine angular velocity vector is ω ( t u r b ) = μ t u r b T ( t u r b ) . (9)

The simulations were performed on the local computing cluster of the Max Planck Institute for the Physics of Complex Systems (Dresden, Germany), consisting of x86-64 GNU/Linux systems. The code was written in C++. We use an Euler algorithm to solve the equations of motion with a time step δt. While this is one of the most simple numerical schemes to solve the equations of our computational model, it often leads to numerical errors and instabilities when modelling molecular dynamics systems for long times [30]. We do not expect that this is a problem in our system due to the stochastic nature of the pili binding and unbinding, which basically represents our system as a series of many short time events, continuously interrupted by rearrangements in the pili network. Hence, we do not expect any differences in the qualitative outcome of the simulations.

If not stated otherwise, we use the parameters provided in Table 1. Most parameters we use are based on previous studies. The excluded volume spring constants k ^(c−c), k ^(turb) and k ^(wall) have no effect on the simulation outcome as long as they are chosen large enough to be able to compete with the pili forces. For the remaining parameters, e.g., the turbine mobility, we do not expect a qualitative difference in the results of the simulation.

We initialize the simulation by randomly distributing the cells in the box and only analyse the turbine rotation after an initialization period of 1,000 s. To calculate the angular velocity of the turbine rotation, we compute ω α t = α t + Δ t − α t Δ t , (10)

with Δt = 0.5 s and the turbine orientation α (see Figure 1B).

To investigate how the adhesion of cells to a turbine affects the turbine rotation, we first simulate the computational model for three different cases: 1) pili detach with detachment times of τ d e t ( t u r b ) = 50 s (see Figure 2A and Supplementary Movie S1), motivated by binding times in the order of a few minutes inferred from cell trajectories on plastic surfaces [17], 2) pili detach with τ d e t ( t u r b ) = 2 s (see Supplementary Movie S2), motivated by considerably smaller binding times experimentally measured for BSA coated beads [18] and 3) pili not bind to the turbine at all (see Supplementary Movie S3) and the only way how cells interact with the turbine is by excluded volume forces. Here, we only vary the turbine detachment time τ d e t ( t u r b ) and not the detachment force F d e t ( t u r b ) . We expect that there will be no qualitative difference in the turbine rotation for both cases, an increase (decrease) of τ d e t ( t u r b ) and a increase (decrease) of F d e t ( t u r b ) , pili will bind more (less) strongly to the turbine. For simplicity, we assume that F d e t ( t u r b ) ≔ F s with the pilus stalling force F _s, see Eq. 1. In Figure 2B, we show the dynamics of the turbine angle α for the three cases. We find that for the largest detachment time, τ d e t ( t u r b ) = 50 s , the rotation is the strongest. If pili are not permitted to attach to the turbine, the rotation appears to be the weakest. Additionally, the distribution of velocities shows that if pili cannot bind to the turbine, the turbine will move with smaller angular velocities ω, while the distributions seem to not depend on the detachment rate τ d e t ( t u r b ) , as long as pili are capable of binding to the turbine. We also find that the turbine does not undergo extended phases in which it is not moving, but instead seems to be continuously rotating (see Supplementary Figure S1).

To quantify how strong the rotation is, we compute the angular mean squared displacement (angular MSD), given by δ α Δ t = ⟨ α t + Δ t − α t 2 ⟩ t , (11)

(see Figure 2C). Indeed, the angular MSD is the highest for τ d e t ( t u r b ) = 50 s and the lowest for the case where pili do not attach to turbine disks. For τ d e t ( t u r b ) = 2 s , the angular MSD is identical to the case τ d e t ( t u r b ) = 50 s for small time differences (Δt). In both cases, τ d e t ( t u r b ) = 2 s and non-attaching pili, the angular MSD follows a linear scaling, δ _α ∝Δt, corresponding to a diffusive regime. For τ d e t ( t u r b ) = 50 s , we find a superdiffusive regime for intermediate time differences Δt ≈ 10, −,10³ s. For particles undergoing translational motion, such superdiffusive regimes usually emerge from ballistic motion. For the rotatable turbine, this corresponds to a regime where the turbine moves persistently in a random direction. Indeed, such a behaviour is observed in Figure 2B, where trajectories of the angle α move persistently in one direction over periods of hundreds to thousands of seconds before the rotation turns towards the opposite direction. For larger time differences, the angular MSD becomes again diffusive, implying that while the rotation is persistent over a certain time scale, for larger time intervals it becomes random again.

To learn more about how persistent the turbine rotation is, we compute the angular velocity autocorrelation function (angular velocity ACF), given by ν Δ t = ⟨ ω α t + Δ t ω α t ⟩ t . (12)

Independently of whether pili can bind to the turbine or not, we find that the correlation function ν(Δt) is decaying with time (see Figure 2D). This confirms our previous observation for the angular MSD where we found that for very large time differences the turbine rotation is diffusive and no longer shows signs of persistence. We also find that for τ d e t ( t u r b ) = 50 s , the normalised angular velocity ACF ν(Δt)/ν(0) has the slowest decay. This behaviour correlates with the persistent rotation for time differences between Δt ≈ 10, −,10³. Surprisingly, we find that the decay of the angular velocity ACF is faster for τ d e t ( t u r b ) = 2 s than for the case where pili do not bind to the turbine. We will provide an explanation for this behaviour in Section 3.3. Before doing that, we will have a closer look at how the turbine rotation depends on the cell adhesion to the turbine disks and the number of cells in the system.

We begin with a systematic analysis of how the angular MSD δ _α (Δt) depends on the detachment time of pili from the turbine disks τ d e t ( t u r b ) , shown in Figure 3A. For small time differences Δt < 10 s, we find that δ _α is independent of τ d e t ( t u r b ) and exhibits a linear scaling, δ _α ∝Δt, corresponding to a diffusive regime. For high enough values of the detachment time τ d e t ( t u r b ) and Δt > 10 s, the angular MSD becomes superdiffusive, δ _α ∝Δt ². The duration of this superdiffusive regime increases with increasing τ d e t ( t u r b ) and is barely observable for detachment times τ d e t ( t u r b ) < 10 s . For even longer time differences Δt, the MSD becomes diffusive again, δ _α ∝Δt.

Next, we investigate the normalised angular velocity ACF ν(Δt)/ν(0) (see Figure 3B) and find that for Δt < 2 s, it rapidly decreases. For Δt ≥ 2 s, it exponentially decreases with the Δt. By fitting a function of the form ν ( Δ t ) ∝ exp − Δ t / τ to this later regime, we identify the characteristic time τ of the exponential decay of ν(Δt) and can investigate its dependence on τ d e t ( t u r b ) . We find that it is increasing with increasing detachment times τ d e t ( t u r b ) . This is in accordance with the increasing duration of the superdiffusive regime of the angular MSD. The stronger pili of a cell bind to the turbine, the more persistently the turbine rotates in one direction before it inverts its direction.

Additionally, we investigated how the rotation of the turbine is affected by the density ρ of cells in the system. We find that with increasing density, the angular MSD δ _α becomes smaller (see Figure 3C). While we observe the diffusive regime of δ _α (Δt) for Δt < 10 s, a pronounced superdiffusive regime for larger Δt vanishes if the cell density becomes too large. For the normalised angular velocity ACF ν(Δt)/ν(0) (see Figure 3D), the characteristic decay time τ is initially increasing with cell density ρ and then decreases when ρ increases further.

In the following section, we provide qualitative arguments that explain the observed behaviours.

To understand how the persistent rotation of the turbine in an otherwise symmetric system can emerge, we first consider how strongly cells bind to the turbine disks with the help of their pili and how pili binding affects substrate attachment of the remaining pili. Our hypothesis is the following: if cells stay attached to the same position on the turbine surface for an extended time, which is considerably larger than average attachment times of individual pili, the cells will continuously pull the turbine in the same direction. This mechanism is dramatically different to the previously reported rotary microdevices driven by swimming bacteria [6, 7], where cells collide with the microdevice and hence, push it. Instead, in our simulations, cells pull on the turbine. In Figure 4A, we show sketch of a cell which is attached to the turbine disk with some of its pili, while the other pili bind or unbind from the substrate. Only pili that protrude towards the turbine disks can attach, while pili protruding away from the turbine can only attach to the substrate. Thus, cells preferentially bind to the substrate in the direction away from the turbine. Since the cell is also attached with some of its pili to the turbine disks, it is pulling the turbine along the direction pointing from the turbine towards the cell. Thus, the cell pulls on the turbine, leading to a rotation of the turbine in the direction of the attached cell. This process is enhanced by the cooperation of multiple pili of a cell. To clarify this, we investigate the durations of cells attaching to the turbine disk before they detach again. To this aim, we solve a stochastic model of pili binding and unbinding, where a cell, possessing in total N ₀ pili that are all growing in the direction of the turbine, has n pili attached to the turbine. With a rate k _att (N ₀ − n) a non-attached pilus attaches to the turbine, while with the rate k _det n a pilus detaches from the turbine (see Figure 4B). By describing the attachment and detachment as Poisson processes, we can numerically solve this system with the help of a Gillespie algorithm and investigate the mean time for a complete unbinding of all pili of a cell from the turbine (see Figure 4C). We find that the mean unbinding time is increasing rapidly with the pili detachment time τ off = exp ( F s / F d e t ( t u r b ) ) / k det , reaching around 10 s for τ _off = 2 s and 10⁷ s for τ _off = 50 s. For simplicity, we assume that all pili pull with their stalling force and that the detachment force is identical to the stalling force, F s = F d e t ( t u r b ) . The duration of how long cells bind to the turbine corresponds to the duration of superdiffusive behaviour in Figure 3A and also exhibits the same qualitative behaviour as the decay time of the angular velocity autocorrelation (see Figure 3B). We do not expect a perfect quantitative agreement between the two time scales (Figure 3B and Figure 4C) as the simplified model of pili binding and unbinding of a cell to the turbine provided here ignores that cells can indeed move even if its pili are bound to the turbine, e.g., parallel to the turbine surface. Additionally, when multiple cells bind to the same region of the turbine, they will interact via excluded volume effects. This will lead to additional forces acting on the involved pili and thus, might enhance their detachment.

We can now provide an explanation of the three time regimes observed in the angular MSD in Figure 2C and angular velocity ACF in Figure 2D. For short times, Δt < 10 s, individual pili stochastically bind and unbind to the substrate, pulling on the turbine, and additionally, cells randomly collide with the turbine. This causes small displacements of the turbine angle α in one direction, also explaining the very sudden drop of the angular velocity ACF for short times since the life time of such displacements is very small. This leads to random fluctuations of the turbine, and as a result, to a diffusive scaling. For larger times, 10 s < Δt < 10³ s, asymmetries in the distribution of cells on the surface of the turbine lead to a unidirectional rotation of the turbine and as a result, to a superdiffusive regime. The persistence in the turbine rotation also explains the positive angular velocity ACF. The duration of this regime depends on how strong pili bind to the turbine disk (see Figure 4C). Finally, for Δt > 10³ s, the distribution of cells on the turbine surface is re-arranged up to such a degree that the turbine can change its direction of motion. This will then lead to a diffusive regime again. This time also corresponds to the characteristic time of the angular velocity ACF decay, confirming that for large enough times, the turbine forgets its initial direction of motion.

For cases where pili bind stronger to the turbine, we expect to not only have more pili to bind to the turbine, but also more cells to be bound to the turbine with its pili. To this aim, we check the number of cells and pili bound to the turbine as a function of time (see Supplementary Figure S2). We find that after an initial phase of pili attachment, for larger pili-turbine detachment times τ d e t ( t u r b ) more cells and more pili are bound to the turbine on average. This is in agreement with our previously introduced hypothesis that the turbine exhibits persistent rotational motion due to cells that are bound to the turbine surface.

Next, we consider the origin of the cell density dependent rotation of the turbine, see Figure 3C and Figure 3D. If the cell density is very small, often no cell will bind to the turbine most of the time and thus, the angular MSD is initially increasing with cell density. For moderate cell densities the difference between cells attaching to different sides of the turbine (see Figure 4D) will be significant and due to the random asymmetry, the turbine will rotate persistently in one direction until cells randomly unbind from the turbine. If the cell density gets too large, more and more cells are pulling in the opposite direction of the turbine rotation, reducing its persistence.

There are additional processes that can have a significant impact on the turbine rotation: 1) Due to the rotation of the turbine, it will constantly “pick up” cells it collides with. For the traditional rotary microdevices that are driven by swimming bacteria that collide with the device and push it [6–8], this would lead to a torque that acts in the opposite direction of the original rotation. When cells pull on the turbine instead, the persistent rotation is being enhanced since asymmetric distribution of cells bound to the turbine disks is getting even stronger (see Figure 4E). 2) In Section 3.1, we found that the angular velocity ACF seems to decay faster for the case with pili attachment and τ d e t ( t u r b ) = 2 s , than for the case without pili attachment (see Figure 2D). A possible explanation for this is that even without pili-turbine attachment, cells can still rotate the turbine due to collisions and resulting excluded volume forces. That way, we reproduce a system where cells push the turbine, instead of pulling on it. This process possesses its own characteristic time scales that are linked to the substrate motility of the cells.

To summarise, the persistent rotation of the turbine originates from persistent attachment of cells to the turbine over times scales much longer than the characteristic detachment times of individual pili. While the resulting rotation is indeed persistent and, depending on the binding properties, can go on over multiple revolutions, it is not permanently unidirectional. There is no asymmetry in the initial rotational direction of the turbine and even though the turbine can move in one direction over an extended time, this rotation will at some point reverse direction. Next, we will provide an example of how cells using type IV pili can drive a permanent unidirectional rotation of a turbine.

In order to produce a permanent unidirectional rotation of the turbine, an asymmetry of the turbine is required. Typically, such assymetries are created by altering the turbine geometry [7, 8, 12], but here we have a chance to exploit asymmetries in the adhesion of cells to the turbine instead. In Figure 5A, Figure 5B and Supplementary Movie S4, we provide an example of how altered pili binding properties on one side of the turbine wings can lead to a unidirectional rotation. Here, we consider different cases: 1) pili bind stronger to the manipulated turbine disks, corresponding to a larger detachment time τ d e t ( t u r b ) or 2) pili bind weaker to the manipulated turbine disks, with a lower value of τ d e t ( t u r b ) . Depending on how pili bind to these regions, the direction of the turbine is affected as, on average, cells are more or less strongly adhered to the manipulated region and thus, more or less cells are bound to these sides of the turbine and mediate the turbine rotation. In case (i), the turbine rotates in the direction of manipulated region because more cells are attached there, for 2) less cells are attached in the manipulated region and the turbine rotates in the opposite direction.

We also estimate the average angular velocity ω of the turbine and find values around 2–4°/s. For the used rotational mobility of the turbine μ _turb (see Table 1), this corresponds to a turbine torque around 10 pN μm, comparable to previously published values for rotary devices in bacterial baths of hundreds of swimming cells [7]. Here, however, the system consists of 40 cells only. Note that this is only the lower limit of the turbine torque and the same turbine might be able to generate larger torques if it experiences an opposing torque. In the absence of a counteracting torque, the rotation speed of the turbine is limited by the retraction velocity of pili, around 2 μm/s [22].

We see evidence that the rotation is permanently unidirectional in the angular MSD (see Figure 5C) which is superdiffusive (ballistic) for arbitrary time differences Δt, equivalent to a rotation in one direction. Additionally, the normalised angular velocity ACF (see Figure 5D) is no longer decaying with time Δt, suggesting a constant rotation in the same direction.

In experiments, such manipulated regions on the surface of the turbine could be generated by coating it with a chemical such as BSA, which has been shown to alter the detachment time of the pili [17, 18, 32].