Cluster and conquer: The morphodynamics of invasion of a compliant substrate by active rods
Mohammad Imaran, Mandar Inamdar, Ranganathan Prabhakar, Raghunath Chelakkot
CCluster and conquer: The morphodynamics of invasion of a compliant substrate byactive rods
Mohammad Imaran , , , Mandar M Inamdar , Ranganathan Prabhakar , Raghunath Chelakkot IITB-Monash Research Academy, Mumbai, India Department of Mechanical & Aerospace Engineering,Monash University, Clayton, Australia Department of Physics,Indian Institute of Technology Bombay, Mumbai, India Department of Civil Engineering,Indian Institute of Technology Bombay, Mumbai, India (Dated: January 12, 2021)The colonisation of a soft passive material by motile cells such as bacteria is common in biology.The resulting colonies of the invading cells are often observed to exhibit intricate patterns whosemorphology and dynamics can depend on a number of factors, particularly the mechanical propertiesof the substrate and the motility of the individual cells. We use simulations of a minimal 2D model ofself-propelled rods moving through with a passive compliant medium consisting of particles that offerelastic resistance before being plastically displaced from their equilibrium positions. It is observedthat the motility-induced clustering of active (self-propelled) particles is crucial for understandingthe morphodynamics of colonisation. Clustering enables motile colonies to spread faster than theywould have as isolated particles. The colonisation rate depends non-monotonically on substratestiffness with a distinct maximum at a non-zero value of substrate stiffness. This is observed tobe due to a change in the morphology of clusters. Furrow networks created by the active particleshave a fractal-like structure whose dimension varies systematically with substrate stiffness but isless sensitive to particle activity. The power-law growth exponent of the furrowed area is smallerthan unity, suggesting that, to sustain such extensive furrow networks, colonies must regulate theiroverall growth rate.
Large-scale collective pattern formation by self-motileelements is a widely studied phenomenon in physics andbiology [1–5] . A prominent class of such collective be-haviour such as flocking and swarming are caused by di-rect local interaction among individual entities withina group, influencing their relative motion [6–8]. Thesurrounding medium can also play a significant part inmediating interactions between active particles. In wetsystems, particularly at low Reynolds numbers, hydro-dynamic interactions are nearly instantaneous. Inhomo-geneous and complex environments qualitatively changethe individual, as well as collective dynamics of such sys-tems [9–16]. In dry systems, particles may leave behindchemical or mechanical cues for other particles to fol-low, thereby increasing the level of medium-induced com-plexity. This phenomenon of stigmergy has been stud-ied widely in the context of ants and termites leavingpheromone trails [17; 18].Recently, a mechanical form of stigmergy has been ob-served in motile bacteria
Pseudomonas aeruginosa and
Myxococcus xanthus when cultured on soft hydrogel sub-strates [19–21]. Under favourable conditions, these bac-teria form extensive networks of permanent furrows asthey move collectively as a monolayer across the soft sur-face in the initial stages of biofilm formation. The rateof expansion of the colony is further observed to be inti-mately related to the morphology of the network and thecellular traffic within the furrows [19–21]. Several recent studies have explored collective phenom-ena in populations of rod-shaped particles using com-puter simulations. The expansion rate and morphologyof colonies of growing and sub-dividing non-motile havebeen shown to strongly depend on the mechanical proper-ties of the passive medium or the extracellular materialsecreted by the rods [22–28]. Studies of self-propelledrods in the absence of coupled interactions with a sub-strate [29–40] reveal self-organized formation of dynamicsmall or large clusters, polar lanes, and nematic defects.In an attempt to understand the behaviour observed inexperiments in colonies of motile rod-shaped cells of
P.aeruginosa growing on agar [19], Zachreson et al. [41; 42]used simulations of self-propelled spherocylinders inter-acting with a continuum model of a deformable substrate.While their results demonstrate that substrate stiffnessand its viscoelastic relaxation time can strongly influencethe morphodynamics of active furrowing, their parame-ters were chosen specifically for the experimental systemat hand. Moreover, the simulations also included celldivision and population growth.Here, we take a closer look by using 2D simulations tostudy the furrow structure that emerges as a single rowof active self-propelled rods advances through a plasticsubstrate. We show that it is motility-induced cluster-ing that generically causes furrow networks to emergewhose fractal dimension depend on substrate plasticity.Clustering further enhances the rate at which the colony a r X i v : . [ c ond - m a t . s o f t ] J a n FIG. 1. (A) Furrow network (in white) formation in a substrate of plasticity, P = . = V ∗ c (filled symbols), and the normalised speed of isolated rods, V ∗ (open symbols). The speeds are normalised by theobserved speed of a single rod in a substrate of zero stiffness. When thermal fluctuations are weak, the average speed of anisolated rod in a fluid-like medium consisting of fully plasticised substrate is V ≈ F a /( γ s + γ r ) = Pe /[ N ( γ s / γ r + )] . edge advances, and this speed gain depends on clustermorphology. Our results also suggest that colonies mustregulate their overall cell growth rate in order to sustainextended furrow networks.In our simulations, each rod of length L = N b σ b is arigid linear array of N b beads of nominal diameter, σ b .The substrate is discretised as randomly-packed isotropicparticles (see Supplemental Information for detailed de-scription [43]). The rods interact with other rods andsubstrate particles through repulsive excluded-volumeforces modeled with the Separation-Shifted Lennards-Jones potential (SSLJ) [31]. The instantaneous config-uration of rod i are characterised by its position r i ( t ) and the unit polarity vector ˆ p i ( t ) . The time evolution ofthese variables are governed by Langevin equations thatinclude random Brownian forces with an energy scale of k b T and a self-propulsion force of magnitude F a that actsalong ˆ p i . Free rods experience frictional resistance totheir motion characterised by the bare rod friction coeffi-cient, γ r = γ b N b , where γ b is the bare friction coefficientof a bead on a rod. A minimal model is used to rep-resent the essential mechanical features of a semi-solid,plastic substrate. The p -th substrate particle is bound toits equilibrium coordinate r p, by a harmonic potential, U p = k ( r p − r p, ) /
2, if ∣ r p − r p, ∣ < l max , where k isthe elastic stiffness of the substrate, and l max is the max-imum displacement of substrate particles beyond whichthey are plasticised ( i.e. U p = k l /
2, a constant) anda permanent furrow is formed (see Fig. S1 [43]). Oncebroken off, the substrate particles remain as free particlesin the system and continue to offer a frictional resistanceto the rods. To model the high viscous resistance of semi-solid substrate to displacements, the friction coefficient,of a substrate particle γ s = γ b in our simulations. Forthe qualitative exploration here, we have neglected di-rect substrate-substrate interactions. The rod-substrateinteraction forces and the binding force derived from the potential U p determine the time-evolution of the sub-strate. The simulations are performed using LAMMPS.The masses of the particles are chosen to be small enoughfor the simulations to be in the overdamped regime (seeSupplemental Material [43]).Each simulation starts with two rows of verticallyaligned rods, introduced at the top and bottom bound-aries with propulsion forces directed towards the undis-turbed substrate that fills rest of the box. The heightof the simulation box , H , is twice its width, W . Peri-odic boundary conditions are imposed on all the sides.The simulation is stopped when any rod from either sidecrosses the middle of the box to the other side. Eachside of the box is then treated as a separate simulation in-stance. We set σ b as the length scale, and τ b = σ γ b / k b T and k b T as the time and energy scales, respectively. Weuse rods with N b = H = σ b and W = σ b . The key dimensionless pa-rameters are the P´eclet number, Pe = F a L / k b T and theplasticity ratio, P = k l max / F a , with l max = . σ b . Wefocus here on collective behaviour of the rods and its ef-fect on furrow formation at high P´eclet numbers at whichnoise does not destroy structure formation. Our choiceof parameters correspond to values of P ≪ i.e. we con-sider plastic substrates for which the stiffness is chosensuch that even single rods can displace the substrate tocreate permanent furrows.Initially, all the rods begin individually creating rela-tively straight and narrow furrows (see Fig. S3; Movie §
1, [43]). Small orientational fluctuations cause the rodsto begin colliding with neighbouring rods to begin form-ing dynamic clusters. As the clusters propel through thesubstrate, they permanently displace substrate particles,forming a complex network of permanent furrows (Fig 1-A). To analyse the furrow formation more quantitatively,we define the colonised area C , which is the y -coordinateof the outermost leader rod at any instant of time (see P ¯ v y , ¯ v R y , ¯ v T y ° º/ º/ ¡ p ( ¡ ) P f T Pe80100150200 Pe = 100 A Train
Raft Pe = 100 P = 0.01 B C ¯ υ T y ¯ υ R y ¯ υ y FIG. 2. (A) Fraction of train clusters, f T , out of all the clusters ( the fraction of rafts, f R = − f T ). (B) Distribution of rodorientation angle p ( φ ) with respect to the vertical axis ( y axis) for rafts (blue curve) and trains (orange curve) for Pe = P = .
01. (C) Individual averages of the y -component of cluster velocity of rafts and trains, and the overall populationaverage, for Pe = Fig. S2 [43]). This area is observed to grow nearly linearlyin time in all our simulations until they are terminated(see Fig. S4 [43]). The colonisation speed V c is estimatedby a linear least-squares fit through the C -vs- t data, andis compared in Fig. 1-B with V , the speed of an isolatedrod through the same substrate. The speeds are nor-malised by observed speed of a single rod in a substrateof zero stiffness. We observe that, at any given Pe and P , V c is systematically larger than V , the speed at whichan isolated rod moves through the substrate. Further,while V decreases nearly linearly with P at a fixed Pe, V c varies non-monotonically with P , displaying a peakcolonisation rate at a non-zero value of P .These interesting differences in colonisation speed fromisolated rod speed clearly arise from collective effects andappear to be related to the behaviour of clusters of activerods that drive the formation of furrows. We used a clus-tering algorithm to identify separate clusters of beads.Individual clusters were then tracked in a frontier regionof depth 8 L from the leading edge of the colony until theybroke up or were joined by new members. Frequency dis-tributions of quantities such as the number of rods in acluster ( i.e. the cluster size), average cluster speed andvertical velocity component, and orientation angle rela-tive to the vertical axis were determined.Based on the relative configurations of the rods in acluster, each cluster is categorised as one of two mutuallyexclusive types: end-to-end trains and non-train rafts that are usually arrow-head shaped (Fig. 1-A; see Supple-mentary Material for mathematical definitions [43]). Wefind that trains become more prevalent as the substratebecomes stiffer (Fig. 2-A; Movie §
2, [43]).Both rafts and trains move more quickly through thesubstrate than isolated rods. Since the propulsive forceon each rod is distributed uniformly over each bead, thetotal propulsive force on a tightly-packed cluster is pro-portional to its area. The resistance it experiences ishowever proportional to its frontal perimeter exposed tothe substrate. Consequently, the resistance experiencedper rod is smaller in a cluster. We therefore find that the average speed of clusters of either kind grows lin-early with cluster size (Fig. S5 in SI [43]). This effectis greater in trains, where the whole propulsive force ofthe train is only resisted by about one or two substrateparticles.In addition, the orientational distribution of clustersshows that trains in the frontier region close to leadingedge of the colony tend to be strongly aligned in theoutward direction (Figs. 2-B and S7 in SI [43]). Rodsin rafts have more diverse orientations, which can causerafts to have a broader range of orientations relative tothe outward direction. The propulsive force on trains,on the other hand, is almost entirely directed along thetrain axis, which keeps them on course.The mean y -component of the velocity ( v y ) of clus-ters within the frontier region is computed as v y = v R y ( − f T ) + v T y f T , where f T is the fraction of trains inthe clusters. In Fig. 2-C, the individual means, v R y and v T y , decrease non-monotonically with P , the overall av-erage v y shows a maximum around the same P at whichthe peak colonisation speed occurs in Fig. 1-B. There-fore, the increase of the colonisation speed with stiffnessat small values of P , appears to be because of the growthin the fraction of faster-moving trains that are also morealigned along the outward direction. When the fractionof trains begins to saturate at high stiffnesses, the de-crease of their speeds with stiffness takes over, and theoverall colonisation speed decreases with stiffness.The furrow network is created by clusters at the headof furrows repeatedly intersecting with previously cre-ated, empty, furrows. If a previous furrow a cluster en-counters is narrow, the cluster can plough through andcontinue unhindered. On the other hand, when the pre-vious furrow has a width comparable to L or greater, thecluster quickly breaks up, as the small orientational noisepresent causes rods at the edge of the cluster to escapethe cluster easily before the cluster passes through thefurrow (see Movie §
3, [43]). This is consistent with ob-servations elsewhere that motility-induced clustering can P . . . D t F Pe 10020080 A B
FIG. 3. (A) Box-counting fractal dimension as a function ofplasticity number P , measured at termination of simulation.(B) The growth of total furrowed area with time, for different P e and P . furrowed area growth shows a non-trivial powerlaw behaviour. The symbols are as follows: P = . ◻ )and P = . ∎ ) at Pe = P = . ▽ ) and P = .
075 ( ▼ )at Pe = P = . ◯ ) and P = . ● ) at Pe = lead to large clusters and lanes in systems of free activerods at sufficiently high densities [31], but in an emptyfurrow, a cluster can quickly “evaporate”. We observe vi-sually that arrow-head shaped rafts are stable for longertimes. The rods at the head of such pointy rafts haveconvergent orientations. This leads to a rectification oftheir propulsion forces along the longitudinal axis of thecluster. Any remaining transverse component of the netpropulsion force causes pointy rafts to swerve and takecurved trajectories. Large rafts can thus create wide andlong furrows that serve as arterial conduits in the network(Fig. S8-D; Movie § § § § D , of the furrow-substrateinterface using a box-counting algorithm [44]. As ex- pected, D has a non-integer value between 1 and 2, andis observed to systematically decrease with P (Fig. 3-A). This is associated with a change from networks createdmostly by interactions of curved raft trajectories ( e.g. Fig. 1-A) to those dominated by the thin and straightfurrows created by trains ( e.g.
Fig. 1-A)This observation suggests a simple model for thegrowth of the area of the furrow network, F , with time.Furrow heads can form randomly at any point on thefurrow interface within the network due to the constantscattering of free rods through the network by clusterbreakup. The furrowing rate is then expected to be pro-portional to the length, Z , of the furrow interface. Thefractal nature of the interface implies that Z ∼ F D / [45].The rate of furrow growth, dF / dt ∼ Z R /F , where R isthe total area occupied by rods in the furrows (a constantin our simulations). Integrating, we obtain F ∼ t /( − D ) at large times. Since the fractal dimension, 1 < D <
2, weobserve in Fig. 3-B that the power-law exponent for thegrowth of F with t is in the range 2 / < /( − D ) <
1. Ifthis can be sustained, while the colonized area C growslinearly with time, the colony morphology will becomeincreasingly sparse.Based on the results above, one can expect similar fur-row networks to generically form when colonies of motilecells spread through soft, plastic environments and inthree dimensions as well. The observations here are inline with many experiments with motile bacteria. Cellrafts similar to those discussed above have been observedto lead the formation of furrow networks in P. aerug-inosa and
M. xanthus [19; 21] cultured under confine-ment on semi-solid agar. These furrow networks extendacross distances around two orders of magnitude greaterthan the length of a single cell. The networks are fur-ther lined with extracellular DNA (eDNA) and extra-cellular polysaccharides (EPS). Other experiments havefurther demonstrated that the three-dimensional latticesof eDNA are necessary for the mechanical integrity ofbiofilms of bacterial pathogens in complex biofluids suchas sputum and otorrhea (ear-discharge fluid)[46]. Thecurrent study provides a natural mechanism for such net-work templates to spontaneously emerge in such systems.Our observation of the power-law growth of F furthersuggests that, in order to create extended furrow net-works such as those observed in experiments, there mustbe mechanisms that limit the growth rate of cells to notexceed the rate at which the furrows are created. Un-constrained exponential growth or even linear growth ofcells quickly obliterates the furrow network, leading toa dense colony (Fig. S9; Movie § P. aeruginosa along with associated release of eDNAand other “public goods” required for biofilm formation[49]. Beyond bacterial collectives, substrate interactionsalso occur during cancer metastasis in which groups ofmigrating cancer cells exploit the plasticity of the extra-cellular matrix (ECM) and actively remodel the ECMfibers to enhance their migration speeds[50–52].In conclusion, the cluster-and-conquer mechanism ob-served in our simulations may be crucial in biofilmformation by motile bacteria in plastic environments,and experimental evidence suggests that bacteria mayhave adaptive mechanisms to exploit the fractal-like spa-tial structure of the furrow network in building robustbiofilms. The results here suggest future theoretical andexperimental avenues for exploring the role played by me-chanical stigmergy in these biological phenomena.
Acknowledgments : RC acknowledges the financial sup-port from SERB, India via the grants SB/S2/RJN-051/2015 and ECR/2017/000744. We acknowledge fi-nancial support from IITB-Monash Research Academyand computational-time grants from the National Com-putational Infrastructure, Canberra, the MonArch facil-ity at Monash University and the SpaceTime facility (IITBombay). [1] C. Dombrowski, L. Cisneros, S. Chatkaew, R. E. Gold-stein, and J. O. Kessler, Physical review letters ,098103 (2004).[2] M. Ballerini, N. Cabibbo, R. Candelier, A. Cavagna,E. Cisbani, I. Giardina, V. Lecomte, A. Orlandi,G. Parisi, A. Procaccini, et al. , Proceedings of the na-tional academy of sciences , 1232 (2008).[3] Y. Katz, K. Tunstrøm, C. C. Ioannou, C. Huepe, andI. D. Couzin, Proceedings of the National Academy ofSciences , 18720 (2011).[4] M. C. Marchetti, J.-F. Joanny, S. Ramaswamy, T. B.Liverpool, J. Prost, M. Rao, and R. A. Simha, Reviewsof Modern Physics , 1143 (2013).[5] S. Ramaswamy, Annu. Rev. Condens. Matter Phys. ,323 (2010).[6] T. Vicsek, A. Czir´ok, E. Ben-Jacob, I. Cohen, andO. Shochet, Physical review letters , 1226 (1995).[7] G. Gr´egoire and H. Chat´e, Physical review letters ,025702 (2004).[8] H. Chat´e, F. Ginelli, G. Gr´egoire, F. Peruani, and F. Ray-naud, The European Physical Journal B , 451 (2008).[9] C. Bechinger, R. Di Leonardo, H. L¨owen, C. Reichhardt,G. Volpe, and G. Volpe, Reviews of Modern Physics ,045006 (2016).[10] J. Elgeti, R. G. Winkler, and G. Gompper, Reports onprogress in physics , 056601 (2015).[11] K. Qi, E. Westphal, G. Gompper, and R. G. Winkler,Physical Review Letters , 068001 (2020).[12] S. M. Mousavi, G. Gompper, and R. G. Winkler, SoftMatter , 4866 (2020).[13] S. Das and R. Chelakkot, Soft Matter , 7250 (2020).[14] A. E. Patteson, A. Gopinath, and P. E. Arratia, Naturecommunications , 1 (2018). [15] L. Theeyancheri, S. Chaki, N. Samanta, R. Goswami,R. Chelakkot, and R. Chakrabarti, Soft Matter , 8482(2020).[16] S. Das, S. Ghosh, and R. Chelakkot, Physical Review E , 032619 (2020).[17] A. B. Attygalle and E. D. Morgan, in Advances in insectphysiology , Vol. 18 (Elsevier, 1985) pp. 1–30.[18] D. E. Jackson, M. Holcombe, and F. L. Ratnieks, Nature , 907 (2004).[19] E. S. Gloag, L. Turnbull, A. Huang, P. Vallotton,H. Wang, L. M. Nolan, L. Mililli, C. Hunt, J. Lu, S. R.Osvath, et al. , Proceedings of the National Academy ofSciences , 11541 (2013).[20] E. S. Gloag, L. Turnbull, and C. B. Whitchurch, Scien-tifica (2015).[21] E. S. Gloag, L. Turnbull, M. A. Javed, H. Wang, M. L.Gee, S. A. Wade, and C. B. Whitchurch, Scientific reports , 26005 (2016).[22] F. Farrell, O. Hallatschek, D. Marenmanzzo, and B. Wa-claw, Physical review letters , 168101 (2013).[23] P. Ghosh, J. Mondal, E. Ben-Jacob, and H. Levine, Pro-ceedings of the National Academy of Sciences , E2166(2015).[24] J. Tchoufag, P. Ghosh, C. B. Pogue, B. Nan, and K. K.Mandadapu, Proceedings of the National Academy ofSciences , 25087 (2019).[25] Z. You, D. J. Pearce, A. Sengupta, and L. Giomi, PhysicalReview X , 031065 (2018).[26] R. D. Acemel, F. Govantes, and A. Cuetos, Scientificreports , 1 (2018).[27] F. D. Farrell, M. Gralka, O. Hallatschek, and B. Waclaw,Journal of The Royal Society Interface , 20170073(2017).[28] Z. You, D. J. Pearce, A. Sengupta, and L. Giomi, Physicalreview letters , 178001 (2019).[29] F. Peruani, A. Deutsch, and M. B¨ar, Physical Review E , 030904 (2006).[30] S. R. McCandlish, A. Baskaran, and M. F. Hagan, SoftMatter , 2527 (2012).[31] M. Abkenar, K. Marx, T. Auth, and G. Gompper, Phys-ical Review E , 062314 (2013).[32] A. Baskaran and M. C. Marchetti, Physical Review E ,011920 (2008).[33] S. Weitz, A. Deutsch, and F. Peruani, Physical ReviewE , 012322 (2015).[34] K. Prathyusha, S. Henkes, and R. Sknepnek, PhysicalReview E , 022606 (2018).[35] C. A. Velasco, M. Abkenar, G. Gompper, and T. Auth,Physical Review E , 022605 (2018).[36] H.-S. Kuan, R. Blackwell, L. E. Hough, M. A. Glaser, andM. D. Betterton, Physical Review E , 060501 (2015).[37] M. B¨ar, R. Großmann, S. Heidenreich, and F. Peruani,Annual Review of Condensed Matter Physics , 441(2020).[38] A. Be’er, B. Ilkanaiv, R. Gross, D. B. Kearns, S. Heiden-reich, M. B¨ar, and G. Ariel, Communications Physics ,1 (2020).[39] G. A. Vliegenthart, A. Ravichandran, M. Ripoll,T. Auth, and G. Gompper, Science advances , eaaw9975(2020).[40] M. C. Bott, F. Winterhalter, M. Marechal, A. Sharma,J. M. Brader, and R. Wittmann, Physical Review E ,012601 (2018).[41] C. Zachreson, C. Wolff, C. B. Whitchurch, and M. Toth, Physical Review E , 012408 (2017).[42] C. Zachreson, X. Yap, E. S. Gloag, R. Shimoni, C. B.Whitchurch, and M. Toth, Physical Review E , 042401(2017).[43] M. Imaran, R. Prabhakar, R. Chelakkot, and M. Inam-dar, (2020).[44] J. Gagnepain and C. Roques-Carmes, wear , 119(1986).[45] B. J. Florio, P. D. Fawell, and M. Small, Powder tech-nology , 551 (2019).[46] A. Devaraj, J. R. Buzzo, L. Mashburn-Warren, E. S.Gloag, L. A. Novotny, P. Stoodley, L. O. Bakaletz, andS. D. Goodman, Proceedings of the National Academy ofSciences , 25068 (2019).[47] A. Nagilla, R. Prabhakar, and S. Jadhav, Phys. Fluids , 022109 (2018). [48] C. Giverso, M. Verani, and P. Ciarletta, Journal of TheRoyal Society Interface , 20141290 (2015).[49] L. Turnbull, M. Toyofuku, A. L. Hynen, M. Kurosawa,G. Pessi, N. K. Petty, S. R. Osvath, G. C´arcamo-Oyarce,E. S. Gloag, R. Shimoni, U. Omasits, S. Ito, X. Yap,L. G. Monahan, R. Cavaliere, C. H. Ahrens, I. G. Charles,N. Nomura, L. Eberl, and C. B. Whitchurch, Nat. Com-mun. , 11220 (2016).[50] B. Lee, J. Konen, S. Wilkinson, A. I. Marcus, andY. Jiang, Scientific reports , 39498 (2017).[51] K. M. Wisdom, K. Adebowale, J. Chang, J. Y. Lee,S. Nam, R. Desai, N. S. Rossen, M. Rafat, R. B. West,L. Hodgson, et al. , Nature communications , 1 (2018).[52] J. Winkler, A. Abisoye-Ogunniyan, K. J. Metcalf, andZ. Werb, Nature communications11