Choristoneura Fumiferana is perilous defoliators of forest lands in North America and many countries in Europe. In this study, we consider mathematical models in ecology, epidemiology and bifurcation studies; the spruce budworm model and the population model with harvesting. The study is designed based on bifurcation analysis. In particular, the results support population thresholds necessary for survival in certain cases. In a series of numerical examples, the outcomes are presented graphically to compare with bifurcation results.

Keywords: Bifurcation analysis,; competition, stability analysis, Spruce budworm, numerical illustration, AMS Subject Classification 2010: 37N25, 49J15, 92D30

Bifurcation analysis is a powerful tool as it allows one to identify in a systematic way where dynamics of interest exist in parameter space. Regions of stable solutions can be determined by computing stability boundaries directly as functions of the relevant system parameters. Sudden qualitative behavior changes in dynamical systems arise as a parameter is changed, which is known as bifurcation analysis.

Spruce budworms (Choristoneura fumiferana) are dangerous defoliators of forest lands in North America.¹ These forests known as coniferous forests have pines, spruces, firs, and larches and the spruce budworm at larvae stage start feeding on needles of these trees. Fully grown insects can go through a tree in a time span of 4-5 years since they eat photosynthetic tissue which is crucial for the development of trees. The population usually stays at a low level in and less forestland because of their avian predators. But every three-four decades there is an outbreak of budworms and it results in a completely destroyed forest. It is assumed by scientists that this phenomenon has been occurring for more than thousands of years. Some papers that have studied bifurcation on various ecological models.^1-4

The dynamics of one dimensional ordinary differential equations are not very diverse, like solutions either converge to a stable equilibrium or diverge. Geometric approaches give better insight on these models. As control parameter is varied bifurcations like saddle node bifurcations, transcritical bifurcation and pitchfork bifurcation occur. When two fixed points annihilate each other or reappear as a parameter is changed it is known as saddle node bifurcation. Transcirital bifurcation is when a steady equilibrium changes its stability. And pitchfork bifurcation works with symmetry, fixed points disappear or appear in symmetric pairs.

The dispersal dynamics of prey-predator, competition, and mutualism is explored in.⁵ A very recent study of diffusion reaction in competition model is.⁶ There are also other studies^7-9 related to diffusion where the effect of dispersion procedures on competitive populations are explored. In another study, the relation between growth rates, carrying capacity, and diffusion strategies is discussed, it is analyzed in⁹ whether coexistence is promoted when diffusion strategies are diverse. Here bifurcation analysis will be used to understand why this happens thus giving a mathematical point of view on this phenomenon.

In this paper, we have an elementary discussion on bifurcations in one-dimensional ordinary differential equations. Graphical methods are used to identify equilibrium points, determine their stability, and to see the behaviors of the numerical solutions near fixed points.^10-12 In section 2, the general forms of bifurcations in autonomous differential equations such as Saddle node bifurcation, transcritical bifurcation, and pitchfork bifurcation will be analyzed. Following this discussion, in section 3 the mathematical model for spruce budworm for single species pray is explored. We talk thoroughly about the function for the rate of predation. In section 4 a dimensionless model is going to be introduced for further analysis of the model. Further into the paper, in section 5 we look into bifurcations that occur in the spruce budworm model when change of parameter values occur. In section 6 we studied harvesting in population models, we took different levels of harvesting into consideration, leading us to the conclusion.

Bifurcation Preliminaries

In this section we have discussed about some common bifurcation behavior.

Pitchfork Bifurcation

Consider the following equation for the classic form of the Pitchfork Bifurcation in one dimensional
form:

$\frac{du}{dt}=\alpha u-u^3$      (1.1)

Pitchfork bifurcation is a local type of bifurcation that deals with symmetry, where equilibrium points show up or dissipate in pairs.^14,15 In Figure 1 we see the vector field of equation (1.1) for multiple values of α. First, when $\alpha =-4$ and $\alpha =0$ there is only one fixed point, and it stable for both times. But when the value of α is increased to  two new equilibrium points appear, and the stability at $u = 0$ changes. The new two fixed points are stable, and this will be true for all $\alpha > 0$  as the fixed points will be symmetrically positioned at $u=\pm \sqrt{\alpha }$ .

Figure 1 Pitchfork bifurcation.

The stable lines in Figure 2 are given by $u = 0$  when $\alpha <0$ and $u={\alpha }^3$ . The origin becomes unstable for α>0 and that is shown using the dotted line. Some other bifurcations are discussed and can be found in Appendix.

Figure 2 Bifurcation diagram for Pitchfork bifurcation.

Spruce Budworm mathematical model

In this study, we propose the following spruce budworm model^3,10 for a single species prey:

$\frac{du}{dt}=ru\left(1-\frac{u}{K}\right)-\frac{\beta u^2}{{\alpha }^2+u^2}$      (2.1)

where r  is the intrinsic growth rate of the spruce budworm. $K$  is the carrying capacity, which
depends on the availability of the foliage of the balsam firs. The rate of predation

$h\left(u\right)=\frac{\beta u^2}{{\alpha }^2+u^2}$      (2.2)

is from avian predators, which associate reward with prey, a learned behavior. They focus on their best sources of prey, which allows a low density prey to escape notice. The functional form of this type of predation follows the Holling’s Type III S-shaped response as defined by $h\left(u\right)$ .^3,10,14  represents the saturation level of the predator, e.g., constant population can eat so much prey, this saturation level works like an upper limit on the predation. α determines the densities of spruce budworm that cause avian predation. Now if taking the limiting value of the predation function,

$\underset{u\to \infty }{\lim }h\left(u\right)=\underset{u\to \infty }{\text{lim}}\frac{\beta u^2}{{\alpha }^2+u^2}=\underset{u\to \infty }{\text{lim}}\frac{2\beta u}{2u}=\beta$

From Figure 3 it is seen that the function is monotonically increasing until it reaches its asymptote at the value β.³ This indicates that no matter how much the budworm population increases the birds can eat a limited quantity. Using basic calculus, the first and second derivative is calculated to find the point of inflection. This point determines when the budworm population increases rapidly.

Figure 3 Rate of Predation.

Differentiating equation (2.2) with respect to u,

$h^{\prime }\left(u\right)=\frac{\left({\alpha }^2+u^2\right)+2\beta u-\beta u^2*\left(2u\right)}{{\left({\alpha }^2+u^2\right)}^2}$

$=\frac{2{\alpha }^2\beta u}{{\left({\alpha }^2+u^2\right)}^2}$      (2.3)

It is obvious that h^' (u)>0  since all the constants are positive. Now again differentiating with
respect to ,

$h\text{'}\text{'}\left(u\right)=\frac{{\left({\alpha }^2+u^2\right)}^2*2{\alpha }^2\beta -2{\alpha }^2\beta u*2\left({\alpha }^2+u^2\right)*\left(2u\right)}{{\left({\alpha }^2+u^2\right)}^4}$

$=\frac{2{\alpha }^2\beta \left({\alpha }^2-3u^2\right)}{{\left({\alpha }^2+u^2\right)}^3}$      (2.4)

Substituting equation (2.4) to zero, yields,

$\frac{2{\alpha }^2\beta \left({\alpha }^2-3u^2\right)}{{\left({\alpha }^2+u^2\right)}^3}=0$

$\Rightarrow u=\pm \frac{\alpha }{\sqrt{3}}$      (2.5)

Population cannot be negative so the negative is omitted. Thus predation dramatically climbs upwards to $\frac{\alpha }{\sqrt{3}}$  from below. After crossing the point of inflection, the function gradually approaches its asymptote at $h=\beta$ .

Dimensionless Model

In order to simplify the differential equation and analyze the behavior better we reduce the number of parameters. Let,

$\rho =mU,\lambda =nt$

Substituting the above values in (2.1),

$\frac{d\rho }{d\lambda }=\frac{m}{n}\left[ru\left(1-\frac{u}{K}\right)-\frac{\beta u^2}{{\alpha }^2+u^2}\right]$

$\Rightarrow \frac{d\rho }{d\lambda }=\frac{m}{n}\left[\frac{r\rho }{s}\left(1-\frac{\rho }{Ks}\right)-\frac{\beta {\rho }^2{s}^2}{{\alpha }^2+{\rho }^2{s}^2}\right]$      (2.6)

If $m=\frac{1}{\alpha }$ , from equation we have,

$\frac{d\rho }{d\lambda }=\frac{r\rho }{n}\left[\left(1-\frac{\rho \alpha }{K}\right)-\frac{\beta }{\alpha n}\left(\frac{{\rho }^2{\alpha }^2}{{\alpha }^2+{\rho }^2{\alpha }^2}\right)\right]$

$=\frac{r\rho }{n}\left[\left(1-\frac{\rho \alpha }{K}\right)-\frac{\beta }{\alpha n}\left(\frac{{\rho }^2}{1+{\rho }^2}\right)\right]$      (2.7)

Assume that $n=\frac{\beta }{\alpha }$ , then from equation (2.7) we get,

$\frac{d\rho }{d\lambda }=\frac{r\rho \alpha }{\beta }\left[\left(1-\frac{\rho \alpha }{K}\right)-\left(\frac{{\rho }^2}{1+{\rho }^2}\right)\right]$

Defining the new scaled parameters $R=\frac{\alpha r}{\beta }$ and $C=\frac{K}{\alpha }$ . Then the dimensionless model is given by,

$\frac{d\rho }{d\lambda }=R\rho \left[\left(1-\frac{\rho }{C}\right)-\left(\frac{{\rho }^2}{1+{\rho }^2}\right)\right]$      (2.8)

To find the equilibrium points, let $\frac{d\rho }{d\lambda }=0$ , then,

$R\rho \left[\left(1-\frac{\rho }{C}\right)-\left(\frac{{\rho }^2}{1+{\rho }^2}\right)\right]=0$      (2.9)

Consider the function, $f\left(\rho \right)=R\rho \left(1-\frac{\rho }{C}\right)$  and $g\left(\rho \right)=\left(\frac{{\rho }^2}{1+{\rho }^2}\right)$ . Thus $\rho =0$ is a trivial equilibrium point. In order to find the other equilibrium points we will graph $f\left(\rho \right)$ and $g\left(\rho \right)$ and look at the intersection points. Here, $f\left(\rho \right)$ =Per capita growth rate of the scaled population variable $\rho$ with respect to scaled time $\lambda = \frac{\beta t}{\alpha }$ =α.$g\left(\rho \right)$ =Per capita death rate of the spruce budworm due to avian predation.

Bifurcation Analysis

Qualitative changes in dynamical systems, like fixed points being created and destroyed or their stability changing is called bifurcation. In one dimensional systems these alterations occur when a parameter value is changed. Here we will see a practical example of saddle node bifurcation. We will plot $f\left(\rho \right)$ and $g\left(\rho \right)$ for different values of the parameters involved in the equation and look for fixed points. The stability of a fixed point is determined by reading the direction of the flow on $x-$ axis. If $f\left(\rho \right)>g(\rho )$  then the flow is to the right otherwise it is directed to the left. Now if the flow is inwards then the point is stable and if the flow is outwards then the point is unstable.¹⁰ We will also plot the numerical solutions of the autonomous equation. First we start with small values of $R = 0.22$  and $C = 3$ . From the Figure 4 it is seen that there is only one fixed point at ${\rho }^{*}= 0.21366$ and it is stable. On the right figure the numerical solutions for equation  from various initial values and all of the solution gradually converge towards equilibrium solution. At this stage the budworm population stays limited.

Figure 4 For R=0.22 and C=3.

Here in Figure 5 the value of  is increased 0.8 to and $C =8$ . Despite this huge increment it is seen from the Figure 5 that the number of fixed points remain the same. This fixed point is also stable, the value is very low and the budworm population remains endemic.

Figure 5 For R=0.8 and C=8.

Figure 6 For R=0:22 and C=12.

We will increase C  to a higher value of  and change R  values. Since the parameter C  is acquainted with the carrying capacity K, holding C  at a higher value will help us analyze the model when the forest is large. Starting at $R=0.22$  we see that the results are quite similar to the ones obtained at Figure 4.

Next in Figure 7 we can see that when we raise R to 0.32  there is one fixed point at $\rho *=0.34843$ and the function is $f\left(\rho \right)$  tangent to $g\left(\rho \right)$ at 5.82579. A saddle-node bifurcation is thus occurring.

Figure 7 For R=0.32 and C=12.

As R  continues to escalate to 0.45 we see from Figure 8 that two new fixed points have emerged at 2.32714 and 9.10661. Here the point 2.32714 is unstable since at the left of the point $f\left(\rho \right) < g\left(\rho \right)$ and the flow is outwards. On the right $f\left(\rho \right)>g(\rho )$  and the flow is inwards.¹⁰ The following fixed point is of course stable. Since two stable points or unstable points cannot occur side by side. As for the numerical solutions the solution curves remain close to the stable equilibria. Values that are larger than the unstable fixed point approaches the larger stable fixed point at ${\rho }^{*}=9.10661$  and values less than it approach the stable point occurring at a low value.

Figure 8 For R=0.45 and C=12.

In figure 9 we notice that a new saddle-node bifurcation arises. The two points 0.56624 and 2.32714 gradually approach each other and collide into a new fixed point at 9.78531, as R is raised to 0.548. The fixed point at 0.49743 is stable, but at a very high value.¹⁵ Which causes an outbreak in the spruce budworm population.

Figure 9 For R=0.548 and C=12.

Figure 10 For R=0.6 and C=12.

Now at $R=0.6$ and the saddle node bifurcation disappears and there is only one stable fixed point at 10.02456.

Population models with harvesting

In this section we will analyze the behavior of the logistic population model when population is removed at some defined rate, more commonly known as harvesting.¹⁶ The harvesting function is defined as , dependent on time $t$ .

Calibration of the model with constant harvesting
Let us consider the following model,

$\frac{du}{dt}=au\left(1-\frac{u}{K}\right)-H$      (3.1)

Equation (3.1) is a population model with constant-yield harvesting. We would like to study this model and see how different harvesting levels affect the model. Some real-life examples are harvesting animals like deers, birds from sanctuaries, or even fishes, for example whales. In these sort of scenarios harvesting is allowed to some limited extent.

In Figure 11, we have plotted $F \left(u, H\right)=1.6u\left(1-\frac{u}{K}\right)-H$  for different values of H. Now solving eqaution (3.1) analytically yields a very complicated structure, so we approach a graphical method to obtain the steady states. Solving the quadratic equation,

Figure 11 F(u,H) for different values of H.

$au\left(1-\frac{u}{K}\right)-H=0$

$\Rightarrow a{u}^2-au+kH=0$      (3.2)

 We get $u=\frac{ka\pm \sqrt{{\left(ka\right)}^2-4aHk}}{2a}$ . Now when ${\left(ka\right)}^2-4aHk>0$ that is, $H<\frac{ka}{4}$ the quadratic equation (3.2) has two real roots, giving two fixed points for the model. H=0  and H=2.6 represents this case. When H=0 the model is constituted as the logistic growth model. Now when H=2.6  there are two fixed points, the one at $u=1.94$ is unstable and the latter fixed point at $u=10.06$ is stable (Figure 12). For initial values, when $1.94<u<10.06$ the solutions move away from the unstable equilibrium, known as the extinction equilibrium and approach the solution curve near $u=10.06$ . But for initial values less than $u=1.94$ the solutions approach negative infinity, thus the population declines.

Figure 12 Numerical solutions for different values of H.

Next when H=4.8 a saddle node bifurcation occurs, two fixed points collide into one at $u=6$ . This fixed point is half stable. The numerical solutions are observed in the second plot Figure 12 shows us that for initial values less than  the solutions diverge, but when they are greater than  they move towards the equilibrium. Gradually as H is increased, that is when $H >\frac{Ka}{4}$ all the fixed points disappear and the solutions diverge. Thus for harvesting, $H>4.8$ the population fails to survive.

Calibration of the model with proportional harvesting

Here, we now consider the following equation,

$\frac{du}{dt}=au\left(1-\frac{u}{K}\right)-uH=G\left(u,H\right)$      (3.3)

In Figure 13 we have plotted $G\left(u,H\right)=1.6u\left(1-\frac{u}{K}\right)-uH$  for various values of .

Figure 13 G(u; H) for different values of H.

Solving equation (3.3) for equilibrium points we get,

$au\left(1-\frac{u}{K}\right)-uH=0$

$\Rightarrow u\left[au+\left(aK+KH\right)\right]=0$

From equation (6.4) we have u = 0 is always an equilibrium point and the other one is $u=\frac{K\left(a-H\right)}{a}$ . When $H=0.3$  we have two fixed points. The extinction equilibrium at $u=0$ and $u=9.75$ . From the first plot in Figure 14 we observe that $u=9.75$  is the stable
equilibrium. This case is true for all $\frac{K\left(a-H\right)}{a}>0$ , that is when $0<H<a$  here $a=1.6$ .

Figure 14 Numerical solutions for different values of H.

When H  is increased to 1.6 transcritical bifurcation occurs and $u=0$  is the only fixed point which is half stable. For $H>a$  or in this case $H>1.6$ , the extinction equilibrium becomes stable and the other fixed point is negative. Thus when $H>a$ , harvesting is too intense and it surpasses population growth.

In this paper, we have studied the local behavior of bifurcations of one- dimensional ordinary differential equations to develop a proper understanding of the topic. We obtained the stable and unstable regions of equilibrium for saddle node bifurcation, transcritical bifurcation, and pitchfork bifurcation. Alongside we analyzed the bifurcation diagrams, where the fixed points dependence on parameter values is represented.

The spruce budworm problem has been analyzed here using bifurcation and numerical illustrations have been shown to justify the results. When the values of R and C are low the population remains endemic. But for higher C values saddle node bifurcation occurs and there are changes noticed in the stability of the system. For R values above 0.45 fixed stable points are found at higher values which lead to outbreaks in the budworm population.

We also studied the dynamics of population models with constant and proportional harvesting. These two cases exhibit examples of saddle node bifurcation and transcritical bifurcation. But the model for constant harvesting is not suitable for lower levels of population. When harvesting exceeds a certain level the population fails to survive. These changes in the qualitative structure of the dynamics of the models shows us that there is much scope for further study.

The author M. Kamrujjaman research was partially supported by University Grants Commission (UGC), 2019-2020, Bangladesh.

The authors declare no conflict of interest exists.

1. C Colona, D Claessen, M Ghil. Bifurcation analysis of an agent-based model for predator-prey interactions. Ecological Modelling. 2014;317:93–106.
2. D Wu, H Zhang. Bifurcation analysis of a two-species competitive discrete model of plankton allelopathy. Advances in Difference Equations. 2014;70:1–10.
3. White A, Huang. Spruce Budworm Population Model in ODE. 2012.
4. S Saha, G Samanta. Modeling of Insect-Pathogen Dynamics with Biological Control. Mathematical Biology and Bioinformatics. 2020;15(2):268–294.
5. L Korobenko, Md Kamrujjaman, E Braverman. Persistence and extinction in spatial models with a carrying capacity riven diffusion and harvesting. Journal of Mathematical Analysis and Applications. 2013;399:352–368.
6. Md Kamrujjaman. Dispersal Dynamics: Competitive Symbiotic and Predator-Prey Interactions. Journal of Advanced Mathematics and Applications. 2017;6:1–11.
7. Md. Kamrujjaman. Interplay of resource distributions and diffusion strategies for spatially heterogeneous populations. Journal of Mathematical Modeling. 2019;7(2):175–198.
8. E Braverman, Md Kamrujjaman. Competitive –cooperative models with various diffusion strategies. Computers and Mathematics with Applications. 2016;72:653–662.
9. E Braverman, Md Kamrujjaman, L Korobenko. Competitive spatially distributed population dynamics models: Does diversity in diffusion strategies promote coexistence?, Mathematical Biosciences. 2015;264:63–73.
10. SH Strogatz. Nonlinear Dynamics and Chaos. First edition, CRC Press, 2000.
11. DG Zill, A First Course in Differential Equations with Modeling Applications. Tenth Edition, Richard Stratton, 2013.
12. L Edelstein-Keshet. Mathematical Models in Biology. First edition, SIAM: Society for Industrial and Applied Mathematics, 2005.
13. LJS Allen. An introduction to mathematical biology. Pearson, 2006.
14. JD Murray. Mathematical Biology. I and II, third edition, Springer, New York, 2002.
15. CS Chou, A Friedman. Introduction to Mathematical Biology. Modeling, Analysis, and Simulations (Springer Undergraduate Texts in Mathematics and Technology, First edition, Springer, 2016.
16. F Brauer, C Castillo-Chavez. Mathematical Models in Population Biology and Epidemiology. Second Edition, Springer, 2010.