Analysis of a model for bacteriophage infections and bacteria defense: a synergetics perspective

A model for bacteriophage infections and bacteria defense is analyzed using the concepts of synergetics. The model order parameter is determined and the corresponding amplitude equations are derived. Within this framework it is shown how the order parameter defines a multi-species building block that captures the organization of infection outbreaks and the initial defense reaction and how the order parameter amplitude determines the corresponding temporal characteristics. Two approximative models with different domains of application are derived as well. In doing so, a supplementary perspective of bacteriophage infections that provides insights beyond the classical state space perspective is provided.

1 Introduction

In general, epidemiological systems are given by complex networks of interacting populations of different species (Pastor-Satorras et al., 1998). While the isolated populations typically exhibit a relatively simple dynamics, a challenge in the field of network physiology is to understand how the interactions between the different types of species shape the overall network dynamics. In this context, an important first step is to consider mean field approximations in terms of ODE and coupled differential equation models (Pastor-Satorras et al., 1998; Granger et al., 2024), which due to their relative simplicity frequently allow for analytical solution methods. For virus infections the ODE three-species TIV model captures the basic dynamics of interacting target cells, infected cells, and virus particles (Nowak and May 2000). Likewise, in the context of bacteriophage infections, we are dealing with target bacteria, infected bacteria, and bacteriophages–the latter act as viruses. Studying bacterial infections and the role of bacteriophages is an important task (Bloch and Wegrzyn, 2024; Geng et al., 2022; Zborowsky et al., 2025) and is an indispensable step when attempting to use phages in modern medicine to cure certain diseases in humans (Li et al., 2021; Zborowsky et al., 2025). To this end, both simplified and generalized three-species models have been studied in the literature (Li et al., 2021; Weitz et al., 2005; Zborowsky et al., 2025). In particular, as part of their comprehensive study, Skanata and Kussell (2021) considered bacteria that can exist as resistant and non-resistant phenotypes with respect to a given invading phage. Increasing the concentration of resistant phenotypes is a defense mechanism against phage attack because this mechanism decreases the effective contact rate between phage and non-resistant bacteria such that under appropriate conditions the infection dies out. While the dynamics systems perspective, in general, is an indispensable tool to analyze bacteriophage infection models, little attention has been paid to utilize the more specific dynamical systems concepts of synergetics (Haken, 2004; Hutt and Haken, 2020; Uhl, 1999; Wunner and Pelster, 2016) in this regard. However, in the wake of the COVID-19 pandemic it has been shown that synergetics can supplement existing dynamic systems approaches to understand epidemiological and virus dynamic models (Frank, 2022). In this Brief Report a parsimonious four-species model for bacteriophage infection and bacteria defense will be considered that involves resistant phenotypes as in Skanata and Kussell (2021). The relative simplicity of the model will allow for an analytical approach. The aim of the study is to derive explicitly the order parameter of the model and to show how it determines the initial organization of a phage attack and the corresponding defense reaction. Moreover, the aim is to derive the amplitude equations that determine the evolution of the system along the order parameter and the remaining directions. Two approximative models for the system dynamics in this context will be derived as well. An exemplary simulation will illustrate some of the analytical results.

2 Materials and methods

As reviewed above, our starting point is the three-species model that involves susceptible $(S)$ and infected $(I)$ bacteria and phage load $(P)$. In line with the TIV model of virus dynamics (Frank, 2022) the model equations read

$\frac{\text{d}}{\text{d}t}S=-k_0\left(P\right)S+\mu S\left(1-\frac{S}{K}\right),\frac{\text{d}}{\text{d}t}I=k_0\left(P\right)S-k_{1}I,\frac{\text{d}}{\text{d}t}P=qI-k_{2}P,(1)$

with $k_1,k_2,q>0$, where $t$ denotes time, $k_0$ describes the $P$-dependent transition rate of $S\to I$ transitions, $k_1$ and $k_2$ denote the decay rates of infected cells and phages, respectively, and $q$ describes the production rate of phages per infected bacteria. The evolution equation for $S$ involves a logistic growth term with the growth rate $\mu$ and the capacity $K>0$. Below we will consider both the case $\mu =0$ when the growth term can be neglected and the more general case $\mu >0$. By comparing these two cases, we will see that the bacterial growth dynamics actually has no effect on the initial outbreaks dynamics captured by the order parameter. Therefore, the growth term may be neglected when (i) changes in $S$ are primarily due to $S\to I$ transitions capture by the $k_{0}S$-term or (ii) the focus is on the initial phase of the outbreak dynamics. The transition rate $k_0(P)$ in Equation 1 depends on the infecting species like $k_0={\beta }_{0}P$, where ${\beta }_0>0$ denotes the effective contact rate (Frank, 2022). In order to take the active defense mechanism mentioned in the introduction into account, the model (1) was modified in two ways. First, it was assumed that when infected bacteria $I$ emerge in the bacteria population then resistant bacteria mutations $(R)$ are grown like

$\frac{\text{d}}{\text{d}t}R=\alpha I\left(1-\frac{R}{R_m}\right),(2)$

where $\alpha \ge 0$ and $R_m>0$ denote the growth rate and the maximal concentration of resistant bacteria, respectively. From a mechanistic point of view, Equation 2 captures the adaptive defense mechanism of bacteria via the so-called CRISPR system (Abedon, 2012; Skanata and Kussell, 2021). The CRISPR system allows bacteria to memorize attacking phages such that they become immune against future attacks. In doing so, in the presence of invading phages phage-resistant bacteria emerge. Second, in general, there are several mechanism with which resistant bacteria $R$ may slow down or stop a bacteriophage infection (Skanata and Kussell, 2021). Again, for sake of brevity, only the effect of $R$ on the $S\to I$ transition rate was considered. By doing so, the transition rate $k_0$ becomes a function of $R$and $P$ and reads (Skanata and Kussell, 2021)

$k_0\left(P,R\right)=\frac{{\beta }_0}{1+\gamma R}P,(3)$

where $\gamma \ge 0$ measures the effectivity of the active defense mechanism. Basically, Equation 3 states that the presence of resistant cells lowers the chance of an effective contact between phages and susceptible bacteria, such that the $R$-dependent effective contact rate reads $\beta ={\beta }_0/(1+\gamma R)$.

Taking a synergetics perspective (Haken, 2004; Frank, 2022), for the model (1–3) the order parameter and the amplitude equations were derived. To this end, using the state vector $\mathbf{X}=(S,I,P,R)$, bacteriophage infections were considered that start close to an initial fixed point ${\mathbf{X}}_{st,0}=(S_{st,0},I_{st,0},P_{st,0},R_{st,0})$ (see Results and Discussions section) of the model. Subsequently, with the help of the eigenvectors ${\mathbf{v}}_j$ obtained from a linear stability analysis the amplitudes $A_j$ were implicitly defined by

$\mathbf{X}={\mathbf{X}}_{st,0}+{\displaystyle \sum _{j=1}^4}{A}_j{\mathbf{v}}_j.(4)$

By constructing a bi-orthogonal basis (Haken, 2004; Frank, 2022) spanned by the vectors ${\mathbf{w}}_j$ with ${\mathbf{w}}_i{\mathbf{v}}_k={\delta }_{ik}$ (Kronecker symbol), the amplitudes were explicitly expressed like

$A_j={\mathbf{w}}_j\mathbf{u}={\mathbf{w}}_j\left(\mathbf{X}-{\mathbf{X}}_{st,0}\right),(5)$

where $\mathbf{u}$ denotes the difference vector $\mathbf{u}=\mathbf{X}-{\mathbf{X}}_{st,0}$. From the model Equations 1–3 and Equation 5, eventually the model amplitude equations of the form

$\frac{\text{d}}{\text{d}t}{A}_j={\lambda }_j{A}_j+N_j\left(\mathbf{A}\right)(6)$

were derived for $j=\mathrm{1,2,3,4}$ with $\mathbf{A}$ constituting the amplitude vector $\mathbf{A}=(A_1,A_2,A_3,A_4)$. In Equation 6 ${\lambda }_j$ denote the eigenvalues of the system and $N_j$ are nonlinear functions in the amplitudes. The order parameter and its order parameter amplitude were identified as the eigenvector ${\mathbf{v}}_j$ and its amplitude $A_j$ corresponding to the potentially positive eigenvalue ${\lambda }_j$ of the model (Haken, 2004; Frank, 2022).

3 Results and discussions

3.1 Amplitude equation perspective

The fixed-point analysis showed that the model (1–3) exhibits the phage-free fixed points defined by

$S_{st}\ge 0,I_{st}=0,P_{st}=0,R_{st}\in \left[0,R_m\right](7)$

for $\mu =0$. For $\mu >0$ Equation 7 holds with $S_{st}=K$. As mentioned in the Methods section, it is assumed that at time $t<0$, i.e., before the infection takes place, the system stays in a fixed point (7). The fixed point is referred to as initial fixed point and denoted by ${\mathbf{X}}_{st,0}$. At time $t=0$ the bacteria population is infected by phages of concentration $P(0)>0$ such that the state is shifted out of its fixed point. Consequently, the model describes infection outbreaks that begin with an initial phage infection of $P(0)>0$ at time $t=0$ and end in a phage-free state defined by Equation 7 or an endemic state if it exists (see below).

The linear stability analysis at ${\mathbf{X}}_{st,0}$ showed that the model for $\mu \ge 0$ exhibits the eigenvalues ${\lambda }_1=-\mu$, ${\lambda }_4=0$, as well as

${\lambda }_{\mathrm{2,3}}=-\frac{k_1+k_2}{2}\pm \sqrt{\frac{{\left(k_1+k_2\right)}^2}{4}-k_1{k}_2+q{\beta }_0{f}_{st,0}{S}_{st,0}}(8)$

with $f_{st,0}=1/[1+\gamma R_{st,0}]$, where the upper (lower) sign holds for ${\lambda }_2$ $({\lambda }_3)$. For $\mu >0$ in Equation 8 and in what follows we must substitute $S_{st,0}=K$. It can be shown that for arbitrary model parameters ${\lambda }_3<0$ holds. In contrast, ${\lambda }_2$ can assume positive or negative values. In this context, note that using the next-generation method, the basic reproduction number $R_0$ of the model (Frank, 2022) can be obtained as $R_0=q{\beta }_0{f}_{st,0}{S}_{st,0}/(k_1{k}_2)$. Case I is defined by $q{\beta }_0{f}_{st,0}{S}_{st,0}<k_1{k}_2\iff {\lambda }_2<0$, which is equivalent to $R_0<1$, such that the fixed point ${\mathbf{X}}_{st,}$ is a neutrally stable/asymptotically stable fixed point for $\mu =0$ and $\mu >0$, respectively. There is no infection outbreak. Tn contrast case II is characterized by $q{\beta }_0{f}_{st,0}{S}_{st,0}>k_1{k}_2\iff {\lambda }_2>0$, which is tantamount to say that $R_0>1$ holds. The fixed point is unstable. The infection dynamics describes an infection outbreak. The inequality means that the infection outbreak scenario, i.e., case II, occurs when the system parameters $q$ and ${\beta }_0$ are relatively large, the initial value $S_{st,0}$ is relatively large, while the initial concentration $R_{st,0}$ is relatively small. For $\mu =0$ the model exhibits only phage-free fixed points. A detailed calculation shows that for $\mu >0$ an endemic fixed point exists if the defense mechanism via the $R$ dynamics cannot stabilized the phage-free fixed point with $S_{st,0}=K$. Mathematically speaking, if ${\lambda }_2>0$ holds for $S_{st,0}=K,R_{st,0}=R_m$, which is equivalent to say that $k_1{k}_2<q{\beta }_0{f}_{st,R(max)}K$ holds (where $f_{st,R(max)}=1/(1+\gamma R_m)$), then an endemic fixed point with $I_{st}\in (0,K)$ and $P_{st}>0$ exists. Having said that since the objective of the study is examine initial infection outbreaks from the phage-free state, we will not dwell on the endemic state.

The linear stability analysis of the phage-free fixed point produced the eigenvectors ${\mathbf{v}}_1=(\mathrm{1,0,0,0})$ and ${\mathbf{v}}_4=(\mathrm{0,0,0,1})$ associated with the zero eigenvalues ${\lambda }_1$ and ${\lambda }_4$. For $j=\mathrm{2,3}$ the eigenvectors read as shown in Equation 9

${\mathbf{v}}_j=\frac{1}{Z_j}\left(\begin{array}{c}\hfill -F_0\left({\lambda }_j+k_1\right)\frac{{\lambda }_j}{{\lambda }_j+\mu }\hfill \\ \hfill F_0{\lambda }_j\hfill \\ \hfill \left({\lambda }_j+k_1\right){\lambda }_j\hfill \\ \hfill F_0\alpha \left(1-R_{st,0}/R_m\right)\hfill \end{array}\right)(9)$

with $F_0={\beta }_0{f}_{st,0}{S}_{st,0}$, where $Z_j$ is a normalization constant such that $|{\mathbf{v}}_j|=1$. It was found that the bi-orthogonal vectors ${\mathbf{w}}_{\mathrm{2,3}}$ of the model exhibit the well-known structure from other epidemiological models (Frank, 2022): ${\mathbf{w}}_2=B^{-1}(0,v_{3,P},-v_{3,I},0)$ and ${\mathbf{w}}_3=B^{-1}(0,-v_{2,P},v_{2,I},0)$, where $v_{j,I}$ and $v_{j,P}$ denote the $I$ and $P$ coordinates of the eigenvectors ${\mathbf{v}}_j$, respectively. Here $B=v_{2,I}{v}_{3,P}-v_{2,P}{v}_{3,I}$. A detailed calculation showed that ${\mathbf{w}}_1$ and ${\mathbf{w}}_4$ associated with ${\mathbf{v}}_1=(\mathrm{1,0,0,0})$ and ${\mathbf{v}}_4=(\mathrm{0,0,0,1})$, respectively, read as shown in Equation 10

${\mathbf{w}}_1=\frac{1}{{\lambda }_3-{\lambda }_2}\left(\begin{array}{c}\hfill {\lambda }_3-{\lambda }_2\hfill \\ \hfill \left({\lambda }_2+k_1\right)\left({\lambda }_3+k_1\right)\left[\frac{1}{{\lambda }_2+\mu }-\frac{1}{{\lambda }_3+\mu }\right]\hfill \\ \hfill -F_0\left[\frac{{\lambda }_2+k_1}{{\lambda }_2+\mu }-\frac{{\lambda }_3+k_1}{{\lambda }_3+\mu }\right]\hfill \\ \hfill 0\hfill \end{array}\right),{\mathbf{w}}_4=\frac{1}{{\lambda }_2{\lambda }_3}\left(\begin{array}{c}\hfill 0\hfill \\ \hfill -{\alpha }^{*}\left({\lambda }_2+{\lambda }_3+k_1\right)\hfill \\ \hfill {\alpha }^{*}{F}_0\hfill \\ \hfill {\lambda }_2{\lambda }_3\hfill \end{array}\right)(10)$

with ${\alpha }^{*}=\alpha (1-R_{st,0}/R_m)$. As in other virus dynamics models (Frank, 2022), the nonlinear functions $N_j$ occurring in the amplitude Equation 6 can be expressed as projections of a nonlinear vector-valued function $\mathbf{G}$ on the biorthogonal vectors ${\mathbf{w}}_j$ like

$N_j\left(\mathbf{A}\right)={\mathbf{w}}_j\mathbf{G}\left(\delta \left(\mathbf{A}\right),I\left(\mathbf{A}\right),P\left(\mathbf{A}\right),\omega \left(\mathbf{A}\right)\right)(11)$

with $\omega =R-R_{st,0}$. That is, $\delta ,I,P,\omega$ are the coordinates of the difference vector $\mathbf{u}$ introduced in the Methods section. A detailed calculation showed the results shown in Equations 12, 13 that

$\mathbf{G}=\left(-G_I-\frac{\mu }{K}{\delta }^2,G_I,0,G_R\right)(12)$

and

$G_I={\beta }_{0}P\left(\frac{S_{st,0}+\delta }{1+\gamma \left[R_{st,0}+\omega \right]}-\frac{S_{st,0}}{1+\gamma R_{st,0}}\right),G_R=-\alpha I\frac{\omega }{R_m}(13)$

As indicated in Equation 11, the variables $\delta ,I,P,\omega$ are expressed in terms of $A_j$. Explicitly, we have $\mathbf{u}=(\delta ,I,P,\omega )={\sum }_j{A}_j{\mathbf{v}}_j$, see Equation 4. Consequently, the amplitude equations defined by Equation 6 and (11–13) form a closed set of coupled differential equations.

3.2 Implications

3.2.1 Order parameter: essential building-block and initial organization

The model exhibits maximally one positive eigenvalue. Consequently, under the case II scenario with ${\lambda }_2>0(R_0>1)$ the system exhibits an order parameter given by ${\mathbf{v}}_2$ and the order parameter amplitude $A_2$ (Haken, 2004; Frank, 2022). Let us split the state dynamics into two parts $\mathbf{X}={\mathbf{X}}_{\mathit{out}}+{\mathbf{X}}_s$, where ${\mathbf{X}}_s=A_3(t){\mathbf{v}}_3$for $\mu =0$ and ${\mathbf{X}}_s=A_1(t){\mathbf{v}}_1+A_3(t){\mathbf{v}}_3$ for $\mu >0$ describes the dynamics along the stable direction(s) and ${\mathbf{X}}_{\mathit{out}}$ captures the remaining dynamics. Initially, i.e., for $t\approx 0$, we have

${\mathbf{X}}_{\mathit{out}}\approx \mathbf{K}+{\mathbf{v}}_2{A}_2\left(0\right)\exp \left({\lambda }_{2}t\right)(14)$

with $\mathbf{K}=$ constant and $\mathbf{K}={\mathbf{X}}_{st,0}+A_1(0){\mathbf{v}}_1+A_4(0){\mathbf{v}}_4=$for $\mu =0$, whereas $\mathbf{K}={\mathbf{X}}_{st,0}+A_4(0){\mathbf{v}}_4=$ for $\mu >0$. Equation 14 describes the dynamics along the unstable direction away from the initial fixed point (i.e., the outwards dynamics). In contrast, ${\mathbf{X}}_s$ describes the dynamics towards the unstable direction, i.e., towards the order parameter. Consequently, the order parameter ${\mathbf{v}}_2$ describes the emergent organization of the multi-species physiological network and its amplitude $A_2$ describes how this organization evolves over time.

Since ${\mathbf{X}}_s$ initially decays in magnitude over time, when considering the initial infection dynamics we may neglect its contribution to the state dynamics. If so, then any state change defined by $\mathrm{\Delta }\mathbf{X}=\mathbf{X}(t)-{\mathbf{X}}_{st,0}$ approximately is given by

$\mathrm{\Delta }\mathbf{X}\approx {\mathbf{v}}_2\mathrm{\Delta }{A}_2\approx {\mathbf{v}}_2\left(\exp \left({\lambda }_{2}t\right)-1\right).(15)$

Equation 15 illustrates again that the order parameter describes the essential building-block that shapes an infection outbreak in the multi-species network under consideration including the defense reaction (see component ${\mathbf{v}}_{2,R}$).

3.2.2 Stopping mechanisms

The exponential increase along ${\mathbf{v}}_2$ as described by Equation 15 is slowed down and eventually stopped at some point in time. In line with the stability analysis let us assume that $\delta ,I,P,\omega$ are small quantities of the order $ϵ$. Then $N_2$ when interpreting $N_2$ as a function of the difference variables can be expanded such that the amplitude equation for $A_2$ reads

$\frac{\text{d}}{\text{d}t}{A}_2={\lambda }_2{A}_2-\underset{>0}{\underbrace{v_{3,p}{\beta }_0{f}_{st,0}}}\underset{\ge 0}{\underbrace{\left(\gamma S_{st,0}{f}_{st,0}P\omega -P\delta \right)}}+O\left({ϵ}^3\right).(16)$

Note that $\delta <0$ for any $t>0$. Consequently, the network physiology produces two mechanisms that slow down the exponential increase of $A_2$: the decay in susceptibles in the presence of phages as measured by the interaction term $-P\delta >0$ and the increase of the number of resistant bacteria again in the presence of phages as measured by the interaction term $P\omega >0$. The former mechanism is a passive mechanism that simply states that the exponential infection outbreak slows down due to a decay of the resources (i.e., susceptible bacteria). The latter mechanism is an active mechanism that states that the introduction of resistant bacteria mutations into the bacteria colony has the desired effect of slowing down the phages invasion.

3.2.3 Linear predictor equations

Equation 15 implies that all species initially satisfy linear regression equations of the form as shown in Equation 17

$X_i\approx a_{i,j}+r_{i,j}{X}_j,r_{ij}=v_{2,i}/v_{2,j}.(17)$

Accordingly, any species of the network can be used to predict any other network species (assuming $r_{ij}\ne 0$ for all $i,j$). The network components are coupled by linear order parameter links. For example, the bacteriophage population size $P$ may be used to construct regression models like

$S=a_{S,P}-\frac{F_0}{{\lambda }_2}P,I=\frac{F_0}{{\lambda }_2+k_1}P,R=a_{R,P}+\frac{\alpha F_0}{\left({\lambda }_2+k_1\right){\lambda }_2}P,(18)$

where $a_{X(i),X(j)}$ are intercept parameters depending on ${\mathbf{X}}_{st,0}$. As indicated in Equation 18, $a_{I,P}=0$ because of $I_{st,0}=P_{st,0}=0$.

3.2.4 Limited impact of bacterial growth term

Clearly, the bacterial growth term may affect the $S$ dynamics. However, it does not affect the order parameter eigenvalue ${\lambda }_2$ and it does not affect the orientation of the order parameter ${\mathbf{v}}_2$ in the 3D space $(I,P,R)$, which is of primary concern. Consequently, the initial outbreak dynamics in the $(I,P,R)$ space as determined by the order parameter dynamics (see Section 3.2.1) is completely unaffected by the bacterial growth term.

3.3 2D infected/infectious species dynamics and double exponential dynamics approximation

The dynamics of the infected and infectious species $I$ and $P$ is completely described by the amplitudes $A_2$ and $A_3$. The reason for this is that the eigenvectors ${\mathbf{v}}_1$ and ${\mathbf{v}}_4$ do not exhibit components in the $I-P$ subspace. The mapping from $A_2$ and $A_3$ to $I$ and $P$ reads

$\left(\begin{array}{c}\hfill I\hfill \\ \hfill P\hfill \end{array}\right)={\mathbf{v}}_2^{\prime }{A}_2+{\mathbf{v}}_3^{\prime }{A}_3,(19)$

where ${\mathbf{v}}_2^{\prime }$ and ${\mathbf{v}}_3^{\prime }$ denote the projections of ${\mathbf{v}}_2$ and ${\mathbf{v}}_3$ into the $I$-$P$ subspace. From Equation 19 it follows that the initial evolution of $I$ and $P$ satisfies a double-exponential function as shown in Equation 20

$\left(\begin{array}{c}\hfill I\hfill \\ \hfill P\hfill \end{array}\right)\approx {\mathbf{v}}_2^{\prime }{A}_2\left(0\right)\exp \left({\lambda }_{2}t\right)+{\mathbf{v}}_3^{\prime }{A}_3\left(0\right)\exp \left({\lambda }_{3}t\right),(20)$

3.4 Scaled model

Using the variable transformations $s=S/S_{st,0}$, $i=I/S_{st,0}$, $p=k_{1}P/(q{S}_{st,0})$, $r=R/R_m$, the model (1–3) becomes

$\frac{\text{d}}{\text{d}t}s=-k_0^{\prime }s+\mu s\left(1-s\right),\frac{\text{d}}{\text{d}t}i=k_0^{\prime }s-k_{1}i,\frac{\text{d}}{\text{d}t}p=k_{1}i-k_{2}p,\frac{\text{d}}{\text{d}t}r={\alpha }^{\prime }I\left(1-r\right)(21)$

with

$k_0^{\prime }\left(r,p\right)=\frac{{\beta }_0^{\prime }}{1+{\gamma }^{\prime }r}p(22)$

and ${\beta }_0^{\prime }={\beta }_0{S}_{st,0}q/k_1$, ${\gamma }^{\prime }=\gamma R_m$, and ${\alpha }^{\prime }=\alpha S_{st,0}/R_m$. Among other things, the scaled model exhibits the following two properties. First, the variables $s,i,p,r$ are dimensionless. Second, the variable transformation $P\to p$ turns the phage variable $P$ into a bacteria-like variable (Frank, 2022). That is, just as the model describes that a susceptible bacteria turns into an infected bacteria, the scaled model describes that an infected bacteria turns into a phage unit in a 1:1 manner when expressing phages in the variable $p$ (rather than in $P$). Mathematically speaking, from Equation 21 it follows that $p$ increases due to the term $k_{1}I$ at the same rate as the number of infected bacteria decays due to the term $-k_{1}I$, which means that the model describes the aforementioned 1:1 transition.

3.5 Simulation

An Euler forward simulation scheme was used to solve Equations 21 and 22. In a first simulation, see Figure 1, only the passive defense mechanism was considered with $\alpha =0$ and $\gamma =0$. In a second simulation, see Figure 2, the active defense mechanism was taken into account with $\alpha =1$ and $\gamma =5$. The remaining parameters and initial conditions were $k_1=0.2$/days (Li et al., 2021), $k_2=1$/day, ${\beta }_0^{\prime }=5$/day, $s_{0,st}=1,i_{st,0}=0,p_{st,0}=0,r_{st,0}=0.01$, and $p(t=0)=0.001$. For sake of brevity, only the most relevant phase curves will be presented and only the first week of the initial outbreak stage will be considered. In this context note that in line with our discussion in Section 3.2.4. In both simulations $\mu =0$ was used.

Figure 1

Figure 1. Comparison of simulated infection dynamics (solid black lines) with the order parameter dynamics (panels) (a,b) and the double exponential approximation (panels) (c,d) shown as gray dotted lines. Phase curves in 2D $s-p$ state spaces (panels (a,c) and $i-p$ state spaces (panels) (b,d) are shown.

Figure 2

Figure 2. As for Figure 1 but for a simulation that takes the active defense mechanism via resistant bacteria into account. Comparison of simulated infection dynamics (solid black lines) with the order parameter dynamics (panels) (a,b) and the double exponential approximation (panels) (c,d) shown as gray dotted lines. Phase curves in 2D s – p state spaces (panels (a,c) and i – p state spaces (panels) (b,d) are shown.

As can be seen in panels (a) and (b) of Figures 1, 2, the phase curves quickly converge towards the order parameter ${\mathbf{v}}_2$ and, subsequently, evolve along ${\mathbf{v}}_2$. By definition, the order parameter ${\mathbf{v}}_2$ does not capture the dynamics along the stable direction ${\mathbf{v}}_3$. As can be seen in panels (c) and (d) of Figures 1, 2, the double exponential approximations can describe the transient initial dynamics towards the order parameter (i.e., the dynamics along ${\mathbf{v}}_3$) as well as the subsequent dynamics along the order parameter ${\mathbf{v}}_2$.

Comparing Figures 1, 2, it can be seen that due to the impact of the phage resistant bacteria the actual dynamics differs from the order parameter dynamics to a greater extent. Likewise, the actual dynamics departs earlier from the double exponential approximations. These observations do not come as a surprise since the active defense mechanism slows down and eventually stops the infection outbreak, see Equation 16. Therefore, the actual dynamics will deviate earlier from the order parameter dynamics, on the one hand, and the double exponential dynamics, on the other hand. Roughly speaking, the active mechanisms weakens the linear order parameter link between the network components.

4 Conclusions and limitations

We conclude that under appropriate conditions as specified in the Methods and Simulation sections the order parameter and its amplitude characterize the (self-)organization of a bacteriophage infection and the corresponding bacterial defense. In physics various experiments have been conducted to test specifically predictions of the theory of self-organization (and synergetics) as presented above. Therefore, just as in physics, the results presented above may serve as a basis for conducting laboratory experiments on bacteriophage infections to test the order parameter hypothesis. Moreover, we conclude that linear regression models as derived above may be used to estimate species that are difficult to observe on the basis of species that can be measured more conveniently. For sake of brevity, in the current study, properties of the endemic fixed point as studied, e.g., by Li et al. (2021) have not been examined in detail. Likewise, the current study was limited to consider one possible defense mechanism while alternative mechanisms (Skanata and Kussell, 2021; Weitz et al., 2005) were ignored. A more comprehensive study (which is beyond the scope of this Brief Report) may overcome such limitations by generalizing the results presented above.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author contributions

TF: Writing – original draft, Writing – review and editing.

Funding

The author(s) declare that no financial support was received for the research and/or publication of this article.

Conflict of interest

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Abedon, S. T. (2012). Bacterial immunity against bacteriophages. Bacteriophage 2, 50–54. doi:10.4161/bact.18609

PubMed Abstract | CrossRef Full Text | Google Scholar

Bloch, S., and Wegrzyn, A. (2024). Editorial: bacteriophage and host interactions. Front. Microbiol. 15, 1422076. doi:10.3389/fmicb.2024.1422076

PubMed Abstract | CrossRef Full Text | Google Scholar

Frank, T. D. (2022). COVID-19 epidemiology and virus dynamics: nonlinear physics and mathematical modeling. Berlin: Springer. doi:10.1007/978-3-030-97178-6

CrossRef Full Text | Google Scholar

Geng, P., Flint, E., and Bernsmeier, C. (2022). Plasticity of monocytes and macrophages in cirrhosis of the liver. Front. Netw. Physiol. 2, 937739. doi:10.3389/fnetp.2022.937739

PubMed Abstract | CrossRef Full Text | Google Scholar

Granger, T., Michelitsch, T. M., Bestehorn, M., Riascos, A. P., and Collet, B. A. (2024). Stochastic compartment model with mortality and its application to epidemic spreading in complex networks. Entropy 26, 362. doi:10.3390/e26050362

PubMed Abstract | CrossRef Full Text | Google Scholar

Haken, H. (2004). Synergetics: introduction and advanced topics. Berlin: Springer. doi:10.1007/978-3-662-10184-1

CrossRef Full Text | Google Scholar

Hutt, A., and Haken, H. (2020). Synergetics. New York: Springer. doi:10.1007/978-1-0716-0421-2

CrossRef Full Text | Google Scholar

Li, X., Huang, R., and He, M. (2021). Dynamics model analysis of bacteriophage infection of bacteria. Adv. Differ. Equations 2021, 488. doi:10.1186/s13662-021-03466-x

CrossRef Full Text | Google Scholar

Nowak, M. A., and May, R. M. (2000). Viral dynamics: mathematical principles of immunology and virology. New York: Oxford University Press.

Google Scholar

Pastor-Satorras, R., Castellano, C., Van Mieghem, P., and Vespignani, A. (1998). Epidemic processes in complex networks. Rev. Mod. Phys. 87, 925–979. doi:10.1103/RevModPhys.87.925

CrossRef Full Text | Google Scholar

Skanata, A., and Kussell, E. (2021). Ecological memory preserves phage resistance mechanisms in bacteria. Nat. Commun. 12, 6817. doi:10.1038/s41467-021-26609-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Uhl, C. (1999). Analysis of neurophysiological brain functioning. Berlin: Springer. doi:10.1007/978-3-642-60007-4

CrossRef Full Text | Google Scholar

Weitz, J. S., Hartman, H., and Levin, S. A. (2005). Coevolutionary arms races between bacteria and bacteriophage. Proc. Natl. Acad. Sci. U. S. A. 102, 9535–9540. doi:10.1073/pnas.0504062102

PubMed Abstract | CrossRef Full Text | Google Scholar

Wunner, G., and Pelster, A. (2016). Self-organization in complex systems: the past, present, and future of synergetics. Berlin: Springer. doi:10.1007/978-3-319-27635-9

CrossRef Full Text | Google Scholar

Zborowsky, S., Seurat, J., Balacheff, Q., Ecomard, S., Mulet, C., Minh, C. N. N., et al. (2025). Macrophage-induced reduction of bacteriophage density limits the efficacy of in vivo pulmonary phage therapy. Nat. Commun. 16, 5725. doi:10.1038/s41467-025-61268-1

PubMed Abstract | CrossRef Full Text | Google Scholar