Multiscale modeling of bacterial colonies: how pili mediate the dynamics of single cells and cellular aggregates

Neisseria gonorrhoeae is the causative agent of one of the most common sexually transmitted diseases, gonorrhea. Over the past two decades there has been an alarming increase of reported gonorrhea cases where the bacteria were resistant to the most commonly used antibiotics thus prompting for alternative antimicrobial treatment strategies. The crucial step in this and many other bacterial infections is the formation of microcolonies, agglomerates consisting of up to several thousands of cells. The attachment and motility of cells on solid substrates as well as the cell–cell interactions are primarily mediated by type IV pili, long polymeric filaments protruding from the surface of cells. While the crucial role of pili in the assembly of microcolonies has been well recognized, the exact mechanisms of how they govern the formation and dynamics of microcolonies are still poorly understood. Here, we present a computational model of individual cells with explicit pili dynamics, force generation and pili–pili interactions. We employ the model to study a wide range of biological processes, such as the motility of individual cells on a surface, the heterogeneous cell motility within the large cell aggregates, and the merging dynamics and the self-assembly of microcolonies. The results of numerical simulations highlight the central role of pili generated forces in the formation of bacterial colonies and are in agreement with the available experimental observations. The model can quantify the behavior of multicellular bacterial colonies on biologically relevant temporal and spatial scales and can be easily adjusted to include the geometry and pili characteristics of various bacterial species. Ultimately, the combination of the microbiological experimental approach with the in silico model of bacterial colonies might provide new qualitative and quantitative insights on the development of bacterial infections and thus pave the way to new antimicrobial treatments.

) of cells within a colony. Here we use the orientation of a cell which is defined its long axis and compute the nematic order relative to the radial vector pointing towards the cell position. The case S = 0 indicates a random distribution of orientations, while S < 0 corresponds to a tangential orientation of the cell relative to the colony surface. This effect is more pronounced for the strong parameter set, characterized by a larger detachment force F d,pil . (B,C) Mean of the normal and tangential net forces acting on a cell relative to the surface of the colony. The attractive pili-pili interactions and the repulsive excluded volume interactions balance each other such that the mean net force acting on a cell vanishes. The weak (blue) and strong (red) parameter sets are defined in table 3.  Figure S4. Dependence of the colony radius on the cell number (blue dots) for a colony without substrate interactions. Here the colony radius corresponds to the radius of a circle having the same area as the envelope of the two-dimensional projection of a colony. The green curve is R = 0.7 · N 1/3 µm. Thus the volume is proportional to the cell number.  Figure S5. Internal properties of individual colonies for F stall = 180 pN, F d,pil = 50 pN and τ d,pil = 50 s. Vertical dashed lines represent the corresponding radii of the colonies R col by fitting the density profile of cells within the microcolony according to equation (10). For such a small detachment force the cells are highly motile. (A) Mean-Squared displacements of cells as a function of their distance d com from the center of the colony. We see that cells do not stay longer than 30-60 seconds in their shell, thus it is not possible to estimate a long time limit diffusion coefficient.   (12)) of cells within a colony. Here we use the orientation of a cell which is defined its long axis and compute the nematic order relative to the radial vector pointing towards the cell position. The case S = 0 indicates a random distribution of orientations, while S < 0 corresponds to a tangential orientation of the cell relative to the colony surface.

Geometry of an individual cell and free pili dynamics
A pilus k is characterized by two points: its start x k . The contour length of a pilus is defined to be l k . New pili are produced with a rate λ p until a cell has a maximal number of pili, N p . In our simulations we do not check for the production of new pili for every time step, but only every ∆t d = 0.02 s. This simplification allows to reduce the numerical cost of our simulations and is valid as long as the time corresponding to the attachment rates is much larger. The start point of the pilus k is randomly distributed on the surface of one of the two cocci (r i ) of the cell i, with the constraints that the start point cannot be inside of any cocci, ensured by the condition r A new pilus k protrudes from the surface of coccus of cell i with a velocity v pro . For the position of the pilus end point it follows where l (fr) k is the free length of a pilus and j = a, b, giving the position of the coccus from which the pilus emerges. The free length is defined as the length of the pilus if there is no force acting on the polymer. The change of the free length of a protruding pilus k is given by The pili switch to the retraction state with a probability that corresponds to a rate λ ret . Each pilus retracts with the velocity v ret , which is related to the free pili length by l (fr) If the end point of a perpendicularly protruding or retracting pilus is inside of the substrate (see figure 1, main text), the end point will slide along the substrate. While the contour length of the pilus keeps its value, the pilus elongates towards the same direction in the x − y−plane as before, but with its z-component positioned exactly at the surface. If the length l (fr) k of a free pilus (in this case the free length has the same value as the contour length) is smaller than 0 the pilus is removed. Furthermore a pilus is not able to switch back to the protrusion state after retraction started.

Attachment to a substrate or other pili
The binding to the substrate is controlled by a stochastic process characterized by the rate λ sub . A pilus will always bind with its end point to a point x at the surface. The binding of two non-attached pili is described in a similar manner, however, in this case it is necessary to describe the binding of two pili in three-dimensional space. We account for a larger binding probability due to thermal fluctuations of a pilus by using the beam-equation for a semi-flexible rod [2]. Solving the beam equation for pilus k gives a mean diameter d beam of a cone as a function of the distance l from the start point x with l p denoting the persistence length. In our model two pili bind stochastically with the rate λ pil if one free pilus, described as the cone (see equation (S4)), and a neighboring free pilus, described as a finite line, intersect. An intersection point x . If this length is smaller than the free length of the pilus l k , binding will take place. Otherwise the binding will not take place because the pilus would need to be longer than its actual length. In case of binding the length of the tail l needs to be saved and will be added to the pilus after detachment. The new contour and free length l

Pilus forces
The pilus is modeled as a Hookean spring with a spring constant k pull [3], which allows to compute the pulling force acting on the cell. The force of a pilus k attached to the substrate is given by Similarly the pulling force of a pilus k attached to another pilus j is given by During our simulations we set the end point of an attached pilus to the start point of the associated pilus in order to easily compute the direction vector x (se) k , thus the contour length of the single pilus l (con) k is now the contour length of the bundle. Since the direction of pili forces between two bound pili only depend on x (se) k , it is sufficient to compute the distance of the start points of the two coupled pili and compare it to the sum of their free lengths. While the contour length solely depends on the motion of the cells and the position of the pili start points, the free length of pilus k is changed due to its retraction. The retraction velocity depends on the pulling force: Here F is either the absolute value of a the pili-pili forces F = F (pp) k , or the absolute value of a force resulting from an attachment to the substrate F = F (ps) k . F stall is the stalling force and determines the characteristic pulling force of a pilus, because the retraction of pili is stalled for larger forces F > F stall [4]. The characteristic pulling force enables us to neglect the length-dependence of the spring constant of the pili. By using a very large spring constant k pull (relative to the F stall ), this force governs the dynamics of the system and not the spring constant itself. While pili that fulfill the condition l (con) k > F stall /k pull are able to generate forces in the order of the characteristic force F stall , the pulling force of shorter pili is no longer able to exceed. While a high spring constant guarantees that the resulting pili need to be extremely short, we set the pulling force of a short pilus equal to the stalling force. This is motivated by the fact that the time to reach the stalling force for a pilus is proportional to its length. The force-dependent binding of pili, similarly to the creation of new pili or their binding, is not checked for every time step, but only after ∆t d = 0.02 s.

Cell forces and motility
An intersection of two cocci of two different cells causes a repulsive force. For the sake of simplicity we call the first coccus i (position vector r (a) i ) and the second coccus j (position vector r with r ij = r i , j = a, b) intersects with the substrate and overlaps by a distance ∆d ov , a repulsive force is acting on the coccus: Here, k cs denotes the excluded volume spring constant and e z is the unit vector pointing perpendicular to the substrate. The force does not act on the surface point, but at the point r i · e z of the coccus. This is realized by a very high spring constant of the surface and a small overlap. In this case, the torque will only be weakly affected.

Simulation details and parameters
Our simulations were performed on the local computing cluster consisting of x86-64 GNU/Linux systems of the MPI-PKS. All machines possess Intel Xeon processors with a clock rate of 2.2 to 3.0 GHz and have between 2 to 4 CPUs. The code was written in C++ and parallelized on CPU by using the library OpenMP. We used the GCC-compiler (version 4.8.1) and were running the simulations on 8 cores in parallel.

S3.1. Velocity autocorrelation of single cells on a surface
We fitted a double-exponential function of the form to the velocity autocorrelation data. The results of the fitting for strong and weak pili-substrate interactions are given in table S1.

S3.2. Profile of diffusion coefficient of cells inside of colony
We fitted an exponential function of the form to the diffusion coefficient D of cells as a function of their distance from the center of mass d com of the colony. The results of the fitting for strong and weak pili-substrate interactions are given in table S2.

S4. Simulation Details
If not mentioned otherwise we used the parameter set, given in table S4.

S4.1. Single Cells
We place an individual cell on top of substrate with a distance similar to the cell radius and random cell orientation. For the first 2 seconds of the simulations the cells only interacted via the repulsive excluded volume interactions to remove overlaps of the randomly distributed cells. After this time we turned on the production and dynamics of pili. To reduce the impact of the initial condition on our results we simulated a time span of 30 minutes and analyzed the last 20 minutes. For every parameter set we analyzed the trajectories of 250 cells and sampled over the parameter set given in table S5.

S4.2. Single Colony on Substrate
We initialized a colony by randomly distributing N cells (with random orientations) in a sphere with radius 0.85 · N 1/3 µm at a position far enough from the substrate so that every cell will have a distance of at least 2.5 µm from the substrate, which is larger  than the characteristic pili length but not too large such that the pili would not be able to attach to the substrate. We allowed initial overlap of the cells and only activated excluded volume forces for the first 2 seconds of the simulation. Afterwards we activated the pili and simulated for 20 minutes, but only analyzed the last 10 minutes. For every parameter set (given in table S6) we analyzed the trajectories of 10 colonies.

S4.3. Free Single Colony
We initialized a colony by randomly distributing 1700 cells (with random orientations) in a sphere with radius of approximately 10µm. For the first 2 seconds we only activated the excluded volume forces in order to remove overlaps of the cells. After this time we switched on the production and the dynamics as described in the main text (see section 2). We simulated in silico colonies for 30 minutes and analyzed the properties of the colonies for the last 20 minutes. For every parameter set (see table S7) we analyzed a single colony. For those sets that are presented in the paper we analyzed 10 different realizations.

S4.4. Colony Coalescence
We initialized two colonies by randomly distributing N cells (with random orientations) per colonies in two sphere with radius 0.85 · N 1/3 µm and distances slightly larger than the sum of the radii of the colonies. For the first 2 seconds we only activated the excluded volume forces in order to remove overlaps of the cells. For the following 100 s we only allowed pili interactions between the cells of the individual colonies in order to create stable colonies. After this time we activated the interactions of pili of both colonies so that the coalescence could start. We simulated the in silico colonies for 30 minutes and analyzed the colonies starting from 60 seconds. For every parameter set we analyzed a single coalescence event (see table S8). For those sets that are presented in the paper we analyzed 14 different realizations (i.e. for the size-dependence of the coalescence).

S4.5. Assembly on a substrate
We initialized 1200 randomly oriented cells on the substrate with a relative distance similar to the cell radius. A fraction γ ∆ of these cells have pili that are not able to retract. For the first two seconds cells only interact via excluded volume forces until all overlaps of the cells with each other and the substrate have vanished. After this time we switched on all pili-mediated interactions and all related interactions of the cells. We simulated a total duration of 60 minutes. The system is characterized by periodic boundary conditions along both axes parallel to the substrate plane with a size 88 × 88 µm. For every parameter set (see table S9) we computed a single assembly realization.

S5. Data Analysis
We used Matlab R2015b and its internal functions to analyze the data.

S5.1. Single Cells
In order to compute the velocity at time t of a cell moving along the trajectory r(t) we computed the displacement ∆r(t) = r(t + τ ) − r(t) and set the velocity v(t) = ∆r(t)/τ with τ = 0.1 s.

S5.2. Single Colony on Substrate
The radii of microcolonies moving over a substrate were computed by projecting the shape of the individual cells onto the plane of the substrate. This two-dimensional projection exhibits to a good approximation a circular shape from which we extracted the radius R(t) of the colony. By averaging over both the time and different realisations of colonies we calculated the mean colony radius as a function of the number of individual cells within a colony. Additionally, we checked that the mean radius was not affected by the sampling time interval of the time average, for example due to the loss of individual cells from the colony. The diffusion coefficient D of the colonies were determined from the mean-squared-displacement where r col (t) denotes the center of the projected two-dimensional colony r col (t). The projection onto the plane of the substrate as discussed above allows us to define the radius and the position of a colony similar to experiments, in which one would detect the edges of a colony and analyze the corresponding binary image [9]. Similar to the computation of the projected colony shape on the substrate we defined the shape of the side-view of the colony by projecting the colony onto a plane perpendicular to the surface and computing its average shape by thresholding the sum of the individual areas of the different colonies.

S5.3. Free Single Colony
In order to compute the diffusion coefficients of cells within a microcolony as a function of their distance from the COM of the microcolony we introduced multiple shells defined by a minimal distance d min and d max from the the center of the microcolony. Afterwards, we picked those parts of the trajectory r i (t) for which cell i solely moved inside of a single shell for at least 10s. From these trajectories we computed the time-averaged mean squared displacement (which was additionally averaged over all cells i) and estimated the diffusion coefficients D from the relation (r i (t + τ ) − r i (t)) 2 i,t = 6Dτ. (S14)

S5.4. Colony Coalescence
In order to compute the height of the bridge forming between two colonies and the properties of an ellipse fitted to their shape, we projected the diplococcus shape of the individual cells onto a plane parallel to the axis between the centers of the two colonies. From the projection of the individual cells we can compute the envelope of the two colonies.
In order to compute the bridge height we defined a line of length 1µm centered around the COM of the projected envelope and parallel to the axis connecting the centers of the two colonies. We moved this line perpendicular to the axis between the colony centers and defined the height of the bridge as the range for which the whole line is inside of the colony envelope. The length of the line was chosen to be small enough so that it was not affected by the circular shape of the projection and large enough so that the effects of the fluctuations of single cells close to the surface of the colonies were reduced. For the late coalescence the results were compared to the short axis of the ellipse. By computing the central moments of the projected two-dimensional area of the microcolonies we were able to compute the properties of the ellipse as explained in [10]. This method allowed us to measure the bridge height and the ellipse properties in a similar manner to the analysis of experimental data.

S5.5. Assembly on a substrate
To define which cells are part of the surface and which cells are part of the bulk of a microcolony moving over a surface, we computed the alpha shape of the collection of cocci points of all cells [11]. The underlying idea of this triangulation technique is the motion of a sensor sphere over the surface of the colonies and the individual cells. Such a sphere is not allowed to intersect the positions of the cocci. If any coccus of a cell is touched by such a sphere (thus laying on the surface of the sphere), the corresponding cell is defined to be a surface cell. The only free parameter is the radius R α = 1 µm which is defined to be in the order of the cell sizes. More details of this method are given in reference [11].