Bifurcation Analysis in a Delay Diﬀerential Equations, which Confers a Strong Allee Eﬀect in Escherichia Coli

. The paper addresses the bifurcations for a delay diﬀerential model with parameters which confers a strong Allee eﬀect in Escherichia coli. Stability and local Hopf bifurcations are analyzed when the delay τ or σ as parameter. It is also found that there is a non-resonant double Hopf bifurcation occur due to the vanishing of the real parts of two pairs of characteristic roots. We transform the original system into a ﬁnite dimensional system by the center manifold theory and simplify the system further by the normal form method. Then, we obtain a complete bifurcation diagram of the system. Finally, we provide numerical results to illustrate our conclusions. There are many interesting phenomena, such as attractive quasi-periodic solution and three-dimensional invariant torus.


Introduction
Allee was initially stimulated by an example only loosely linked to the current interpretation of the Allee effect: he showed that goldfish grew faster in water which had previously contained other goldfish, than in water that had not (see [1]). Allee saw these phenomena as "automatic cooperation", believing that the beneficial effects of numbers of animals present in a population represented a fundamental biological principle. By 1953, E. P. Odum was referring to "Allee's principle" as the concept that "undercrowding (or lack of aggregation) may be limiting" (see [16]). Reduction in fitness or population growth at low abundance has received considerable attention in conservation genetics, under such guises as the "50/500 rule" (see [21]), and is also widely debated in fisheries science, where it is usually referred to as depensation (see [15]). Depensation is principally a population level phenomenon, which may or may not arise from changes in individual fitness, and thus need not be directly analogous to the Allee effect.
Although the Allee effect is reasonably well known, the concept has a range of meanings, not all of which are acknowledged by contemporary use. Allee did not provide a definition but he clearly considered "certain aspects of survival values" rather than total fitness. Until 1999, Stephens P. A. (see [22]) defined the Allee effect as: a positive relationship between any component of individual fitness and either numbers or density of conspecifics. After that, the Allee effect be thought as a phenomenon in biology characterized by a correlation between population size or density and the mean individual fitness of a population or species (see [3,5,18,20]). Differential equation is a very powerful tool to explore problems, so a lot of differential equations models about Allee effect appear. Holmes E.E.(see [10]) studied a partial differential equations(PDE) in ecological which add Allee effect.Čiegis R. and Bugajev A. (see [4]) analysed a model of bacterial self-organization and gave the specific numerical approximation. Painter K.J. and Hillen T. (see [17]) obtained spatio-temporal chaos in a chemotaxis model. The stability analysis of such systems is a very important topic in theoretical and applied mathematics. For a long time, it has been recognized that delays can have very complicated impact on the system (see [2,19,23,29]). Yan J. (see [29]) discussed an impulsive delay differential equation with Allee effect, and obtained a periodic solution and its stability. Lingchong You, Robert Smith Randy et al. (for details, see [20]) engineered a gene circuit to confer a strong Allee effect in Escherichia coli and examined its impact on spread and survival. They used the LuxR/LuxI quorum-sensing (QS) system from Vibrio fischeri (see [14]) and the CcdA/CcdB toxinCantitoxin module to control population survival. The detailed operational principle of they circuit was described in [20], so we will leave it at that. Their circuit can be modeled by the following delayed differential equations: where C and [A] represent the bacterial density and the concentration of AHL (µM ),respectively. µ represents the maximum specific growth rate (h −1 ), k A represents the synthesis rate constant of AHL (µM h −1 ), k dA represents the degradation rate constant of AHL (h −1 ), τ represents the time delay of the activation of gene expression by the LuxRAHL complex (h), t represents simulation time (h), γ is a lumped term that represents the killing rate of CcdB (µM h −1 ), and β is a lumped term that represents the amount of CcdA leading to half-maximal killing rate of CcdB (µM ). Three cases were considered, with the circuit OFF ,ON and ON+ rescue. They declared that their protocol may be amenable for the study of Allee effects in natural systems, including Drosophila melanogaster and One could dispersed an established population of flies to new medium in a separate culture and examine reproductive success.
But in the whole paper, they did not discuss the effects of the time delay τ on the bacterial density. In our recent works [24,25], we discussed the stability and the Hopf bifurcation of some models with delay-dependent parameters. Moreover, it is well known that the delayed logistic differential equation is used to model the evolution of a single species x(t)(see [11,12,30]). Along this line, we established a new model as follows: where σ represents the time of bacterial split (h), K represents the maximum capacity of container for bacterial. In our paper, we investigate the dynamics near the equilibriums and hope to derive the perfect theoretical results. In fact, we not only discuss the Hopf bifurcation when τ or σ as a parameter, but also consider a codimension-2 bifurcation on a two-parameter plane, i.e. double Hopf bifurcation. For double Hopf bifurcation, we usually can observe lots of interesting and complicated dynamic behavior. Especially, strong resonant double Hopf bifurcation can present more complex nonlinear behavior because it may be a codimension-3 bifurcation.
In order to get a deeper insight into the double Hopf bifurcation analysis, we may apply and extend the methods and results in [6,7] to the case of two bifurcation parameters, and obtain the normal forms for double Hopf bifurcation without strong resonance. For strong resonant case, we will consider in another paper.
The rest of this paper is organized as follows. In Section 2, we discuss nonnegativeness and boundedness of solutions for the system (1.1). In Section 3, we consider the stability and the local Hopf bifurcation of the equilibriums when the time delay σ and τ as parameters respectively. In Section 4, the existences of double Hopf bifurcation is discussed. The normal form method and the center manifold theory are used to analyze the double Hopf bifurcation for system. In Section 5, we give some numerical examples to illustrate our results at the previous sections. In Section 6, we summarize our results.

Nonnegativeness and boundedness of solutions
We denote by X = C([−κ, 0], R 2 + ) the Banach space of continuous functions mapping the internal [−κ, 0] into R 2 + equipped with the sup-norm, where κ = max(σ, τ ). By the standard theory of functional differential equations, we know that for any φ ∈ X there exists a unique solution T (t, φ) = (C(t, φ), [A(t, φ)]) of the system (1.1), where T 0 = φ. The initial conditions are given by Theorem 1. Solution of system (1.1) corresponding to the above initial conditions remains nonnegative and ultimately bounded for all t ≥ 0.
Proof. Nonnegative: The proof of nonnegative for C(t) is similar to that of Lemma 2.1 of [28], and thus details are omitted here for the sake of brevity. For [A(t)], we assume the contrary, and let t > 0 be the first time such that [A( t)] = 0, then by the second equation of system (1.1), we have [A( t)] = k A C( t) ≥ 0. If C( t) = 0, then C(t) = 0, [A(t)] = 0 for all t ≥ t. We omit this trivial case. Hence, [A(t)] < 0 for t ∈ ( t − , t), where > 0 is sufficiently small, which contradicts [A(t)] > 0 for t < t. Ultimately bounded: From the first equation of (1.1), we obtain Hence, for t > τ we get that C(t) ≤ C(t − τ )e µτ , i.e., On substituting the above into (1.1), we obtain which clearly implies that, We complete the proof.
3 Stability of the equilibrium and existence of the Hopf bifurcation For the system (1.1), it is obvious that it has a trivial fixed point E 0 (0, 0) and two possible fixed points The characteristic equation of its corresponding linear system around the E 0 is given by Obviously, if γ β > µ, then E 0 is stable for all σ ≥ 0, τ ≥ 0. If γ β < µ, then E 0 (0, 0) is unstable for all σ ≥ 0, τ ≥ 0. By the translation u 1 (t) = The linearization of (3.1) at (u 1 , u 2 ) = (0, 0) is Thus, the characteristic equation of (3.2) is We can obtain the following results about the fixed point bifurcation.
Proof. By the hypothesis Kk A = k dA β, when µ = γ β , the (3.3) with σ = τ = 0 has a zero root and another negative root. Note that, if µ ≤ γ β , then there is only equilibrium E 0 , it is stable for µ ≤ γ β and unstable for µ > γ β . At the same time, on the side of µ > γ β , two new stable equilibria are created and given by E 1 , E 2 . Therefore the system (1.1) undergoes a pitchfork bifurcation at µ = γ β .
Here, we further explored the system for bifurcation analysis. This enabled us to find the system also undergoes a Hopf bifurcation for some σ or τ . By the above theorem, the system maybe undergoes a high co-dimensional bifurcations, i.e., Hopf-pitchfork bifurcation. For this bifurcation, the system may exhibit chaos via period-doubling bifurcations as the unfolding parameter values are far away from the critical point. That is, you can find the existence of highly irregular and chaotic-like dynamics in the system. We will discuss this situation in the next paper. In this paper, we only focus on Hopf and Double Hopf bifurcation.
In the rest of this section, we show the Hopf bifurcation of the system when the time delay is used as parameter.

Case σ = 0
We first consider the situation with no delay σ. If iω(ω > 0) is a root of (3.3) Separating the real and imaginary parts, we have It is easily to see that if either H > 0, then (3.5) has a unique positive root From (3.4), we find 3) and taking the derivative with respect to τ , we have which, together with (3.4), leads to

So we have
Notice that when τ = 0, the equation (3.3) becomes If b < ak dA , then the equation (3.8) has two roots with negative real parts. Otherwise, if b > ak dA , there is at least one root with positive real parts. Therefore, we can obtain the following lemma.
1. If H < 0, then the equation (3.3) has no root with positive real parts for all τ ≥ 0.
From the above lemma and (3.7) and the Hopf bifurcation theorem for functional differential equaitons (Theorem 1.1 in Chapter 10 of Hale and Lunel [8]), we have the following results on stability and bifurcation to system (1.1).
Remark 1. When τ = 0, the model is the same as [20]. But in [20], there is a lot of directions about the experimental phenomena, but no theoretical analysis.
In this paper, we analysis the dynamics behavior of the model in detailed, and give the reasons which cause the experimental phenomena from a mathematical point.
Proof. By direct calculation, one can easily obtain that .
The proof of that three conclusions is similarly.
Remark 3. By the widely used method which is based on the normal form theory and the center manifold theorem introduced by Hassard et al. [9].The direction, stability, and the period of the bifurcating periodic solutions can be obtain. On the other hand, from [13,26,27], we known that the properties of Hopf bifurcation of discrete schemes are the same as that of the corresponding delay differential equations, such as Runge-Kutta or strictly stable linear multistep method. In the last section, we will use some numerical simulations to illustrate the analytical results we obtained in previous sections and to show the properties of Hopf bifurcation, such as the direction of bifurcation and stability of periodic solutions.

The existence and normal form of the double Hopf bifurcation
In this section, for system (3.1), the existence of the double Hopf bifurcation from origin (0, 0) is obtained and the normal form will be also described. In the above section, we find that the system (3.1) undergoes a Hopf bifurcation at the origin (0, 0) when σ = σ ± j , j = 0, 1, . . .. Thus, a possible double Hopf bifurcation occurs when some σ 0 = σ + j = σ − l , where j, l = 0, 1, . . .. The equality σ = σ + j = σ − l implies that the linearized system on the trivial equilibrium has two pairs of purely imaginary eigenvalues ±iω − σ 0 and ±iω + σ 0 . In this paper, we assume ω − : ω + = n 1 : n 2 ∈ Z + . That is, we only consider the non-resonant double Hopf bifurcation.
We consider the value a and time delay σ as bifurcation parameters of system (3.1). Suppose system (3.1) undergoes a double Hopf bifurcation from the trivial equilibrium at the critical point: a = a 0 , σ = σ 0 . From Theorem 4, If H > 0, ∆ > 0 and (H 1 ) is satisfied. Then the equation (3.9) has two roots with positive and negative real parts when σ = 0, respectively. Therefore, we have reason to believe the existence of such possible situation: its roots except ±iω − σ 0 and ±iω + σ 0 have negative real parts at the first intersection point by the transverse condition. Thus, the center manifold can be used to describe the dynamics of the whole system in this case.
The system (3.1) may be transformed into a functional differential equation (FDE) in C = C([−1, 0], R 2 ), where the phase space C = C([−1, 0], R 2 ) as the banach space of continuous functions from [−1, 0] to R 2 with the supremum norm (see [8]) as where U (t) = (u 1 (t), u 2 (t)) T ∈ R 2 , U t ∈ C is defined by When α = (0, 0) T , then we have L 0 (φ) and F (0, φ) respectively by We set 4 . Then when α = (0, 0) T , the linearization equation at the trivial equilibrium of (3.1) is and the bilinear form on C * × C is Then, the phase space C is decomposed by ±iω − σ 0 , ±iω + σ 0 as C = P Q, where P is the center subspace spanned by the basis vectors associated with the imaginary characteristic roots and Q is the complement subspace of P .
Then the normal form of equation (4.1) on the center manifold of the origin near α = 0 up to quadratic order terms is given bẏ where g 1 2 (z, 0, α) is determined by g 1 2 (z, 0, α) = P roj (Im(M 1 2 )) c ×f 1 2 (z, 0, α) with f 1 2 (z, 0, α) is the function giving the quadratic terms in (z, α) for v = 0 defined by the first equation of (4.2). Thus, Neglecting the conjugate situation, we only consider z 1 and z 3 , the normal form (4.4) becomes that Next, we need calculate higher order terms of the normal form. Similarly, let M 1 3 denote the operator defined in V 4 3 (C 4 × Kerπ) with where V 4 3 denotes the linear space of the third order homogeneous polynomials in four variables (z 1 , z 2 , z 3 , z 4 ). When there is no strong resonance, it is easy to check that one may choose the decomposition V 4 3 = (ImM 1 3 ) (ImM 1 3 ) c with complementary space (ImM 1 3 ) c spanned by the elements z 2 1 z 2 e 1 , z 1 z 3 z 4 e 1 , z 1 z 2 2 e 2 , z 2 z 3 z 4 e 2 , z 1 z 2 z 3 e 3 , z 2 3 z 4 e 3 , z 1 z 2 z 4 e 4 , z 3 z 2 4 e 4 .

Numerical Simulations
In this section, we use some numerical simulations to illustrate the analytical results we obtained in previous sections. We first consider system (1.1) with µ = 0.4, K = 0.3, γ = 0.3, k dA = 0.09, β = 0.6, and k A = 0.5, σ = 0. Under these parameter values, system undergoes a Hopf bifurcation when τ . = 8.93635. By Theorem 3, the equilibrium E * is unstable for all τ ≥ 0.  As shown in Figure 1, when τ < τ 0 , the equilibrium E * is unstable and the other equilibrium E 0 is asymptotically stable. That is, when the killing module and rescue module are both on, for small specific growth rate µ < γ β , the the bacterial will become extinct at last. We choose we choose τ = 9.5 > τ 0 , Figure 2 shows that there a periodic solution with small amplitude exists. A strong Allee effect is observed, the specific growth rate is negative if the initial C is below the Allee threshold (see Figure 1). When the initial C is above the Allee threshold, although the time of rescue delay is large, the density of E. Coli. has oscillation near the threshold. Next, we change the values of γ to 0.05, i.e., the killing rate of the circuit is reduced. Hence, it follows that H < 0 and µ > γ β . Now, we know that E 0 becomes a saddle from sink and E * is a sink. Figure 3 shows that the stability of E * .  At last, we can give some analyses and examples to observe the behavior near the double Hopf bifurcation of system (4.1). For convenience, we assume that a is a independent parameter. Choosing µ = 0.4, K = 0.3, γ = 0.3, k dA = 0.01, β = 0.6, k A = 0.5, τ = 0, then b = 0.0008147257349. We can obtain σ 0 = 24.8210724, a 0 = 0.0714267. The unfolding parameters are Re(b 1 ) = 1.4246752743α 1 + 0.3245334987α 2 and Re(b 3 ) = 2.2157433961α 1 − 0.1126549270α 2 . we also obtain that Re(b 112 ) = 6.4350970712, Re(b 134 ) = 12.87019414, Re(b 123 ) = −7.300691906, Re(b 334 ) = −3.6503459527. Thus, the pdifficultq case happens under the parameters above. Let ρ 1 = Re(b 112 )r 2 1 , ρ 2 = −Re(b 334 )r 2 2 ,t = 2t, then we have the following planar system in terms of ρ 1 and ρ 2 : Note that N 0 = (0, 0) is always an equilibrium of (5.1). The other two semi-trivial equilibria given in terms of perturbation parameters are N 1 = (−Re(b 1 ), 0) and N 2 = (0, Re(b 3 )). There also exist a nontrivial equilibrium ). Since there does not exist unstable manifold containing the equilibrium, according to the center manifold theory, the solutions on the center manifold determine the asymptotic behavior of the solutions of the system (4.1). Therefore, if equation (5.1) has stable(unstable) semi-trivial equilibria, then (4.1) has stable(unstable) periodic solutions in the neighborhood of the trivial equilibrium. Since equation (5.1) has a nontrivial equilibrium N 3 , then (4.1) has a quasi-periodic solution in the neighborhood of (0, 0). So, we shall call the periodic solution or the quasi-periodic solution the source (respectively, saddle, sink) periodic solution of (4.1) if the semi-trivial equilibrium or the nontrivial equilibrium of (5.1) is a source (respectively, saddle, sink), respectively. The corresponding bifurcation diagram and phase portraits are shown in the parameter plane (α 1 , α 2 ) by Figure 4, where the critical bifurcation lines L 1 : α 2 = −4.389917466α 1 , L 2 : α 2 = 19.66841092α 1 , L 3 : α 2 = 127.1253084α 1 , L 4 : α 2 = −14.99637164α 1 , L 5 : α 2 = −71.46428697α 1 . We explain the bifurcations in the anticlockwise direction, starting from region 1 . First, in region 1 , there is only one trivial equilibrium which is a saddle. When the parameters are varied across the line L 1 from region 1 to 2 , N 0 becomes a sink, and anther semi-trivial equilibria N 1 (saddle) appears. That is, for the system (4.1), an periodic solution(saddle) appears due to Hopf bifurcation. In the region 3 , a more sink N 2 exists, i.e., another periodic solution(sink) appears, and N 0 becomes a saddle. When the parameters are varied from region 3 to region 4 , N 2 becomes a saddle, a new focus N 3 appears which corresponding a quasi-periodic solution of the system (4.1). When the parameters are varied from region 4 to region 5 , a limit cycle is present in the system (5.1). At this time, for system (4.1), a torus may be appears. When the parameters are further changed from region 5 to region 6 , the limit cycle disappears, N 3 collides with N 1 and then also disappears, and N 1 becomes a source. From region 6 to region 7 , N 1 disappears, and the N 0 becomes a source. When the parameters to region 1 , N 2 disappears at last. Moreover, a heteroclinic cycle may be formed between N 1 and N 2 , which is the main reason of the disappears of limit cycle.  To demonstrate the analytic results obtained in above, here we present some interest numerical simulation results. We choose there groups of perturbation parameter values: (α 1 , α 2 ) = (−0.001.0.021), (0.001, −0.075), (0.001, −0.050), belonging to the regions 3 , 4 and 5 , corresponding to a stable periodic solution as depicted in Figure 5, a stable quasi-periodic solution, see Figure 6, and a torus displayed in the three dimensional u 1 − u 2 − u 1 space, see Figure 7. It is clear that the numerical simulations agree very well with the analytical predictions.

Conclusions
In this paper, we have discussed the dynamical behaviors of a delay differential model which confers a strong Allee effect in Escherichia coli. Firstly, we have shown that the Hopf bifurcation exists as the time delay τ or σ crosses some critical values.
Secondly, we have considered the double Hopf bifurcation in system (4.1) with delay σ. We have obtained the vector field reduced to the center manifold, and derived the normal forms and their unfolding with perturbation parameters. Furthermore, we have given the bifurcation diagram as the unfolding parameter values are near the critical point. We improved the model in [20], and the research presented is a new step about our system in studied from the view point of high co-dimensional bifurcations. The study will help understanding the interpreting biological phenomena in theory.