PNAS Nexus, 2022, 1, 1–11
https://doi.org/10.1093/pnasnexus/pgac269

Research Report

Takuro Shimaya

and Kazumasa A. Takeuchi*

Department of Physics, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, 113-0033 Tokyo, Japan∗
To whom correspondence should be addressed: Email: kat@kaztake.org
Edited By: Peter Fratzl

AbstractRod-shaped bacteria, such as Escherichia coli, commonly live forming mounded colonies. They initially grow two-dimensionally ona surface and finally achieve three-dimensional growth. While it was recently reported that three-dimensional growth is promoted
by topological defects of winding number +1/2 in populations of motile bacteria, how cellular alignment plays a role in nonmotile
cases is largely unknown. Here, we investigate the relevance of topological defects in colony formation processes of nonmotile E. coli
populations, and found that both ±1/2 topological defects contribute to the three-dimensional growth. Analyzing the cell flow in thebottom layer of the colony, we observe that +1/2 defects attract cells and −1/2 defects repel cells, in agreement with previous studies on
motile cells, in the initial stage of the colony growth. However, later, cells gradually flow toward −1/2 defects as well, exhibiting a sharp
contrast to the existing knowledge. By investigating three-dimensional cell orientations by confocal microscopy, we find that verticaltilting of cells is promoted near the defects. Crucially, this leads to the emergence of a polar order in the otherwise nematic twodimensional cell orientation. We extend the theory of active nematics by incorporating this polar order and the vertical tilting, which
successfully explains the influx toward −1/2 defects in terms of a polarity-induced force. Our work reveals that three-dimensional
cell orientations may result in qualitative changes in properties of active nematics, especially those of topological defects, which may
be generically relevant in active matter systems driven by cellular growth instead of self-propulsion.
Keywords: bacteria, active nematics, topological defect, colony growth

Significance Statement:
Biofilms are three-dimensional dense aggregates of bacteria that cause various problems from industry to daily life. Using nonmotile growing Escherichia coli, we discovered that topological defects, i.e., locations of misalignment in cell orientation, attract
cells, and promote local 3D growth. Regarding the two types of defects labeled “±1/2,” while previous studies showed that only
+1/2 defects attract cells, we found that −1/2 defects also become attractive in growing nonmotile bacteria. We revealed that thisis realized by 3D cell tilting around defects, which leads to the emergence of a polar order and generates a polarity-induced force.
These may be characteristic of growing active matter without self-propulsion of constituent particles, with possible implications
for the morphogenesis of various cell populations and tissues.

Introduction
Numerous species of bacteria live in dense populations, which
often take the form of biofilms (1). Besides being a challenging
subject for biologists and physicists, because biofilms cause a variety of problems in medicine, industry, and our daily life (2, 3),
understanding the mechanism of biofilm formation is a crucial
mission across diverse disciplines. In the early stage of biofilm
formation processes, two-dimensional colonies are first formed,
then a three-dimensional structure is eventually constructed (1).
Because mechanical interactions between cells are important at
this stage, many studies have attempted to understand structureformation dynamics from a physical perspective (4).
In particular, rod-shaped bacteria, irrespective of whether they
are motile or not, are aligned with each other and behave like anactive nematic liquid crystal in a dense two-dimensional space(5–13). For motile bacteria, it has recently been reported that +1/2
topological defects promote three-dimensional growth of Myxococcus xanthus populations (13). Besides bacteria, it is known that
topological defects also play decisive roles in various kinds of cell
populations (14, 15), such as epithelial cells (16), neural stem cells(17), fibroblasts (18, 19), and actin fibers in Hydra (20). However, forgrowing but nonmotile bacteria, while some studies investigated
how nonmotile cells initiate three-dimensional growth (21–30),
the relevance of local cell alignment to three-dimensional growth,in particular, that of topological defects, remains unknown.
Here, by observing colony formation processes of nonmotile Escherichia coli between a coverslip and a nutrient agar pad [Fig. 1(a)],we find an indication that both +1/2 and −1/2 topological defects promote three-dimensional growth of colonies. This finding is put on solid ground by analyses of the two-dimensional

Competing Interest:The authors declare no competing interest.
Received: July 18, 2022. Accepted: November 30, 2022

C The Author(s) 2022. Published by Oxford University Press on behalf of National Academy of Sciences. This is an Open Access article distributed under the terms of the Creative Commons Attribution-NonCommercial-NoDerivs licence (https://creativecommons.
org/licenses/by-nc-nd/4.0/), which permits non-commercial reproduction and distribution of the work, in any medium, provided theoriginal work is not altered or transformed in any way, and that the work is properly cited. For commercial re-use, please contact
journals.permissions@oup.com

Downloaded from https://academic.oup.com/pnasnexus/article/1/5/pgac269/6948003 by University of Tokyo Library user on 22 December 2022

Tilt-induced polar order and topological defects ingrowing bacterial populations

2 |PNAS Nexus, 2022, Vol. 1, No. 5

(a)

(b)

Plastic cover

Frame
seal

Agar

SpacerCoverslip
Stage
Objective lens

10μm

Fig. 1. Morphology of three-dimensional colonies formed from numerous cells observed by end-point confocal microscopy. The confocal data weretaken 14 hours after the cells had filled the bottom plane. (a) Experimental setup. Bacterial cells were between a coverslip and a nutrient agar pad. (b) A
three-dimensional image of colonies (184.52 × 184.52 μm2 ), where cells were stacked three-dimensionally. The region surrounded by the red rectangle
is enlarged and displayed in the upper right side. (c) A two-dimensional cross-section showing the bottom layer. Red comets and blue trefoils indicate+1/2 and −1/2 defects, respectively. The arms of the symbols reflect the structure of the director field as illustrated in the insets of (b). (d) Mean colony
height at the location of the defects in the bottom layer and that far from defects. At each position in the xy-plane, the height was evaluated by thelength of the profile along the z-axis whose intensity is higher than 20% of the maximum. The heights at the defect positions were extracted from allof the hundreds of defects that were sufficiently far (>9 μm) from each other. For the colony heights far from defects, we randomly picked up 1000
points, which were sufficiently far (>9 μm) from any defect (see the experimental methods discussed in Appendix A). The error bars indicate thestandard error from the ensemble averaging. The optical resolution was about 250 nm in the vertical direction; see also Fig. S2(a) to (c) for the colony
height distributions. Here, we show the result of a single measurement (uniform colony, end-point confocal #1); see Fig. S2(e) for the result of anotherbiological replicate (uniform colony, end-point confocal #2), where we confirmed that the colony height was again higher at the locations of the defects.

velocity field around topological defects, which reveal that cells
are transported toward both +1/2 and −1/2 defects, albeit weakly,
implying upward growth there. Remarkably, this influx towardboth types of defects is contrary to the existing knowledge thatcells escape from −1/2 defects (5, 6, 13, 16, 17, 19), and cannot be
explained by the conventional active nematic theory. Combining
confocal observations and theoretical modeling, we find that thethree-dimensional tilting of cells is promoted around topologicaldefects, which can induce additional force around defects. Crucially, we uncover the formation of a polar order due to threedimensional asymmetric tilting of cells around defects, whichturns out to be the key to theoretically account for the emergenceof the influx toward −1/2 defects.

Results
Topological defects promote three-dimensional
growth of bacterial colonies
First, we studied the relation between cell orientation and colony
structure, using nonmotile E. coli placed between a coverslip anda nutrient agar pad [Fig. 1(a); see the experimental methods discussed in Appendix A]. We put cell suspension on the coverslip sothat cells are initially distributed densely and uniformly. Then, we
cultured it for 14 hours after cells had filled the bottom plane andobserved the resulting three-dimensional colony, which consistedof multiple layers of tilted cells, by confocal microscopy. Here, wetook only a single confocal image at this end point, to take a highquality image without photobleaching.
To test the relevance of cell alignment to the three-dimensional
growth, we investigated whether the presence of topological defects influenced the colony height. First, we noticed that the orientation of cells in the bottom layer was nearly horizontal [Fig. 1(b)and (c)], tilted only weakly [typically ∼10◦ in this end-point observation; see Fig. 3(a) and descriptions thereof], and many topological defects were observed. We also observed subsequent layers
and found that the cell orientation therein was essentially guided
by that of the bottom layer [Fig. S1(c) to (e)]. This led us to focus
on the bottom layer and regard it as a quasi-two-dimensional active nematic system. We measured the two-dimensional orientation of cells, n(R) at position R in the bottom layer, from theimage intensity using the structure tensor method [see Fig. S1(a)and the experimental methods discussed in Appendix A). We thendetected topological defects [Figs. 1(c) and S1], and measured thecolony height at the positions of the defects (see the experimental
methods discussed in Appendix A). For comparison, we also measured the colony height at randomly selected locations that are
sufficiently far from topological defects. We found that the meancolony height is slightly higher at the positions of the defects[Fig. 1(d)], both +1/2 and −1/2, than in the regions far from thedefects. The statistical significance is confirmed by the Wilcoxon
rank-sum test (Fig. S2). For the null hypothesis that the median of
the height distribution at the positions of the defects is identicalto that far from defects, the P-value was 0.018 for the +1/2 defects
and 0.043 for the −1/2 defects. While these results were obtained
from observations of 20 separate regions in a single experiment,
the reproducibility was also confirmed by another biological replicate using a different substrate and agar pad [Fig. S2(e)]. These results suggest that topological defects promote the vertical growthof colonies, slightly but statistically significantly.

Two-dimensional velocity fields around
topological defects
To clarify the origin of the promoted three-dimensional growth,
we investigate how cells in the bottom layer were displaced near

Downloaded from https://academic.oup.com/pnasnexus/article/1/5/pgac269/6948003 by University of Tokyo Library user on 22 December 2022

(d)

(c)

Shimaya and Takeuchi | 3(a)phase contrast

corresponding
side view
10μm-20 min

0 min

(c)

(b)

20 min(d)defectdefect(+or-)

(e)

(f)

or
defect(g)(h)

Phase contrast

Fig. 2. Phase-contrast images and results of two-dimensional velocity analyses for uniform colonies formed from numerous cells observed byphase-contrast microscopy. The time interval of the time-lapse observation was 1 min. The double sign ± corresponds to the sign of the defects. These
data were obtained by a single measurement (uniform colony, phase contrast #1), but the reproducibility was confirmed in Fig. S5 by using anotherbiological replicate with a different substrate and agar pad (uniform colony, phase contrast #2). (a) Phase-contrast images of cells taken at t = −20, 0,
20 min, from left to right. t = 0 is the moment when cells filled the bottom plane; see also Videos 1 and 2. The sketches on the right side illustrate how
three-dimensionally tilted cells appear in the phase-contrast images. (b and c) Velocity field v± (r) (black arrows) and its divergence (color map) around+1/2 defects (b) and −1/2(c). (d) Schematic illustration of the definition of the radial velocity v±
r (r, t). With this, the mean radial velocity is defects
1
defined by v̄±dφv±
r (r, t) ≡ 2π
r (r, t), which corresponds to the average of the radial components of the velocity over a circumference of radius r centered±at the defect. The sign of v̄±r indicates the direction of net flow, including the existence of vertical growth if v̄r < 0. (e and f) Time evolution of the mean
radial velocity v̄±r (r, t) around +1/2 defects (e) and −1/2 defects (f). Here, we used the velocity field averaged over 0 min ≤ t ≤ 15 min, 15 min ≤ t ≤30 min, and 30 min ≤ t ≤ 45 min. The error bars indicate the time average of the standard error evaluated from each frame. (g) Time evolution of theminimum of v̄±r (r, t) in the region r < r0 = 2 μm near the defect. The moving average taken from (t − 5 min) to (t + 5 min) is shown with thecorresponding error bar. (h) Comparison of the mean radial velocity v̄±r (r) between the experimental data (symbols) and theoretical curves based on
the conventional equation (Eq. 1) (lines). The displayed experimental data (symbols) are identical to those shown in (e and f) for 15 min ≤ t ≤ 30 min(for which the influx toward −1/2 defects was strongest). The vertical dashed lines indicate the defect core radius [see Fig. S6(a) and (b)]. The dottedlines represent the results for  =  0 and an = a0n ( = 0.25, an /ξ 0 = 0.055 μm2 /min, rS = 1.2 μm; see Appendix B), which correspond to the conventionalcase of extensile active nematics. The solid lines are the results for Eq. (1) with three-dimensional nematic tilting, i.e., an = a0n cos θn± (r). Specifically, weused θn± (r, φ) = θn∞ + (θn0 − θn∞ ) exp(−r2 /rθ2 ) with θn∞ = 0.3 and θn0 = 0.75, with the other parameters left unchanged (see Appendix B). The inset is a
close-up of the results for −1/2 defects.

topological defects. We conducted a time-lapse phase-contrast
observation of the bottom layer of cells, cultured from densely
and uniformly distributed populations as in the confocal observation (see the experimental methods discussed in Appendix
A). Cells then filled the two-dimensional plane rather homoge-

neously, without forming visible microcolonies, and after a shortwhile, cells started to tilt upward, almost simultaneously [Fig. 2(a)and Video 1; the appearance of a dark spot in the cell body indicates the tilting of the cell, as sketched in Fig. 2(a)]. Based on theuniformity of this initial two-dimensional growth, as compared

Downloaded from https://academic.oup.com/pnasnexus/article/1/5/pgac269/6948003 by University of Tokyo Library user on 22 December 2022

Bottom layer

4 |PNAS Nexus, 2022, Vol. 1, No. 5

Theoretical analyses and relevance of
three-dimensional tilting of cells
To account for the observed influx toward −1/2 defects, we developed a theory based on two-dimensional extensile active nematics, extended to incorporate characteristics of growing nonmotile colonies we observed. Following earlier studies (7, 9, 13,17), we describe the cell alignment by the nematic order tensor
Q(r, t) ≡ S(2n ⊗ n − 1), with the scalar nematic order parameter

S(r, t), the director field n(r, t), and the identity matrix 1 (see Appendix B). As a result of cell growth along the long axis of the cellbody, interacting with nearby cells, cells exert the extensile active stress σ = −an Q with the active stress coefficient an ( > 0) even
without the motility (7, 9). This stress induces the force f = ∇ · σ
and drives the velocity field v(r, t). In the overdamped and low
Reynolds number limit, this active force is balanced by the friction originating from cell–substrate interaction, giving the following linearized equation:
ξv = ∇ · (−an Q ),(1)with the friction tensor ξ. We assume that the friction is
anisotropic with respect to the cell alignment: ξ = ξ0 (1 − Q ) withthe friction anisotropy parameter . As suggested in ref. (32),
we may reasonably assume that it is easier for E. coli cells toslide along their longitudinal axis, hence,  > 0. Setting Q withthe theoretical director configuration for ±1/2 defects, n± (r) =(cos(±φ/2), sin(±φ/2)) with azimuth φ of the coordinate r, and using the experimentally determined core radius [see Appendix B
and Fig. S6(a) and (b)], we calculated the mean radial velocity v̄±r (r)(Supplementary Information, Section I) and display the curves inFig. 2(h) (dotted lines). Being consistent with the literature (13, 17),this shows influx only for +1/2 defects (v̄+r (r) ≤ 0) and outflux onlyfor −1/2 defects (v̄−r (r) ≥ 0). Therefore, to explain the experimentally observed influx toward −1/2 defects, we need to extend theexisting theoretical framework described so far. At this point, it
dv̄−
dv̄+
is useful to remark that drr (0) = 0, while drr (0) < 0 [see Eq. (S7)in Supplementary Information], so that the curve v̄−r (r) for −1/2
defects stays near zero for small r [blue dotted line in Fig. 2(h)].
Therefore, even a minor effect may be sufficient to change thesign of v̄−r (r) near the defects. One may consider that it might be
the density heterogeneity, in particular small voids observed at−1/2 defects in the early stage of the process (Video 1). However,this is unlikely to explain the observed influx, because more voids
existed at earlier times whereas the influx developed later [compare Video 1 and Fig. 2(g)]. We also examined the possibility thatthe cell growth may generate an influx toward −1/2 defects, byadding a growth term to the hydrodynamic equation, but this did
not yield the influx (see Supplementary Information).
Instead of the growth and the density heterogeneity, here, we
focus on the three-dimensional orientations of the cells, because
the influx toward −1/2 defects became strong when cells began
to tilt three-dimensionally [Fig. 2(f) and (g) and Videos 1 and 2) despite the decay of the overall flow speed. We experimentally measured the tilt angle θ n of cells from the horizontal plane [see the illustration in Fig. 3(a)], at a late time t from the end-point confocal
data, by the structure tensor method for the three-dimensional
space, applied to the bottom layer (see the experimental methods
discussed in Appendix A). Taking average over the regions around±1/2 defects, we obtained a field of the tilt angle, θn± (r) (Fig. 3).
While almost all cells were already tilted (hence, θn± (r) > 0 everywhere) at the moment of the end-point observation, we found thatthree-dimensional tilting was strongest at the core of both defects [Fig. 3; see also Fig. S7(a) and (b) for the results of anotherbiological replicate]. The peak of θn± (r) is well approximated by
a Gaussian function centered at the defect core plus a constant[Fig. 3(a)], and θn± (r) turned out to be essentially isotropic [Fig. 3(b)and (c)].
We considered that this tilting may have weakened, in our twodimensional description, the local active stress and the friction
anisotropy around the defects. More quantitatively, we assume
that the local active stress coefficient and the friction anisotropy

Downloaded from https://academic.oup.com/pnasnexus/article/1/5/pgac269/6948003 by University of Tokyo Library user on 22 December 2022

to the growth of a single circular colony discussed later, we shallrefer to the present case as “uniform colony” in the following. Using the images after the two-dimensional plane was filled, we detected topological defects from the two-dimensional cell orientation n(R, t) of the bottom layer, where t ≥ 0 is the time elapsed
since the bottom plane was filled. The density of defects initially
increased slightly, then stayed approximately constant from t ≈30 min [Fig. S3(a)]. As expected from the absence of cell motility,
the defects moved only a little in our system, typical displacements being a few microns over the observation time [Fig. S3(b)].
We then measured the velocity field around defects by particle
image velocimetry (PIV) (see the experimental methods discussed
in Appendix A). In Fig. 2(b) and (c), the arrows show the velocity
field v± (r, t) around ±1/2 defects, time-averaged over 30 min ≤ t≤ 105 min, where r indicates the position relative to the defectand the double sign corresponds to the sign of the defect [see also
Fig. S4(a) and (b)]. While the structure of v± (r, t) resembles thosearound defects in typical extensile active nematic systems (31),their divergence ∇ · v± (r, t) [Fig. 2(b) and (c); see also Fig. S4(e)] reveals a distinguished character of our system: we found negative
divergence around both types of defects, not only around +1/2 defects [Fig. 2(b)] as previously reported for systems of motile cell
populations (5, 6, 13, 16, 17, 19), but even around −1/2 defects[Fig. 2(c)], as opposed to those earlier studies. Since negative divergence indicates influx of cells, this implies that cells are moving
toward both types of defects in the bottom layer and pushed out
upward. This is consistent with the result of the confocal observation that the colony height was higher at the positions of the±1/2 defects. To inspect the time evolution of this influx, we examined the mean radial velocity at a distance r from +1/2 or −1/2

1
±
dφv±
defect, v̄±r (r, t) ≡ 2π
r (r, t), where vr (r, t) is the radial component of the velocity v± (r, t) at polar coordinates r = (r, φ) centeredat the defect [Fig. 2(d)]. For the +1/2 defects [Fig. 2(e)], we find thatv̄+r (r) is essentially negative all the time, but the depth of the minimum decreased with increasing time. This may be because of
decay of the overall flow speed throughout the colony [Fig. S4(f)],possibly due to nutrient starvation, pressure increase, and/or quorum sensing. In contrast, for the −1/2 defects [Fig. 2(f)], v̄−r (r) wasinitially positive for all r, but decreases near the defect and eventually becomes negative. To see the time-dependent influx toward
the defects more clearly, we plotted minr<r0 v̄±r (r) with r0 = 2 μm
in Fig. 2(g). While the strength of the influx toward the +1/2 defectmonotonically decreased, that toward the −1/2 defect increased
until t 25 min. These suggest an intrinsic change in the dynamics around the −1/2 defect that cannot be explained by the decayof the overall flow speed. The reproducibility was confirmed by anindependent biological replicate (Fig. S5). Although the strength of
the influx toward −1/2 defects is weak [roughly 4 × 10−3 μm/min
lasting over 15 min; see Fig. 2(f) and S5(b)], it is consistent withthe amount of the increment in the colony height we observed[Fig. 1(d)]. As it is reproducible and not understandable with thecurrent knowledge of active nematics, this led us to seek for a
possible mechanism of the influx toward −1/2 defects.

Shimaya and Takeuchi | 5

(a)

(b)

(c)

are given by an (r, t) = a0n cos θn (r, t) and (r, t) = 0 cos θn (r, t), respectively, with constants a0n and  0 . Using this, we solved Eq. (1)and found that the influx toward −1/2 defects can emerge[Fig. 2(h) blue solid line and inset; see also Fig. S6(c)] within thereasonable range of parameter values. However, the strength of
the influx was too small to account for the experimental result[Fig. 2(h) blue symbols, to be compared with the blue solid line].
This led us to seek for another key factor for the influx toward−1/2 defects.
Here, we propose a key mechanism for the observed influx toward −1/2 defects. So far, we assumed that active force is induced
only by nematic alignment. However, when cells are tilted threedimensionally, the sign of the tilt angle θ p may break the nematic
symmetry and make it possible to develop a polar order [Fig. 4(a)].
If this happens, the violation of the nematic symmetry may result in the generation of an additional force term that is otherwise
forbidden, which needs to be included in the force balance Eq. (1).
Such a polarity-induced force is expected to be proportional to thestrength of the polar order, i.e., θ p , in its lowest order, and act in the
direction of the director. In this context, it is interesting to refer topast experiments on densely packed vibrated granular rods (33,34), which indeed showed the formation of the polar order due
to rod tilting and the resulting horizontal transport of the rods,
driven by the polarity-induced force. This suggests that a similar
polarity-induced force may arise in our growing bacterial populations, resulting from the extensile active force of cells, if the polar order is formed. Inspired by this possibility, we measured θp (r)around both types of defects by end-point confocal microscopy.
Note that the single-cell tilt angles fluctuate largely from cell tocell (see Videos 1 and 2), and this is why θ p and θ n , i.e., the signed
and unsigned averages of the tilt angles, respectively, differ. The
sign of θ p is determined by choosing the direction of the head of
the nematic director n [see Fig. 4(a) and Appendix A]: here, we set
n± (r) = (cos(±φ/2), sin(±φ/2)) for the director field around ±1/2
defects. Figure 4(b) displays the result on the +x-axis. This shows
nonvanishing θp (r) for both ±1/2 defects, specifically θp (r) > 0 (upper end oriented outward) for +1/2 defects and θp (r) < 0 (upper
end oriented inward) for −1/2 defects on the +x-axis, demonstrating the emergence of the polar order in our growing bacterial
populations. Consequently, the above-mentioned symmetry argument predicts the polarity-induced force f p to arise, which satisfies f p ∝ θp n for small θ p . In Fig. 4(c) and (d), we show θp± (r)n± (r)around ±1/2 defects, which represent the strength and the direc-

tion of the polarity-induced force f p (r). What contributes to the±mean radial velocity is its radial component fp,r(r), proportional to
(r),
which
is
shown
in
Fig.
4(e)
and
(f),
with n±θp± (r)n±r
r (r) being theradial component of the director n± (r). These results show that,while the polarity-induced force around +1/2 defects drives thedefects toward their comet tail, that around −1/2 defects acts inward, leading to the influx toward the defects. We confirmed the
reproducibility of the main structure of the polarity-induced force
by taking a biological replicate (Fig. S7), as well as by a time-lapse
confocal observation of yet another biological replicate [Fig. S8(d)and (e)].
To quantitatively deal with the effect of the polar order upon
the mean radial velocity, we incorporate the polarity-induced
force f p into Eq. (1). With an = a0n cos θn and ξ = ξ0 (1 − 0 cos θn Q ),we obtain the following equation:

ξ0 (1 − 0 cos θn Q )v = ∇ · (−a0n cos θn Q ) + a0p θp n.

(2)Then, we experimentally determined θn (r, t) and θp (r, t) for ±1/2
defects, by combining end-point and time-lapse confocal observations (see Appendix B for details). We are to determine three
unknown parameters, a0n S0 /ξ0 , a0p /ξ0 , and  0 S0 , where S0 is thescalar nematic order parameter sufficiently far from defects. In
the leading term,  0 appears in the form of a product with a0n S0 /ξ0[see Supplementary Information, Eq. (S23) for example], so that
we only need to tune a0n S0 /ξ0 and a0p /ξ0 . Then, while the nematic
contribution solely could not reproduce the experimental result
as we described above, we found, remarkably, that the addition of
the polar contribution a0p /ξ0 strengthened the influx toward −1/2
defects significantly [Fig. 4(g) solid curves]. In particular, we wereable to find such values of a0n S0 /ξ0 and a0p /ξ0 that satisfactorily−
reproduced the experimental data of both v̄+r (r) and v̄r (r) simultaneously (see the experimental methods discussed in Appendix
A). We tested this for all of the four pairs of the phase-contrast
data sets (#1 or #2; see Table A1) and the end-point confocal data
sets (#1 or #2) and found agreement of similar quality (Fig. S9).
This demonstrates that the three-dimensional tilting and resulting polar order were the keys to understand the unusual
influx toward −1/2 defects we observed in our growing nonmotile
bacterial populations.

Downloaded from https://academic.oup.com/pnasnexus/article/1/5/pgac269/6948003 by University of Tokyo Library user on 22 December 2022

Fig. 3. Results on the three-dimensional nematic tilting obtained by the end-point confocal observation. To obtain them, we first measured thethree-dimensional cell orientations by the structure tensor method and obtained the nonnegative tilt angle with respect to the xy-plane (see theexperimental methods discussed in Appendix A). We then took the ensemble average over all defects. These data were obtained by a single
measurement (uniform colony, end-point confocal #1), but the reproducibility was confirmed in Fig. S7(a) and (b) by using another biological replicate
with a different substrate and agar pad (uniform colony, end-point confocal #2). (a) Nematic tilt angle θn± (r) around ±1/2 defects. The results on thex-axis, i.e., r = (x, 0), are displayed. The error bars indicate the standard error from the ensemble averaging. The two insets illustrate the definition of
the axes for +1/2 (right inset) and −1/2 (left inset) defects. (b and c) Spatial profiles of the nematic tilt angle θn± (r) for the +1/2 defects (b) and the −1/2
defects (c). The white rods represent the nematic director field.

6 |PNAS Nexus, 2022, Vol. 1, No. 5

(a)

(b)

(d)
(g)

-1/2
(e)

(f)
(outward)(outward)+1/2(inward)(inward)Fig. 4. Results on the polar order obtained by the single end-point confocal measurement (uniform colony, end-point confocal #1) [see also Fig. S7(c) to(g) for the results of another biological replicate (uniform colony, end-point confocal #2)] and theoretical calculations of the mean radial velocity basedon our theory with polar order. (a) Illustrations of the polar tilt angle θ p and the polar director np = nθp /|θp |. By choosing the direction of the head of
the nematic director n, we can uniquely determine the sign of θ p . The polar tilt angle θ p takes a positive (negative) value if the cell end at the head
(tail) of the director n is lifted above the substrate. Note that θ p changes its sign if n is reversed, but np remains unchanged, always pointing thedirection of the upper end of the cell. The polarity-induced force f p ∝ θp n is oriented toward np (see the experimental methods discussed in Appendix
A). In the following, we use n± = (cos(±φ/2), sin(±φ/2)) for the director field around ±1/2 defects, with azimuth φ of the coordinate r. (b) The polar tilt
angle θp± measured on the +x-axis of ±1/2 defects [see Fig. 3(a) for the definition of the axis]. The ensemble average over all defects is shown. The error
bars indicate the standard error from the ensemble averaging. (c and d) θp± n± , representing the strength and direction of the polarity-induced force f ±paround the +1/2 defect (c) and the −1/2 defect (d). This direction indicates which ends of cells are lifted up on average. (e and f) θp± n±
r , which is±
proportional to the radial component of the polarity-induced force, fp,r , around the +1/2 defect (e) and the −1/2 defect (f). The negative radial
component indicates that the polarity-induced force is directed toward the defect. The black rods represent the nematic director field. The outlined
arrows illustrate the direction of the radial component of the polarity-induced force near the defect. (g) Comparison of the mean radial velocity
between the experimental data (symbols) and our theoretical curves with nematic tilting and polar order (lines) (Eq. 2). The theoretical curves (lines)
use the experimental data of θ p from the end-point confocal observation #1, and the direct experimental data of the mean radial velocity (symbols)
are from the phase-contrast observation #1 [those shown in Fig. 2(e) and (f)]; see Fig. S9 for the results using other data sets. The shaded bands indicate
the range of uncertainty, evaluated from the standard error of the experimental data of θp± (r, t), which were directly used in the theoretical evaluation.
The parameters on the nematic tilt angle were set to be θn∞ = 0.2 and θn0 = 0.25 on the basis of the time-lapse confocal observation data [Fig. S8(b); see
Appendix B]. The other parameter values were  0 = 0.25, a0n /ξ0 = 0.055 μm2 /min, and a0p /ξ0 = 0.8 μm/min (see Appendix B). The vertical dashed lines
indicate the defect core radius.

Relation to circular colonies formed fromisolated cells
Although many earlier studies have already investigated how
nonmotile bacteria construct three-dimensional structures, mostof them have focused on the process where isolated cells grow andform circular colonies (21, 23, 24, 26, 30). In this situation, it has
been reported that the in-plane stress derived from cell growth ismaximized at the center of the colony (26, 35–37), which causes
a few cells to be verticalized first, locally, near the center (21, 23,
26). This is contrasted to the case of our experiments starting fromdensely and uniformly distributed cells, in which cells were verticalized almost homogeneously and simultaneously (Video 1). We
checked if cell alignment plays any role in such circular colonies(Video 3), but detected no significant correlation between the position of the first verticalization and the strength of the local orientational order (Fig. S10; see Supplementary Information, Sec-tion IV for details). Instead, we confirmed that shorter cells tend
to be verticalized first [Fig. S10(e)], in agreement with the recent
theory based on the torque balance (26). These suggest that, insuch isolated circular colonies, the spatially nonuniform stressindeed constitutes a major contribution to the start of the threedimensional transition, as reported earlier (21, 23, 26), regardless of topological defects. Conversely, by using uniform colonies,
we reduced the effect of nonuniform stress and thereby revealed
the intriguing role of topological defects in the three-dimensionaltransition.
Note also that a previous study (23) on circular colonies reported that collisions between colonies also triggered the cell
extrusion in that case. In the case of uniform colonies fromnumerous cells we studied, groups of cells merged and filledholes to complete the formation of the two-dimensional bottom layer [Video 1 and Fig. S11(a)]. We tested the possible influence of such collisions upon our defect analyses, by exam-

Downloaded from https://academic.oup.com/pnasnexus/article/1/5/pgac269/6948003 by University of Tokyo Library user on 22 December 2022(c)

Shimaya and Takeuchi | 7

Concluding remarks
In summary, we showed the relevance of topological defects tothe three-dimensional growth of growing nonmotile E. coli populations, unveiling the emergence of polar order and resulting novel
properties endowed with this active nematic system. When cultured from densely and uniformly distributed populations, cells
started to construct the three-dimensional structure a short while
after they filled the bottom plane. Since then, the net influx toward both +1/2 and −1/2 defects appeared, which may have promoted the vertical growth of colonies. The influx toward −1/2 defects was unexpected from the existing theory of active nematics, but we revealed that this resulted from the three-dimensional
tilting of cells around defects and the polar order induced
thereby. We extended the active nematics theory to incorporate
these effects and successfully accounted for the experimental
observation.
Our results suggest the role of −1/2 defects in the formationof three-dimensional structures of nonmotile cell populations,which has been overlooked compared to that of +1/2 defects supported by many recent studies on motile cells (5, 6, 13, 16, 17,19). Although it is fair to recall that the strength of the effects
reported here, i.e., the height increase at the defects and the cell
influx thereto, is not large in our setup, it is possible that vertical growth may have been prevented by the presence of the agar.
Further investigation is needed to see whether the colony height
above the defects can grow further, by alternative methods thatcan stably measure sessile E. coli populations for a longer periodof time, and whether the orientation and topological defects inintermediate layers may also affect the colony height. Besides, it
is important to contemplate the possibility of physiological significance that topological defects may ultimately have. In Bacillus
subtilis colonies, it has been found that the roughness of the colony
surface can change the wettability of the biofilm, making it moreresistant to droplets that may contain toxic substances (38–42). It
is therefore tempting to investigate the possibility that the local
vertical growth mediated by topological defects might be involvedin such surface morphology.
Finally, the emerging polar order and the influx toward −1/2
defects reported in this work may provide a novel characterization of nonmotile but growing active matter, contrasted with thestandard active matter for self-propelled particles. As such, theseresults may also shed a new light on other cellular systems withthree-dimensional structures. In this context, it is of great importance to elucidate how the polar order is formed when cells
start to tilt. Our observations show that the direction of the polar order [Fig. 4(c) and (d)] and that of the velocity field [Fig. 2(b)and (c) arrows, typical of extensile active nematics (31)] tend to
be oriented oppositely. This suggests that the polar order may
be driven by the active stress originating from the nematic orientation. The recently reported instability of the in-plane orientation in extensile active nematics (43) may also be a hint. It isalso important to understand how the absence of motility is in-volved in this mechanism; qualitatively, we may argue that thelack of cell motility would help maintain the cell tilting. Further,elucidation of the mechanism of the polar order formation andquantitative prediction of the resulting polar angle as well as thepolarity-induced force are key tasks left for future studies, which
will also clarify the relevance of our findings to other cellular
populations.

Acknowledgments
We are grateful to S. Ramaswamy for motivating us to investigate polar order in three-dimensional orientations. We thank Y.T.
Maeda and H. Salman for sharing the plasmid DNA pZA3R-EYFP
and K. Inoue for producing the strain MG1655-pZA3R-EYFP. We
also acknowledge discussions with K. Kawaguchi, D. Nishiguchi,
M. Sano, and Y. Zushi.

Supplementary Material
Supplementary material is available at PNAS Nexus online.

Funding
This work is supported by KAKENHI from the Japan Society for
the Promotion of Science (JSPS) (No. 19H05800, 20H00128), by KAKENHI for the JSPS Fellows (No. 20J10682), and by JST, PRESTO Grant
No. JPMJPR18L6, Japan.

Authors’ Contributions
T.S. and K.A.T. designed research. T.S. performed all bacterial experiments. T.S. and K.A.T. did the modeling, wrote the codes forthe numerical evaluation, analyzed data, and performed the theoretical calculations. T.S. and K.A.T. wrote the manuscript.

PreprintsA
preprintof
thisarticleDOI:10.48550/arXiv.2106.10954.

ispublishedatData Availability
The data that support the findings of this study, as well as relevant
microscope images and scripts, have been deposited in Zenodo athttps://doi.org/10.5281/zenodo.7388260.

References
1.
2.
3.
4.
5.

6.

Flemming H, et al. 2016. Biofilms: an emergent form of bacterial
life. Nat Rev Microbiol. 14: 563–575.
Mattila-Sandholm T, Wirtanen G. 1992. Biofilm formation in theindustry: a review. Food Rev Int. 8: 573–603.
Shirtliff M, Leid JG. 2009. The role of biofilms in device-related
infections. Berlin, Heidelberg: Springer.
Allen RJ, Waclaw B. 2018. Bacterial growth: a statistical physicist’s guide. Rep Prog Phys. 82: 016601.
Peng C, Turiv T, Guo Y, Wei Q, Lavrentovich OD. 2016. Commandof active matter by topological defects and patterns. Science.
354: 882–885.
Genkin MM, Sokolov A, Lavrentovich OD, Aranson IS. 2017. Topological defects in a living nematic ensnare swimming bacteria.
Phys Rev X. 7: 011029.

Downloaded from https://academic.oup.com/pnasnexus/article/1/5/pgac269/6948003 by University of Tokyo Library user on 22 December 2022

ining whether hole filling events affected the defect formation.
Detecting the locations of the holes at t = −5 min [Fig. S11(a)]and those of the defects at t = 0 [Fig. S11(b)], we confirmedthat hole filling events did not promote the formation of defects [Fig. S11(c) and (d)]. Therefore, we conclude that our resultson the relevance of topological defects to the cell flow and thethree-dimensional growth are not significantly affected by cellcollisions that preceded the formation of the complete bottom
layer.

8 |
7.

8.

9.

11.
12.

13.

14.
15.
16.
17.

18.

19.
20.

21.
22.

23.

24.

25.
26.

27.
28.
29.

30.

You Z, Pearce DJG, Sengupta A, Giomi L. 2018. Geometry and mechanics of microdomains in growing bacterial colonies. Phys Rev
X. 8: 031065.
Yaman YI, Demir E, Vetter R, Kocabas A. 2019. Emergence of active nematics in chaining bacterial biofilms. Nat Commun. 10:
2285.
Dell’Arciprete D, et al. 2018. A growing bacterial colony in two
dimensions as an active nematic. Nat Commun. 9: 4190.
Doostmohammadi A, Thampi SP, Yeomans JM. 2016. Defectmediated morphologies in growing cell colonies. Phys Rev Lett.
117: 048102.
Sengupta A. 2020. Microbial active matter: a topological framework. Front Phys. 8: 184.
Meacock OJ, Doostmohammadi A, Foster KR, Yeomans J.M,
Durham WM. 2021. Bacteria solve the problem of crowding by
moving slowly. Nat Phys. 17: 205–210.
Copenhagen K, Alert R, Wingreen NS, Shaevitz JW. 2021. Topological defects promote layer formation in Myxococcus xanthus
colonies. Nat Phys. 17: 211–215.
Saw TB, Xi W, Ladoux B, Lim CT. 2018. Biological tissues as active
nematic liquid crystals. Adv Mater. 30: 1802579.
Doostmohammadi A, Ladoux B. 2022. Physics of liquid crystals
in cell biology. Trends Cell Biol. 32(2): 140–150.
Saw TB, et al. 2017. Topological defects in epithelia govern cell
death and extrusion. Nature. 544: 212–216.
Kawaguchi K, Kageyama R, Sano M. 2017. Topological defects
control collective dynamics in neural progenitor cell cultures.
Nature. 545: 327–331.
Duclos G, Erlenkämper C, Joanny J, Silberzan P. 2017. Topological defects in confined populations of spindle-shaped cells. Nat
Phys. 13: 58–62.
Turiv T, et al. 2020. Topology control of human fibroblast cells
monolayer by liquid crystal elastomer. Sci Adv. 6(20): eaaz6485.
Maroudas-Sacks Y, et al. 2021. Topological defects in the nematicorder of actin fibres as organization centres of Hydra morphogenesis. Nat Phys. 17: 251–259.
Su P, et al. 2012. Bacterial colony from two-dimensional division
to three-dimensional development. PLoS One. 7: e48098.
Farrell FDC, Hallatschek O, Marenduzzo D, Waclaw B. 2013. Mechanically driven growth of quasi-two-dimensional microbialcolonies. Phys Rev Lett. 111: 168101.
Grant MAA, Wacław B, Allen RJ, Cicuta P. 2014. The role of mechanical forces in the planar-to-bulk transition in growing Escherichia coli microcolonies. J R Soc Interface. 11: 20140400.
Duvernoy M, et al. 2018. Asymmetric adhesion of rod-shaped
bacteria controls microcolony morphogenesis. Nat Commun. 9:1120.
Beroz F, et al. 2018. Verticalization of bacterial biofilms. Nat Phys.
14: 954–960.
You Z, Pearce DJG, Sengupta A, Giomi L. 2019. Mono- to multilayer transition in growing bacterial colonies. Phys Rev Lett. 123:
178001.
Warren MR, et al. 2019. Spatiotemporal establishment of dense
bacterial colonies growing on hard agar. eLife. 8: e41093.
Hartmann R, et al. 2019. Emergence of three-dimensional orderand structure in growing biofilms. Nat Phys. 15: 251–256.
Takatori SC, Mandadapu KK. 2020. Motility-induced bucklingand glassy dynamics regulate three-dimensional transitions ofbacterial monolayers. arXiv:2003.05618.
Dhar J, Thai ALP, Ghoshal A, Giomi L, Sengupta A. 2022. Selfregulation of phenotypic noise synchronizes emergent organization and active transport in confluent microbial environments.
Nat Phys. 8(8): 945–951.

31. Doostmohammadi A, Ignés-Mullol J, Yeomans JM, Sagués F.
2018. Active nematics. Nat Commun. 9: 3246.
32. Doumic M, Hecht S, Peurichard D. 2020. A purely mechanicalmodel with asymmetric features for early morphogenesis ofrod-shaped bacteria micro-colony. Math Biosci Eng. 17: 6873–
6908.
33. Blair DL, Neicu T, Kudrolli A. 2003. Vortices in vibrated granularrods. Phys Rev E. 67: 031303.
34. Volfson D, Kudrolli A, Tsimring LS. 2004. Anisotropy-driven dynamics in vibrated granular rods. Phys Rev E. 70:051312.
35. Volfson D, Cookson S, Hasty J, Tsimring LS. 2008. Biomechanicalordering of dense cell populations. Proc Natl Acad Sci USA. 105:
15346–15351.
36. Boyer D, et al. 2011. Buckling instability in ordered bacterialcolonies. Phys Biol. 8: 026008.
37. Los R, et al. 2020. Defect dynamics in growing bacterial colonies.
arXiv:2003.10509.
38. Epstein AK, Pokroy B, Seminara A, Aizenberg J. 2011. Bacterial
biofilm shows persistent resistance to liquid wetting and gas
penetration. Proc Natl Acad Sci USA. 108: 995–1000.
39. Trejo M, et al. 2013. Elasticity and wrinkled morphology of Bacillus
subtilis pellicles. Proc Natl Acad Sci USA. 110: 2011–2016.
40. Werb M, et al. 2017. Surface topology affects wetting behavior of
Bacillus subtilis biofilms. NPJ Biofilms Microbiomes. 3: 11.
41. Hayta EN, Rickert CA, Lieleg O. 2021. Topography quantifications
allow for identifying the contribution of parental strains to physical properties of co-cultured biofilms. Biofilm. 3: 100044.
42. Zabiegaj D, et al. 2021. Wetting/spreading on porous media andon deformable, soluble structured substrates as a model system
for studying the effect of morphology on biofilms wetting and forassessing anti-biofilm methods. Curr Opin Colloid Interface Sci.
53: 101426.
43. Nejad MR, Yeomans JM. 2022. Active extensile stress promotes
3D director orientations and flows. Phys Rev Lett. 128: 048001.
44. Sveen JK. 2004. An introduction to MatPIV 1.6.1. Eprint no. 2, ISSN
0809-4403. Dept. of Mathematics, University of Oslo.
45. Pismen LM. 2006 Patterns and interfaces in dissipative dynamics. Berlin, Heidelberg: Springer.
46. Pismen LM. 1999. Vortices in nonlinear fields. From liquid crystals to superfluids. from non-equilibrium patterns to cosmic
strings. Oxford: Oxford University Press.

Appendix A: experimental methods
Strains, culture media, and sample setup
We used a wild-type Escherichia coli strain MG1655 and its mutant MG1655-pZA3R-EYFP that contains a plasmid pZA3R-EYFP
expressing enhanced yellow fluorescent proteins. We used LB
broth (tryptone 1 wt%, sodium chloride 1 wt%, and Yeast extract
0.5 wt%) and TB+Cm medium (tryptone 1 wt%, sodium chloride
1 wt%, and chloramphenicol 165 μg/ml). To prepare nutrient agar
pads, we added agar powder to medium, solidified it by a microwave oven, then cut it into squares of size 13 × 13 mm. For
each observation, we inoculated bacterial suspension on a coverslip and put an agar pad on the suspension. We then attached
the following on the coverslip, surrounding the agar pad, to prevent the agar from drying out [Fig. 1(a)]: a frame seal (SLF0601,
Bio-Rad), a 3D printed PLA spacer (5 mm height, hollow square,inner dimensions 14 × 14 mm and outer dimensions 22 × 22 mm),another frame seal, then a plastic cover that enclosed the inner
region. Details on the strain and the culture condition in each experiment are provided below and in Table A1. The E. coli strains weDownloaded from https://academic.oup.com/pnasnexus/article/1/5/pgac269/6948003 by University of Tokyo Library user on 22 December 2022

10.

PNAS Nexus, 2022, Vol. 1, No. 5

Shimaya and Takeuchi | 9
used did not swim at all in our experimental conditions (Videos 1,
2, and 3).

Confocal observations of uniform colonies
formed from numerous cells⎡

Analysis of confocal images
For each region, we chose the plane corresponding to the bottom
layer and measured the two-dimensional cell orientation n(R) by
the structure tensor method. The image pixel size was ≈0.18 μm.
After sharpening the images by a high-pass filter, we calculated
the structure tensor J(R) at a given pixel R = (X, Y ) by
J(R) =[X I, X I]R , [Y I, X I]R[X I, Y I]R , [Y I, Y I]R


,

(A1)

with the image intensity I(X, Y), X I ≡ I(X + 1, Y) − I(X− 1, Y), Y I ≡ I(X, Y + 1) − I(X, Y − 1), and [g, h]R ≡





σ


(X  ,Y  )∈ROIR g(X , Y )h(X , Y ) f R (X , Y ). Here, the summation istaken over a region of interest ROIR , which is a square of size ≈ 7.2 μm (40 pixels) centered at R, and fRσ (X  , Y  ) is the Gaussian(X  −X)2 +(Y  −Y )2] with σ ≈ 1.8 μm
kernel defined by fRσ (X  , Y  ) ≡ exp[−
2σ 2(10 pixels). Then the cell orientation n(R) is given by the eigenvector of J(R) associated with the smallest eigenvalue λmin (R). The
orientation n(R) can also be represented by angle ψ (R) such thatn = ±(cos ψ, sin ψ ) with −π /2 ≤ ψ < π /2.
To detect topological defects, we first calculated the nematicorder parameter by
S(R) = sin 2ψ2ROI + cos 2ψ2ROI ,
R

R

(A2)

where ·ROIR denotes the spatial average within ROIR . Then, we located the positions of local minima of S(R) as candidates of topological defect cores. For each candidate point, we calculated the1
topological charge q = 2π
C dψ, where C is a square closed path⎤[x I, x I]r y I, x I r [z I, x I]r⎢⎥
J(r) = ⎣ x I, y I r y I, y I r z I, y I r ⎦ ,[x I, z I]r y I, z I r [z I, z I]r(A3)

where x I ≡ I(x + 1, y, z) − I(x − 1, y, z), y I and z I are defined like
wise, [g, h]r ≡ r ∈ROIx ,y ,z g(r )h(r ) frσ (r ). Here, the summation isr , ,

taken over a three-dimensional region of interest ROIrx y z , which
is a cuboid of size x × y × z centered at r, with x = y ≈4.3 μm (24 pixels) and z ≈ 3.8 μm (24 pixels). The Gaussian ker(r −r)2
nel frσ (r ) is defined by frσ (r ) ≡ exp[− 2σ 2 ] with σ = 2.2 μm. Then,the three-dimensional cell orientation n3 (r) is given by the eigenvector of J(r) associated with the smallest eigenvalue. The orientation n3 (r) is then represented by angles ψ (r) and θ (r) such thatn3 = (cos θ cos ψ, cos θ sin ψ, sin θ ) with 0 ≤ ψ < 2π and −π /2 ≤ θ <π /2. As is clear from the definition, the angle ψ (r) specifies the twodimensional cell orientation n(r) by n = (cos ψ, sin ψ ) and θ (r) indicates the angle between the three-dimensional orientation andthe xy-plane. Note that n3 (r) and −n3 (r) are equivalent, so that thesign of n(r) and θ (r) can be changed simultaneously.
To investigate statistical properties of the cell tilt angle around±1/2 topological defects, we need to define tilt angles whose sign
can be determined unambiguously. The simplest choice is to takethe ensemble average of |θ (r)|, which can be used to detect thepresence of the three-dimensional tilting. We took this averageover isolated defects of each sign, separated by a distance longer
than 9 μm from the nearest defect, and this defines our θn± (r). To
characterize the polar order, we need an angle that can take both
positive and negative values. To define such an angle, we use the(r) ≡ (cos(±φ/2), sin(±φ/2)) around ±1/2 defects,
director field n±
refDownloaded from https://academic.oup.com/pnasnexus/article/1/5/pgac269/6948003 by University of Tokyo Library user on 22 December 2022

We used the mutant strain MG1655-pZA3R-EYFP that expresses
enhanced yellow fluorescent proteins. Before the observations,
we inoculated the strain from a glycerol stock into 2 ml TB+Cm
medium in a test tube. After shaking it overnight at 37◦ C, we transferred 20 μl of the incubated suspension to 2 ml fresh TB+Cm
medium and cultured it until the optical density (OD) at 600 nm
wavelength reached 0.1–0.5. The bacterial suspension was finallyconcentrated to OD = 5 by a centrifuge, and 1 μl of the suspension
was inoculated between the coverslip and the agar pad (1.5 wt%agar).
The sample was placed on the microscope stage, in a stage-top
incubator maintained at 37◦ C. The microscope we used was Leica SP8, equipped with a 63× (N.A. 1.40) oil immersion objective
and operated by Leica LasX. The data shown in Figs. 1, 3, 4, S1, andS2 were obtained by a single end-point observation, in which we
cultured the colonies without excitation light until 14 hours after
the cells had filled the observation area. We also show data obtained by another biological replicate with a different substrate
and agar pad in Fig. S7. For each set of these data, we captured
three-dimensional images of size 184.52 × 184.52 × 16 μm from
20 separate regions. The optical resolution, as evaluated by theformula of the point-spread function, was about 140 nm in thehorizontal plane and 250 nm in the vertical direction. The confocal pinhole size was 0.21 Airy unit. For the data shown in Fig. S8,
we carried out a single time-lapse observation and obtained images of size 184.52 × 184.52 × 6.4 μm from four separate regionswith the time interval 15 min. The image pixel size was ≈0.18 μm
in the xy-plane and 0.16 μm along the z-axis.

with a side of about 3.6 μm (20 pixels) centered at the candidate
point. The candidate point is regarded as a topological defect if q
= ±1/2, and dismissed otherwise. To determine the angle of thearm of each defect [Fig. 1(d) inset], we used the profile of |ψ − φ|on C, where φ is the azimuth with respect to the defect core. A
single minimum of |ψ − φ| exists for each +1/2 defect, while there
are three local minima for each −1/2 defect. Each minimum point
corresponds to an arm of the defect. Blue trefoils indicating −1/2
defects in Fig. 1(c) were drawn by setting one of the arms of the trefoil at the angle of the global minimum, with the other two arms
added by rotating the first arm by 120◦ . We thereby obtained thetwo-dimensional locations of all defects and their signs.
To investigate the dependence of the colony height on topological defects, we picked up hundreds of isolated defects, separated
by a distance longer than 9 μm from the nearest defect. For comparison, we also randomly selected 1000 points, which are separated more than 9 μm from defects. For a given position in thexy-plane, we obtained the image intensity profile along the z-axis,
with the interval of z-slices being 0.16 μm. The height was then determined by the length of the region whose intensity was higher
than 20% of the maximum intensity in this profile.
The three-dimensional tilting of the cells around defects wascharacterized as follows. First, for each defect, we rotated the confocal image horizontally so that the defect arm was orientated in
the positive direction of the x-axis. For −1/2 defects, we did this rotation for each of their three arms and obtained a set of three images from each defect. Then, for each rotated confocal image I(r),where r is the coordinate relative to the defect, we obtained thethree-dimensional cell orientation n3 (r) by the three-dimensionalversion of the structure tensor method. For each pixel r = (x, y, z),which was chosen from the plane corresponding to the bottom
layer in each region, we calculated the three-dimensional structure tensor:

10 |PNAS Nexus, 2022, Vol. 1, No. 5

Table A1. List of experimental measurements in this study.
Initial celldensity

Strainuniform colony, end-point confocal #1
uniform colony, end-point confocal #2
uniform colony, phase contrast #1
uniform colony, phase contrast #2
uniform colony, time-lapse confocal #1
circular colony, phase contrast #1
circular colony, phase contrast #2

a mutant MG1655-pZA3R-EYFP
a mutant MG1655-pZA3R-EYFP
a wile-type MG1655
a wile-type MG1655
a mutant MG1655-pZA3R-EYFP
a wile-type MG1655
a wile-type MG1655

with the azimuth φ of the position r in the xy-plane, and took the(r) · n(r)] over isolated defects (withaverage of the field θ (r)sign[n±refthe same criterion on the distance from other defects). This is ourθp± (r) which characterized the polar order. The polarity-induced±±
force is then f ±p (r) ∝ θp (r)nref (r). This right-hand side is shown in
Fig. 4(c) and (d), and its radial component in Fig. 4(e) and (f).

Phase-contrast observation of uniform colonies
formed from numerous cells
We used the wild-type strain MG1655. Before the time-lapse observation, we inoculated the strain from a glycerol stock into 2 ml LB
broth in a test tube. After shaking it overnight at 37◦ C, we transferred 20 μl of the incubated suspension to 2 ml fresh LB brothand cultured it until OD at 600 nm wavelength reached 0.1–0.3.
The bacterial suspension was finally concentrated to OD = 5 bya centrifuge, and 1 μl of the suspension was inoculated between
the coverslip and the LB agar pad (1.5 wt% agar).
The sample was placed on the microscope stage, in an incubation box maintained at 37◦ C. The microscope we used was Leica
DMi8, equipped with a 63× (N.A. 1.30) oil immersion objective anda CCD camera (Leica DFC3000G), and operated by Leica LasX. The
image pixel size was ≈0.17 μm. For the data shown in Figs. 2, S3,
and S10, we used a single substrate and carried out a time-lapse
observation with the time interval 1 min for 30 separate regions ofdimensions 110.03 × 81.97 μm. We also obtained a biological replicate using another substrate for the data shown in Fig. S5. For eachregion, we determined the frame at t = 0, i.e., the frame in whichcells filled the observation area for the first time. We then measured the cell orientation n(R) and detected topological defects
in all frames, by the method described below. We used isolated
topological defects only, each separated by a distance longer than
9.5 μm from the nearest defect. As a result, we obtained hundreds
of defects for each time.

Phase-contrast observation of circular colonies
formed from a few cells
We used the wild-type strain MG1655. We cultured bacteria in thesame way as for the observation of uniform colonies. The bacterial suspension was finally diluted to OD = 0.01, and 1 μl of thesuspension was inoculated between the coverslip and the LB agar
pad (2.0 wt% agar).
The imaging process and the condition during the observation
were the same as those for the observation of uniform colonies.
We carried out time-lapse observations with the time interval
1 min for 30 isolated colonies, which started to form from a few
cells. We repeated the experiments twice using different substrates and acquired data from 60 colonies in total. From eachcolony, we chose the frame right before the first extrusion ofa cell from the bottom layer took place. We used 60 such im-

high
high
high
high
high
low
low

Data
Figs. 1, 3, 4, S1, and S2
Fig. S7
Fig. 2, S3, and S10
Fig. S5
Fig. S8
Fig. S10
Fig. S10

ages from the 60 colonies for analysis. For each colony, we binarized the image, and obtained the area A by the total numberof pixels, the center position by the center of mass, and the radius Rmax by π R2max = A, using the regionprops function of MATLAB. The first extruded cell was detected manually, by using a
black spot that a tilted cell exhibits in the phase-contrast image(see Video 3). We manually labeled pixels contained in each extruded cell, and obtained the position as well as the mean and thestandard deviation of the coherency over the labeled pixels (seethe section “Analysis of phase-contrast images” for the method
to evaluate the coherency). To obtain the spatial dependence ofthe coherency shown by boxplots in Fig. S10(b), we divided thespace into regions bordered by concentric circles, with the radiithat increased by R/Rmax = 0.1. The length of cells was evaluatedmanually from the major axis of each cell by using a painting
software.

Analysis of phase-contrast images
Using phase-contrast images from uniform and circular colonies,
we measured the two-dimensional cell orientation n(R) and detected topological defects, in the same manner as those for confocal observations. The image pixel size was ≈0.17 μm. The structure tensor was calculated with the ROI size  ≈ 6.8 μm (40 pixels)
and the characteristic length of the Gaussian filter, σ ≈ 1.7 μm (10pixels). The detection of topological defects was carried out withthe closed path C with a side of about 3.4 μm (20 pixels), as inFig. 2(a) and Video 2.
In addition to the cell orientation n(R), we also obtained thecoherency parameter C(R) defined by

C(R) =

λmax (R) − λmin (R),λmax (R) + λmin (R)(A4)

with the largest eigenvalue λmax (R). This quantifies the degree ofthe local nematic order.
For uniform colonies, we also measured the velocity field of
the cells around the detected defects, by particle image velocimetry (PIV). For this, we used MatPIV (44) (open source PIV toolbox for MATLAB), with the PIV window set to be a square of
size ≈2.7 μm (16 pixels). To take averages over defects, for eachdefect, we rotated the image so that the defect arm was oriented in the positive direction of the x-axis. For −1/2 defects, we
did this rotation for each of their three arms, and all of the resulting velocity fields were used for the ensemble average. We
thereby obtained the ensemble-averaged velocity field v(r, t), asa function of the coordinate r = (x, y) relative to the defect, andtime t.
The divergence of v(r) = (u(r), v(r)) was calculated as follows(here, we omit t from the argument for simplicity). We first obu(x+1,y)−u(x−1,y)v(x,y+1)−v(x,y−1)+
with the pixel size δ ≈tained D(r) =2δ
2δ

Downloaded from https://academic.oup.com/pnasnexus/article/1/5/pgac269/6948003 by University of Tokyo Library user on 22 December 2022

Measurement

Shimaya and Takeuchi | 11
0.17 μm. We then calculated the divergence field by



σ 

(x ,y )∈ROIr D(r ) f r (x , y ),σ (x ,y )∈ROIr f r (x , y )(∇ · v)(r) =(A5)where ROIr and the Gaussian kernel frσ (x , y ) were defined asabove, but with  ≈ 2.7 μm (16 pixels) and σ ≈ 0.68 μm (5 pixels).

To theoretically account for the experimental result of the mean
radial velocity v̄±r (r), in particular the influx toward −1/2 defects shown in Fig. 2(e) and (f), we solved the force balance equations Eqs. (1) and (2). While detailed descriptions on the solutions
are given in Supplementary Material, here, we outline the theoretical assumptions and the methods to obtain the theoretical results shown in Fig. 4(g), which satisfactorily reproduced the experimental data when the influx toward −1/2 defects was strongest.
First, we assume the director field winding uniformly arounda +1/2 or −1/2 defect, n± (r, φ) = (cos(±φ/2), sin(±φ/2)), where (r,φ) is the two-dimensional polar coordinate, centered at the defectcore. The nematic order tensor Q ± (r, φ) is then given by
Q ± (r, φ) = S(r)cos(±φ) sin(±φ),
sin(±φ) − cos(±φ)(B1)

with the scalar nematic order parameter S(r) left as a free parameter. Based on the assumption that Q ± minimizes the nematic freeenergy, S(r) can be theoretically expressed by the following Padé
approximant (13, 45, 46):
S(r) = S0 F(r/rS ),F(x) ≈ x

0.34 + 0.07x2,
1 + 0.41x2 + 0.07x4(B2)

with the defect core radius rS and S0 = S(∞). To determine thevalue of rS , we fitted Eq. (B2) to the experimental data of the coherency C(r) [Fig. S6(a) and (b)] and obtained rS = 1.2 μm. Note that,
because the angle field ψ (r) does not contain information of thedefect core, the nematic order parameter evaluated by Eq. (A2) is
not suitable for estimating rS . Concerning S0 , it always appears as
a product with either  or an , so that we fix S0 = 1 without loss ofgenerality.
The case without three-dimensional cell tilting, described by
Eq. (1), was already dealt with by earlier studies (13, 17). Since
Eq. (1) is linear, we can readily solve it and obtain, for the mean
radial velocity:v̄±r (r) = −anS (r) ± S(r)/r
S(r)
.
ξ0
1 −  2 S(r)2(B3)

Then, we can show, with Eq. (B2), that it is negative for +1/2 defects
and positive for −1/2 defects, for all r > 0 (see Supplementary Ma-

v± (r, φ) = ξ0−1 (1 − 0 cos θn± (r, φ)Q ± )−1

× ∇ · (−a0n cos θn± (r, φ)Q ± ) + a0p θp± (r, φ)n± .

(B4)Regarding the first term that describes the contribution by
nonuniform nematic tilting, we determined θn± (r, φ) by time-lapse
and end-point confocal observations. Because we could not obtainclear spatial profile of θn± (r, φ) from the time-lapse observation due
to photobleaching, we used high-quality, end-point confocal images to determine the spatial profile, then calibrated its amplitude
by the time-lapse observation to account for the time period of interest. First, on the spatial profile, our end-point confocal observation (Fig. 3) suggests that θn± (r, φ) = θn∞ + (θn0 − θn∞ ) exp(−r2 /rθ2 ) withconstants θn∞ , θn0 , rθ , regardless of φ and the sign of the defect.
From the spatial profile, we obtained rθ = 1 μm. For the peakheight, we used time-lapse observations for 200 min ≤ t ≤ 250 min,during which the influx toward −1/2 defects was strongest
for this strain (Fig. S8a), and estimated θn∞ = 0.2 and θn0 = 0.25
[Fig. S8(b)].
To see the influence of the nematic tilting, we numerically∞0
calculated v̄±r (r) with θn = 0.3 and θn = 0.75, which were estimated from the end-point confocal observation, without polar order [Fig. 2(h), the solid lines]. The other parameters were  0 = 0.25,
an /ξ 0 = 0.055 μm2 /min, and rS = 1.2 μm. The strength of the influx
toward −1/2 defects obtained thereby was smaller than the experimental result, indicating that the nematic tilting is insufficnent
to quantitatively explain the influx toward −1/2 defects.
For the polar contribution to Eq. (B4), we determined thespatial structure of θp± (r, φ) by the end-point confocal observation [Fig. 4(c) and (d)]. Then, we calibrated the amplitude by multiplying the ratio of θp± 0<x<10 μm,y=0 fromthe time-lapse observation for 200 min ≤ t ≤ 250 min[Fig. S8(c)] to that from the end-point observation [Fig. 4(c)and (d)].
We are finally left to determine the following parameters:  0 ,a0n /ξ0 , and a0p /ξ0 . First, we note that, in the leading term,  0 appears
in the form of a product with a0n S0 /ξ0 (see Supplementary Information, Eq. (S23) for example), so that we chose  0 = 0.25. Then,
we tuned a0n /ξ0 and a0p /ξ0 to reproduce the experimental data ofv̄±r (r). When the end-point confocal observation data #1 were used,
we obtained a0n /ξ0 = 0.055 μm2 /min and a0p /ξ0 = 0.8 μm/min, withthe results shown in Fig. 4(g) [see also Fig. S9(a) to (c)]. When theend-point confocal observation data #2 were used instead, we obtained a0n /ξ0 = 0.055 μm2 /min and a0p /ξ0 = 1.4 μm/min, with theresults shown in Fig. S9(b) to (d).

Downloaded from https://academic.oup.com/pnasnexus/article/1/5/pgac269/6948003 by University of Tokyo Library user on 22 December 2022

Appendix B: theoretical calculations

terial). In Fig. 2(h), by the dotted lines, we showed v̄±r (r) for  = 0.25,
an /ξ 0 = 0.055 μm2 /min, rS = 1.2 μm.
In fact, even in the presence of three-dimensional cell tilting
and polar order, i.e., in the case of Eq. (2), it is linear in v and thesolution for the case of ±1/2 defects is given by