Biomechanical modeling of spatiotemporal bacteria-phage competition

To investigate various aspects of bacteria-phage competition, we performed four major computational experiments. Experiment #1 examines the peripheral collision between a growing bacterial colony and a single phage, a common scenario observed in Petri dish assays involving colonies encountering phages⁴⁹. Experiment #2 models predation and growth within a monolayer colony containing 20 infected cells near its center. It highlights the organized flow of cells toward the phage-affected region while others evade predation. Experiment #3 investigates the influence of key phage parameters (latent period, burst size, and adsorption coefficient) on cell movement into phage-affected regions. Finally, Experiment #4 focuses on the lysis event in a monolayer colony that triggers 3D growth. The resulting 3D phage-plaque exhibits a less symmetric, more disorganized void region, accelerating the colony’s destruction. In addition, Experiment #4 shows fundamental differences between well-mixed models for phage-cell dynamics and our spatially structured simulator. In these simulations, the same initial monolayer colony develops in 3D, surviving (or not) the attack of lytic phages originating exclusively from different relative positions.

Computational experiment #1. A growing colony meets a lytic phage

We examined interactions between E. coli-like cells and lytic phages, focusing on how variations in initial conditions can affect the morphology of phage-induced plaques. The model is described in Methods, and the parameters of phages are summarized in Table 1. Comparison with the well-mixed case is provided in Supplementary Note 1 and Supplementary Fig. 1. By controlling the initial placement of cells and phages, we can observe distinct patterns in bacteria-plaque co-development. Figure 1 illustrates a bacterial monolayer encountering a lytic phage. Over five hours of continuous cell growth and phage-induced lysis, a noticeable wedge-shaped sector emerges within the colony. This sectoral formation is facilitated by the specific alignment of the sensitive E. coli cells, turning the main axis of the rod-shaped cells perpendicular to the wedge’s border.

After 15 h of growth, the colony encounters a phage that infects a single cell near its border. As the infection spreads, lysed cells create a wedge-shaped (missing) sector within the colony, which remains consistent throughout the simulation. The growing colony actively pushes cells into the phage-affected region, effectively containing the infection within a localized area. This dynamic alignment and compression ensure a steady supply of cells to sustain the phage region while allowing cells in unaffected areas to continue their normal growth and life cycles. This observation explains most of the shapes reported in the Petri dishes in ref. ⁴⁹.

Computational experiment #2. Bi-directional flow induced by phage inside a colony

In a typical numerical experiment, the initial setup consists of 1500 E. coli cells forming a confluent monolayer (bacterial lawn). To simulate infection dynamics, 20 of these cells, located near the colony’s center, are initially infected. This initial configuration will be consistently applied across all our subsequent simulations to maintain uniformity and comparability in our results. However, our code is versatile and allows arbitrary initial conditions. This flexibility enables the exploration of a wide range of infection scenarios and colony morphologies beyond the specific case presented here. For this computational study, we employed the T5-like phage; see Table 1. This phage has the longest latent period compared to other phages in this family (i.e., the least virulent T-phage). In our simulation parameters, the latent period is set to 44 minutes. Following lysis, each infected cell vacates the space and releases approximately 290 new phage particles. Additionally, we set the adsorption coefficient, which represents the rate of phage attachment to host cells, to 2 × 10⁻¹⁰ (ml/min), reflecting the probability of successful attachment over time.

Proliferating bacteria exert forces and torques on each other. Within the colony, the pressure field exhibits minimal values near the outer border and reaches its maximal values closer to the colony’s center; see Fig. 2. It creates a high-pressure core and lower-pressure periphery, influencing the physical dynamics of cellular movement and colony expansion. For simplicity, we did not include the dependence of bacterial growth rate on the pressure⁵⁰. Figure 2 corroborates that when phages are introduced near the colony center, they shift this high-pressure location due to the localized damage they inflict by lysing cells. As the infection progresses, with phages killing more cells, the colony’s pressure field tends to equilibrate across annular regions that encircle the initial infection site. This outward shift in the pressure field’s maximum region generates a dual influence on the surrounding cells. First, it compels cells to move inward, effectively filling the void left by the phage-affected regions. Second, it promotes an outward expansion of the colony, extending its radial spread. Together, these effects create a dynamic response in which the colony compensates for the phage-induced voids by inward cell migration while simultaneously pushing the colony’s outer boundary outward. This process results in a complex interaction between colony expansion and infection dynamics, as pressure distribution drives both containment of the phage impact and growth in unaffected areas.

The pressure on any cell p_i is determined by evaluating the external work induced by interactions with neighboring cells and the substrate:

where F_ij is the force exerted on the i-th cell due to the j-th cell (or substrate), δ(r_ij) represents the overlap between the cells (or substrate), and v_i denotes the volume of the i-th cell. Figure 3 illustrates the angle-averaged pressure as a function of the radial coordinate at two instances (17.5 and 20 hours). Without phage, maximal pressure is consistently observed near the colony’s center, driving purely radial spreading via a negative pressure gradient that induces outward flow. A simple estimate for a confluent monolayer colony follows from the momentum balance ηv = − ∇ p and the mass balance ∇ ⋅ v = λρ³⁷. Here \(\rho ={{{\rm{const}}}}\) is the bacterial number density, λ is the growth rate, and η is the effective viscosity. For pure radial growth, we readily obtain linear radial velocity and parabolic pressure profiles

where R₀ is the colony’s external radius. In contrast, the phage predation shifts the location of maximum pressure, resulting in bi-directional flow dynamics: outward movement in regions with negative pressure gradients and inward migration toward phage-affected zones in regions with positive gradients.

Spatial organization of growing colonies of bacterial and eukaryotes cells is a rapidly expanding research field, see., e.g., refs. ^{36,37,38,51,52,53,54,55,56,57}. Mature colonies exhibit topological defects altering the organization and motion of cells. The topological defects arise due to the collective motion of cells. Since the bacterial cells under study are apolar (symmetric to reflection by 180 degrees), the topological defects are  ± 1/2 point disclinations (comet-like for 1/2 type and triangular for −1/2 type), similar to that in nematic liquid crystals⁵⁸. While it is unclear how the defects affect long-term colony growth, they may influence local phage concentration and bacteria flow around them. We compare the evolution of two microcolonies, (i) the no phage case and (ii) the center phage infection; Figure 4. While the density of defects fluctuates, it remains roughly constant in the process of colony growth, with or without phages. However, the effects of phages’ infection result in large-scale phage transport through the growing colony beyond the local motion of phages. As the phage-plaque grows, the surviving cells align to the radial direction, maximizing the cells’ probability of encountering phages. Increasing the radial order decreases the local number of topological defects since the defects close to the hole in the colony center are “flushed” by the bacterial flow.

Computational experiment # 3. Assessing the conditions to observe the bi-directional flow induced by phages inside a colony

To explore the robustness of the observed cell inflow towards the phage affected region, we perturbed key parameters of the T5 phage, like the latent period, burst size, and adsorption coefficient, allowing uniform variations towards the default value of 50%, for each phage parameter to observe variations in phage-plaque evolution while maintaining the same initial conditions; see Supplementary Fig. 1. We analyzed the corresponding probability density functions (PDFs) linked to each parameter, highlighting our two extreme “virulent” cases. To derive general conclusions, we perturbed each phage parameter and generated 20 samples, using each as input for our simulator to produce computational plaque assay solutions. Due to the inherent non-linearities in our model (some explicitly defined), first-order methods like a Taylor expansion could be insufficient for a reliable “input-output” link. Instead, we adopted a conservative approach, directly correlating the model outputs with the sampled parameters^59,60, see Supplementary Fig. 2, offering a broader perspective than traditional local sensitivity analysis.

In the following, we explored the phage effect on the radial order within the colony. We define the in-plane radial order parameter S_r in colonies with and without cells’ predation by phages:

Here, the order parameter is evaluated in the region of interest (ROI)_i, N_i is the number of bacteria in the ROI, φ_j is the orientation of bacteria at the radial vector r_j and polar angle θ_j with respect to the center of the colony. The radial order parameter⁵³(S_r) quantifies the misalignment of any cell with the radial direction. Ranging from  −1 (tangential orientation, φ_j − θ_j = ± π/2) to 1 (aligned to the radial direction, φ_j − θ_j = ± π).

Figure 5 a presents a comparative analysis of the radial order among surviving cells, highlighting colonies that generate monolayer structures with approximately centered void regions, similar to those in Fig. 2. Without phage, the colony does not establish any preferential alignment, with the radial order fluctuating around 0. In some instances, it achieves negative values (tangential alignment). When phages are introduced, the surviving cells increase radial order as the phage-plaque spreads. In the most virulent case, the radial order reaches a peak of 15% after 20 h, indicating a mild but consistent radial alignment. In contrast, the colony remains disordered under less virulent phage conditions, with cells exhibiting random alignment. An increase in the radial order corresponds to a reduction in topological defect density. Figure 5b illustrates how topological defects (point singularities of the orientational order) emerge in regions where the colony displays significant radial disorder. When cell orientation is inconsistent with a radial alignment, the structural disruptions lead to increased occurrences of these defects. To quantify this, we assessed the orientation of each cell, generating a coarse-grained orientation field over the colony. Within this field, we applied the criteria established in ref. ⁵⁶ to identify singular points in the orientation field that correspond to topological defects. This approach allowed us to systematically evaluate how phage-induced radial alignment impacts the defect landscape across the colony, revealing that a more ordered radial configuration effectively minimizes defect density. Comparative analysis of the density of topological defects within the simulated colonies, Fig. 5a and b highlight differences between colonies with and without phage presence. As observed in prior simulations, the introduction of phages reduces the radial spreading speed of the colony and encourages a reorganization of cell orientation that favors radial alignment. This shift in cellular alignment correlates with a reduction in topological defects.

Figure 6a and c describe the average radial order 〈S_r〉 of all the simulated colonies at 17.5 hours (a) and 20 hours (c), respectively. As expected, the radial alignment increases thanks to the void space generated by the lysed cells. Figure 6b and d illustrate the variability of the averaged radial velocity 〈V_r〉 of several monolayer colonies facing phages with different parameters at two different instants, 17.5 h (b) and 20 h (d). As a result of the phage-induced bacteria killing, the surviving cells increase the radial order. This alignment effect is stronger as the phage-plaque matures. In the absence of phages, the colony shows no consistent order at 17.5 h (as the radial coordinate increases, the violet curve has tangential order, followed by radial alignment and developing tangential order at the colony’s border see Figs. 2 and 7), or no relevant order at 20 h (here the violet curve fluctuates with small amplitude around 0 indicating disorder). In the most virulent cases, the radial order stabilizes around 15% (see Fig. 6a and c), whereas the least virulent behave as the no phage case.

Without phages, as the colony grows larger, its expansion rate also increases proportionally. This trend suggests a uniform outward progression as the colony cells proliferate and extend radially (see Fig. 6b and d). When phages are introduced into the simulation, this outward radial velocity is reduced. The presence of phages inhibits the natural expansion of the colony by targeting and lysing cells within infected regions. This reduction in averaged radial speed reflects the colony’s struggle to grow outward while battling phage-induced cell death. In the most virulent cases, the averaged radial velocity of the colony becomes negative, indicating an inward movement of cells (see Fig. 6b and d, most virulent case). This negative velocity occurs when the colony shifts its growth dynamics to supply cells toward the infected regions rather than expanding outward. In these high-virulence scenarios, the colony’s cells are mobilized to “feed” the phage-affected zone, filling in areas of cell loss due to phage lysis. At the same time, the flow of bacteria toward the infected region could be interpreted as a colony “defense” against the infection: the inflow pushes infected bacteria to reopen the hole. It reduces the infection probability in the bulk of the colony.

Computational experiment # 4. Bacteria-phage competition in three-dimensional colonies

Phage-induced cell removal patterns differ significantly in 3D colonies compared to 2D monolayers. In a phage-free 3D colony, cells display a more random alignment due to the increased spatial freedom in three dimensions. It contrasts with the relatively uniform, planar arrangement observed in 2D monolayers. This random bacteria orientation in 3D structures contributes to a more complex, less predictable response to phage attack, Fig. 7, at 20 h. When phages are introduced to a 3D colony, the infection begins near the center and disrupts the typical radial symmetry seen in 2D cases, Fig. 2. After five hours of continuous cell predation, the affected region expands outward, eventually reaching the colony’s border. At this stage, the extensive lysis causes structural disconnections within the colony, fragmenting it into isolated cell clusters. These new cell clusters form as the phage infection erodes the colony’s integrity, severing intercellular connections and producing spatially distinct groupings of surviving cells. Figure 7 shows a damaged colony characterized by disrupted continuity and newly formed clusters that may be more susceptible to environmental stresses and further phage predation.

Obtaining the best killing efficiency with the same amount of phages and cells, changing the initial location

To corroborate the dependence between cell lysis and initial phages’ placement, we fixed the number of cells and phages (5000 cells and 500 T7 phages), changing the phage-colony relative localization. We analyze three different cases (i) placing the phages in the center of the microcolony, (ii) placing the phages within the border of the microcolony, and (iii) placing the phages randomly on top of the microcolony. Figure 8 displays different phage-plaque shapes. On surfaces, the initial phage/cell placement plays the most relevant role in the colony’s predation. In particular, Fig. 8e and f maintain circular shapes, being Fig. 8e the no phage case. Figure 8g and h show severe cases of cell predation, notoriously altering the shape of the microcolonies. Such structural differences can be recognized in the population plots depicted in Fig. 8i and J. Randomly placing phages maximizes cell killing and phages’ production, returning compatible population trends as well-mixed models. Similar results were reported in the experimental work of Vidakovic et al.⁵¹, where the mixing of phages and E. coli cells removes all spatial dependence. Placing phages on the colony’s border and colony’s center illustrates the resilience that the colony gains from the cells’ collective behavior. Placing phages on the colony center is the less effective technique to eradicate cells (see Fig. 8i). Eventually, a few cells escape (by pressure due to cell division) the phage zone and create new microcolonies. Figure 8j presents center and border infections as efficient ways to produce phages continuously.