Steady states and outbreaks of two-phase nonlinear age-structured model of population dynamics with discrete time delay

ABSTRACT In this paper we study the nonlinear age-structured model of a polycyclic two-phase population dynamics including delayed effect of population density growth on the mortality. Both phases are modelled as a system of initial boundary values problem for semi-linear transport equation with delay and initial problem for nonlinear delay ODE. The obtained system is studied both theoretically and numerically. Three different regimes of population dynamics for asymptotically stable states of autonomous systems are obtained in numerical experiments for the different initial values of population density. The quasi-periodical travelling wave solutions are studied numerically for the autonomous system with the different values of time delays and for the system with oscillating death rate and birth modulus. In both cases it is observed three types of travelling wave solutions: harmonic oscillations, pulse sequence and single pulse.


Introduction
Age-structured models are an important and useful mathematical tool for describing and studying the processes of population life cycles: birth, aging, proliferation and death. This class of models was presented at first in the mathematical theory of epidemics in the pioneer works of Kermack and McKendrick [35][36][37], theory of population of Von Foerster [51] and nonlinear theory of population dynamics in more late papers of Gurtin and Mac-Camy [27,28]. The main feature of age-structured model is a non-local boundary condition which describes the fertility process in population in explicit form. As the reproduction of individuals in populations plays an important role in nature the age-structured models are widely studied and used in theoretical and applied researches of biological, demographical and ecological systems [9,12,13,16,17,22,23,26,29,[31][32][33][42][43][44]53]. This class of models is well adapted to simulation of the dynamical properties of experimental cultures of cell aggregations including the cell-division cycle with variable intermitotic time. The total cell population in age-structured models is partitioned into proliferating and quiescent phases [11,14,25,30]. In [25] authors define the age of a biological cell as the time elapsed from its last mitosis. Starting from this definition, they derive a linear age-structured model giving the dynamics of populations of proliferating and quiescent cells. The comparison of numerical solutions in [25] to the experimental observations presented in [50] approves the efficiency of age-structured models in simulation of cell population dynamics. In this paper we consider the nonlinear age-structured population model which generalizes the model mentioned above in the following ways: (i) Nonlinear death rate with density independent term and density-dependent term (density feedback loop, according to the common Gurtin and MacCamy's approach [27,28]. (ii) Discrete time delay of the density-dependent death rate [20,39,40,47]. (iii) Proliferation phase and quiescent phase running in parallel [25].
These three aspects provide for more realistic and accurate modelling of cell population and populations of other biological organisms, particularly regarding their long-term dynamics. For instance, linear models lead to either unlimited growth or unlimited decay. Due to the density feedback loop used here the fertility and the mortality can balance out leading to predictable long-term behaviour, for example, asymptotically stable solutions of autonomous system, quasi-periodic solutions (including population outbreaks) of system with big time delay or with oscillating coefficients of equations.
Section 2 provides derivation of the model. The age distribution of the polycyclic population in proliferation phase is modelled via nonlinear transport equation with delay and integral boundary condition. The age of the individuals in the quiescent phase does not play a major role since there is no proliferation in this phase. Hence, the size of the population in this phase rather than its age distribution is considered and represented via a delay ordinary differential equation. This system is subjected to theoretical analysis and numerical investigation in the subsequent sections. Since we consider the well-known class of nonlinear age-structured model and use the verified mathematical methods of logical deduction and deriving of numerical recurrent algorithm for solving the above problem, a thorough study of the existence and uniqueness of continuous travelling wave solution are beyond the scope of this paper. We refer readers to the corresponding works for the age-structured problems without and with time delay [3,5,6,18,20,24,27,40,52]. Our objective here is to focus the attention of specialists in life science and applied mathematics to the obtained in the series of simulations new dynamical properties of nonlinear agestructured models with time delay for the different scenarios and regimes of population dynamics.
There are two main approaches in the numerical simulations of population dynamics for age-and size-structured models with time delay. The first one is numerical approach. It uses the theory of difference schemes or finite element method and based on the difference approximation of initial-boundary value problems for ODE, PDE [2,4,10,49]. As a result, the initial differential problem is reduced to the difference problem and can be solved by means of the methods of computational mathematics. The strength of such approach is that it provides the possibility of simulation of complex nonlinear models with non-local boundary conditions and complex domains, for which an analytical analysis is difficult or even impossible. The weakness of such approach is an error of difference approximation of the problem, which can lead to significant loss of accuracy in the simulation of evolutionary problems. The second approach uses the method of characteristics of the theory of hyperbolic PDE which allows us to reduce an initial-boundary value problem for nonlinear transport equation with time delay to the initial Cauchy problem for the nonlinear delay ODE and renewal integral equation [3,6,10,17,32,33,34,41]. Using the method of steps from the theory of delay ODE [24,39,41,48] and methods of solution of some particular classes of nonlinear ODE we can write the explicit recurrent algorithm for solving the initial problem for the nonlinear delay ODE and computing the travelling wave solution for the initial-boundary value problem for semi-linear transport equation. The strengths of this approach are much more high accuracy of simulation in comparison with the numerical approach, because the numerical solution in this case does not contain the error of difference approximation. The numerical error in this approach appears when we apply at the each time step the quadrature formulas for computing integrals of solution obtained at the previous steps of time. The weakness of such approach is higher complexity of computational algorithm in comparison with first approach and possibility of constructing of explicit algorithm only for some particular classes of nonlinear age-or size-structured models. The comparative analysis of these two approaches for nonlinear monocyclic age-structured model is performed in our papers [3,4].
In this paper we derive the original numerical algorithm for computing the travelling wave solutions of age-structured model of polycyclic population dynamics which is presented in Section 3 and used in Sections 4-6 for simulation of different scenarios and dynamical regimes. The age-structured models with the density-dependent death rate are based on the semi-linear transport equations which belong to the class of hyperbolic partial differential equations of the first order known also as wave equations. The solutions of such equations have a form of travelling wave which propagates along the characteristics. A feature of the travelling wave solution is that it depends from characteristic variable and time. It does not depend explicitly from the age. The numerical algorithm uses the explicit form of recurrent formulae in contrast to numerical algorithm obtained in [6] for the linear age-structured model of polycyclic population dynamics with travelling wave solution in the form of infinite series. Explicit form of new algorithm allows us to compute the travelling wave solution for finite number of steps with high accuracy.
Travelling wave solution is represented at any given time in terms of time-independent parameterized class of algebraic functions. Sections 4 and 5 deal with numerical simulations in the case of small and large values of discrete time delays respectively. In the case of small value of time delay, the asymptotic behaviour of the solutions is qualitatively similar to the case of no delay with the solutions approaching an asymptotically stable steady state in increasing, decreasing or oscillatory manner depending on the initial state. Increasing the time delay leads to the outbreak solutions of two types: pulse sequence and single pulse, which describe the populations with cyclical eruption dynamics and population with pulse eruptive dynamics respectively [15]. In general, the coefficients of the model may depend on time, so that we have an non-autonomous system which can accommodate environmental changes in time. Simulations illustrating the properties of the solutions in the case of periodic coefficients are presented in Section 6. Some concluding remarks are given in Section 7.

Two phase nonlinear age-structured model with delay
We consider a nonlinear model of one-sex population with proliferating and quiescent phases. The model of vital dynamics is developed on the basis of Quiescence-Growth Model of cell population dynamics which is presented in the experimental work [50] and is schematically shown in Figure 1. Each individual in the first, proliferation phase age and can proliferate with some probability within a fixed age period beginning from minimal age x p up to maximal age x d . It gives birth to γ new individuals (in cell populations γ = 2) with zero age which can remain in the proliferating phase with some time-dependent probability f (t) or join the second, quiescent phase with probability 1 − f (t). The individuals in the first phase who did not divide nor die and reach the maximum age x d move to the quiescent phase. In the second, quiescent phase, the individuals age, die with some probability, that is they are the same as the individuals in the first phase except that they do not proliferate. It consists of two age groups, younger than x d and older than x d . The first group is replenished only by newly born individuals from proliferating phase. The second group is replenished by individuals from the first group in this phase and from the individuals in the proliferating phase who reach the age x d . The maximum age of this second group in the quiescent phase these individuals is theoretically unlimited but is clear that there will be practically insignificant number beyond certain age.
Let p(x, t), q (1) (x, t) denote the population densities of the proliferating phase and of the first part of the quiescent phase at age x and time t in the domain¯ (1) (2) (x, t) denote the population density of the second part of the quiescent phase for x and t in the domain (2) The respective numbers for the first and second parts of quiescent phase are Q (1) (2) (x, t) dx, so that the total population at time t is W(t) = P(t) + Q (1) In constructing the mathematical model we use the linear age-structured model of cell population presented in [25] with additional assumption that the death rate depends on the total population, thus providing for population size feedback reflecting environmental constraints. The considered two-phase population is represented via a system of initial-boundary value problems for nonlinear delay transport equations. For the density of proliferating phase: whereθ is a birth modulus, is a natural death rate of individuals; l = dx/dt = const > 0 is a speed of individual's 'maturing' (biological clock); 0 < r(x, t) is a coefficients of proportionality of nonlinear death rates; 0 ≤ p 0 (x, t) is an initial value of individuals density; γ > 0 is an individual's reproduction rate; 0 < f (t) ≤ 1 is a probability of those that new born individuals rest in the proliferating phase; 0 < τ is a time delay or time lag -some given positive parameter. Nonlinear death rate in Equation (1) characterizes the features of the interaction of growing population with nutritional environment, the sustainability and survival of individuals of overcrowded population in the bounded area, etc. In particular, the growth of resource consumers in the habitat of population leads to the starvation and increasing of mortality in population. Further, it is assumed that this feedback is delayed in time. Time lag τ in this case may characterize the features of the interaction of growing population with nutritional environment, the sustainability and survival of individuals of overcrowded population in the habitat with limited food resource. In all cases time delay is a time interval of feedback reaction of all subpopulations growth on the mortality process of individuals in quiescent subpopulation. In this model we do not specify this parameter because it can take different values depending from the particularities of biological population and complexity of their habitat (food chains). In the following sections we carry out the numerical study of dynamical regimes of population with the different values of τ . For the individuals of first and second parts of quiescent phase we have: where 0 ≤ q (1) 0 (x, t), 0 ≤ q (2) 0 (x, t) are the initial values of population density in first and second parts of quiescence phase. Considered system (1)-(9) may be used efficiently for study and modelling of applied problems in biological systems. But, since our objective in this paper is only to study the principal patterns and asymptotically stable states of considered system, we can make some insignificant simplification. Let us assume, that functions α 0 (x, t) and r(x, t) in Equations (4) and (7) are independent from age and may be substituted by new functionsᾱ(t) andr(t) (e.g. average by age values of functions).
Then, instead of systems (4)-(6) and (7)- (9) we can consider the dynamical system for the number of all individuals in the quiescent phase Q(t) = Q (1) (t) + Q (2) (t) which is can be obtained from Equations (4)-(9) in the form Cauchy problem for the nonlinear delay ODE:Q Thus the subject of our analysis is a system (1)-(3), (10), (11). We will analyse both the behaviour of solutions of this system and dynamics of number of all individuals:

The travelling wave solution and recurrent steps method
For the problem (1)-(3), (10), (11) we develop a hybrid numerical algorithm combining the well-known approaches of biomathematics and bioinformatics -the method of steps [24,39,41,48] from the theory of time-delay nonlinear ODE, method of characteristics [3,6,10,17,32,33,34,41] ('travelling wave' solution) from the theory of hyperbolic equations and numerical integration methods. Our objective here is to create a numerical toolkit which allows us to study the main features of processes of two-phase population dynamics with high accuracy. Without loss of generality we will assume that time cut [0, T] may be divided on the set of which create the consequence of sets: Let us define also the supplementary set of cuts for each natural number from the range k = 1, . . . , K:ˆ where [x] is an integer part of number x. Using the method of characteristics for linear hyperbolic equations we reduce the problem (1)-(3) to the Cauchy problem for nonlinear delay ODE. Integrating this equation and exploiting the method of steps we can write the travelling wave solution of initial-boundary value problem (1)-(3) and explicit solution of initial problem (10), (11): where α( is given: The values of functions F (k) (v) for the other numbers from range k = 2, . . . , K are given: where the set of supplementary functions F (k) m (v) is defined for k = 2, . . . , K as: We have to add to these expressions the compatibility condition for the first and second parts of travelling wave solution (19). This condition has a form of restriction for the initial distribution of population in an initial point of time: The solution of problem (10), (11) is given by: Using the steps method for delay nonlinear ODE we can write the recurrent algorithm for calculating the travelling wave solution (19)- (30) and solution of Cauchy problem (31), that is, we compute the numerical solution of system (1)-(3), (10), (11). The general algorithm is defined as follows. First, we calculate functions p(x, t), P(t) on the first time cut t ∈ [0, τ ]using the expressions (19)-(30), (12) and the known values of P(t − τ ), Q(t − τ ) from the period t ∈ [−τ , 0]. Then, we calculate Q(t) (32), S(t) (12) using function p(x, t) and known values of P(t − τ ), Q(t − τ ) from the previous time cuts. Afterwards we consider the next time cut t ∈ [τ , 2τ ]and repeat this calculating procedure. Finally, we stop this process and output result when time reaches the maximum value T. The recurrent algorithm is used in the following numerical simulation of two-phase population dynamics for different scenarios.

Simulation of population dynamics which evolves to the asymptotically stable states for small values of time delays
Asymptotically stable states of age-structured population dynamics model means that the values of population density in all age points evolves to some fixed values and does not change for a long time. In the major part of biological applications the asymptotically stable behaviour of reproductive population is a crucial question. In fact, steady states mean achieving of a balance between proliferation and mortality processes for a sufficiently long time period. The linear age-structured models considered in papers [5,6], do not involve the dynamic mechanism of population self-regulation. In work [14] authors discuss the different lows of proliferating cell population growth in tumour and normal (healthy) tissue. While the size of cell population including proliferating and quiescent cells in a normal tissue remains bounded so as to maintain tissue homeostasis the number of proliferating cells in a tumour tissue increases continuously. The models of cell population of healthy tissue have to include some mechanisms of self-regulation which provide the existence of steady states of dynamical systems within some fixed time interval. Nonlinear model of age-structured population dynamics which was proposed by Gurtin and MacCamy [27] considers more complex mortality process in population. According to their approach the growth of population density leads to the individual's mortality growth by means of a selfregulating feedback. This dynamic mechanism of self-control was introduced in the death rate of the evolutionary equation as an increasing function of the number of all individuals at the fixed time. In this paper we develop a nonlinear age-structured model of population dynamics and consider the feedback of population growth in death rate as a function of number of all individuals given in some previous moment of time. This means that population dynamics is described by delay nonlinear transport equation. The numerical study of asymptotically stable states for autonomous system is presented for an academical test. In all formulae with integrals we used trapezoid quadrature formula with the order of approximation O(h 2 ) (h is a step of difference mesh). We consider a sufficiently small value for the time delay τ = 0, 02x d for autonomous system (1)-(3), (10), (11). Such approach allows us to eliminate the impact of dynamical properties of all coefficients and delay process features in the dynamical properties of evolution system. In system (1)-(3), (10), (11) the set of stationary coefficients is given by: where θ(x) is a probability density function of the continuous linear increasing distribution by age for the adult individuals, α 0 (x, t) is a monotone increasing by age function. All parameters are correct and are chosen according to the conditions of the theorems of existence and uniqueness of travelling wave solution of age-structured models. We do not specify here any biological population and define all parameters (33)- (36) only for the illustration of existence of specific dynamical regimes of age-structured population model. The constant ρ 0 is defined from the restriction (31) which led to the following nonlinear algebraic equation The set of constants for all functions in (33)-(37) is given in Table 1.
The numerical experiments are carried out on the uniform mesh with discrete points The universal step by age and time is less or equal to the time delay h = 0, 5. Analysis of asymptotically stable states of autonomous system is performed only for reproductive population when the reproduction process dominates the mortality process. In this case the total value of individuals in population tends at a long time to some fixed value S(t) → S * max withṠ(t) > 0. The dynamics of non-reproductive (vegetative) population is expected and predictable: the number of all individuals in population decreases exponentially (Ṡ(t) < 0) and evolves to some fixed value S(t) → S * min for a long time. In Figures 2-9 we exhibit the results of the numerical simulations for three different regimes of population dynamics. In each experiment the solution evolves to asymptotically stable states at the final moments of time. In the first case (Figures 2 and 3) the number of proliferating and quiescent individuals (like cells of normal tissue) P(t), Q(t) and number of all cell S(t) increases on the first stage and evolves then to some fixed value. The population density p(x, t) in each age point is also increasing function which evolves to some stationary state for a long time. In the second regime (Figures 4 and 5) solutions P(t), Q(t), S(t) evolve to some fixed point after negligibly small oscillations in the neighbourhood of initial value. The population density p(x, t) starts with little wiggles (oscillations with negligibly small amplitude) in the neighbourhood of the equilibrium initial state and after that evolves to the stationary distribution which is very close to its starting, initial value p 0 (x, 0). The third regime of population dynamics (Figures 6 and 7) is characterized by decreasing of the number of proliferating and quiescent individuals P(t), Q(t) and number of all cell S(t) with following attracting to some fixed values for a long time. The population density p(x, t) decreases in the first time period and then evolves Table 1. Constants of system (1)-(3), (10)-(11) and functions (33)- (37).  to some fixed distribution for a long time. The third regime can be considered as opposed to the first one.
In all numerical experiments each type of asymptotically stable state of reproductive population is obtained for different values of the parameter of proportionality r 0 (35) of nonlinear death rate. The other functions (33)-(36) are fixed. In fact, r 0 plays the role of a self-control population regulator, which changes the balance between two opposite processes -proliferation and mortality -thorough the nonlinear death rate. The results of numerical study exhibit that for the set of stationary functions (33)-(36) there is a so-called 'equilibrium' value r * 0 for which the density of population and number of all individuals evolve to the second type of asymptotically stable state. The first and third types of asymptotically stable states of autonomous system correspond to the values r 0 < r * 0 and r 0 > r * 0 , respectively.  The graphs of density of proliferating phase p(x, t), number of individuals of proliferating phase P(t), quiescent phase Q(t) and number of all individuals in the population S(t) obtained in the numerical simulations provide insight into the temporal and aging evolution and asymptotical behaviour for a long time of two-phase population for the autonomous system.
The results of numerical experiments presented in Figures 2-7 motivate us to study one important property of asymptotically stable states of autonomous system. It is well known that the different solutions which correspond to the different initial values of autonomous system with fixed coefficients of equations evolve to the same stationary values for a long time. We consider two opposite regimes of population dynamics with increasing population density (first type, Figures 2 and 3) and decreasing one (third type) which corresponds  to the bigger initial value (key parameter ρ 1 = 40, 0) and the same parameters of dynamical system. The result of simulation and comparison of two opposite population dynamics are shown in Figures 8 and 9. We can conclude that for the fixed set of coefficients of Equations (33)-(36) of autonomous system there exists an asymptotically stable state which can be considered as attracting final distributionp(x) for different positive initial values of population density p(x, t), for a long time.

Numerical simulations of quasi-periodic oscillations and pulse sequences of population dynamics
Delay differential equations play a very important role in theoretical population dynamics, biology, demography, ecology, etc. More realistic mathematical models of complex dynamical systems and processes should include some of their previous states. We consider the  two-phase age-structured model of population dynamics which involves a time-delayed feedback influence, of the density growth, on their death rate. Note that nonlinear delay differential equations in the most cases have quasi-periodic solutions which can evolve to some quasi-periodic or even to chaotic random functions (for example, as a result of bifurcation cascade) [38,39,41]. We use in the simulation the functions (33)-(37) and constants given in Tables 1 and 2. We consider three different regimes of population dynamics with delay feedback for middle time lag τ = 5, 0x d , large time lag τ = 10, 0x d and very large time lag τ = 20, 0x d .   The results of simulations are shown in Figures 10-15. For the first value of time delay we obtained harmonic oscillations of density of the proliferating phase, numbers of individuals in each phase and number of all individuals in population. The dynamics of individual's density in the proliferating phase is described by the increasing function in the initial time period due to the dominated reproduction process. Then this function reaches a finite maximum value which is a result of the balance between reproduction and mortality processes. And finally the population density of the proliferating phase takes the form of damping harmonically oscillating function that converges to an asymptotically stable state. In fact, we have obtained a regime of population dynamics which is very close to the first dynamic regime presented in Section 4.    For instance, the time series of NIH 3T3 cell population in the form of pulse sequence were obtained by using the FUCCI method in work [16] (Figures 1 and 2). Such extraordinary changes of population density, so called 'outbreak dynamics', are often observed in populations of insects such as desert locusts [19], herbivorous insects [1], forest Lepidoptera [45], Cardiaspina albitextura Taylor [15], etc. According to the classification of theory of outbreaks pulse sequence is inherent in the populations with cyclical eruption dynamics [15].
The third regime of population dynamics deals with the very large values of time lag τ = 20x d . We obtain the density of proliferating phase, number of individuals in each phases and number of all individuals in the form of single pulse, shown in Figures 14 and 15. Large time delay of the nonlinear death rate at the initial time period leads to the domination of reproduction process over mortality and as a result we observe the tremendous growth of population density ( ∼ 10,000%) with growth of number of all individuals ( ∼ 200,000%) up to some maximum values. Solution reaches the maximum value at the time moment close to the time delay. As the nonlinear death rate of proliferating phase is directly proportional to its density the tremendous growths of death rate lead to the dominance of mortality in the next time period with tremendous decreasing of population density and number of all individuals in population. This fall of population density is so great that, in fact, the population cannot reproduce itself. The maximum value of p(x, t) in the subsequent time period is less than the minimum floating-point number in a PC's memory ≈ 10 −30 . As a result, the population stops the reproduction and disappears. Solutions of system (1)-(3), (10), (11) take the form of single pulses (Figures 14 and 15).
Single pulse solution describes the process of population elimination within a fixed interval of time due to the nonlinear death rate with time delay. Such behaviour of cell population was observed in laboratory experiments [50]. According to the theory of outbreaks [15] these populations are classified as population with pulse eruptive dynamics. The dynamics of such populations is divided into four stages: non-outbreak phase (initial phase), building phase (increasing part of pulse), outbreak phase (top of the pulse) and decline phase (decreasing part of pulse), Figures 14 and 15.

Numerical simulations of population dynamics for oscillating death rate and birth modulus
For the non-isolated biological systems the influences of environment often lead to the significant changes in the biological parameters of individuals. For example, in the works [6,7] we applied age-structured models of cell population for the plant hop whose mean lifespan is three years. Coefficients of equations were considered in the set of quasi-periodic algebraic functions with the time period of one year. Such approach allowed us to fit theoretical curves to the experimental data with acceptable accuracy. The dynamical properties of biological parameters in mathematical models play a very important role in the simulation and solving of different applied problems like parameter estimation and optimal control of biological systems which include the impact of environmental factors [7,8]. According to the theory of outbreaks [15] the general causes of the explosive increase in the abundance of a particular species are numerous external factors such as dramatic changes in the physical environment, environmental stresses, influence of natural enemies and predators. The study of impact of environmental factors in outbreak population dynamics is presented in work [1,15,19,45,46]. Theoretical study of quasi-periodical regimes of population dynamics is carried out for the following set of oscillating functions: where the constant ρ 0 is a solution of nonlinear algebraic equation (37) which is related the compatibility condition (31). The set of constants for functions (38)-(42) is given in Table 3. We consider three different regimes of population dynamics for functions (38)-(42) with three types of quasi-periodic birth modulus (38): with (I) small, (II) middle and (III) large magnitudes (see parameters θ 0 and θ 1 in Table 4). The birth modulus θ(x, t) and natural death rates α 0 (x, t),ᾱ(t) has two oscillation periods defined by the parameter σ (42). In order to eliminate the impact of delay factor in this simulation we use small value of time lag τ = 0, 02x d . In all dynamic regimes we fixed the values of death rates α 0 (x, t), α(t), parameter of proportionality r 0 and change the value of supplementary empirical parameter R rp (t) which is considered as an indicator of reproduction process: In the first experiment (first dynamic regime) the proliferating phase of population contains the reproductive individuals. The value of reproductive parameter is always more than critical threshold R rp (t) > 1, 25 (plot 1 in Figure 16). In the second experiment (second regime) the proliferating phase of population contains the reproductive individuals for almost all time intervals with large value of reproduction parameter R rp (t) ≥ 1, 25 (plot 2 in Figure 16). In the third experiment (third regime) the proliferating phase of population contains the reproductive individuals at the less than half of time intervals and individuals of this phase are considered as semi-reproductive (plot 3 in Figure 16).
The result of simulation for the first regime of population dynamics is shown in Figures 17 and 18. For the always reproductive population with the slightly oscillated death rate and birth modulus we obtain the quasi-periodic oscillations of density of proliferating phase, numbers of individuals in quiescent phase and number of all individuals in population. It looks like a resultant function of dynamic balance between reproduction and mortality processes with harmonic oscillations.
The result of simulation for the second regime of population dynamics is shown in Figures 19 and 20. For the almost always reproductive system with middle values of oscillating Table 3. The set of constants for the experiments with oscillating biological parameters.   Figure 16.     death rate and birth modulus we obtain the pulse sequence of population density of proliferating phase, numbers of individuals in quiescent phase and number of all individuals in population.
Although the behaviour of obtained solution in this case is very similar to this one for the population dynamics regime with time delay shown in Figures 12 and 13, there are some other explanation of such behaviour. The magnitude of oscillations for all obtained in this experiment functions is significantly lower than magnitude of solution for second regime of population dynamics in the previous section (Figures 12 and 13). We consider almost always reproductive population with increasing and decreasing birth modulus. On the descending section of birth modulus curve the mortality process dominates on the reproduction process and population density decreases to the minimum but not critical value. Latest allows population to renew the reproduction process with the small value of natural death rates α 0 (x, t),ᾱ(t). Due to the oscillations of death rate and birth modulus this process is repeated two times and results in the pulse sequence of population density p(x, t), number of individuals in each phases S(t), P(t) and number of all individuals Q(t).  The result of simulation of the third dynamic regime -semi-reproductive population is shown in Figures 21 and 22. We obtain the population density of proliferating phase, number of individuals in each phases and number of all individuals in the form of single pulse. The forms of these plots are very close to these ones shown in Figures 14 and 15 for the third regime of population dynamics with time delay in previous section. It is worthwhile to note that the amplitudes of plots in this case are not so huge in comparison with the amplitudes of plots in Figures 14 and 15. The fall of population density p(x, t), number of individuals in each phases S(t), P(t) and number of all individuals Q(t) is a result of decreasing of oscillating birth modulus less than some critical value. This decreasing takes place for sufficiently long time period (plot III in Figure 16). Total domination of mortality over the reproduction process is so long, that the maximum value of p(x, t) in the next time period is extremely small, less than the minimum floating-point number in a PC's memory ≈ 10 −30 . As a result, the population stops the reproduction and disappears, solutions of system (1)-(3), (10), (11) takes the form of single pulses (Figures 21 and 22).
The results of numerical study of quasi-periodical regimes of population dynamics carried out in this section allow us to conclude that all types of dynamical regimes obtained in the previous section for age-structured model with time delay can be reproduced by the specific oscillating death rate and birth modulus without time delay. The actual results of our study of nonlinear age-structured model with time delay are close to the conclusions of the outbreak theory [15] formulated for the single-species population model for nonlinear ODE. Two main factors in the model of nonlinear age-structured population dynamics lead to the outbreak solutions: density-dependent death rate with time delay and/or quasiperiodical oscillations of birth modulus and natural death rate due to the influences of external environmental factors.

Concluding remarks
In this paper we have developed a two-phase polycyclic age-structured model of population dynamics for proliferating and quiescent phases used in [25] for simulation of biological cell population dynamics. This model is based on the system of initial-boundary value problem for semi-linear delay transport equation which describes the evolution of density of proliferating phase and initial problem for nonlinear delay ODE which describes the dynamics of number of individuals in quiescent phase. The method of characteristics and method of steps enabled us to obtain the explicit recurrent algorithm for computing the solution of this system and construct the effective numerical method for study the features of two-phase nonlinear age-structured model with time delay.
First crucial question is the existence of asymptotically stable states of autonomous system. Existence of asymptotically stable states of nonlinear age-structured model was proved for particular systems in [18,27] and was studied also for nonlinear age-structured monocyclic population models in [3,4]. Numerical study of asymptotically stable states of two-phase population for autonomous system was carried out with small value of time delay and stationary set of coefficients of equations. We observed that in all numerical experiments dynamic system always evolves to some asymptotically stable state independently from the their initial state, for a long time. There exist only three different types of solutions which describe three different dynamic regimes -increasing, quasi-equilibrium and decreasing behaviour of density of proliferating phase and total numbers of individuals in two phases of population. The main results of this study are very close to that presented in the papers [3,4] and obtained for the similar nonlinear age-structured monocyclic model of population dynamics. We can conclude that the feedback of nonlinear death rate in agestructured model plays part of self-regulating factor in population dynamics and always provides the asymptotically stable states of autonomous systems.
Next question is the existence and possible forms and types of quasi-periodic solutions of the autonomous system with different values of time delay. The quasi-periodic solutions of the systems of nonlinear ODE with discrete time delay were obtained a long time ago and studied in numerous works (see [21,38,39] and references therein). Our objective here was to obtain and study all possible dynamic regimes of age-structured polycyclic population dynamics with nonlinear death rate and discrete time delay. Performed simulations indicated existence of the only three regimes of population dynamics which are described by the following solutions: (1) harmonic oscillations which evolve to some asymptotically stable states for a middle value of time delay; (2) pulse sequence for a large value of time delay; (3) single pulse for a very large value of time delay. Last two outbreak dynamical regimes describe the populations with cyclical eruption dynamics and population with pulse eruptive dynamics respectively [15]. The tremendous growth of population density and number of individuals in outbreak dynamical regimes is a result of dominance of reproduction process in population over the mortality due to the lack of feedback of the nonlinear death rate with large value of time delay. Following decreasing of population density within the subsequent time period exhibits the huge impact of feedback of nonlinear death rate on the mortality of individuals. This impact is so huge that the maximum values of population density and number of individuals become less than minimum floating-point number in a PC's memory ≈ 10 −30 . Population is not able to reproduce itself and vanishes. In our opinion, this solution reflects the dynamical properties of nonlinear age-structured model with time delay. Interpretation of this solution may differ for different real-world applications and may be of interest for scientists who work in the field of life science.
Final question of our paper is study of different regimes of population dynamics for oscillating with different amplitudes death rate and birth modulus for the fixed negligibly small value of time delay. Our goal here is to reproduce the same dynamical regimes which were observed in the previous simulations for some other scenarios. The oscillations of death rate and birth modulus may be caused by impact of external environmental factors such as season changes of climate, periodic treatment of population by electromagnetic fields or radiation, periodic chemical treatment etc. We consider three different dynamical regimes for always reproductive, almost always reproductive and semi-reproductive populations. In these simulations we observed three different types of solution: (1) harmonic oscillations for always reproductive population; (2) pulse sequence for almost always reproductive population; (3) single pulse for semi-reproductive population. So, we have reproduced all obtained in the previous simulations regimes of population dynamics including two types of population outbreaks. In contrast with previous simulations we used the oscillating coefficients of equations and negligibly small value of time lag. It means that both scenarios of quasi-periodic dynamic regimes -harmonic oscillations, cyclical eruption dynamics and pulse eruptive dynamics may be exploited for the identification of biological populations. The applied mathematical problems such as parameter identification and optimal control of the biological, ecological systems with oscillating characteristics are often considered in the sets of parameterized quasi-periodic algebraic solutions. In the context of these problems it is worthwhile to know if we can identify the biological system with quasi-periodic oscillated characteristics using the models with nonlinear death rate with time delay or model with oscillating coefficients of equations.
Obtained in this paper result of theoretical analysis and numerical simulations for the two-phase nonlinear age-structured polycyclic population dynamics model with density depended death rate and discrete time delay can be exploited in a variety of complex biological and ecological systems for solving of different problems in the real-world applications.

Disclosure statement
No potential conflict of interest was reported by the authors.