## Abstract

This paper provides a mathematical analysis of a virus–marine bacteria interaction model. The model is a simplified case of the model published and used by Middelboe (Middelboe, M. 2000 *Microb. Ecol*. **40**, 114–124). It takes account of the virus, the susceptible bacteria, the infected bacteria and the substrate in a chemostat. We show that the numerical values of the parameters given by Middelboe allow two different time scales to be considered. We then use the geometrical singular perturbation theory to study the model. We show that there are two invariant submanifolds of dimension two in the four-dimensional phase space and that these manifolds cross themselves on the boundary of the domain of biological relevance. We then perform a rescaling to understand the dynamics in the vicinity of the intersection of the manifolds. Our results are discussed in the marine ecological context.

## 1. Introduction

In the early 1990s, the very large abundance of viruses in marine ecosystems (Bergh *et al*. 1989) and their potential impact on the bacterial biomass (2.9) were demonstrated. Suttle (2005) considers that the biomass of viruses in the world’s oceans could be compared to that of 75 million blue whales and is the ‘second largest component of biomass after prokaryotes’ (see Wilhelm & Suttle (1999) and Suttle (2005) and references therein). Furthermore, viruses increase the mortality of bacteria and primary producers. This lysis affects the nutrient cycle in the sea by increasing the amount of limiting nutrients and also by modifying the quality of the nutrient sources (Wilhelm & Suttle 1999). In order to understand the biogeochemical cycles, it is thus necessary to integrate the viral component into the ecosystem study. Moreover, the functioning of ecosystems is greatly affected (Fuhrman 1999). For instance, it has been demonstrated that the microbial community structure can be modified by viral activities (Fuhrman 1999; Middelboe *et al.* 2001). Though some models have been proposed in the literature to quantify some aspects of the role of viruses in ecosystems (Thingstad & Lignell 1997; Middelboe 2000; Middelboe *et al.*2001, 2003), only a few of them deal with the mathematical analysis of the dynamics (Beretta & Kuang 1998). These last authors consider a three-dimensional system describing the interactions between a bacterial population and a bacteriophage. The bacterial population is subdivided into susceptible and infected subpopulations. The number of viruses produced by a lysed bacterium (the burst coefficient) is chosen as the bifurcation parameter. The authors show that, if this number is low, then the virus disappears. Above a critical value of this parameter, there exists a positive stable equilibrium. Finally, there is another critical value above which a Hopf bifurcation occurs and fluctuations appear. Carletti (2002) and Carletti *et al.* (2004) analyse the noise-induced effects on the Beretta–Kuang model dynamics, but the main result shows that, in the investigated situations, the deterministic model properties are robust.

In the present paper, we consider a simplified version of the model described by Middelboe (2000). In the model on which we focus here, four variables are described, namely the susceptible and infected bacteria, the virus and the substrate in a chemostat. The introduction of an explicit variable for the substrate allows us to consider mass conservation. The previous work of Beretta & Kuang (1998) is based on the logistic growth of the susceptible bacterial subpopulation. Taking explicit account of the substrate allows us to understand how the virus affects the interaction between the bacteria and their resource. It is one of the questions addressed in the paper of Fuhrman (1999) for instance. We thus have to deal with a more realistic but more complex model than the Beretta–Kuang model. The values of the parameters of the model are provided by Middelboe (2000) and it turns out that the burst coefficient (which is a dimensionless parameter) is very high. This fact seems to be in accordance with field evidence (Wilhelm & Suttle 1999). We shall take advantage of this high value to use the geometrical singular perturbation theory. Consequently, we are able to analyse the mathematical dynamics of the five-dimensional system. This paper also aims to illustrate a method that allows one to deal with more realistic and complex models. For instance, in order to understand the role of a virus in a microbial community structure and in nutrient cycles, we need to include the nutrients, the virus and the bacterial populations in the model. This leads to more complex systems. The method used in this paper belongs to the so-called aggregation of variables methods (Auger & Roussarie 1994; Auger & Poggiale 1998; Auger *et al*. 2008*a*,*b*). The aggregation methods aim to reduce the complexity of the system in order to make the mathematical study easier and to gain insight into the biological knowledge. Our approach consists in using a time-scale separation argument and taking advantage of this separation to split the mathematical study into fast and slow dynamics. First, we describe the model analysed in this paper and prove that the dynamics is bounded. Next, the time-scale separation is performed. Then, we use the geometrical singular perturbation (GSP) theory (see Fenichel 1971, 1979; Sakamoto 1990; Jones 1994; Wiggins 1994; Jones & Khibnik 2000) in order to analyse the model. The paper ends with a discussion and conclusions section where the perspectives of this work are announced.

## 2. The virus–bacteria model

The model can be seen as an extension of the Beretta & Kuang (1998) model to the chemostat environment. It reads
2.1
where *s*, *x*, *y* and *z* are the respective amounts (concentration or whatever) of substrate, susceptible bacteria, infected bacteria and virus; *k* denotes the contact rate, i.e. the proportion of bacteria that are infected per unit virus and per unit time; and *λ* is the lysis rate induced by infection. Lysis leads to release of the cellular content and of the produced virus; *δ* is the amount of virus released by each bacterial lysis, and is often called the burst coefficient. We first proceed to a non-dimensionalization in order to eliminate some parameters and to simplify the mathematical analysis. We thus introduce the new variables , , , and , and the new parameters , , and . Replacing the variables and parameters and omitting the bars over the symbols (for ease of notation) provides the following system:
2.2a
2.2b
2.2c
2.2d
Let us consider the domain *Ω*={(*s*,*x*,*y*,*z*);*s*≥0,*x*≥0,*y*≥0,*z*≥0} and the compact set ; we now prove the positive invariance of these sets.

### Theorem 2.1.

*Let us consider (s(0),x(0),y(0),z(0))∈ Ω. Then the solution (s(t),x(t),y(t),z(t)) of system (2.2) remains in Ω for all t≥0. Furthermore, the ω-limit sets of trajectories initiated in Ω lie in*

*.*

### Proof.

The proof of the invariance of *Ω* is straightforward and is based on the fact that, for each variable, the corresponding derivative is positive when the variable is null. Let us prove the second result. We consider *v*=*s*+*x*+*y*+*z*/*δ*. It follows that
Thus, the *ω*-limits are contained in the set {*v*=1}, which is included in . ▪

In this paper, we prove that, with system (2.2), a Hopf bifurcation can occur: either a globally stable equilibrium exists (figure 1) or a stable limit cycle appears (figure 2). The proof is based on perturbation theory arguments. We prove that a two-dimensional invariant manifold exists and we reduce the analysis on this manifold. Then we prove that the dynamics on this manifold is a perturbation of a centre. Melnikov techniques allow us to conclude. Middelboe (2000) provides an estimation of the model parameters, which can be seen in table 1. The burst coefficient *δ* is usually high (Beretta & Kuang 1998; Middelboe 2000). This is also the case for the parameter *λ*. In order to analyse system (2.2), we consider a small dimensionless parameter *ε*, which allows us to represent the different orders of magnitudes of the different rates involved in the system. Let and . We get the following system, in which we again omit the bars over the symbols in order to simplify the notation:

### Theorem 2.2.

*Let us consider the parameter μ=kδ/λ and let us assume that a>b+1. If μ<ab/(1+b)(a−b−1)+O(ε), then system (2.2 ) has a globally stable equilibrium; otherwise, there are fluctuations.*

In order to prove this theorem, we consider an equivalent model, written as a slow–fast system, and prove the theorem on this new model. We thus perform the following change of variables:
2.3
With these new coordinates, the system (2.2) reads
2.4a
2.4b
2.4c
2.4d
2.4e
In the new coordinates setting, the sets *Ω* and 𝒦 are, respectively, transformed into
and
From theorem 2.1, these sets are invariant under the flow defined by system (2.4). In this system, the variables *x* and *y* are fast with respect to variables *u* and *v*, since the derivatives of the latter are proportional to *ε*. This system is called a slow–fast system. GSP theory provides geometrical tools for the study of such systems. The main ones are the Fenichel (1971, 1979) theorems on the persistence of normally hyperbolic manifolds.

### Lemma 2.3.

*For ε=0 and u≠0, the system (2.4) admits two normally hyperbolic manifolds 𝒫*_{i,0}*, with i=1,2. Furthermore 𝒫*_{1,0} *admits a two-dimensional normally stable manifold, while 𝒫*_{2,0} *admits a one-dimensional normally stable manifold and a one-dimensional normally unstable manifold.*

### Proof.

We consider the non-perturbed system obtained from equation (2.4) by letting *ε*=0 and we get
2.5a
2.5b
2.5c
2.5d
2.5e
Let us consider the invariant sets, which are the sets of equilibria of system (2.5):
and
We call *X*_{ε} the vector field associated to system (2.4). The linearizations of *X*_{0} along the previous invariant sets 𝒫_{0,1} and 𝒫_{0,2} are given by the matrices
and
respectively, where the symbol * denotes entries that are unimportant for our study. The matrix *D**X*_{0,1}(*u*) has three null eigenvalues and two negative eigenvalues, −*k**δ**u* and −*λ*. It is also straightforward to see that the matrix *D**X*_{0,2} admits one positive eigenvalue, one negative eigenvalue and three null eigenvalues. The product of the non-vanishing eigenvalues of *D**X*_{0,2} is −*k**δ**λ**u*. It follows that, for each *u*>0, the invariant manifolds 𝒫_{0,i}, *i*=1,2, are normally hyperbolic. Moreover, the invariant manifold 𝒫_{0,1} has a normally stable manifold of dimension two, and the second invariant manifold 𝒫_{0,2} has a normally stable manifold of dimension one and a normally unstable manifold of dimension one. ▪

These invariant manifolds 𝒫_{0,i} cross themselves on the set {*u*=0}, where they lose the normal hyperbolicity property. Let us denote by *p*_{c} the point of coordinates (0,0,0,1) and by *B*_{pc}(*r*) the ball centred at *p*_{c} with radius *r* for the usual Euclidean norm. We now prove the following lemma.

### Lemma 2.4.

*Let α be a positive real number and* *be an initial condition for system (2.4) with strictly positive coordinates. If ε is small enough, there exists T>0, such that, for some t>T, (x(t),y(t),u(t),v(t))∈B*_{pc}*(α).*

### Proof.

For any initial condition in , the trajectory enters an *α*-neighbourhood of the compact set ; this follows from the *v* equation of system (2.5). We then consider the compact set . We thus just have to consider initial conditions in the set . We denote by *p*_{0}=(*x*(0),*y*(0),*u*(0),*v*(0)) the initial condition in . If *p*_{0}∈*B*_{pc}(*α*), then the result is already obtained. Let us assume that *p*_{0}∉*B*_{pc}(*α*). It follows that *u*(0)≥*α*. Indeed, if *u*(0)<*α*, since *x*(0)≤*u*(0) and *y*(0)≤*u*(0) and as we assumed that ∥1−*v*(0)∥<*α*, then *p*_{0}∈*B*_{pc}(*α*). Thus *u*(0)≥*α*. Now, let us consider a positive number *q*<1 and the compact set . According to Fenichel’s theorem, there exists a real number *ε*_{0}>0, such that, for all *ε*<*ε*_{0}, the invariant manifold 𝒫_{0,1} persists and there exists an invariant manifold 𝒫_{ε,1} for system (2.4). On this invariant manifold, system (2.4) reduces to
2.6a
2.6b
Moreover, after a transient time, there exists a solution of system (2.6) that, for the solution (*x*(*t*),*y*(*t*),*z*(*t*),*v*(*t*)) of equation (2.4), satisfies
2.7
and
2.8
as long as the trajectories remain in Δ_{q} and *ε* is small enough. Now, there exists *T*>0 such that and we consider the lowest *T*>0 for which this occurs. It follows that *u*(*t*)=*q**α*+*O*(*ε*)<*α* if *ε* is small enough. Since *u*(*T*)<*α*, (*x*(*t*),*y*(*t*),*u*(*t*),*v*(*t*))∈*B*_{pc}(*α*). This ends the proof. ▪

The previous result claims that *u* decreases slowly and that, after a transient time, the trajectories of system (2.4) reach a neighbourhood of *p*_{c}. In order to understand the dynamics generated by the system (2.4), we thus need to understand the dynamics around the point *p*_{c}, which is in the region where the above-mentioned normally hyperbolic manifolds have lost their normal hyperbolicity. Let us perform a rescaling around this set:
The system (2.4) then reads
2.9a
2.9b
2.9c
2.9d

### Lemma 2.5.

*We assume that a>1+b. For ε small enough, the system (2.9) admits a two-dimensional invariant manifold ℳ*_{2,ε} *on which it can be reduced to a perturbation of a Lotka–Volterra predation model. Moreover, let μ*_{0}*= ab/(1+b)(a−b−1)+O(ε). The reduction on ℳ*_{2,ε} *of system (2.9) has a positive equilibrium, which is stable when the parameter μ=kδ/λ is smaller than μ*_{0}*, while it is unstable when μ>μ*_{0}*: μ*_{0} *is a bifurcation value of the parameter μ and a Hopf bifurcation takes place when μ crosses the bifurcation value.*

### Proof.

System (2.9) is a slow–fast sytem. The invariant manifold when *ε*=0 is the set *Y* =0. According to Fenichel’s theorem, there exists an invariant manifold ℳ_{3,ε}, for *ε*>0 small enough, on which the previous system can be reduced to
Let *R*=*U*−*X* and *t*=*ε**τ*. Then we obtain
2.10a
2.10b
2.10c
Notice that system (2.10) can be split into two subsystems: the first one corresponds to equations (2.10a,*b*) and the second subsystem is just equation (2.10c), which is independent of the variables *X* and *R*. As a consequence, we can see the first subsystem (2.10a,*b*) as a non-autonomous but asymptotically autonomous system (e.g. Markus 1956; Thieme 1992). The dependence of equations (2.10a,*b*) on the variable *V* is in fact in the *ε*-order terms and this may have important consequences since the case with *ε*=0 is not structurally stable. From a geometrical point of view, we can say that the invariant manifold ℳ_{3,ε} is the graph of the map
2.11
The previous remark means that there is a two-dimensional invariant manifold ℳ_{2,ε}, which is the graph of and on which the system (2.9) is reduced to the first two equations of system (2.10) in which we set *V* =0. Since, in general, *V* ≠0 (but of course will tend to 0), we will need the Thieme (1994*b*) theorem to conclude. As a consequence, we now focus on equations (2.10a,*b*) in which we set *V* =0 and we will prove later that this permits us to conclude for the whole system (2.10). We thus analyse the dynamics on ℳ_{2,ε}. If *a*<1+*b* then (0,0) is globally asymptotically stable. However, if *a*>1+*b*, the subsystem formed by the first two equations is an *ε*-perturbation of the Lotka–Volterra system, for which one can show that the following application is a first integral:
where *H*_{0} is a constant. For *ε*=0, the positive equilibrium is
The constant *H*_{0} is chosen such that *H*(*X*_{E},*R*_{E})=0. The Lotka–Volterra model is not a structurally stable differential system. As a consequence, the *ε*-order terms in the system (2.10) may play a key role in the dynamics. In order to determine the dynamics on the invariant manifold ℳ_{ε}, we perform an expansion of this manifold with respect to the small parameter *ε*. Fenichel’s theorem claims that the invariant manifold ℳ_{2,ε} is the graph of the map
with *Y* (*X*,*R*,0)=0. Thus, *Y* (*X*,*R*,*ε*)=*ε**w*_{1}(*X*,*R*)+*o*(*ε*). Let us determine *w*_{1}:
2.12a
2.12b
Thus we get
This expression allows us to obtain a first-order approximation with respect to *ε* of the reduction of system (2.4) on the invariant manifold ℳ_{2,ε}. We get
2.13a
2.13b
This system admits an equilibrium (*X*_{E,ε},*R*_{E,ε})*ε*-close to (*X*_{E},*R*_{E}). Let us consider the half straight line from this equilibrium and parallel to the *X*-axis, and let us parametrize this line by the levels of *H*. On this half straight line, there exists a first return map (Poincaré map), denoted by *P*_{ε}. Since, for *ε*=0, all the trajectories in the positive quadrant are close curves (except the equilibrium (*X*_{E},*R*_{E})), it follows that *P*_{0}≡id. Let us define the displacement application:
The sign of this application provides us with the existence and the stability of periodic orbits and of the positive equilibrium. There is a periodic orbit by the point *h*_{0} if and only if *P*_{ε}(*h*_{0})=*h*_{0}. In order to determine the displacement application, we consider the dual differential form associated to the vector field defined by the differential equations (2.13a) and (2.13b) divided by the positive function (*X*,*R*)↦*X**R*:
where
From the Poincaré lemma, we get
We thus need to determine the sign of . Stokes’ theorem implies that
In the vicinity of the equilibrium, *R*≃*R*_{E}=[*a*−(1+*b*)]/(*b*+1)*k**δ*. Thus the term in the above double integral equals
If *k**δ*/*λ*>*a**b*/(*b*+1)[*a*−(1+*b*)] and if *h* is small enough, then , and thus the equilibrium is unstable. When *h* is large enough, the contribution of the term *a**b*/(*b*+1)^{2}*R* in the integral is larger since large *h* implies that trajectories visit regions with small *R*. It follows that the integral is thus negative. As a consequence, if *h* is large enough and if *ε* is small enough, then *P*_{ε}(*h*)<*h*. If the inequality *k**δ*/*λ*>*a**b*/(*b*+1)[*a*−(1+*b*)] is satisfied, then there exists a periodic trajectory. This ends the analysis that we needed on ℳ_{2,ε}. We now complete the proof to show that system (2.9) has the same dynamics as its reduction to the two-dimensional manifold ℳ_{2,ε}. We already know from Fenichel’s theorem that system (2.9) has the same dynamics as its reduction on ℳ_{3,ε}. In other words, the *ω*-limits of system (2.9) are those of system (2.10). Moreover, they are on ℳ_{2,ε} because of the form of the third equation in system (2.10). Furthermore, as we explained earlier, system (2.10) can be seen as a non-autonomous but asymptotically autonomous system, for which the asymptotic limit system is given by the reduction on ℳ_{2,ε}. Markus (1956) provided some results on such asymptotically autonomous systems, which have been extended by Thieme (1992, 1994*a*,*b*). The subsequent articles provide results *à la Poincaré–Bendixon*. We apply theorem 1.2 in Thieme (1994*b*) to conclude that the *ω*-limit sets of system (2.10) are either the positive equilibrium of the reduction on ℳ_{2,ε} or a limit cycle of this reduction. This completes the proof of lemma 2.5. ▪

We can now prove theorem 2.2.

### Proof of theorem 2.2.

If *μ*<*μ*_{0}, then the result is straightforward. If *μ*>*μ*_{0}, the system (2.9) has no stable equilibrium. Let *p* be an initial condition. If the *ω*-limit set of *p* is in *B*_{pc}(*α*), then it is the *ω*-limit of a trajectory of (2.9). Since *μ*>*μ*_{0}, we proved in lemma 2.5 that the *ω*-limit sets of system (2.9) are limit cycles. Finally, if the *ω*-limit set of *p* is not in *B*_{pc}(*α*), then the trajectory fluctuates. Indeed, according to the latter assumption, *ω*(*p*) have some points outside of *B*_{pc}(*α*). Furthermore, according to lemma 2.4, it follows that *ω*(*p*) intersects *B*_{pc}(*α*). As a consequence, the trajectory of *p* is fluctuating. ▪

## 3. Discussion and conclusions

We have proved that the presence of the virus may lead to fluctuations of the susceptible bacterial population. Even though this result has already been obtained by Beretta & Kuang (1998), we show it on a more realistic representation of the system, since the model we analysed here explicitly accounts for the substrate. However, the underlying mechanism of these fluctuations is the same in both models and is due to the large value of the burst coefficient, that is, the amount of virus liberated in the environment by one lysed bacterium. Furthermore, this result is important in the context of biogeochemical cycles since the presence of various species of bacteria in the environment leads to the degradation of a wide spectrum of organic compounds. It will be the aim of future work to understand the effect of the virus on the speed of substrate degradation. Here we used time-scale separation and aggregation techniques based on GSP theory to simplify the mathematical analysis of the model. We aim to continue this approach in order to study the effect of virus on microbial food webs. Such food webs usually involve time-scale separation, as reported for instance by Muratori & Rinaldi (1992), Kooi *et al.* (1998) or Deng & Hines (2002). We could even have more than two different scales. Thus, even if the number of variables is increasing, the reduction to invariant submanifolds in the phase space helps us to perform the mathematical treatment of the models. Finally, we did not prove in the present paper that the fluctuations of bacterial abundances correspond to periodic solutions of the model nor that a generic Hopf bifurcation occurs to generate these fluctuations from a steady-state situation. This will be the aim of our future work.

## Footnotes

One contribution of 17 to a Theme Issue ‘From biological and clinical experiments to mathematical models’.

- © 2009 The Royal Society