Co-dimension two bifurcations analysis of a delayed tumor model with Allee effect

The mathematical model has become an important means to study tumor treatment and has developed with the discovery of medical phenomena. In this paper, we establish a delayed tumor model, in which the Allee effect is considered. Different from the previous similar tumor models, this model is mainly studied from the point of view of stability and co-dimension two bifurcations, and some nontrivial phenomena and conclusions are obtained. By calculation, there are at most two positive equilibria in the system, and their stability is investigated. Based on these, we find that the system undergoes Bautin bifurcation, zero-Hopf bifurcation, and Hopf–Hopf bifurcation with time delay and tumor growth rate as bifurcation parameters. The interesting thing is that there is a Zero-Hopf bifurcation, which is not common in tumor models, making abundant dynamic phenomena appear in the system. By using the bifurcation theory of functional differential equations, we calculate the normal form of these Co-dimension two bifurcations. Finally, with the aid of MATLAB package DDE-BIFTOOL, some numerical simulations have been performed to support our theoretical results. In particular, we obtain the bifurcation diagram of the system in the two parameter plane and divide its regions according to the bifurcation curves. Meanwhile, the phenomena of multistability and periodic coexistence of some regions can be also demonstrated. Combined with the simulation results, we can know that when the tumor growth rate and the delay of immune cell apoptosis are small, the tumor may tend to be stable, and vice versa.


Introduction
Cancer, also known as malignant tumor, grows fast and is very difficult to cure, which seriously affects people's life and health [1][2][3][4]. In recent years, with the vigorous development of biomathematics and interdisciplinary, it has become a new research trend to solve the related problems in biomedical field by establishing mathematical models. Therefore, many scholars have devoted themselves to the study of tumor models and made remarkable achievements, including tumor-immune and various types of tumor models, such as melanoma model, autoimmune disease model, hepatitis B model, and so on [5][6][7]. Most of them study the model from the perspective of stability and bifurcation, and consider the delay [8,9] and other effects in the model. For example, the authors in reference [10] stud-ied the symmetric-breaking bifurcation properties of a tumor model with delay and free boundary. In [11], a three-dimensional ordinary differential equation model is established to explore the growth law of tumor cells, and its bifurcation and multistability are studied. These models have different definitions of tumor growth, most of which are exponential growth and logistic growth. However, many scholars have pointed out that the Allee effect exists not only in the macro population entities but also in the growth of tumors [12][13][14][15][16][17]. Therefore, considering the delay of apoptosis of immune cells and the Allee effect of tumor cell growth, the following functional differential equation model is established: where x, y represent the number of human immune cells and tumor cells at time t respectively. α, b, and d stand for the growth rate, semi-saturation constant, and apoptosis rate of immune cells respectively, and the delay of apoptosis is expressed by τ . r, β represent the growth rate and apoptosis rate of tumor cells respectively. The maximum carrying capacity of tumor cells in human body is defined as k, and m is Allee threshold. In particular, inspired by [18][19][20], we use α y y+b instead of the usual immune cells growth term αy, whose growth is independent of tumor cells, and for the apoptosis of tumor cells, we follow the models in [19,21].
Although many scholars have studied models similar to (1), most of them stop at codimension one bifurcations, and there is less research on co-dimension two bifurcations. Co-dimension two bifurcations are usually accompanied by abundant dynamic phenomena, especially when there are multiple codimension two bifurcations in a model. Therefore, this paper focuses on some interesting dynamic phenomena of model (1) from the perspective of high co-dimension bifurcations.
The rest of this paper is organized as follows. In Sect. 2, we analyze the existence and stability of positive equilibria of system (1). Hopf bifurcation and Bautin bifurcation are investigated in Sect. 3, and the key values required for Bautin bifurcation analysis are calculated. We study the existence and normal form of zero-Hopf bifurcation in Sect. 4. In Sect. 5, we analyze another co-dimension two bifurcation, namely Hopf-Hopf bifurcation. To support our theoretical analysis, the two-parameter bifurcation diagram is obtained in Sect. 6, and given the bifurcation diagram, we give some numerical simulations for different parameters. Finally, a brief conclusion and expectation section completes the paper.

Existence and stability of positive equilibria
It is obvious that system (1) has a zero equilibrium E 0 = (0, 0) and no boundary equilibrium. Considering the positive equilibrium E * (x * , y * ), one satisfies After sorting out equation (2), we can get the following result: and for convenience, we give the derivative of (3) The two zeros of function (4) are Using some basic knowledge of quadratic function, we have the following.
Without loss of generality, this paper focuses on the stability of positive equilibrium E 2 . The linearization of system (1) at this equilibrium is actually where Thus, the characteristic equation of the system can be easily obtained as follows: When τ = 0, equation (6)  When τ > 0, suppose that λ 1,2 = ±iω are the roots of equation (6), and take it into the equation Separating the real and imaginary parts of (7), we can get By squaring both sides of (8) and adding the two equations, the following result can be obtained: Let v = ω 2 , and equation (9) is now and the corresponding two roots are expressed as It is easy to see that when η 1 -4η 2 > 0 and -η 1 2 > 0 hold, equation (10) has positive roots, and without loss of generality, suppose that v 1 and v 2 are positive roots. Accordingly, equation (9) has two positive solutions ω 1 = √ v 1 and ω 2 = √ v 2 , and the corresponding critical value of time delay is n }, the corresponding ω n is ω 0 . Let λ = ξ (τ ) + iω(τ ) be the root of equation (6), where ξ and ω satisfy ξ (τ 0 ) = 0 and ω(τ 0 ) = ω 0 , respectively. Differentiating both sides of the characteristic equation ( .

Bautin bifurcation
When the first Lyapunov coefficient is equal to 0, the Hopf bifurcation may degenerate and Bautin bifurcation occurs. In this section, Bautin bifurcation is investigated at positive equilibrium E 2 with τ and r as bifurcation parameters. Next, according to the research in [22][23][24], the properties of Hopf bifurcation and Bautin bifurcation are studied. Let (τ 0 , r 0 ) be the Bautin bifurcation point for bifurcation analysis. After scaling t → t τ and introducing μ = ττ 0 , ς = rr 0 , u 1 = x(t)x 2 , u 2 (t) = y(t)y 2 , system (1) can be written as where The above system (11) is transformed into a functional differential equation on the phase where For ϕ(θ ) = (ϕ 1 (θ ), ϕ 2 (θ )) ∈ C, we have

Zero-Hopf bifurcation
When the characteristic equation (6) has a simple zero root and a pair of pure imaginary roots and the other roots have negative real parts, the zero-Hopf bifurcation [25][26][27]will appear in system (1). With the emergence of zero-Hopf bifurcation, the dynamic phenomena of coexistence of periodic solutions and multistability are produced. Then, we study the existence of zero-Hopf bifurcation and calculate its normal form.

Existence of zero-Hopf bifurcation
To find the existence of zero-Hopf bifurcation, we need to study equation (6). Let For (15), we have the following.

Hopf-Hopf bifurcation
From the point of view of bifurcation diagram, when two Hopf bifurcation curves intersect, a Hopf-Hopf bifurcation may occur at these intersections. Theoretically, in this situation, the characteristic equation (6) has two pairs of simple pure imaginary eigenvalues, denoted as ±iω 1 and ±iω 2 , and the other eigenvalues have strictly negative real parts [28][29][30][31]. In the following, we use the methods in [32] to analyze the Hopf-Hopf bifurcation. Suppose that the eigenvectors corresponding to A(0, 0) and A * are p 1 , p 2 , p * 1 , p * 2 , respectively, and satisfy Continuing to use the disturbance μ = ττ 0 , ς = rr 0 , then, on the center manifold, we havė where g 1 and g 2 are in the form of The third-order normal form near a Hopf-Hopf point is given by Similarly, the key quantities to determine the bifurcation are Re e 2100 .
Remark 5.1 Here, we do not give the specific expressions of e 2100 , e 1011 and e 1100 , e 0021 in (21), and the calculation process can refer to the introduction of Hopf-Hopf bifurcation in reference [32].

Simulations
In this section, we perform some simulations to support the theoretical results about codimension two bifurcations in system (1). The simulations are divided into two parts, one for Bautin and zero-Hopf bifurcations, and the other for the Hopf-Hopf bifurcation. Moreover, the phenomena of multistability and periodic coexistence are also simulated.
From the diagram, we can see that there are two Bautin points and one zero-Hopf point of system (1), corresponding to GH 1 , GH 2 , and ZH respectively. In fact, for GH 1 , some calculations demonstrate τ = 4.84628, r = 0.10548, and the first Lyapunov coefficient Re l 1 (0, 0) ≈ 0, the second Lyapunov coefficient Re l 2 (0, 0) = 0.019362 > 0. Similarly, for GH 2 , we have τ = 5.65348, r = 0.10914 and Re l 1 (0, 0) ≈ 0, Re l 2 (0, 0) = -0.0091855 > 0, and at the top right of the point, the Hopf bifurcation is supercritical, and the bottom left of this point, the Hopf bifurcation is subcritical. At ZH, τ = 5.0408, r = 0.09749, it is easy to verify that the content of Lemma 4.2 can be satisfied; moreover, we have s = 0.00062874, θ = -0.6416. Meanwhile the limit point of cycles, marked green, originating from Bautin bifurcation and the torus bifurcation, marked cyan, originating from the zero-Hopf bifurcation are also drawn. Due to the appearance of these bifurcation curves, the τr plane is divided into 11 regions, as shown in Fig. 1, labeled I-XI. Next, we take different values of (τ , r) in these regions for numerical simulations. Using Matlab function DDE23, the simulation results are shown in Fig. 2, Fig. 3, and Fig. 4, and the black dot in the graph represents the initial value of the solution. Specifically, there are a stable periodic solution for the system in region I and coexistence of two unstable periodic solutions in region IX(see Fig. 2). Figure 3 shows the coexistence of two unstable periodic solutions in region II, one of which is eventually stable to the origin, which reveals that immune cells and tumor cells go to apoptosis together, an unstable periodic solution and an asymptotically stable solution in region III, and shows two opposite phenomena, one of which is asymptotically stable, and the other is unstable, corresponding to regions IV and X, respectively. An asymptotically stable solution and an unstable periodic solution are obtained in region V, an asymptotically stable solution appearing in region VII, and an unstable periodic solution and an asymptotically stable solution showing in region VII, two stable periodic solutions coexisting in region VIII(see Fig. 4). Finally, for region XI, after calculation, we find that there is no positive equilibrium in this region, so the simulation is not carried out.
In a conclusion, in general, the larger the values of τ and r, that is, the closer the values of (τ , r) in Fig. 1 to the upper right regions, the more likely the solutions to be unstable, such as regions III, IV and V, VII(Although there may be stable solutions in a small range of positive equilibrium). On the contrary, when τ and r are small, the solutions are asymptotically stable or periodically stable, such as regions I, II, and VIII. Therefore, if the tumor cell proliferation rate and the delay of immune cell apoptosis are small, it is conducive to the tumor in a stable state.

Simulation of Hopf-Hopf bifurcation
For Hopf-Hopf bifurcation, we simply give a numerical example to illustrate the existence of Hopf-Hopf bifurcation and draw a bifurcation diagram near the bifurcation point, as In fact, system (1) has more than one Hopf-Hopf bifurcation. With the increase of time delay, more Hopf-Hopf points appear in the two-parameter plane τr, so more complex phenomena of periodic solutions and chaos will exist in the system.

Conclusion and discussion
Co-dimension two bifurcations are a common bifurcation phenomenon in delay differential equations, which leads to complex dynamic behavior of the system. In this paper, a delayed tumor model is established, and the Allee effect is considered in this model. It is found that there are at most two positive equilibria in the system, and the stability condition of the equilibria is given. Due to the time delay and Allee effect, the system has abundant co-dimension two bifurcations, accompanied by multistability and periodic coexistence. Particularly, after adding Allee effect, the system changes from a single equilibrium to two different positive equilibria, which is an important reason for the appearance of coexistence, and the time-delay also makes the original two-dimensional ordinary differential equation more complex, such as the emergence of zero-Hopf bifurcation and Hopf-Hopf bifurcation. We calculate the normal form of Bautin bifurcation, zero-Hopf bifurcation, and Hopf-Hopf bifurcation of the system. For Bautin bifurcation and zero-Hopf bifurcation, the two-parameter bifurcation diagram is also given in numerical experiments in Sect. 6. By the simulation of different regions, the multistability and periodic coexistence are demonstrated. For Hopf-Hopf bifurcation, we only perform a numerical example. In fact, there are stable periodic solutions and chaos near the Hopf-Hopf bifurcation point.
To conclude, this paper analyzes the stability and bifurcations of a tumor model. The theoretical and simulation results are of biological significance, which is helpful for the future tumor research at the theoretical level. More specifically, the time-delay and growth rate of tumor cells play a key role in the development of tumor, which is conducive to our further research in the future.