{'fraction_non_alphanumeric': 0.04152087707021227, 'fraction_numerical': 0.016736645672964776, 'mean_word_length': 4.411486273272326, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 10, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Numerical Continuation and Bifurcation Analysis in a Harvested Predator-Prey Model with Time Delay using DDE-Biftool

Time delay has been incorporated in models to reflect certain physical or biological meaning. The theory of delay differential equations (DDEs), which has seen extensive growth in the last seventy years or so, can be used to examine the effects of time delay in the dynamical behavior of systems being considered. Numerical tools to study DDEs have played a significant role not only in illustrating theoretical results but also in discovering interesting dynamics of the model. DDE-Biftool, which is a Matlab package for numerical continuation and numerical bifurcation analysis of DDEs, is one of the most utilized and popular numerical tools for DDEs. In this paper, we present a guide to using the latest version of DDE-Biftool targeted to researchers who are new to the study of time delay systems. A short discussion of an example application, which is a harvested predator-prey model with a single discrete time delay, will be presented first. We then implement this example model in DDE-Biftool, pointing out features where beginners need to be cautious. We end with a comparison of our theoretical and numerical results.


Introduction
Time delay has been incorporated in models to reflect certain physical or biological meaning. Examples include optical feedback in laser systems [1,2,3], maturation age in stage structured population models [4], and delayed information in queueing models [10] just to name a few. The theory of delay differential equations (DDEs) [9,12], which has seen extensive growth in the last seventy years or so, can be used to examine the effects of time delay in the dynamical behavior of systems being considered. Numerical tools to study DDEs have played a significant role not only in illustrating theoretical results but also in discovering interesting dynamics of the model. DDE-Biftool [8,11], which is a Matlab package for numerical continuation and numerical bifurcation analysis of DDEs, is one of the most utilized and popular numerical tools for DDEs.
In this paper, we present a guide to using the latest version of DDE-Biftool targeted to researchers who are new to the study of time delay systems. A short theoretical discussion of the example, which is a harvested predator-prey model with a single discrete time delay, will be given in the next section. In section 3, we implement this example model in DDE-Biftool and compare the theoretical and numerical results. We conclude the paper with a summary and thoughts on using DDE-Biftool in studying time delay systems.

Harvested Predator-Prey Model with Time Delay
We first discuss the model that we are going to use as an example for numerical continuation and numerical bifurcation analysis. The theoretical results presented here will be compared to the numerical results obtained in the succeeding section.

The Model
We consider the following model, studied in [14], with a single discrete time delay parameter τ > 0 Here, x(t) and y(t) are the state variables denoting, respectively, the densities of the prey and predator populations at time t. The parameter a is the ratio of the intrinsic growth rate r of the prey and the carrying capacity K for the prey population in the absence of the predation. The rate of consumption of prey by the predator is given by the parameter b while c measures the conversion of prey consumed into the predator reproduction rate. The death rate of the predator is represented by d. Both species are assumed to have economic value and are harvested. The parameters h and k denote the harvesting rates of the prey and predator populations, respectively. All parameters in system (1) are positive real numbers.
The equilibrium solutions of system (1) are solutions (x(t), y(t)) of system (1) satisfying dx(t)/dt = 0 and dy(t)/dt = 0, and hence are obtained by solving for constant values x and y in the following system of nonlinear equations Since h > 0 and k > 0, the value of x nor the value of y can not be zero. From the equations in system (2), we get y = (rx − ax 2 − h)/bx and y = k/(cx − d).
The characteristic equation corresponding to the linear system (6) is where I is the 2 × 2 identity matrix. Equation (9) is obtained by using the ansatz X(t) = e λt x 0 y 0 to the linear system (6). If all roots of the characteristic equation (9) lie in the open left-half plane, i.e. Re λ < 0 for all roots λ of equation (9), then the equilibrium E * is locally asymptotically stable.
Using equations (7) and (8) and the functions f and g in (5), we obtain Thus, we can write the characteristic equation (9) as where When the time delay τ = 0, equation (10) reduces to the quadratic equation Both roots of equation (11) have negative real part if and only if Hence, for the case τ = 0, the equilibrium E * is locally asymptotically stable whenever conditions in (12) are satisfied.
We wanted to know if E * , under the conditions in (12) will become unstable as we vary the time delay parameter. Suppose that the conditions in (12) are satisfied and consider now the case where τ > 0. Initially, E * is locally asymptotically stable, i.e. all roots of the characteristic equation (10) with τ = 0 lie in the open left-half plane. If one or more roots of equation (10) cross the imaginary axis and move towards the open right-half plane as τ is increased, then E * will switch stability and becomes unstable. We have two possibilities: either a real root of equation (10) will cross the imaginary axis, i.e. λ = 0 is a root of equation (10) at some critical delay value, or a pair of complex conjugate roots of equation (10) cross the imaginary axis, i.e. λ = ±iω is a root of equation (10) at some critical delay value where ω is a nonzero real number.

Lemma 1. The number of positive roots of equation
where (u, h(u)) is the vertex of the parabola given by the graph of the function h(u) in equation (17).
Remark 1. The conditions in (3.13) of [14] are in fact the three inequalities in the last row of the table in Lemma 1, i.e.
where equation (17) has exactly two positive roots.

Transversality Conditions
We saw earlier that at τ = τ − k (resp. τ = τ + k ) for k = 0, 1, 2, 3, . . . , the characteristic equation (10) has purely imaginary roots ±iω − (resp. ±iω + ). We wanted to know if these roots along the imaginary axis will move towards the open righthalf plane or towards the open left-half plane. We address this by determining if the rate of change of the Re λ with respect to τ at the critical time delay values τ ± is positive or negative, where λ = λ(τ ) is a root of the characteristic equation (10).

Numerical Continuation and Bifurcation Analysis
Numerical tools to study DDEs have played a significant role not only in illustrating theoretical results but also in discovering interesting dynamics of the model. DDE-Biftool [8,11], which is a Matlab package for numerical continuation and numerical bifurcation analysis of DDEs, is one of the most utilized and popular numerical tools for DDEs. It was originally developed by K. Engelborghs [8] as part of his PhD work at the KU Leuven under supervision of D. Roose. DDE-Biftool provides a set of capabilities that is similar to what a range of alternative tools do for ordinary differential equations (ODEs) and maps, such as Matcont [6], COCO [5] and AUTO [7]. Another tool performing a similar set of tasks for DDEs, particularly for time-dependent DDEs with time-dependent delays, is Knut [13].
Aside from continuation of steady-state and periodic-orbit solutions which are typically done by varying a single parameter, DDE-Biftool can also continue bifurcations in two parameters. This includes steady-state folds, Hopf bifurcations, folds of periodic orbits, period doublings, and torus bifurcations. It can also perform normal form analysis for equilibria. DDE-Biftool is GNU Octave compatible and has a BSD licence such that it can be run completely as free software. The most recent version, DDE-Biftool v3.1.1, is maintained by J. Sieber and can be downloaded from https://sourceforge.net/projects/ddebiftool. The manual for this newest version of DDE-Biftool is provided at [11].
In this section, we illustrate the use of DDE-Biftool to perform numerical continuation and bifurcation analysis of system (1) varying the delay parameter τ . The boxed commands are the required commands and can be saved in a single m-file for convenience.  x(1,1).*x(2,1) -p(6); p(4)*x(1,1).*x(2,1) -p(5)*x(2,1)-p(7)]; funcs = set_funcs('sys_rhs',pp_sys,'sys_tau',@() [8]); We start with addpath which identifies the location of the folders ddebiftool and ddebiftool utilities containing the functions that we need for continuation and bifurcation analysis. This should be adjusted depending on where the user intends to do their computations and where the folders ddebiftool and ddebiftool utilities were copied.
We then set-up the function structure funcs identifying the previously defined system right-hand side pp sys as 'sys rhs'. The time delay τ will be the eighth parameter in our parameter list. Thus, assigning [8] in 'sys tau'. In parbd, we set the minimum bound, maximum bound, and maximum stepsize for the time delay τ , which is our main continuation and bifurcation parameter. Here, we choose τ to be from 0 to 15 because τ is non-negative and since we want to see the dynamics as we vary the time delay τ up until τ + 2 = 12.582660 (see Table 2).
Next, we set-up the branch of equilibria which we denote here by br. For 'parameter', we use the parameter values from Table 1. Note that the intial value set for τ , which is eighth in the parameter list, is 0.00. For 'x', we use the initial guess [6.00; 3.00] targeting the equilibrium E 1 with values given in equation (3). Our continuation parameter is τ , so 'contpar' is 8. After running the commands above, we get the following results. This means that our attempt to set-up a branch of equilibria is successful. The equilibrium branch br now contains two points. The first branch point in br has τ = 0 while the second branch point in br has τ = 0.02 since 'step' is assigned a value 0.02. The corrected value for E 1 can be obtained by typing the following in the Command Window.
>> format long; br.point(1).x ans = 6.061458241811056 3.254242134274988 The value of the time delay parameter τ for the first and second points in the equilibrium branch br are obtained by typing the following in the Command Window. We now continue the equilibrium branch br and then determine the stability of the continued branch. The above commands yield two plots of the equilibrium branch br as shown in Figure 1. The plot on the left panel of Figure 1 shows br with additional 300 points. So now br contains a total of 302 points which is achieved by using the function br contn. Meanwhile, the plot on the right panel of Figure 1 shows the same branch br but with stability information. This is obtained using the function br stabl. The stable and unstable parts of the branch are in green and red, respectively, while the Hopf bifurcation points are marked with asterisks ( * ). It is worth noting that in this particular example the stability switches occur at the Hopf bifurcation points. Here, there are five stability switches.
The function df measr gives the default xm and ym for the equilibrium branch br. We can check what these are by typing xm and ym in the Command Window.  Hence, for the plots in Figure 1, the horizontal axis is terms of the time delay parameter τ while the vertical axis is in terms of the state variable x(t). If you want the vertical axis to be in terms of y(t) instead of x(t), you need to type in ym.row=2 in the Command Window to make the desired change. Here, we keep the vertical axis in terms x(t) for future plots.
br.method.stability.minimal_real_part = -2; nunst = GetStability(br); ind_hopf = find(abs(diff(nunst))==2) The function GetStability yields nunst which is the number of characteristic roots in the open right-half plane. Hopf bifurcation occurs when a complex conjugate pair of simple characteristic roots crosses the imaginary axis. The above set of commands gives ind hopf which is the list of points in the equilibrium branch br where Hopf bifucation occurs. The value of the time delay parameter τ at the branch points given in ind hopf are obtained as follows.
>> for i = 1:length(ind_hopf) critical_tau(i) = br.point(ind_hopf(i)).parameter (8) Comparing with values in Table 2, the first three values above are correct up the first decimal digit while the last two are not. We have to keep in mind that the above values are mere approximations since the maximum stepsize that we assigned for the time delay parameter τ in parbd is just 0.05.  The above value is a correction to the initial guess 1.348598400000001 obtained from br.point(ind hopf (1)).parameter (8). Moreover, this corrected value matches the value of τ + 0 given in Table 2 Similarly, the correct value of τ for the next four Hopf bifurcation points are obtained by setting up Hopf branches. Here, we denote the next four Hopf branches as br hopf2, br hopf3, br hopf4, and br hopf5. The correct values of τ for the next four Hopf bifurcation points are obtained as follows.
>> [br_hopf2.point (1).parameter (8) These values match the values in Table 2 of τ − 1 , τ + 1 , and τ − 2 , but not that of τ + 2 . We remedy this by adding 'excludefreqs',br hopf4.point (1).omega in the previous code which removes undesired eigenvalues from consideration. The correct value for τ + 2 is now obtained by typing in the following in the Command Window. Branches of periodic solutions can be obtained from the identified Hopf bifurcations. Among these branches of periodic solutions, the branch emanating from the third Hopf bifurcation shows some interesting dynamics. We focus on this branch for the rest of this section.
[br_psol3,success] = SetupPsol(funcs,br,ind_hopf (3)) The function SetupPsol allows us to set-up a branch of periodic solutions. Here, we denote by br psol3 the branch of periodic solutions that emanates from the third Hopf bifurcation which as we saw earlier occurs approximately at the point ind hopf (3)  The above commands yield a plot of the periodic-solution branch br psol3 together with the equilibrium branch br as shown in Figure 2. Plotting the branch br psol3 requires defining xm psol and ym psol. Here, xm psol is the same as xm, i.e. the time delay parameter τ . For ym psol, we use the maximum function taking just the maximum value of the periodic x(t). Observe that br psol3 is initially stable and then it becomes unstable at a point marked with ( ). This stability switch actually occurs at a period-doubling bifurcation. We can get the value of τ where this period-doubling bifurcation occurs by setting up the branch of period-2 solutions per2 using the function DoublePsol. The value of τ where the period-doubling bifurcation occurs is obtained by typing in the following commands in the Command Window. We continue the branch of period-2 solutions per2 emanating from the perioddoubling bifurcation ( ) along the branch of periodic solutions br psol3. Figure  3 shows the plot this continued branch together with the periodic solutions branch br psol3. The period-2 solutions branch per2 also undergoes stability switch at a point marked with ( ). As before, this is also a period-doubling bifurcation. We leave it to the reader to verify that this second period-doubling bifurcation occurs at the value τ = 8.757752002502176. That is, beyond this value, we can expect a period-4 solution. Figure 4 shows the time series plots of x(t) and y(t) for different values of the time delay parameter τ showing period-1 solutions (τ = 7.10), period-2 solutions (τ = 8.60), and period-4 solutions (τ = 8.78).
We end this section with an exercise for the readers to check if the perioddoubling bifurcations will go on and will lead to a cascade of period-doubling bifurcations and eventually to chaos.

Conclusions
In this paper we presented how to use DDE-Biftool to obtain and analyze branches of solutions to a system of delay differential equations. We did this by revisiting the work of Toaha and Hassan [14] on a harvested predatory-prey model with time delay. We first discussed this model theoretically using a more simplified approach. Then, we implement the model in DDE-Biftool to study it numerically. The values of the time delay where Hopf bifurcations occur were numerically obtained and matches the theoretical results. In addition, branches of periodic solutions were also obtained using numerical continuation. At one of the branches of periodic solutions, some interesting dynamics occured. The occurrence of period-doubling bifurcations, which could lead to chaotic behavior of the system, was observed. Here, we emphasize the importantance of both theory and numerics in studying models. We hope that this paper served its purpose of introducing researchers to time delay systems and its implementation in DDE-Biftool which reveals more dynamical behavior of the model being considered.