Mathematical model of interaction Escherichia coli and Coliphages

1.   Introduction

Bacteriophages (phages), which are viruses that infect bacteria, are enormously abundant in the biosphere. Viral particles are estimated to number about $10^{31}$, including bacteriophages ^[1]. Bacteriophages can present lytic cycles or lysogenic cycles. In the lytic cycle, phages infect the host and replicate inside, after that, they destroy its host cell to release virion progeny. While in the lysogenic cycle, phages infect the host but do not replicate, and they integrate into the host or exist as plasmids within their host cell, the genetic material then is transmitted to the next generation of bacteria ^[2,3]. Phages have gained importance in medicine and industry for their application or use in the fight against antibiotic resistance (phage therapy), as vehicles for vaccines delivery, or as display system for many proteins and antibodies, among others ^[4,5]. Phage therapy offers a number of benefits and only a few drawbacks or concerns in the fight against the global problem of multidrug-resistant bacteria. Most of these problems can be solved with a combination of proper phage selection, efficient formulation, and increased clinician knowledge and familiarity with product application ^[6]. Understanding the dynamics of these viruses and the bacterial host from different perspectives is vital to gain insight into the synergy that exists in this interaction. In this sense, many mathematical models using ordinary differential equations have been proposed to understand this dynamics, for example in ^[7], they studied the dynamics of lytic RNA phage MS2 and its host E.coli C-3000, a mathematical model simulating the interaction between bacteriophages and their bacterial hosts considering PH and temperature was developed in ^[8]. In ^[9] they studied the non-linear kinetics between the pathogen bacteria Campylobacter jejuni and a lytic phage considering susceptible and resistant bacteria, infected cells and free phage particles. A model considering the lytic and lysogenic life cycles of phages and the prophage induction in the interaction between phages and bacteria is presented in ^[10]. In ^[11], a SI type model is considered, where the equation for phages is not explicit. In ^[12], the authors formulated a SIR type mathematical model where the infection rate is a law of mass action. In all these articles either the approach is a mathematical development without the contrast with real data ^[10,11,12] or a development with real data without bothering with mathematical analysis ^[7,8,9]. In this sense, we formulated a mathematical model considering a free lytic phage, sensitive bacteria and infected bacteria where the incidence rate deals with a saturation process due to the free phages, where we explored the stability analysis of equilibrium points using the Lyapunov indirect method. In addition, we used global sensitivity analysis to identify the model parameters that are most relevant and we used data growth of bacteria E.Coli in presence of coliphages to fit the parameter of the model using a genetic algorithm to study the interaction between bacteria-phages. The organization of this paper is as follows: In Section 2, the mathematical model is formulated. In Section 3, we found three equilibrium points, the trivial equilibrium, the phages extinction equilibrium and the coexistence equilibrium for E.coli sensitive, E.coli infected and coliphages. In Section 4, the stability analysis of equilibria is proved. In Section 5, sensitivity analysis is performed. In Section 6, the most relevant parameters are estimated. Finally, in Section 7 we discuss our results.

2.   Model formulation

Let $E(t)$ and $I(t)$ denote the populations sizes of sensitive bacteria and infected bacteria at time $t$, respectively and $C(t)$ the concentration of free phages at time $t$. Bacteria reproduce at a constant per capita rate, which depends on the species of bacteria, in this study we consider the E.coli bacteria type for make the simulations. We consider that there exists an intraspecific and interspecific competition between sensitive and infected bacteria as population size increases and resources become more limited, which is modeled by logistical growth with carrying capacity $N$ and reproduction rate $k$. For the interaction between bacteria and phages we consider the Holling Ⅱ functional response because we believe that there is saturation in the infection process, i.e., the more free phages there are, the less sensitive bacteria there will be. In this sense the sensitive bacteria acquire infection at rate $\beta E\dfrac{C}{a+C}$, where $\beta$ is the infection rate of bacteria and $a$ the half-saturation phages density, the mortality rate of this infected bacteria is $\upsilon$. The free phages growth proportionally to the concentration of infected bacteria $I$ at rate $\alpha I$, where $\alpha$ is the release rate of viral particles and decay at rate $\lambda$. From the above suppositions we derive the following system of non-linear differential equations

The summary of the parameters present in the model is shown in Table 1. The set of biological interest is given by

Table 1.  Model parameters.

Parameter	Definition	Range	Value	Reference
$k$	bacterial growth rate	[0.4, 1.8]	$1.11$	This study
$N$	carrying capacity	[2.5, 4]	$3.2$	This study
$\beta$	The infection rate of bacteria	[1, 20]	Estimated for different MOI	-
$a$	Half-saturation phages density	[0.07, 20]	Estimated for different MOI	-
$\upsilon$	Lysis rate of infected bacteria	[0.05, 4]	$1.002$	[7]
$\alpha$	The release rate of viral particles	[1.4, 1.8]	$1.63$	[8]
$\lambda$	decay rate of viral particles	[0.0003, 2]	$1.032\times 10^{-2}$	[9]

---

 | Show Table

DownLoad: CSV

Lemma 1. The set $\Omega$ defined above is positively invariant for the solutions of the system (2.1).

Proof. Let $E$, $I$, $C$ solutions of system (2.1) and $\Phi = E+I$. Now, adding the first two equations of (2.1) we obtain,

If $\Phi(0)\leq N$, the solutions to before equation satisfies $0\leq \Phi(t)\leq N$, for all $t\geq 0$. Now considering the equation for the free phages in the system (2.1), we deduce,

If $C(0) = C_0\in \Omega$, we conclude that

and that if $t\to \infty$ then $C(t)\to \frac{\alpha N}{\lambda}$. Therefore, the set $\Omega$ is positively invariant.

3.   Equilibrium points

The equilibrium points of system (2.1) are given by the solutions of the system of equations

From the last equation of above system we obtain

By substituing $I$ in the first two equations of (3.1), we obtain

We observe that the equilibrium solutions of system (3.3) are totally determined by $C = 0$ and $C\neq 0$. If $C = 0$, from system (3.3) we get the equation

which has the solution $E = 0$ or $E = N$. Therefore, for the case $C = 0$, we obtain the equilibrium points $P_0 = (0, 0, 0)$, $P_1 = (N, 0, 0)$ that always exists.

Finally, we determine the equilibrium solutions for $C\neq 0$. In this case, the second equation of system (3.3) is written as

and substituting it in the first equation of system (3.3), we get

where $R_0 = \frac{\beta N\alpha}{ \lambda\upsilon a}$. Note that, when $R_0 < 1$, the quadratic Eq (3.4) have two negative values of $C$, and there is not a positive equilibrium point from the Eq (3.1), it follows that a necessary and sufficient condition for the biological meaning of $C$ is $R_0 > 1$, and in this way, we get the equilibrium point $P_2 = (E^*, I^*, C^*)$, where

and $b = \left[a(R_0-1)-a\left(1+\frac{\beta}{\upsilon}\right)-\frac{\alpha\beta^2 N}{k\upsilon \lambda}\right]$.

The above results are summarized in the following proposition.

Proposition 1. System (2.1) always has a trivial equilibrium $P_0 = (0, 0, 0)$ and the equilibrium point $P_1 = (N, 0, 0)$. If $R_0 > 1$, there exists an equilibrium $P_2 = (E^*, I^*, C^*)$, in which E.coli sensitive, E.coli infected and bacteriophages co-exist.

Le us interpret the parameter $R_0 = \frac{\beta N\alpha}{ \lambda\upsilon a}$. One phage during its average lifetime, $\frac{1}{\lambda}$, infects one sensitive bacteria, with rate $\beta N$, this infected bacteria releases $\frac{\alpha}{\upsilon a}$ number of phage particles. In this way, $R_0$ is the net number of phages produced by a phage in a lytic cycle, in a concentration of sensitive bacteria.

4.   Stability

In this section, we determine the asymptotic stability of the equilibrium solutions of the system (2.1). Linearization of this system around an equilibrium $P$ is given by $x' = J(P)x$, where $x = (E, I, C)^T$ and the matrix evaluated at $P$ is

By evaluating $J$ at $P_0$, we verify that the eigenvalues of $J(P_0)$ are $k$, $-\upsilon$ and $-\lambda$, in consequence $P_0$ is always unstable. The eingenvalues of $J(P_1)$ are $-k$, and $\frac{-(\upsilon+\lambda)\pm\sqrt{(\upsilon+\lambda)^2+4(R_0-1)}}{2}$, which have negative real part if $R_0\leq1$.

The result is summarized in the following proposition.

Proposition 2. The trivial equilibrium $P_0$ is always unstable and the equilibrium $P_1$ is locally asymptotically stable if $R_0\leq1$, and unstable otherwise.

Proposition 3. If $R_0\leq 1$, the equilibrium point $P_1$ is globally asymptotically stable in $\Omega$.

Proof. Let $L$ the function defined by

It is easy to check that $L(P_1) = 0$ and $L > 0$ in $\Omega-\{P_1\}$. Besides the orbital derivative of $L$ along solutions of the system (2.1) is given by

which is equivalent to

From the inequality $\frac{1}{1+\frac{C}{a}}\leq 1$, we obtain

Finally, from last equation $\dot{L} < 0$, in $\Omega-\{P_1\}$ if and only if $R_0\leq 1$. Therefore, we have that $P_1$ is globally asymptotically stable.

Finally, to study the stability of point $P_2$, we use the following lemma

Lemma 2. Let $A$ be a $3\times 3$ real matrix. If $tr(A)$, $det(A)$ and $det(A^{[2]})$ are all negative, then all of the eigenvalues of $A$ have negative real part.

A proof of this lemma can be found in ^[13]. The matrix $A^{[2]}$ in lemma is the second additive compound ^[14], which is define by

where $A$ is a $3\times 3$ real matrix

Now, we show that the trace and determinant of the Jacobian (4.2) and the determinant of the second additive compound are all negative, which allows us to conclude that the point $E_2$ is locally asymptotically stable. To see this, we evaluated $J$ defined in (4.1) at $P_2$ to get

Evidently $tr(J(P_2)) = -\frac{kE^*}{N}-\upsilon-\lambda < 0$ and

because from the Eq (3.5), we have

The second compound of Jacobian matrix (4.2) is

Whose determinant is

From Eq (4.3), we have

Then, by removing the brackets of the first two terms, the positive term is canceled and the last term is negative if $R_0 < \frac{k}{\lambda}$ and we can get $det(J^{[2]}(P_2)) < 0$. Since $R_0 > 1$, we must impose that $k > \lambda$. We summarize it in the next proposition.

Proposition 4. If $R_0 > 1$ and $R_0 < \frac{k}{\lambda}$ ($k > \lambda$), the point $P_2$ is locally asymptotically stable.

Remark 1. We emphasize that $R_0 < \frac{k}{\lambda}$ and $k > \lambda$ is one of the many restrictions that can be imposed and is not a determining condition to guarantee mathematically that the determinant of matrix (4.4) is negative. In Section 7, we explore the condition $k > \lambda$, $k < \lambda$ with $R_0 > \frac{k}{\lambda}$.

5.   Sensitivity analysis

In order to determine which are the parameters that most affect the competition dynamics between sensitive and infected bacteria with bacteriophages, we use uncertainty and sensitivity analyses, which focus on measuring the effect of the parameters inputs $k$, $N$, $\beta$, $a$, $\upsilon$, $\alpha$, and $\lambda$ on the output variables $E(t)$, $I(t)$ and $C(t)$. Therefore, we use the data in Table 1 and the Latin hypercube sampling (Lhs) and Partial rank correlation coefficients (PRCCs) to assess global uncertainty and sensitivity analyses following the methodology proposed in ^[15].

We begin by describing the sensitivity analysis for the sensitive bacteria $E$ shown in the Figure 1. There are three parameters that have a strong influence on this variable, namely $k$, $\beta$ and $a$. The parameter $k$ is positively correlated and is initially responsible for the sensitive bacterial load increase, immediately this parameter loses relevance, the parameter $a$ becomes positively correlated, helping the variable $E$ to grow. While, the parameter $\beta$ is negatively correlated which means that if we increased this parameter there will be a decrease in the E.coli bacteria.

Figure 1.  PRCCs of the sensibility analysis performed for the sensitive bacteria $E$, plotted over a time course.

DownLoad: Full-Size Img PowerPoint

We now continue to describe the sensibility analysis for infected bacteria variable $I$. There are two influencing parameters in this variable, these are $\beta$ and $a$. The parameter $\beta$ is strongly positively correlated over the time and that is because it is the infection rate of bacteria, therefore we have an increase in the infected bacteria load while we rise this parameter. The parameter $a$ is negatively correlated over the time, which means that the larger this parameter is, there will be a decrease in infected bacteria load because there will be fewer virus particles free to infect bacteria see Figure 2.

Figure 2.  PRCCs of the sensibility analysis performed for the infected bacteria $I$, plotted over a time course.

DownLoad: Full-Size Img PowerPoint

Finally, we discuss the sensibility analysis for Bacteriophages $C$. The parameters that most affect this variable are $\beta$ and $a$. We emphasize that these parameters are not in the equation for $C$. The parameter $\beta$ (the infection rate of bacteria) is strongly positively correlated over the time and this is because the larger this parameter is, the more bacteria will be infected and therefore the more phages will be released. The parameter $a$ (Half-saturation phages density) is negatively correlated over the time, this means that if we increase this parameter the phage growth will be slower because there will be fewer infected bacteria releasing virus due to the presence of fewer free viral particles see Figure 3.

Figure 3.  PRCCs of the sensibility analysis performed for the Coliphages $C$, plotted over a time course.

DownLoad: Full-Size Img PowerPoint

6.   Parameter estimation

From the last section we note that $\beta$ and $a$ are the parameters that most affect the output variables $E(t)$, $I(t)$ and $C(t)$. Therefore, we select these parameters to be estimated using experimental data published in ^[16], in particular those for E.coli ATCC® 11775™ (gram-negative) and T4-like A coliphages, where they elaborate curves for bacterial growth curves in the presence of bacteriophages using the MOI 0.01, 0.1 and 1 to compare and analyze the activity of bacteriophages. The data used here are shown in the Figure 4. Since the data in ^[16] only show the total remaining bacteria, i.e., the data used do not distinguish between infected and susceptible bacteria, nor does it show coliphages growth data, we find the least square between the data and the sum of sensitive and infected bacteria $E(t)+I(t)$. To do this, we use a evolutionary algorithm called genetic algorithms (ga) where the function to minimize was $fmin(\beta, a) = \sum_i\{data_i-[E(i)+I(i)]\}^2$, where $E(t)$ and $I(t)$ are the outputs of the model (2.1), for an explanation of how the ga works see ^[17].

Figure 4.  Growth of E.coli bacteria under different Multiplicities of infection.

DownLoad: Full-Size Img PowerPoint

First, we estimated parameters for the control curve, for this purpose we use the equation

the simulations give the values $k = 1.11$ and $N = 3.2$, and the result is shown in figure 5. Second, we estimated the parameter $\beta$ and $a$ of the model (2.1) when MOI is $0.01$, $0.1$ and $1$. The estimates of these parameters are shown in the Table 2 and the curve fits are shown in the Figures 6(a), 7(a) and 8(a). In the case MOI $0.01$, we observed that the growth curve of the bacteriophages is slow compared to the other two MOI used, this generates a higher growth of the E.coli bacteria and in this sense it would take more time for its extinction. In addition, we observed a higher growth in the coliphages due to a higher amount of sensitive bacteria. There is a change between sensitive bacteria and infected bacteria that occurs faster when the MOI is higher, from which it is observed that the faster this change occurs, the less sensitive bacteria and bacteriophages we will have at the end.

Figure 5.  E.coli control growth without presence of coliphages.

DownLoad: Full-Size Img PowerPoint

Table 2.  Parameter estimation for different MOI.

MOI	$\beta$	a	Initial condition [S(0) I(0) C(0)]
$0.01$	$2.6$	$0.22$	$[0.24\ \ 0\ 0.0024]$
$0.1$	$3.81$	$0.4428$	$[0.19\ \ 0\ 0.019]$
$1$	$18.92$	$6.68$	$[0.18\ \ 0\ \ 0.18]$

---

 | Show Table

DownLoad: CSV

Figure 6.  The Figure (6a) shows the estimated growth curve of susceptible E.coli bacteria and infected E.coli bacteria in the presence of coliphages when the MOI is 0.01 and (6b) shows the growth curves of each population.

DownLoad: Full-Size Img PowerPoint

Figure 7.  The Figure (7a) shows the estimated growth curve of susceptible E.coli bacteria and infected E.coli bacteria in the presence of coliphages when the MOI is 0.1 and (7b) shows the growth curves of each population.

DownLoad: Full-Size Img PowerPoint

Figure 8.  The Figure (8a) shows the estimated growth curve of susceptible E.coli bacteria and infected E.coli bacteria in the presence of coliphages when the MOI is 1 and (8b) shows the growth curves of each population.

DownLoad: Full-Size Img PowerPoint

7.   About stability condition $R_0 < \frac{k}{\lambda}$

In the proposition 4, we impose the condition that $R_0 < \frac{k}{\lambda}$ and $k > \lambda$, and this assure the local stability of coexistence point $P_2$, but what happen if eventually $k < \lambda$ or $R_0 > \frac{k}{\lambda}$. To explore this we use data in table 1 and the values of $\beta$ and $a$ in the table 2 for MOI $0.01$. In this case, we have $k > \lambda$, and $R_0 > \frac{k}{\lambda}$. If we change the value of $\lambda$ to $0.8$ we will have $k > \lambda$, and $R_0 > \frac{k}{\lambda}$, and the simulation that appear in the Figure 9.

Figure 9.  The Figure (9a), $k > \lambda$ and $R_0 > \frac{k}{\lambda}$ $R_0 = 4662.3$ $\frac{k}{\lambda} = 84.091$, the eingenvalues of Jacobian are $-1.01734$, $0.00086\pm 0.09097i$ and (9b) $k > \lambda$, and $R_0 > \frac{k}{\lambda}$, $R0 = 76.927$ $\frac{k}{\lambda} = 1.3875$, the eingenvalues of Jacobian are $-1.77200$, $-0.02718 \pm 0.54719i$.

DownLoad: Full-Size Img PowerPoint

We conclude that the condition in proposition 4, is sufficient but not necessary. We want to note that in all the parameter estimates shown in the Figures 6, 7 and 8, we have that $k > \lambda$ and $R_0 > k/\lambda$, and the eigenvalues of the Jacobian matrix have complex conjugates with positive real part, this means that if we increase the observation window to 600 hours we will see oscillations as in the Figure 9a. Now we explore $k < \lambda$, and $R_0 > \frac{k}{\lambda}$ and $k > \lambda$, and $R_0 < \frac{k}{\lambda}$. If we change the value of $\lambda$ to $1.46$ we will have $k < \lambda$, and $R_0 > \frac{k}{\lambda}$. the results are shown in Figure 10.

Figure 10.  The Figure (10a), $k < \lambda$, and $R_0 > \frac{k}{\lambda}$, $R0 = 42.152$ $\frac{k}{\lambda} = 0.76027$, the eingenvalues of Jacobian are $-2.37630$, $-0.06455\pm 0.64459i$ and (10b) $k > \lambda$ and $R_0 < \frac{k}{\lambda}$ $R_0 = 1.7352$ $\frac{k}{\lambda} = 2.22$, the eingenvalues of Jacobian are $-1.38423$, $-0.63084$, $-0.25147$. For this simulation we use $\lambda = 0.5$, $\beta = 1.2$ and $a = 7.2$.

DownLoad: Full-Size Img PowerPoint

8.   Discussion

The relationship between phages and bacteria has become relevant because it could play an important role in antibiotic resistance, however this dynamic is complex and has been described by a variety of dynamics at different levels. Here we formulated a mathematical model where we consider the populations of sensitive bacteria, infected bacteria and bacteriophages, and for incidence rate we used an nonlinear function with an inhibitory effect. We found three possibilities for Bacteria-phages relationships: uninfected bacteria, infected bacteria and phages become extinct simultaneously (trivial equilibrium); extinction of infected bacteria and phages (bacteriophages extinction equilibrium); sensitive bacteria, infected bacteria and phages coexist (coexistence equilibrium). The local stability of the three equilibrium point was analyzed by using the Lyapunov's indirect method and the second additive compound matrix. The trivial equilibrium is always unstable, the phages extinction equilibrium is locally asymptotically stable when $R_0\leq 1$ and unstable otherwise. The coexistence equilibrium is locally asymptotically stable when $R_0 > 1$ and $R_0 < \frac{k}{\lambda}$, we clarify that the last inequality is based on our calculations and it is a sufficient condition but not necessary to obtain stability, since our simulations suggest that stability may occurs when $R_0 > \frac{k}{\lambda}$. Something interesting happened with the results of the parameter estimation, although the curve fits are good, we have that the stability condition is not satisfied because $k > \lambda$ and $R_0 > k/\lambda$ and moreover the eigenvalues of the Jacobian matrix have complex conjugates with positive real part. The effect of this is that we will have oscillations that will last in time, which means that we will have periodic increases in the populations of coliphages, the infected and sensitive bacteria as shown in Figure 9a.

From the analysis of sensitivity we concluded that the more influential parameter in the model are the infection rate of bacteria ($\beta$) and Half-saturation phages density ($a$). The parameter $\beta$ helps to increase the density of infected bacteria and the density of phages and decrease the density of sensitive bacteria. while the parameter $a$ plays an inhibitory role, the larger it is, the less number of free viral particles will infect sensitive bacteria, which implies a smaller number of infected bacteria.

The model and the parameter estimation suggested there is a higher growth of the E.coli bacteria when the MOI is $0.01$ than $0.1$ and $1$, this has two implication the coliphages curve growth slow and time of extinction of bacteria is greater, besides due to a higher amount of sensitive bacteria we will have a higher density of phages at the end. One limitation of these estimates is that the data set used does not distinguish between sensitive bacteria and infected bacteria, nor is there data for the growth of coliphages. More experiments are needed where the growth curves of the infected E.coli bacteria and that of the coliphages are considered to improve the estimates made for our model. Despite the enormous effectiveness of phages to control pathogenic bacteria in vitro, there are some problems when applying them to clinical situations, we mention some of them: the therapeutic bacteriophages are not easily identifiable ^[18], the stability of therapeutic formulations is entirely dependent on bacteriophage ^[19]. Bacteria have defense mechanisms that can prevent infection by phages, but a phage cocktail can reduce this risk^[6].

Acknowledgments

The authors express their gratitude to the Government of the Department of Nariño, the University of Nariño and the Ministry of Science, Technology and Innovation (MINCIENCIAS) for funding the project entitled: "Study of the diversity of fecal contamination indicators Escherichia coli and somatic coliphages in Lake Guamuez, strategic productive ecosystem of Nariño", code: 110481865200, CT-80740-517-2019.

Conflict of interest

The authors declare there is no conflict of interest.