A SHORT SURVEY ON DELAY DIFFERENTIAL SYSTEMS WITH PERIODIC COEFFICIENTS

In this survey, we collect a few results scatted in the literature covering advances in systems of periodic delay differential equations including both modeling and theory. We will present the (almost) most oldest and (almost) most recent contributions made for this subject.


Introduction
The temporal variation of environment plays an important role in the nonlinear dynamics of life sciences and engineering systems. For example, the effects of periodically varying environment on the evolution of epidemiological systems are significant, as the selective forces on systems in the fluctuating environment differ from those in a stable and temporally constant environment. The periodicity of the environment (such as seasonal weather, food supplies, mating habits and reproduction activities) is often reflected by the periodic coefficients of the modeling dynamical systems. As time lags are also involved in the dynamic feedback of these systems, systems of delay differential equations with periodic coefficients (or periodic systems of delay differential equations) arise most naturally. Delay differential equations have been studied for at least 100 years (see, for example, the work by Schmitt in 1911 [121]) and the study of periodic systems of differential equations must have an even longer history. It is interesting to note that both periods and delays provide the most convincing mechanisms for observed nonlinear oscillations in science and engineering. However, which mechanisms are dominating for specific oscillation and how delay and periodicity interact to generate complicated behaviors remain challenging questions. Answers to these questions, even partially, are expected to find applications in a wide range of important applications in ecosystem management, avian influenza and vector-borne disease prevention and control in template environments. One can describe the above problems in the following mathematical settings.
Let r be a nonnegative constant or infinity. We denote by C the space of continuous functions mapping [−r, 0] into R n with the uniform norm defined by A general form of periodic systems of delay differential equations can be written as where F is a smooth function in (t, φ, α) ∈ R × C × R p for some p ∈ N and satisfies F (t + T, φ, α) = F (t, φ, α) for some T > 0. The solution segment y t ∈ C is defined by y t (θ) = y(t + θ), −r ≤ θ ≤ 0 and t ≥ 0.
The solution operator U : R → C defined, for t ∈ R and φ ∈ C, by U (t)φ := y t (φ) normally does not have a closed form, even for specifically identified periodic solutions (such as globally attractive nontrivial solutions). Thus, some of the theories (such as Floquet multipliers and bifurcations) developed for autonomous delay differential equations which have been extended to periodic systems (theoretically) can be hardly applied (directly) to provide practically tractable qualitative solutions even for the simplest scalar delay differential equations with periodic coefficients. The aim of this survey, therefore, is to revisit some of the old and recent results on the qualitative theory of periodic systems of delay differential equations, to select a few prototype models, to collect relevant results and to identify the theory from the dynamical systems point of view.
In what follows, we try to arrange our presentation by topics and in chronological order. In Section 2, we present some periodic models with delay that arise from life sciences, economics and engineering. Section 3 is devoted to the theoretical results for periodic differential equations with finite, distributed, periodic and statedependent delays. Finally, a brief conclusion will be given in Section 4. However, we must declare clearly that we do not aim for a complete account of the theory and applications of periodic delay differential systems which will require a much more substantial effort.

Modeling
In this section we describe a few differential equations with delay and periodic coefficients that arise from epidemiology, ecology, season migration, biology, engineering and economics.

Epidemiology
Many epidemiological problems can be modeled as periodic systems of nonlinear differential equations with delay. An example is a seasonal model with immunity (for diseases which are not lethal and do not confer immunity) developed by Cooke and Kaplan in 1976 [29]. The resulting model is a periodic scalar delay equation of the form I (t) = a(t)I(t)(1 − I(t)) − a(t − τ )I(t − τ )(1 − I(t − τ )), (2.1) where τ is the length of time an individual is infective, and a(t) is the contact rate defined as the average number of effective contacts with other individuals per infective per time period. This equation has been the subject of intense study, for example, by Nussbaum in 1977 [104], Smith in 1977 [126] and Nussbaum in 1978 [105].
Modeling the recurrent outbreaks of measles, chickenpox and mumps in cities, London and Yorke in 1973 [91] obtained a system of periodic differential equations with two delays, where the first delay represents the time from exposure to infectivity and the second delay represents the duration of infectivity. The period here represents the seasonal variation which is caused by the close contacts made by children, particularly during the coldest months, when school is in session. To describe their model, we denote by E(t) the number of individuals exposed to the virus, S(t) the number of susceptibles and I(t) the number of infectives at time t. Let β(t) be the contact rate which reflects the social behavior of the individuals of the considered community and the ease with which the disease is transmitted. The contact rate is supposed to be periodic since it may vary greatly during a year due to such seasonal factors as the weather and the fact that children are in school only during certain months. Since London and Yorke were interested in diseases that can be acquired only once, the rate of change of susceptibles is given by dS/dt = γ − β(t)S(t)I(t). It is also assumed that all individuals exposed at time t will become infectious after time T 1 and remain infectious only for a time interval of length T 2 . Thus, the rate of change of infectives is given by dI/dt = E(t − T 1 ) − E(t − T 1 − T 2 ). Therefore, the periodic model with delay is given by E (t) = β(t)S(t)I(t), Another similar example of models was presented in 1985 by Dietz and Schenzle [36]. Their model described the fluctuations in childhood diseases induced by the cycle of the school year with interruption of the infections chain during vacations and the inclusion of new susceptible individuals at the beginning of each school year, which gave rise to a periodic system with delay. Based on the assumptions of the model of London and Yorke [91], an epidemiological model with vaccination was developed by Zhang et al. in 2010 [150] as follows: Let S(t), I(t), R(t) and V (t) denote the numbers of susceptible, infective, recovered, and vaccinated individuals, respectively, at time t. The time-dependent parameters are denoted by: A(t) the recruitment rate of the population (into the susceptible class), µ(t) the natural death rate, α(t) the disease-induced death rate, p(t) the fraction of vaccinated susceptible and β(t) is the transmission coefficient of the disease. Let τ be the immunity period provided by the vaccine. Then the term p(t − τ )e − t t−τ µ(x)dx S(t − τ ) represents the individuals vaccinated at time t − τ . Therefore, the model of Zhang et al. is given by 3) The period and the delay then represent the cycle of the school year and the average length of the immune period, respectively. In 2009, Lou and Zhao [94] considered a periodic SIS epidemic model with maturation delay where B(t, N ) and d(t), respectively, are the time-dependent birth and death rates of the population at the adult state, d 1 (t) is the death rate of the population at the juvenile stage, and τ > 0 is the maturation delay. Under the assumption that all parameters are T -periodic in t for some T > 0, they established sufficient conditions for the single population growth equation to admit a globally attractive positive periodic solution. Furthermore, they obtained the threshold dynamics for the system above in terms of the basic reproduction number denoted by R 0 and defined as the number of secondary infections produced by a single infective in a totally susceptible population. In 2010, Lou and Zhao [95] presented a malaria transmission model with periodic birth rate and age structure for the vector population, which is governed by a periodic and time-delayed system of seven differential equations. They first introduced the basic reproduction ratio R 0 for this model and then proved that there exists at least one positive periodic state and the disease persists when R 0 > 1. It was also shown that the disease will die out if R 0 < 1, provided that the invasion intensity is not strong. They further used these analytic results to study the malaria transmission cases in KwaZulu-Natal Province, South Africa. Some sensitivity analysis of R 0 was performed, and in particular, the potential impact of climate change on seasonal transmission and populations at risk of the disease was analyzed.
During the last two decades, the definition and computing methods of the basic reproduction number R 0 for structured populations have been widely developed mainly in the context of epidemic models, and R 0 has become a most important key concept in mathematical epidemiology. Originally it quantifies the threshold condition of population growth in a constant environment, so it cannot be applied to formulate the threshold principle for populations in time-heterogeneous environments, which are described by non-autonomous dynamical systems. Since the 1990s, several authors (Heesterbeek and Roberts [63,64], Bacaër and Guernaoui [7], Bacaër and Ouifki [8], Bacaër [4], Wang and Zhao [144], Thieme [136], Bacaër and Ait Dads [5,6], Inaba [68,69], Zhao [155]) proposed some ideas to extend the definition of R 0 to the case of periodic environment. In particular, the definition of R 0 in periodic environments by Bacaër and Guernaoui in 2006 [7] was the first most important, because as was shown by Bacaër and Ait Dads in 2011 [5] and in 2012 [6], their definition of periodic R 0 can be interpreted as the asymptotic ratio of the size of successive generations of newborns. In 2007, Bacaër and Ouifki [8] obtained a closed approximate formula for the basic reproduction number R 0 of equation (2.1) when a(t) = a 0 (1 + cos ωt) and is small. This formula showed that the periodic factor increased or decreased R 0 , depending on the choice of parameter values. Using the generation evolution operator, Inaba in 2012 [68] showed that the definition given by Diekmann et al. [34] and by Bacaër and Guernaoui [7] completely allows the generational interpretation and the spectral radius of the generation evolution operator coincides with the spectral radius of the next generation operator, so it gives the basic reproduction number and the new definition is an extension of the definition by Diekmann et al. [34] and by Bacaër and Guernaoui [7]. Recently, the classical threshold theory for R 0 and the Malthusian parameter in constant environments was extended to the case of periodic environments without loosing its essential features. In fact, Inaba [69] proposed another approach to establish the sign relation for R 0 and the Malthusian parameter λ 0 in a periodic environment by focusing on the existence of exponential solutions. In what follows, we summarize an important result by Inaba [69]. Let us introduce I(t, τ, a) as the density of infected population at time t and infection-age τ whose chronological age at infection is a. Let S(t, a) be the age density of susceptibles at time t and age a, and R(t, τ, a) be the density of recovered individuals at time t, duration (the time elapsed since recovery) τ and age at recovery a. Then the age-duration-dependent homogeneous SIR epidemic system is formulated as follows: where B is the number of susceptible newborns per unit time, µ(a) is the force of mortality, γ(t, a) is the rate of recovery at infection-age τ and the age at infection a, and κ is the force of infection. Inaba assumed that the force of infection is given by where m(t, a, σ) is θ−periodic with respect to t and can be interpreted as the probability that a susceptible individual with age a encounters with an infected individuals with age σ at time t, and the function f (τ ) is the probability of successful transmission of infective agents from infective individuals with infection-age τ . Assuming that m, γ, f and µ are bounded nonnegative measurable functions for t ∈ R + and a, σ, τ ∈ R + ,we state the following proposition Finally, it should be stressed that since delay equations in population dynamics can usually be considered as particular cases of age-structured systems of partial differential equations, the definition of the basic reproduction number R 0 for periodic versions of the latter introduced by Bacaër and Guernaoui [7] and Inaba [69] can be applied as well to the case of delay models (See also the work by Wu et al. in 2015 [146]). Furthermore, the result in Bacaër and Guernaoui [7] considered a periodic epidemic model with a Γ-distributed delay as application. However, its application is much more difficult and not straightforward for some cases. Recently, Zhao in 2015 [155] established the theory of basic reproduction ratio R 0 for periodic and time-delayed population models with compartmental structure as follows. Consider a linear periodic functional differential system: where F (t) : C → R m is a positive function given by and V is a matrix such that −V is cooperative. Let Φ(t, s), t > s, be the evolution matrice of the linear equation Let v(t), w-periodic in t, be the initial distribution of infectious individuals. Then the basic reproduction ratio associated to system (2.4) is defined by is the spectral radius of L and C w is the set of continuous w−periodic functions. Furthermore, the author proved that R 0 is a threshold value for the stability of the zero solution for periodic system (2.4). This theory is the most recent and gives the most important result for the definition of R 0 in periodic systems with delay.

Ecology
In 1972, Fretwell [49] claimed that "environmental variation is often an important aspect in the dynamical nature of populations with most populations experiencing at least seasonal fluctuations in their 'resource' or 'food' availability, and thus in the carrying capacity of their environment". In 1976, Nisbet and Gurney [103] considered a periodic delay logistic equation of the form where N (t) is the population at time t, r is a sinusoidal growth rate of the species in an unlimited environment, and K is a sinusoidal carrying capacity of the habitat.
They carried out a numerical study of the influence of the periodicity in r and K on the intrinsic oscillations of the equation such as those caused by the time delay.
In 1978, Rosen [113] noted the existence of a relation between the period of the periodic carrying capacity and the delay of the logistic equation. Basing on the quote by Fretwell [49], Gopalsamy et al. [51] proposed in 1990 the following time delayed w-periodic equation which describes the dynamics in a "food-limited" population system. Here, N and r are, respectively, the number and the growth rate of Daphina population, K is the mass of saturation and w is the time delay which incorporates the population fluctuations.
In 1992, Freedman and Wu [48] considered several single-species models with time delays where both the coefficients and the delays are periodic functions. These models are based on laboratory observation of the population growth of rotifers and are of the following general form where D(t, x) is a periodic function in t which represents the death rate and B(t, φ), φ ∈ C is a periodic function which represents the birth rate. This model represents the case where there is a delay in the per capita birth rate, whereas the death rate is instantaneous. A particular case of the general model developed by Freedman and Wu [48] was considered in 2008 by Berezansky and Idels [15] in order to model the Fox production harvesting. Their model is given by where r(t) is an intrinsic factor, F (t) is a variable harvesting rate, K(t) is a varying environmental carrying capacity and τ (t) is the time to develop from newborns to reproductively active adults. Periodic logistic equations with distributed delay were studied in 1977 by Cushing [30]. In 1982, Bardi and Schiaffino [11] considered the following integro-differential equation with periodic coefficients where a, b, c and K are positive periodic continuous functions. The functions a and b describe the growth of a single population whose size is N (t) while c and K describe the hereditary factors. In 2001, Li and Kuang [84] presented a delay Lotka-Volterra system with periodic coefficients for multiple species population growth which is a generalization of the single species population growth model developed by Freedman and Wu. In particular, the following periodic Lotka-Volterra cooperative system with distributed delays was considered See also Li and Xu [87] for a periodic prey-predator system with infinite delay. For other related works on periodic logistic equations with distributed delay, we refer to Bardi [10], Cohen and Rosenblat [26], Karakostas [77] and the references therein.
In a classic study by Nicholson in 1957 [102], the regulation of population size in laboratory colonies of the Australian sheep-blowfly Lucilia cuprina was analyzed. The dynamics of Nicholson's blowflies was described using a delay differential equation with constant coefficients by Gurney et al. [53]. However, during the experiments of Nicholson, the population size showed large-amplitude quasi-periodic oscillations. One of the cultures was maintained for a period of two years and the oscillations continued for this entire period of time. Thus, in 2002, Saker and Agarwal [117] modified the model in Gurney et al. [53] by assuming that all the parameters of this model are w-periodic. Their model is given by where N (t) is the size of the population at time t, P is the maximum per capita daily egg production, 1/a is the size at which the population reproduces at its maximum rate, δ is the pair capita daily adult death rate, and mw is the generation time. Taking spatial distribution of the blowfly population into account, the above periodic Nicholson blowflies model with delay was recently extended by Pang and Wang in 2015 [106] to the following nonlinear diffusive form, where ∆ is the Laplacian operator modeling nonlinear diffusion. The given boundary condition means that the habitat Ω is surrounded by a lethal environment. See also Ruan [116] for a brief survey on periodic single species models with delay.

Seasonal migration and avian influenza
Bird migration is the regular seasonal movement, often north and south along a flyway, between breeding and wintering grounds. Periodic systems of delay differential equations arise naturally from modeling the bird migration. In 2010, Gourley et al. [52] formulated a deterministic mathematical model capturing some features of the spatial dynamics of migratory birds as follows: A typical migration process involves different phases based on the consequences of biological activities and seasonality, as wintering, spring migration, breeding, maturation and autumn migration. Migration routes are interrupted by major nodes, the so called stopovers, which provide the resting locations between the flights for refueling. Therefore, the proposed model took the form of a periodic system of differential equations with multiple delays which account for the time needed for migratory birds to fly from one stopover to the next. In particular, if S i (t) denotes the number of birds in patch i, then one obtains the following model system where the first patch is the summer breeding area, the last patch (the n-th patch) is the winter refuge ground, and patches in between are the major stopovers along the migratory route. b(S 1 (t), t) is the birth rate function and d i,j (t) accounts for migration into or out of each patch. Model parameters except delays are all timeperiodic functions. The framework of Smith in 1987 [126] was used here to develop a threshold dynamics of such a model system. In 2010, Bourouiba et al. [19] developed a periodic model with delay to discuss the effect of repeated epizootic at specific migratory stopovers for bar-headed geese.
To develop some qualitative insights of the role of migratory birds in the global H5N1 spread, the model by Gourley et al. [52] was further expanded by Bourouiba et al. in 2011 [18] to incorporate the disease progress within each patch and the interaction of migratory birds with domestic poultry for a patchy model of avian influenza spread. The aim was to understand the role of the interplay between migratory birds and nonmigratory birds, particularly poultry, in the persistence and recurrence of H5N1 in endemic regions.
In 2012, Wang and Wu [141] developed a model of single species bird migration between the summer breeding ground and the winter refuge site. The obtained model is a simplified version of the models in Gourley et al. [52] and Bourouiba et al. [19]. The model is given by the following periodic differential equations with delay, where x s (t) and x w (t) are the numbers of birds in the summer and winter sites, respectively, the functions m ws (t) and m sw (t) describe the migration rates between the two patches. The subscription w represents winter, s represents summer, ws represents migration from the winter patch to the summer patch and sw represents migration from the summer patch to the winter patch. The time delay τ ws represents the time (days) flying from the winter site to the summer site, while τ sw represents the time (days) flying from the summer site to the winter site. The corresponding rates of in-flight mortality are denoted by µ ws and µ sw , and the death rates in the summer and winter sites, respectively, are denoted by µ s and µ w . The birth rate in the summer site is the intrinsic birth rate γ s (t) regulated by the density through a logistic term. All the coefficients in the system are periodic functions with the period T = 365 days. The authors showed that the periodic delay differential system is completely characterized by a finite-dimensional ordinary differential system. As a consequence, they derived the threshold condition, explicitly in terms of the model parameters, for the extinction and persistence of the considered bird species. This model was then further expanded, in 2014, by the same authors Wang and Wu [142] to incorporate the spread of avian influenza within each patch.

Biology
A few periodic delay models have been developed in biology. In 2002, Saker and Agarwal [117] modified the delay model of dynamic disease involving respiratory disorders, introduced by Mackey and Glass in 1977 [97], by assuming the periodicity of the volume of air passed through the lungs during a single breath times the frequency of respiration. In 2001, Andersen and Mackey [2] developed different versions of models of delayed periodic differential equations to evaluate the effects of periodic chemotherapy administration on normal and leukemic (a cancer of the white blood cells) cell populations. First, Andersen and Mackey assumed that all proliferating cells are affected by the drug which leads to the following system where N and P , respectively are the numbers of nonproliferative cells and proliferative cells, γ is the rate of cell loss from nonproliferative phase, τ is the time spent in the proliferative phase, β is the feedback function, rate of recruitment from non-proliferative phase. The cells in the proliferative phase die due to an increased apoptotic rate γ 0 + γ c (t). Moreover, the parameter γ c (t) is the loss rate due to chemotherapy and is periodic with the period π. They concluded, using this formulation, that in every case the protocols are more destructive to the normal cell population than they are to the leukemic cell population. Thus, they extended this model to show the effect of chemotherapy administration only on the specific phase S (DNA synthesis) to obtain the following: where S is the number of cells in DNA synthesis phase and M is the number of cells in mitosis (Note that the variable S in the S−equation above is written as P in Andersen and Mackey's paper but it is only a typo!). This last formulation leads the following conclusion: S-phase-specific resonance (optimal therapeutic success) is more pronounced than cycle-non-specific resonance, and M -phase-specific resonance is even sharper than the S-phase resonance. This (N ; S; M )-model was further extended to a (N := G0/G1; S; P := G2, M )-model to track the total cell number and fraction of cells in each phase for host and tumor during a course of periodic chemotherapy by Bernard et al. in 2010 [16]. Furthermore, they assumed that transition from one phase to another depends on the circadian time. The model is given by the following periodic differential equation with periodic delays, , where all the time dependent parameters are 24-hours periodic. Concerning In-host models, to the best of our knowledge, there was only one work -up to date-that consider the incorporation of time-periodic combination antiviral drug therapy (reverse transcriptase inhibitors and protease inhibitors) into the in-host model with delay for HIV. In fact, in 2004, Dixit et al. [37] developed a model of HIV infection that combines drug pharmacokinetics and intracellular delay. Here we will only present the model which consider the effect of reverse transcriptase inhibitors drug, given by where the viral replication period (τ ) is divided into pre-(τ 1 ) and post-(τ 2 ) drug action parts, T are uninfected CD4 + T cell lymphocytes which are generated at a rate λ and die at rate d. In the presence of virus V , healthy cells are converted, at a rate k, into productively infected cells, T * an interval τ after they are infected.

Engineering
Several engineering problems give rise to time delays and periodicity at the same time. Examples include transmission lines (Lopes [93]), neural networks (Hagan et al. [54]), parametric control of robotic systems (Insperger et al. [70]), machine tool vibrations in cutting processes (Bayly et al. [12], Zhao et al. [151], Insperger et al. [73], Hartung et al. [62], Bachrathy and Stepan [9]), aeroelastic systems with structural nonlinearities (Trickey et al. [137]), optimal controls (Deshmukh et al. [33]) and reaction-diffusion BAM neural networks (Zhou et al. [156]). Here we choose to briefly describe the model by Insperger et al. in 2000 [70] where a single degree of freedom mechanical model for milling process is proposed. This models a machine tool where a workpiece is rotating, the tool cuts the surface that was formed in previous cut. The chip thickness is determined by the current and the previous position of the tool and the workpiece. In the case of milling, the direction of the cutting force is changing due to the tool rotation, and the cutting is also interrupted as each tooth enters and leaves the workpiece. Consequently, the system can be modeled as a one degree of freedom oscillator that is excited by the cutting force. The resulting equation of motion is a delay differential equation with a time periodic coefficient given bÿ where m is the modal mass, ξ is the damping ratio, w n is the natural angular frequency and F x is the cutting force in the x−direction for a zero helix cutter. The total cutting force acting in the x-direction, F x , is given by where b is the nominal depth of cut, K s and f 0 are τ −periodic functions and characterize the instantaneous chip width and τ is the tooth pass period. Another milling-process model, introduced by Long et al. in 2007 [92], that give rise to a delay system with periodic coefficients is described in detail in Section 2.7.

Economics
A very few periodic delay models have been developed in economics. There is an interesting work of Boucekkine et al. in 1999 [17] where the authors showed that job creation can be modeled by a delayed differential equation with periodic coefficients, where the delay describes the optimal age of capital goods and the period describes the length of an exogenous profitability cycle. Their model is given by where h(t) represents employment associated with vintage t, k 2 (t) = 1−β β e −γT 0 −p(t) , β is a share of the surplus going to the worker, 1 − β is a share going to the firm, e −γT 0 are the produced units (of good from one unit of the non-produced good) from each unit of labor associated with an operating vintage T 0 and p(t) is the periodic exogenous price. The mathematical analysis and numerical simulations lead to some conclusions. Firstly, creation is asymptotically periodic with the same period as the profitability cycle. Secondly, replacement echoes seem to play a fundamental role in the short run dynamics of job creation and job destruction.

Models with periodically time-varying delays
In some cases there are good reasons for allowing the delay to vary with time, and the use of a periodic delay is especially relevant since it can model daily, seasonal or annual fluctuations. In comparison with the aforementioned systems where the delay is a constant, the situation becomes more complex, even for modeling, when the system has periodically time-varying delays such as where F and τ are smooth functions in (t, y, α) ∈ R × R × R p for some p ∈ N and τ satisfies τ (t + w) = τ (t) for some w > 0. In 2000, Schely and Gourley [119] considered the effect of seasonal fluctuations on populations in a predator prey model. Liu  In 2007, Long et al. [92], presented a milling-process model with a variable time delay associated with each cutting tooth as follows, where the tool has two degrees of freedom and the workpiece has two degrees of freedom. The variables q x and q y , respectively, represent the tool dynamic displacements measured along the X and the Y directions. The variables q u and q v represent the workpiece displacements measured respectively along the U and V directions. The quantities m x , m y , m u , and m v are the modal masses, the quantities c x , c y , c u and c v are the modal damping coefficients, and the quantities k x , k y , k u and k v are the modal stiffness parameters associated with the motions along the X, Y, U and V directions, respectively. Furthermore, τ (t, i, z) is a periodic variable time delay introduced in the governing equations through the cutting-force components F.
Recently in 2015, Wu et al. [146] developed a model describing the evolution of the ectothermic insects population in each biologically distinctive stage under seasonally varying periodic environment. The authors used the formulation in Webb [145], for general structured population dynamics models, and a subdivision of the life cycle of the given population into n stages to obtain the following stagestructured population system with temporally periodic delay: for t > A n (t), where x 1 represents the size of the mature subpopulation who are able to produce offspring (egg-laying females), x j , j = 2, . . . , n, is the size of subpopulation at the j th stage, δ is the fixed sex ratio for such population, A i−1 and A i are the timedependent minimum and maximum ages of those individuals who are developing within the specific i th stage, α i (t, t − A i (t)), i = 2, . . . , n, represents the densitydependent survival probability of an egg who was born at time t − A i (t) and is able to live until time t when the egg eventually belongs to the stage x i with full maturation, b(x 1 (t)) is the birth rate and µ i (x i (t)) is the mortality rate for the i th stage.
For more models that incorporate periodic delay, see also the work by Berezansky and Idels [15] described in Subsection 2.2 and the work by Bernard et al. [16] described in Subsection 2.4.

Theory
The theory of periodic systems of delay differential equations has been studied for a long time and the study of linear periodic systems with time lag dates back at least to the work of Hahn in 1961 [55], Stokes in 1962 [128], Shimanov in 1963 [120] and Zver'kin in 1967 [157]. Despite the relatively matured theoretical framework already established, the main difficulty, however, lies in realizing such a theoretical framework for specific systems to provide qualitative insights of the dynamics of these specific systems using the model parameters as explicitly as possible.
Let us look at the following seemingly simple periodic differential equation with delay with a fixed delay r > 0 and constant parameters α, ν ∈ R. Define the period map T (s): C → C as follows: T (s)(φ) = y s+2π/ν (φ, s).

(3.2)
We observe that the linear continuous operator T (s) is given by the following variation of parameters formula Thus, it is obvious that T can not be given explicitly if the delay r is not a multiple of the period 2π/ν. However, if we set, for example, r = 2π and ν = 1, then T becomes which has a close form. The case of periodic delay difference equations with integer time lags has received considerable attention. Efforts have been made to develop the expansion of solutions in terms of a series of Floquet solutions. Earlier results and discussions about the difficulties involved can be found in Hahn [55] and Zver'kin [157].

Floquet theory and stability
The early paper of Halanay in 1961 [56] started the development of the infinite dimensional version of the Floquet theory for delayed periodic systems, but it gave only a theoretical possibility for the stability analysis of the most general case. However, the Floquet theory have been completely developed in Hale and Lunel by [59] in 1993, Diekmann et al. in 1995 [35] and Avdonin and Germanovich in 1995 [3], when the delay is a multiple of the period. See also the recent development, by Ford and Lamb in 2007 [47], of some basic of Floquet theory for small solutions. When the delay is not a multiple of the period, but commensurable with it, then in some cases, the Floquet multipliers can be deduced from the explicit solution of a system of ordinary differential equations. By linearizing along a periodic solution of period 3 of a retarded autonomous equation of delay 1, Skubachevskii and Walther in 2003 [124] have obtained some information on Floquet multipliers of this problem, and they formulated, under some conditions, a criterion for the hyperbolicity of the 3-periodic solution. A generalization of this result for irrartional periods is established, for general linear periodic differential equation with delay 1, by the same authors in 2006 [125]. By constructing a special boundary value problem for ordinary differential equations, Dolgii in 2006 [40] determined conditions of asymptotic stability, of linear periodic delay systems, in terms of the spectrum of the monodromy operator (Floquet transition matrix).
The most difficult case arises if the delay is incommensurable with the period. In 1993, Hale and Lunel [59] showed by example that, in this case, a complete Floquet theory will not exist. There is only a Floquet representation on generalized eigenspaces of characteristic multipliers. For periodic delay differential equations, the difficulty is that, in general, the monodromy operator has no closed form (see example given for the linear equation (3.1)), so no closed form stability conditions can be expected directly except for a few and very special cases. To overcome the difficulty due to the lack of closed forms for periodic differential equations with delay, other approaches and alternatives should be developed and exploited. For example, Mallet-Paret and Sell in 1996 [98] considered a special case of delayed periodic equations with a cyclic feedback structure. They defined a discrete Lyapunov function and related their values to the real parts of the Floquet multipliers.
In particular, they proved that all Floquet subspaces are at most two-dimensional. A characteristic function in the form of a continued fraction expansion is given for a scalar delay equation by Just in 2000 [76]. If the given system is monotone (i.e., the comparison principle holds), we may expect to have some simple but practical criteria of stability. Below we present an illustrative example.

3) is determined by that of the linear ordinary differential equation u (t) = (a(t) + b(t))u(t), which can be obtained by ignoring the time delay in (3.3).
In the case where a(t) = −1, b(t) = αg(t), τ = 1 and g : [0, ∞) → [0, ∞) is continuous and periodic with arbitrary period T > 0, Chen and Wu in 2013 [25] used a discrete Lyapunov functional and the monodromy operator to show that there exists a positive value α * such that the solution is stable if α ∈ (0, α * ) and unstable if α > α * . When a is a negative periodic function and b is a positive periodic function with arbitrary period T > 0, such that b(t + 1) − a(t) does not change its sign, Nah and Röst in 2016 [101] showed that r = 0 is the stability threshold where r := T 0 (a(t) + b(t))dt. Namely, they showed that zero is unstable when r > 0, stable when r = 0 and asymptotically stable when r < 0.
In 1980, Simpson [123] used the two-timing method to state a Floquet theory for a class of linear integro-differential equations with periodic coefficients of the form where x is a column vector in R n , A(t) is an n × n matrix function, periodic in t with period w > 0, K(t, s) is the kernel function, an n × n matrix function defined for s ≥ 0, periodic of period w in its first argument, K(t, s) = O(e −γs ) as s → ∞ for some γ > 0 uniformly in t. It is assumed also that, for each continuous function φ : [0, ∞) → R such that sup θ≤0 | e γ + θ φ(θ) |< ∞ for some γ + < γ, there exists a unique x(t) solving equation (3.4) for t > 0 such that x 0 = φ. Simpson stated and proved the following theorem.
Theorem 3.1. Corresponding to (3.4) there exists a set , called the spectrum of (3.4), consisting of finitely many points in the complex plane, each µ ∈ satisfying | µ |> e −γ + w . Corresponding to each µ ∈ there exist (i) an r × r matrix B (consisting of complex constants) whose sole eigenvalue is λ where λ = (1/w) log µ (Reλ > −γ + ; r depends on µ) and (ii) a set of r column vector functions P µ (t) in C n (i = 1, ..., r) which are each periodic of period w. Then for each solution x(t) of (3.4) generated by initial data satisfying φ there are, for each µ ∈ , unique r × 1 constant complex vectors a µ such that x(t) differs from µ∈ (P µ1 (t), ..., P µr (t)) e Bt a µ (3.5) by a function x R (t) which is O(e γ + t ) as t → +∞. Conversely, sums of the form (3.5) are solutions of (3.4) for arbitrary choices of a µ .
An application of this theory results in a computation of relevant Floquet exponents associated with the linearization of the equations about the periodic solutions. Using a reduction method of integro-differential equations to systems of ordinary differential equations obtained by Domoshnitsky in 2003 [38], a study of Floquet theory for integro-differential equations with periodic coefficients was presented by Agarwal et al. in 2005 [1]. The authors employed this theory to establish a relationship between oscillation and asymptotic properties of solutions of integro-differential equations.
See also the work by Tang et al. in 2006 [132] for periodic n-species Lotka-Volterra competition system with periodic delays and the works by Davidson and Gourley in 2001 [31] and by Pao in 2005 [107] for periodic semilinear parabolic equations with delays and monotonic cases.
A large flow of numerical methods have been proposed to approximate the infinite dimensional monodromy operator with finite dimensional one with application to periodic delay differential equations. Examples of numerical methods include, but are not limited to, the methods based on averaged coefficients method by Minis and Yanushevsky in 1993 [100], semidiscretization method by Stepan in 2002 and2004 [71, 72], full discretization method by Yilmaz et al. in 2002 [149], temporal finite element method by Bayly et al. in 2003 [13] and by Mann et al. in 2004 [99], Chebyshev polynomials by Butcher et al. 2004 [24], by Bueler in 2007 [21] and by Butcher and Bobrenkov in 2011 [23]. A more general result, due to Szalai and its collaborators in 2006 [130], presented a construction of a characteristic matrix for general periodic delay differential equations, but was only useful in the numerical continuation of bifurcations. In 2011, Sieber and Szalai [122] showed that matrices constructed by Szalai et al. [130] can have a discrete set of poles in the complex plane, which may possibly obstruct their use when determining the stability of the linear system. Then they modified and generalized the original construction such that the poles get pushed into a small neighborhood of the origin of the complex plane. In 2012, Tweten et al. [138] provided a comparison of the semidiscretization, spectral element, and Legendre collocation method and showed that the spectral element method has the quickest convergence rate. Thus, Khasawneha and Mann in 2013 [78] developed a general spectral element approach for periodic systems with multiple delays. This approach was recently generalized by Lehotzky et al. in 2016 [83] for periodic systems with multiple delays and distributed delay.

Bifurcation theory
The bifurcation theory studies persistence and exchange of qualitative properties of dynamical systems under continuous perturbations. Very few "rigorous" works have been done for "constructive" bifurcation theory in periodic differential equations with delay and periodic linearized equations. The first bifurcation result, to the best of our knowledge, was established by Busenberg and Cooke in 1978 [22]. Busenberg and Cooke exploited some basic techniques on positive operators and coincidence theorem to show a transcritical (forward) bifurcation in a class of scalar periodic delay differential equation of the form (3.6) Using Floquet theory and spectral projection method, Röst in 2005 [114] obtained a formula to compute the coefficient that determines the direction of the Neimark-Sacker bifurcation for some scalar periodic delay differential equations of the form where γ is a real parameter, a(t+1) = a(t), f : R×R → R is a C 4 −smooth function satisfying f (t + 1, ξ) = f (t, ξ) and f (t, 0) = 0 for all t, ξ ∈ R. When a(t) = 0, Röst in 2006 [115] determined some conditions leading to bifurcation of invariant curves and four-periodic orbits for the equation (3.7). Bifurcations in integro-differential equations were also studied. Simpson in 1980 [123], dealt with problems of bifurcation of periodic solutions for systems of the form (3.4). In his work, Simpson used the two-timing perturbation procedure to determine conditions which lead to different kind of bifurcations and stability of their bifurcating branches, namely, Hopf bifurcation and other multi-codimension bifurcations. Another result of bifurcation in periodic differential equations with delay was given by Britton in 1990 [20] which carried out a linear stability analysis for the positive steady-states of a periodic integro-differential reaction-diffusion equation with delay. Domoshnitsky and Goltser in 2002 [39] developed a method of reduction of nonlinear integro-differential equations (IDE) to a system of ordinary or delay differential equations to study stability and bifurcation for a wide class of IDE. In particular, Domoshnitsky and Goltser considered a Hopf bifurcation result in IDE.
In 2001, Elmer and van Vleck [27] considered traveling wave solutions to delayed spatially discrete reaction-diffusion equations with periodic diffusion. The authors have obtained analytic solutions for the traveling wave problem using a piecewise linear nonlinearity and presented numerical analysis of period doubling bifurcation in the wave speed.
The Center manifold and normal forms theories are one of the most effective tools in the study of nonlinear dynamical systems when analyzing bifurcations of a given type. In 1997, Faria [43] presented the method of normal forms, to calculate the elements of bifurcations, for periodic functional differential equations of the form (3.11) where N is a purely nonlinear term. The reason for imposing such a restriction is due to the difficulty to provide a suitable change of coordinates in the neighborhood of periodic solutions for a general periodic equations of the form (1.1). Faria proved that, if certain nonresonance conditions are satisfied, the normal form on center manifolds of the trivial equilibrium for periodic equations of the form (3.11) coincides with the normal form for the averaged equation where F 0 (φ) = 1 w w 0 F (s, φ)ds for φ ∈ C. However, Hale and Weedermann in 2004 [61] established a local system of coordinates near a periodic orbit of autonomous delay differential equations of the form (3.8) where f : C → R n is C k , k ≥ 1. The key idea to set this system is summarized as follows: Suppose Γ = {p t : 0 ≤ t ≤ w} is a nontrivial periodic orbit of (3.8), and consider the linear variational equation Since 1 is always a Floquet multiplier of (3.9), consider a decomposition of the phase space C with respect to µ = 1. One possible choice for the basis of E 1 (s) is Φ 1 (s) = {ṗ s , φ 2 , ..., φ d1 } , becauseṗ t is a corresponding eigenvector for µ = 1. Let T (t, s), t ≥ s, be the solution operators for (3.9), defined by T (t, s)φ = x t (s, φ) for all t ≥ s and φ ∈ C, where x = x(s, φ) is the solution of (3.9) such that x s = φ.
Assume that Γ is non-degenerate, that is , µ = 1 is a simple Floquet multiplier of (3.9). In this case a decomposition of C is given by where Q 1 (s) is such that T (t, s)Q 1 (s) = Q 1 (t). So, the above decomposition is a natural moving coordinate system in C. Every element in C has a unique representation of the form φ = q s , w sṗ s + L s w, where L s is explicitly defined as L s w = w − q s , w sṗ s . Furthermore, they gave a description of the flow of (3.8) on the center manifold of the periodic orbit. In 2007, Qesmi and Hbid [112] considered a class of nonlinear periodic differential equations with delay of the form y (t) = f (y(t − 1), a ( sin(νt)), cos(νt)) , α) , (3.10) where , ν, α are real parameters, f : R 3 → R and a : R 3 → R are a C ∞ −smooth function satisfying f (0, a(0, 0), α) = 0 and ∂ dy f (0, a(0, 0), α) = 0 for all α ∈ R. Thus, equation (3.10) can be viewed as a perturbation, in a special way, of an autonomous differential equation. They proved the existence of integral manifolds with periodic structure without any restriction on the relationship between period and delay. Their method was based on calculation of center manifolds coefficients for autonomous differential systems with discrete delays. Consequently, they investigated the saddle-node bifurcation for this class of periodic equations with delay. In 2010, Szalai et al. [129] developed an analytic method to compute bifurcation elements (Poincar-Lyapunov constant) for some periodic delay equations where the delay is equal to the period. Their method used center manifolds theory but avoided the computation of their coefficients (Note that despite the publication of this paper in 2010, it was already accepted for publication in 2008). Deshmukh et al. in 2008 [32] used a nonlinear extension of the infinite dimensional monodromy operator defined by Halanay in 1966 [57], for linear periodic delay equations, to reduce the system of T -periodic delay differential equations of the forṁ on center manifolds. Their method does not involve the explicit computation of the adjoint equation and formalize the center manifold theory. However, it is applied, once again, only when the period is the same as the delay.

Oscillations and periodic solutions
A solution of a given dynamical system is called oscillatory if it has arbitrary large zeros. Otherwise it is called non-oscillatory; that is, if it is eventually positive or eventually negative. Oscillation theory has long been rich and attractive areas of mathematical research. This section will be organized following different class of differential equations: Periodic equations with constant delays, Periodic equations with periodic delays, periodic neutral differential equations and periodic semi-linear parabolic equations.
Oscillations in periodic differential equations with constant delay of the form where L(φ) is a linear map in φ, and N (t, φ, µ) is a perturbation term which is small for small φ and µ relative to linear part, has been studied first by Driver in 1963 [41], Bellman and Cooke in 1963 [14], Cooke in 1965 [28] and Perello in 1968 [108]. The authors used the theory of linear systems, given by Hale in 1963 [58], in order to establish procedures for studying the oscillations of perturbed autonomous linear systems by a periodic terms. By appealing to theory of global attractors and steady states for uniformly persistent dynamical systems (see Zhao [152] and Long et al. [92] for more details), Zhao [154], in 2008, established the existence of interior periodic solutions for dissipative periodic functional differential equations. This result is then applied to a multi-species competitive system and a SIS epidemic model. The following result is from Theorems 3.1 of Zhao [154].
Consider a periodic retarded functional differential equation with f : R × C → R m a continuous and bounded map, f (t, φ) a locally Lipschitz function in φ, and f (t + ω, φ) = f (t, φ) for some ω > 0. (ii) Solutions of (3.12) are uniformly persistent with respect to M 0 in the sense that there exists η > 0 such that lim inf Then system (3.12) has at least one ω-periodic solution x * (t) with x * t ∈ M 0 for all t ≥ 0.
Freedman and Wu [48], in the best of our knowledge, was the first who proved the existence of positive periodic solutions of single-species population growth models with periodic delay and periodic coefficients given by (2.5). In particular, using Horn's asymptotic fixed point theorem, they stated and proved the following theorem. Theorem 3.3. Assume that (H1) For any continuous w−periodic function x : R → R, t → B(t, x t ) is w−periodic function.
(H2) There exists a positive w−periodic continuously differentiable function K(t) such that D(t, K(t)) = B(t, K t ) for t ∈ R.
(H3) For all t ∈ [0, w], we have where | B φ (t, φ) | denotes the operator norm of the bounded linear operator B φ (t, φ) : (H4) There exists a constant δ > 0 such that for every δ 0 ∈ (0, δ), and for any Then the model equation (2.5) has a positive w−periodic solution Q(t). Moreover, if for all [0, w], then Q(t) is globally asymptotically stable with respect to positive solutions of (2.5).
Freedman and Wu pointed out that it would be of great interest and difficulty to consider the existence of positive periodic solutions of higher-dimensional systems with periodic delays, such as predator-prey or competitive systems. In 1997, combining the theory of monotone flow generated by FDEs, Horn's asymptotic fixed point theorem and linearized stability analysis, Tang and Kuang [133] answered the Freedman-Wu's question and proved the existence, uniqueness and global asymptotic stability of positive periodic solutions for periodic Lotka-Volterra system with periodic time delay of the forṁ x 1 (t), · · · x n (t), x 1 (t − τ (t)), · · · x n (t − τ (t))) , i = 1 . . . n.
In 1998, Li [86] studied the periodic differential equations with several periodically varying delays of the form y (t) = F (t, y(t − τ 1 (t)), ...., y(t − τ n (t))), . . , n, and w > 0 is a constant. By using some techniques of the Mawhin coincidence degree theory, Li proved the existence of a positive periodic solution of the initial value problem (3.13) as follows. In 2000, Schely and Gourley [119] considered scalar models incorporating a small discrete periodic delay. They shows that the steady state (if it can ever bifurcate) will bifurcate to periodic solutions. Furthermore, two-timing method was used to investigate how the stability of the system differs from that with a constant delay. Using the coincidence degree method, existence of positive periodic solutions for periodic differential equations with periodic state-dependent delays was considered by Li and Kuang in 2001 [85].
By applying a perturbation method to a class of differential equations with periodic delay close to zero, Qesmi and Hbid in 2006 [111] transformed the studied equation into a perturbed periodic ordinary differential equation. This technique can be a powerful tool to show existence of positive periodic solutions for general periodic systems with small periodic delays (See the work by Qesmi [110]).
Very little is known about positive periodic solutions for general non-monotone periodic systems with periodic delays. To the best of our knowledge, the most important work was done by Faria [44] in 2016. The author used the Schauder fixed point theorem and relies on the permanence of the system to prove the existence of periodic solutions for n-dimensional periodic delay differential equations, with patch structure and multiple time-varying delays, of the form where all the coefficients and delay functions are assumed to be continuous, nonnegative and periodic on t, with a common period w > 0, and η ik (t,s) are bounded, nondecreasing on s, locally integrable on t and w-periodic.
Periodic solutions and oscillation theory of periodic neutral differential equations have received much attention. In 1974, Hale and Mawhin [60] applied the Fredholm alternative technique and a generalized Leray-Schauder theory to prove the existence of w−periodic solutions of neutral differential equations with w−periodic coefficients of equations of the form where D(t)φ = φ(0) − A(t)φ and A(t) : C → R n for t ∈ R and φ ∈ C. In 1975, Lopes [93] showed the existence and stability of periodic solutions of a class of periodic neutral differential equations, which include many applications to problems arising in transmission lines area, of the form where g is monotonic nondecreasing function, p is periodic and µ, q satisfy some conditions. The relationship between the oscillation of all solutions of a class of linear neutral differential equations and the corresponding multiplier Floquets was established first, in 1990, by Huang and Chen [66]. In particular, they proved that linear neutral differential equations with periodic coefficients of the form (3.15) where 0 ≤ p ≤ 1, r > 0, τ i = ir and q i (t) ≥ 0, i = 1, ..., n are continuous, r−periodic functions, are oscillatory if and only if all the roots of the corresponding characteristic equation are purely complex. The same results were proved by Ladas in 1984 [82] and Huang in 1988 [65], but for non-neutral delay differential equations, (i.e, p = 0), with periodic coefficients, without sign restrictions on the coefficients q i . In 1991, Ladas, Philos and Sficas [81] generalized the result for neutral equations with two delays of the form where Q is w−periodic, n 1 , n 2 are positive integers. More precisely, they showed that every solution of the neutral equations with periodic delay oscillates if and only if every solution of an associated neutral equation with constant coefficients oscillates. For non-neutral delay differential equations, (i.e, p = 0), but with general relationship between delays and period (i.e τ i = ir), the oscillation phenomena for equation (3.15) was established by Philos in 1992 [109] for constant delays and by Kordonis and Philos in 1999 [79] for periodic delays (i,e. τ i = τ i (t)) and periodic integro-differential equations. For nonlinear neutral equations with two delays, In 2002, Tanaka [131] gave a necessary condition of oscillation of all solutions of the following equation where σ = ±1, γ > 0 and h, q are positive functions which can be periodic or not periodic.
In 1990, using Fredholm alternative technique for periodic solutions of some linear neutral equations and a generalized Leray-Schauder theory, Erbe, Krawcewicz and Wu [42] show the existence of periodic solutions for some periodic boundary value problem (BVP) of neutral functional differential equations of the form where f : R × C → R is a completely continuous mapping w−periodic with respect to the first variable. Their approach consists to reduce this BVP to a BVP of retarded periodic differential equations of the form avoiding the neutral character of the BVP. This reduction principle turned out to be a powerful and useful tool to extend the existence results established for the periodic boundary value problems of certain retarded equations to neutral equations.
Zhao [154], in 2008, generalized Theorem 3.2 of existence of interior periodic solutions to dissipative periodic neutral functional differential equations of the form (3.14) where D(t) is now defined as follows: Suppose that A(t, θ) is continuous m × m matrix function and D 0 : R → L(C, R m ) is a C 1 -continuous map such that both A and D 0 are ω-periodic in t and the zero solution of the linear functional equation D 0 (t)y t = 0, y 0 = φ is uniformly asymptotically stable for φ ∈ C D0 =: {φ ∈ C : D 0 (0)φ = 0}. For a given t ∈ R, we define D(t) : C → R m by We have the following (ii) Solutions of (3.14) are uniformly persistent with respect to M 0 . Then system (3.14) has at least one ω-periodic solution x * (t) with x * t ∈ M 0 for all t ≥ 0.
In applications of the above theorem and Theorem 3.2, one may prove the uniform persistence for a given periodic and time-delayed system with the help of the persistence theory for periodic semiflows, see the subsequent section 3.5 for some references.
In 2015, using topological theory, Gao et al. [50] investigated the existence of periodic solutions of a class of Linard type equations with a deviating argument of the form where p > 1, τ ¿0,

Spreading speeds
The theory of asymptotic speed of spread (in short, spreading speeds) for periodic systems with delay has attracted less attention as well. The intuitive interpretation for the spreading speed v 0 in a spatial epidemic model can be stated as follows: If one runs at a speed v > v 0 , then one will leave the epidemic behind; whereas if one runs at a speed v < v 0 , then one will eventually be surrounded by the epidemic. To the best of our knowledge, this theory has been developed first, for delayed periodic systems by Liang, Yi and Zhao in 2006 [88]. Using the Poincaré(period) map method, the authors extended the results by Liang and Zhao in 2007 [89], for spreading speeds traveling waves for monotone autonomous semiflows, to periodic systems with general delay. In the following, we will state the results obtained by Liang, Yi and Zhao [88]: Define, for r ∈ C (the set of all bounded and continuous functions from [t, 0] × R to R k ), and v ∈ C r . Let w > 0 and r ∈ C r with r ≥ 0 be given. A family of mappings {Q t } t≥0 is said to be an w−periodic semiflow on C r provided Q t has the following properties: The mapping Q w is called the Poincaré map associated with this periodic semiflow. Now we are in the position to state the main result in [88]. However, we will avoid some technical assumptions (A1-A5 in [88]) on the w−periodic semiflow and we refer the reader to the same paper for more details.
As application of this theory, Jin and Zhao in 2009 [75] proved the existence of the asymptotic speed of spread for a periodic reaction-diffusion population model with delay.

Persistence
A basic and important ecological question associated with the study of mathematical population interaction models is the long-term coexistence of the involved populations. Mathematically, this is equivalent to the so-called persistence of the populations. We say a population x(t) is persistent if lim inf t→∞ x(t) > 0, and we say a system is persistent if all its populations are persistent. The persistence theory of periodic systems with delay has also received some attention. An excellent survey of persistence of nonautonomous and infinite dimensional systems was given by Hutson and Schmitt in 1992. Persistence aspects for general nonautonomous delay Kolmogorov-type population models of the form x (t) = x(t)f (t, x t , y t ) y (t) = y(t)g(t, x t , y t ), (3.16) where f, g : R×C × C → R are continuous and Lipschitz in (φ, ψ) ∈ C × C, were presented by Kuang and Tang in 1993 [80]. They proved that system (3.16) is uniformly persistent under the following assumptions: (K1 f ) : There exist positive constants δ 1 = δ 1 (f ), δ 2 = δ 2 (f ), K 1 = K 1 (f ), K 2 = K 2 (f ), with K 1 < K 2 such that for all t ≥ 0, f (t, x t , 0) > δ 1 , for x(t + θ) ∈ [0, K 1 ], θ ∈ [−τ, 0], and f (t, x t , 0) < −δ 2 , for x(t + θ) ∈ [K 2 , ∞], θ ∈ [−τ, 0].
Their method is based on the construction of a set of proper autonomous ordinary differential systems whose solutions can serve as lower or upper bounds for the delayed system (3.16) in certain regions.
In 1993, Thieme [134] proved that, under relaxed point dissipativity, uniform weak persistence implies uniform strong persistence. A more general result of uniform strong persistence with applications to periodic delay differential equations was presented by Thieme in 2000 [135]. Zhao in 1995 [152] and 2003 [153], respectively, proved that the uniform persistence of a periodic semiflow is equivalent to that of its associated Poincaré map. Thus, one can use the theory of uniform persistence for discrete-time semiflows to obtain the uniform persistence for a periodic and time-delayed system. In 2005, Xu and Zhao [147] considered the following periodic system of competing mature populationṡ u i (t) = u i (t − τ i )F i (t, u i (t − τ i )) − u i (t)G i (t, u 1 (t), . . . , u m (t)) = f i (t, u 1 (t), . . . , u m (t), u i (t − τ i )), 1 ≤ i ≤ m, (3.17) where the continuous function f i (t, u 1 , . . . , u m , v i ) is T -periodic in t, and Lipschitzian in (u 1 , . . . , u m , v i ) in any bounded subset of R m+1 + , 1 ≤ i ≤ m. By appealing to theories of monotone dynamical systems, periodic semiflows and uniform persistence, they analyzed the evolutionary behavior of this system and established sufficient conditions for competitive coexistence and exclusion.

Conclusions
We intended to give a short survey of some existing studies on periodic systems of delay differential equations. In principle, all general results of periodic processes and maps on infinite dimensional spaces can be used to study the qualitative theory of periodic systems of delay differential equations. Thus, our choice of topics has to be limited to those addressing the interaction of delay and periodicity. Faced with the wide number of results connected with delay systems with periodic coefficients, we hope that this overview has provided some enlightenment to the matter. Despite significant progress in theories for such systems, there are still many questions that need to be discussed for more general cases. Some open questions that come to mind are, for example, explicit stability threshold for general linear delay equations with periodic coefficients, the influence of the periodicity of the coefficients on the intrinsic oscillations of the equations such as those caused by the time delay. Which mechanisms are dominating for specific oscillation and how delay and periodicity interact to generate complicated behaviors remain challenging questions too. Furthermore, spreading speeds and persistence theories has received little attention in the literature. Although some progress known recently for relating theory to application (such as establishing relationship between the basic reproduction number R 0 and the spectral radius of the generation evolution operator and/or the Malthusian parameter), there is still a great challenge to link general analytic results to specific systems arising from important applications.
Delay equations with periodic coefficients are widely considered in epidemiology, ecology, seasonal migration and engineering, however applications in economics and biology (especially in-host modeling for infectious diseases and cell biology) are lacking. Finally, it is possible to argue that it is essential or at least interesting to include stochastic effects, threshold delays and age structures in biological delay models with periodic parameters. Overall, many challenges are thus opening the way to the new generations of interdisciplinary researchers.