Stenosis triggers spread of helical Pseudomonas biofilms in cylindrical flow systems

Biofilms are multicellular bacterial structures that adhere to surfaces and often endow the bacterial population with tolerance to antibiotics and other environmental insults. Biofilms frequently colonize the tubing of medical devices through mechanisms that are poorly understood. Here we studied the helicoidal spread of Pseudomonas putida biofilms through cylindrical conduits of varied diameters in slow laminar flow regimes. Numerical simulations of such flows reveal vortical motion at stenoses and junctions, which enhances bacterial adhesion and fosters formation of filamentous structures. Formation of long, downstream-flowing bacterial threads that stem from narrowings and connections was detected experimentally, as predicted by our model. Accumulation of bacterial biomass makes the resulting filaments undergo a helical instability. These incipient helices then coarsened until constrained by the tubing walls, and spread along the whole tube length without obstructing the flow. A three-dimensional discrete filament model supports this coarsening mechanism and yields simulations of helix dynamics in accordance with our experimental observations. These findings describe an unanticipated mechanism for bacterial spreading in tubing networks which might be involved in some hospital-acquired infections and bacterial contamination of catheters.


SI text
The mathematical description of an elastic filament as a unidimensional centerline was discussed in References. 62,63 The rod is represented by an adapted framed curve Γ = {γ; t, m 1 , m 2 } (Fig. 9a). γ(s) is the centerline parametrized by the arc length in IR 3 . {t(s), m 1 (s), m 2 (s)} is the orthonormal material frame characterizing the local orientation. t(s) = γ l (s) is tangent to the curve. t l = κ is the curvature (normal) vector. The Bishop frame B(s) = {t, u, v} defines the "rest orientation" at any point of the centerline. It is an adapted frame with zero twist, i.e., uʹ • v = −vʹ • u. Assigning the Bishop frame at s = 0 uniquely determines the remaining values by parallel transport. The twist of the thread at any point is implemented rotating Γ(s) an angle θ with respect to B(s) around the plane defined by {u(s), v(s)} (Fig. 9a): m 1 (θ ) = cos(θ )u + sin(θ )v, m 2 (θ ) = − sin(θ )u + cos(θ )v.
is the unit tangent vector per edge. To construct a Bishop frame we set and define . The frames at other edges are found by parallel transport , . P i are rotation matrices about the curvature binormal satisfying , . If , P i is the identity. The condition must hold during the simulation. It is reestablished by parallel transport in time (instead of space). Let θ i be the angles defining the material frames by rotation of the Bishop frames. The material frame vectors at each edge are:

cos
, sin When the undeformed filament is straight and its elastic behaviour is isotropic, the elastic energy due to twisting and bending is given by the expression: where α and β are the bending and torsion modulus, respectively. ̅ is the length of the segments where is the curvature binormal. The material frame is updated in a quasistatic way, imposing 0, for all segments i not fixed by a boundary condition. 35 The angle configuration minimizing the energy of the rod is determined by this system of equations. When the edges are clamped, we assign the material frame for i = 0 or i = n. Twist at the edges is fixed assigning the values of θ . For stress free ends no boundary condition is assigned. The dynamics of the nodes is governed by Newton's second law . (2) where f represents the external forces and the elastic forces, that can be evaluated once the angles are known. 35 M is the mass matrix, that is set equal to a multiple of the identity M=mI for an isotropic thread. This system is integrated using a Verlet solver combined with a manifold projection method to enforce the inextensibility constraint for each segment, see Ref. 64 The force exerted by the fluid on the filaments can be evaluated following Ref. 65 Where a node hits the tube, a virtual force is applied to impose the spatial constraint generated by the presence of the wall (penalty method), keeping the node inside. Addition of biomass may be represented by increase of segment length. These equations can be coupled to the motion of bodies

3/11
placed at the edges of the thread following. 35 Cells coupled to the filament are expected to be trapped rotating in the vortices formed at narrowings. 48-50 This is a source of twist, that can be introduced in the simulation following different procedures of increasing complexity: fixing a time dependent angle at the edge, coupling to spinning particles as in references. 66,67 The fluid force is modeled using Reference. 65 For numerical purposes we nondimesionalize the equations (2). Choosing λ=1mm and T=1s as reference lengths and times, the change of variable x = λ x' , t = T t' yields: (3) Revising the definition of the energy E (1), this change brings about the controlling parameters 2 3 and 2 3 . We are working with one dimensional filaments, therefore neglecting the cross-sections. A real cylindrical thread of density , small radius r, and length L >> r, would be approximated in this framework by a discrete filament with N nodes and N-1 edges, with mass m =   r 2 L /(N-1). Knowing ranges of values of ', ' that lead to different types of numerical helical structures we might guess ranges for ,  provided that the density  and radius r of the thread are known.
We have performed series of tests in tubes with diameters 0.2, 0.5, 1, 2 mm and lengths 10, 40, 80 mm, and explored parameter ranges for which different helical dynamics are observed. Fig. 6 in the main text and Supplementary figure S4 and movie S2 and S3 are generated setting α' = 1.345, β' = 0.789, N = 100. The tube diameter and length are 2 mm and 80 mm, respectively. Still images of Fig. 6 and S4 show snapshots of the simulations for 3200, 6300, 9600 and 16200 computational iterations. Simulations start from an originally straight thread in which the positions of the nodes have been randomly perturbed, resulting in a small length increase compared to the tube length (an excess around 0.49% of the tube length). The left edge of the thread is subject to smoothly increasing twist, reaching 7 × π at the end. We alternate steps in which we solve (3) with steps in which we increase the length of the edges of the thread (edge length increase rate = 2·10 -5 mm/iteration) and reset the reference edges. The helix pitch reduces, as shown by Supplementary Fig. S4 and Supplementary movie S3. The final helix pitch is constrained by the excess length to one helical loop per 5 mm. Adding together the times during which equation (3) is numerically solved, and the time the biofilm thread needs to increase its length (due to biomass production, or other factors) by the specified amount, we would find the total evolution time. However, we lack experimental values for this later parameter.
When we couple the dynamics of the filament to the flow, the whole structure moves downstream very slowly as the helix develops, in a similar way to Supplementary movie S1. The twist is not really necessary for the thread to evolve into a helix. In tubes of smaller diameter, threads naturally evolve to helices wrapping around the walls provided the initial excess length allows it. However, the presence of twist improves the stability of the evolution as the length of the thread increases, avoiding knots and messy structures.

4/11
Supplementary figures Figure S1. Several experimental setups were designed (a,c,e,g,i) and tested (b,d,f,h,j) to study helix formation. Figure S2. The helical pattern adapts its geometry to the tube radius where it develops. Helical pitches in tubes of 2 mm inner diameter (a) are usually larger than those in 1 mm tubes (b). Brightness and contrast were adjusted to enhance the image. Figure S3. The biofilm expands to fill the length of the tubes they colonize. Colonized zones glow in blue as the result of an autofluorescence phenomenon typical of Pseudomonas putida (see reference 1 provided at the end of supplementary material) when exposed to light at 312 nm. Dark zones in the tubes contain air bubbles. Brightness and contrast were adjusted to enhance the image. Figure S4. Numerical simulation of helix dynamics using our discrete filament model. As more biomass accrues (a-d), the pitch of a helix diminishes. The length of a complete turn of helix of pitch h and circumference 2πr is (h 2 + 4 π 2 r 2 ) 1/2 . A thread of length Lh wraps around a tube of radius r and length Lt forming k loops of pitch Lt/k if Lh 2 ~Lt 2 +4 π 2 r 2 k 2 . As Lt grows, the helix pitch decreases according to this rule. In our dynamical simulation of filament growth, an imposed twist facilitates a stable evolution avoiding knots and messy structures. The helix pitch for the last image is 5mm. Tube showing a helix formed by this coiling effect into an unperturbed tube. The inoculated strain was P. putida mt-2 and the imposed flow rate was 0.15 ml/min. The photograph was taken at the position of the assembly depicted by the yellow circle in (a). Brightness and contrast were adjusted to enhance the image. Figure S7. Poiseuille flow profile evolution. The flow evolves from the entrance of the channel to adopt the final laminar steady state profile. The laminar steady flow of an incompressible viscous fluid driven by a pressure drop Δp established between two ends of a straight tube of uniform circular cross-section is described by the final parabolic Poiseuille profile. Writing down the Navier-Stokes equations in cylindrical coordinates (r,Θ,z), the velocity of the fluid U=(u r ,u Θ ,u z ) is given by u r =u Θ =0 and u z = Δp (R 2 -r 2 )/4μ, where R is the tube radius and μ the fluid viscosity. Additional information can be found in reference 2 provided at the end of supplementary material.

11/11
Movie S1 (Download MovieS1_lowq.avi from supplementary material). Accelerated movie (1 second in the movie corresponds to 1.75 hours of experimental time) showing that for a fixed position in the circuit, the biomass thread changes its geometry slowly with time, in a similar way to numerically simulated threads. Note that flow velocities near the walls are much slower than fluid average velocities, already low (about 0.8 mm/s). For high resolution video, visit: http://www.mat.ucm.es/~acarpio/helix/MovieS1.avi