On Population Models with Delays and Dependence on Past Values

: The current values of many populations depend on the past values of the population. In many cases, this dependence is caused by the time certain processes take. This dependence on the past can be introduced into mathematical models by adding delays. For example, the growth rate of a population depends on the population τ time units ago, where τ is the maturation time. For an epidemic, there is a time τ between the contact of an infected individual and a susceptible one, and the time the susceptible individual actually becomes infected. This time τ is also a delay. So, the number of infected individuals depends on the population at the time τ units ago. A second way of introducing this dependence on past values is to use non-local operators in the description of the model. Fractional derivatives have commonly been used to provide non-local effects. In population growth models, it can also be done by introducing a new compartment, the immature population, and in epidemic models, by introducing an additional exposed population. In this paper, we study and compare these methods of adding dependence on past values. For models of processes that involve delays, all three methods include dependence on past values, but fractional-order models do not justify the form of the dependence. Simulations show that for the models studied, the fractional differential equation method produces similar results to those obtained by explicitly incorporating the delay, but only for specific values of the fractional derivative order, which is an extra parameter. But in all three methods, the results are improved compared to using ordinary differential equations.


Introduction
Many biological processes depend on the population in past times.For a birth process, the number of births depends on the mature population, that is, the population τ time units ago, with τ the maturation time of the population.In an epidemic, after contact between a susceptible individual and an infected one, it takes a certain time for the infection-causing agent (virus, bacteria, etc.) to reproduce inside the body and make the individual sick.So, the number of newly infected individuals depends on the number of contacts that occurred at a certain time in the past.At the cell level, immune cells take time to become active and eliminate the bacteria or debris.Of course, not all biological processes depend on the past.For a population with a fixed death rate, the number of deceased individuals at the present time depends only on the number of individuals at the present time.Mathematical models based on ordinary differential equations (ODEs) do not consider the influence of the past.But for models with processes depending on past values, incorporating this dependence into the model equations will change the solutions for the better since the model will be more realistic.
There are three widely used methods for introducing the dependence of a process in the past.The first one is to use delay differential equations (DDEs).The second one is to use fractional differential equations (FDEs).The third one is to introduce additional compartments in the model, where the individuals spend an average time before moving out of the compartment.Introducing delays is a straightforward way of adding the dependence of the past to a biological process.Adding delays may change the dynamics of the process [1][2][3].Some possible changes are the introduction of oscillations, the introduction of discontinuities in the time derivatives, the non-uniqueness of solutions, and different stability regions.There is a vast amount of literature spanning many years on epidemic or disease transmission models using delay differential equations (DDEs) (see, for example, [1,[4][5][6][7][8]).For epidemic models at the population level, the time between a susceptible having contact with an infective and actually getting the disease (time of infection) is the most commonly considered delay.But there are also other delays like maturation times, the time it takes a vaccine to become effective, the time it takes an infected individual to recover after getting the disease, etc.For a within-host infection, delays include the time it takes the pathogen to replicate or reproduce and spread over the affected organ, the time it takes the immune system to react, and the time it takes cells to reproduce, among others.There are different types of delays.Discrete-constant delays are used when each delay has the same value for the entire population.Distributed or continuous delays are used when the delay varies according to a given distribution, such as normal or gamma, around the mean value.In this paper, we will only consider discrete delays.Fractional calculus started with a letter from L'Hôpital to Leibniz asking about derivatives in non-integer order.There are several different definitions for fractional derivatives and integrals, with the most common ones being the Riemann-Liouville and the Caputo versions.There are many good papers and books about differential equations using fractional derivatives (FDEs), such as [9][10][11][12].In recent years, a vast number of papers have used fractional derivatives in epidemic models, including many for specific epidemics such as influenza, dengue, and COVID-19 [13][14][15][16][17].The Caputo form is the most widely used fractional derivative.It has the advantages that its derivative of a constant is zero and for FDEs, the initial conditions are given in terms of the unknown and, if necessary, integer-order derivatives at the initial point.Exact solutions of FDEs are hard to find.There are several numerical methods [18][19][20][21].Fractional derivatives are non-local operators that involve the solution from the initial time until the current time.So, the solution at the current time depends on the solution at all previous times, so FDEs incorporate memory.ODEs neglect such effects.In addition, when fitting data, fractional models have more degrees of freedom through the orders of the fractional derivatives compared to ODE models.In fact, many papers that use fractional derivative models for epidemics justify their use by demonstrating superior data fitting compared to ordinary differential equation models.This has been done for different epidemics such as dengue, Ebola, and COVID-19 [14,22,23].Discretetime fractional differential equations have also been used [24][25][26] but not as widely as continuous-time fractional differential equations.Possible reasons for this are that popular methods become more popular and the widespread availability of ready-to-use software for FDEs.A third way of introducing dependence on the past in a population mathematical model is to add more compartments.For models including birth processes, this can be achieved by adding an immature population, with individuals staying in this compartment for an average time equal to the maturation period.Examples involving insect populations can be found in [27,28].In epidemic models, this can be achieved by adding an exposed population [29][30][31].The dependence of solutions on past values is a topic of very high interest, especially in epidemics, with some recent references being [32][33][34][35][36].The main objective of this paper is to study how to introduce dependence on past values into ordinary differential models of population growth and disease transmission.This is discussed in Section 2. In Section 3, results and numerical simulations are presented.Finally, in Section 4, a discussion of the results is presented.

Materials and Methods
In this section, we present the addition of dependence on the past using delays, fractional derivatives, and additional compartments to several compartmental models of population growth, epidemics, and within-host infection propagation.
We start with very simple models and continue to more complicated ones.For all examples, dependence on the past is added to the ODE model using three methods.The first method is to introduce discrete delays.While this can be done in several ways, not all of them make biological sense.The second method is to use non-local fractional derivatives instead of time derivatives.The third method is to add more population compartments, such as immature or exposed populations, and incorporate dependence on the past based on the time individuals spend in these compartments.For the simplest model considered, with one population and described by a linear differential equation, it is possible to calculate exact solutions for the original ODE and all models based on the three methods.For the delay differential equation models, the solution requires numerical approximation of the solution of a transcendental equation.For the logistic model, there is an exact solution only for the ODE, but stability can be found for all models.For the epidemic models, there are no exact solutions, and the local asymptotic stability of the disease-free equilibrium can be determined for all methods.But the expressions become more complicated with the complexity of the model.
The pyramid flowchart depicted in Figure 1 illustrates the sequence of models used in this paper.These models progress from linear, with one state given by the Malthus model, to non-linear, with one state given by the logistic model, and finally, to non-linear epidemic models with several states, each increasing in complexity.ODE refers to the corresponding ordinary differential equations model, DDE refers to the delay differential equations model, EXP refers to the model with immature or exposed populations, and FRAC refers to the fractional differential equations model.E indicates that the model has an exact solution, S indicates the local asymptotic stability of the steady states for the Malthus and logistic models, and of the disease-free equilibrium for the infection models.These results are shown later in this section.In Section 3, the simulations show that all methods that incorporate dependence on the past can produce similar results, except for the fractional derivative equations models, where this holds true only for specific values of the fractional order.

Population Growth Models
Many mathematical population models can be derived by starting with a discrete-time population model [37,38].All these models state that the current population is equal to the population at the previous time plus its change during the time step.For example, for a model considering only birth and death processes, the population at time t + ∆t can be given as the population at time t plus the new births minus the new deaths.The discretetime model is where N(t) is the population at time t, b > 0 is the birth rate, and µ > 0 is the death rate.By taking the limit as ∆t tends to 0, we obtain the simplest continuous-time population model known as the Malthus model For many populations, there is a delay τ > 0 due to the maturation time in the birth process.But the number of deaths depends only on the current population.By adding a delay to the discrete-time model, we obtain where the delay is τ = k∆t, with k being a positive integer.Again, taking the limit as ∆t tends to 0 and keeping k∆t constant, we obtain the delayed continuous-time model that takes into consideration the maturation time.The Malthus model and the delayed Malthus model have only one equilibrium point, N(t) = 0.One issue with using the DDE model is that it requires initial values on the history interval [t 0 − τ, t 0 ], where t 0 is the initial time.
Instead of having a single value, the maturation time may take values in an interval [k i ∆t, k f ∆t], where 0 ≤ k i ≤ k f .The model is where G(t − k∆t) is a weighting function with an average of one.By taking the limit as ∆t tends to 0 and keeping the terms k∆t constant, k = k i , . . ., k f , one obtains the distributed delay equation The function G is taken to have compact support and an average of 1. Usually, it has a very small support.For a discrete delay, it is taken as a Dirac delta function.The birth function bN can be replaced with another birth function or with a function F representing other processes happening in the past.Also, the death term µN can be substituted with a function H representing all the processes happening at the present time.A more general integro-differential equation for the growth of one population can be written as A discrete-time model that considers an immature population is where Y(t) is the immature or young population, all new births are immature, and σ = 1/τ is the rate at which the immature population matures.Taking the limit as ∆t tends to 0, we obtain the ODE model that consists of a system of equations for the two age groups, mature, N, and immature, Y [27,28], ( A modification of Malthus's continuous-time model using fractional derivatives is introduced later in this subsection.
A second population growth model that adds intra-species competition for resources is the logistic growth model given by where N(t) is the population, r = b − µ is the growth rate, and K is the carrying capacity.
The maturation time delay can again be introduced into the birth term: The competition and death terms depend on the current time t.One disadvantage of delay differential equations is that they require initial conditions given on the interval [t 0 − τ, t 0 ], where t 0 is the initial time.This interval is usually called the history interval.Following the approach used in the Malthus model, an alternative way of adding dependence on the past is to introduce an immature population where r = b − µ and σ is 1 over the maturation time.We assume that there is no competition between the adult and the immature populations, as is the case for insects, where the larvae feed on different food than the adults.But there are several other ways of introducing the immature population.Before writing the model using FDEs, we introduce some definitions.The Caputo fractional derivative, which is now the most commonly used fractional operator in mathematical biology models, is obtained by modifying the definition of the Riemann-Liouville fractional derivative to regularize the initial conditions, so they are stated only in terms of ordinary derivatives [9,10].In our definitions of Riemann-Liouville and Caputo fractional derivatives, we consider that the initial time is t = 0.The Riemann-Liouville fractional derivative is and the Caputo fractional derivative is where α is the order of the derivative.These two definitions are not equivalent to each other.Their difference is expressed by .
The Caputo operator C 0 D α t has the advantage for differential equations with initial values in that these initial values are given in terms of ordinary derivatives.The initial values for Caputo differential equations are given as In our examples, m = 1.From the definition of the Caputo derivative, C 0 D α t y(t), it can be seen that it is a non-local operator over the interval (0, t) and that it does not depend on values for t < 0. From the kernel of the integral, it is clear that it is weighted toward the values of y(t) at time t.The replacement of the integer-order derivative with a fractional-order one in population models is usually justified by stating that the population has memory, that is, it depends on the past.But it is not stated how it takes into account the actual dependence on the past and whether it is due to maturation or infection times, or to other processes.
Another definition of the fractional derivative D α t y(t), the Grünwald-Letnikov definition, is based on taking finite differences on equidistant time steps in [0 and use the usual notation of finite differences, Then, the Grünwald-Letnikov definition is [9,19] Take the limit α → 1 to obtain the explicit or implicit Euler method.Compared with linear multistep methods for the approximation of the fractional derivative, the sum of divided differences becomes longer and longer [9,19], so it is not usually used in practice.But it is useful for comparing an approximation of a fractional differential equation to that of an ODE.We apply the Grünwald-Letnikov definition to the following fractional differential equation and let N(t) be the exact solution in the interval [0, T].If N k denotes the approximation of the true solution N(t k ), then the explicit or implicit Grünwald-Letnikov method on a uniform grid is given by Here, r α n+1 = (n+1) −α Γ(1−α) , with Γ the Gamma function.Taking n = 0, the equation for the first step is ).This equation does not agree with the first step of any standard finite approximation of Equation ( 3).The Caputo derivative incorporates past values of the solution starting at the initial time t = 0, but it also has terms that do not come from a standard discrete-time conservation equation.The population Malthus growth model in terms of Caputo fractional derivatives is

SIRS Models
A model commonly used for the propagation of many different diseases in a population is the SIRS model with no demographics [37].It consists of three compartments: susceptibles (S(t)), infectives or infected (I(t)), and recovered (R(t)).A susceptible turns into an infective with a given probability after contact with an infective.An infective recovers with a rate γ.A recovered loses immunity at a rate ν and converts into a susceptible.The total population is assumed constant, S(t) + I(t) + R(t) = N, and the disease is assumed to be short enough so that births and deaths need not be considered.The infection coefficient is denoted by β/N.All parameters are positive.The system of ordinary differential equations (ODEs) describing the model is Since the total population is constant, one equation can be eliminated using R(t) = N − S(t) − I(t).The new system of equations is This continuous-time model can also be derived from a discrete-time model [39,40].There are several different methods of introducing delays in the infection into the model.The first one, which we call SIRS Model 1 [41][42][43], is given by where the hypothesis is that an infective has contact with a susceptible and it takes time τ for this susceptible to turn infective.Until that time, it is counted as a susceptible.The second model is SIRS Model 2 [44], which is described by where the assumption is that after contact between a susceptible and an infective, time τ has to pass before the susceptible turns into an infective and is able to infect another susceptible.
So, there is a new infective after contact between the susceptible and an individual infected a time ago equal to the delay.The third model is SIRS Model 3, which is described by where after contact with an infective, the susceptible leaves the susceptible population but takes time τ to be included in the infective population.It conserves the total population, and the infection contact is between individuals at the same time.These three models can all be derived from discrete-time models.As shown in [45], temporary immunity can be included in the SIRS model.The authors used a distributed delay, but by introducing the delay after exposure τ 1 , the delay due to the minimum duration of the disease τ 2 , and the delay due to the minimum time with immunity τ 3 , we obtain the model These are the equations for SIRS Model 1, but the equations for SIRS Models 2 and 3 are similar.Models with delays introduced in different forms may have the same equilibrium solutions but the actual time-dependent solutions will be different.There is also the question of the influence of the initial conditions given in the history interval.For an epidemic, they can always be given as the equilibrium solutions with no disease.That is, the natural assumption that the epidemic starts at t = 0 can be used.
The SIRS fractional model is The same fractional order, α, was used for all three fractional derivatives.This is not necessary but is usually done for fractional epidemic models [13,14,46].Two examples of papers using different orders are [22,47].These fractional derivative models include dependence on all population values between the initial and the current times, but, again from looking at finite difference approximations, they do not arise from standard discretetime models.But simulations may produce good agreement with real data due to the extra parameters, the fractional derivative orders.
A model that takes into account dependence on the past caused by the delay in transmission by including an exposed population, E(t), is The additional parameter is σ = 1/τ.Models of this form are usually called SEIRS models, and the dependence on the past comes from the average time an individual spends in the exposed compartment.

Plant Virus with Vector Transmission Models
Humans, animals, and plants can all be affected by vector-transmitted viruses.The processes involved are the same.For simplicity, we consider a model for cultivated plants, for which it is reasonable to assume a constant plant population since farmers replace dead plants.Two examples of virus-causing diseases in plants are the cassava virus [48] and the cacao swollen shoot virus [49].Again, we introduce the delay in different ways to see the differences in the solutions.We consider a simple model with two populations of plants: susceptible (healthy), S(t), and infective (already infected), I(t).Usually, plants do not recover so we do not consider a recovered class.There are also two populations of vectors: susceptible, X(t), and infective, Y(t).Vectors are only carriers and do not get the disease so they do not recover.A susceptible plant becomes infected when an infected vector feeds on it, and a susceptible vector becomes infected by feeding on an infected plant.This model is a simplified version of the models presented, for example, in [43,50].
The assumptions of the model are as follows: all new plants and vectors are susceptible; the total population of plants is a constant K; the infection terms between vectors and plants, and vice versa, are of mass-action type; the viruses only kill plants; and neither plants nor vectors recover from the disease.
The system of ODEs for the model is The parameters of the model are taken from [43,50] and are all positive: K is the total number of plants, β is the infection rate of plants due to infective vectors feeding on the plant and infecting it, β 1 is the infection rate of vectors due to feeding on an infected plant and becoming infected, µ is the natural death rate of plants, d is the additional death rate of plants due to the disease, m is the natural death rate of vectors, and Λ is the replenishing rate of vectors (due to birth and/or migration).
In this virus transmission via a vector model, there are two delays.The first one is the time it takes the virus to spread in the plant after infected contact.The second one is the time it takes the virus to spread in the vector after becoming infected through feeding.Since the virus usually stays around the feeding organs of the vector and does not replicate inside it, this second delay is much smaller than the first.For simplicity, we consider only the first delay.
Following the approach used in the SIRS model, we introduce the delay in the transmission in two different ways.The first one uses the assumption that after contact with an infective, a susceptible takes the delay time to become infective itself [42,43].That is, a newly exposed susceptible remains a susceptible until a time equal to the delay elapses and only then does it become an infective.This is known as Model A1: The second model with infection delay uses the assumption that after a contagion, the susceptible immediately stops being susceptible but it takes the delay time to become infective.It takes into account the probability of the newly infected plant dying before becoming infective using the term exp(−µτ), which is the average survival percentage in the time period τ.This is known as Model A2: Delay differential equations, while more realistic, have the drawback that it is necessary to specify initial conditions over an interval instead of just at a point, as with ordinary differential equations.These initial conditions may be unknown, and, if so, are usually assumed to be constant.What is more, they are usually assumed to represent disease-free values.A second drawback is that the analysis is much harder due to the infinite number of eigenvalues [1,51].
A common alternative is to use fractional derivatives, which have their own advantages and difficulties, as already mentioned for the population growth models and the SIRS models.The DFE model is A third alternative for accounting for dependence on the past is to introduce another compartment, the exposed population (E).After contact with an infective, a susceptible is converted into an exposed or latent, one that cannot yet infect.The exposed turns into an infective at a rate of ϵ = 1/τ.The model is Epidemic models with an exposed class are very common (for plant virus propagation, see, for example, [30,52]).The advantages of using models with exposed populations are that they are based on ODEs, so they only require the initial condition at the initial time, and, as with all epidemic models based on ODEs, the local asymptotic stability at the beginning of the epidemic can be determined using the next-generation matrix method [53,54].But, from a biological point of view, there is the complication that an individual only stays in the exposed state for a time equal to the delay, whereas some individuals may stay in this state for a much longer time.

Within-Host Virus Infection Models
A basic model of within-host virus infection considers three populations: susceptible cells, infected cells, and virus particles (virions).Susceptible cells are recruited and die naturally, and they can become infected through contact with a virus particle.Cells infected with a virus burst and release a certain number of virions.These free virions can infect healthy cells or die.More realistic models also include the elimination of infected cells and virions by the effector cells of the immune system, as well as cell-to-cell transmission.The equations for such models are [55][56][57][58] Here, x, y, v, and e are the concentrations of susceptible cells, infected cells, virions, and effector cells, respectively.λ is the recruitment rate of new susceptible cells, µ x is their death rate, β is the infection rate, µ y is the killing rate of infected cells by the virus, B is the number of virions produced per infected cell, and µ v is the death rate of virions.Susceptible cells can be infected through direct contact with an infected cell at a rate β xy ; effector cells eliminate infected cells at a rate β ey ; free virus particles are introduced into susceptible cells at a rate β; and effector cells are recruited at a constant rate s and are also recruited depending on the number of infected cells y at a rate α 1 .In addition, as effector cells eliminate infected cells and virus particles, they are also eliminated from the system with rates β ye and β ve , respectively.Finally, effector cells die naturally at a rate µ e .All parameters are non-negative.We add two delays to the model ( 22), τ 1 , which is the time it takes the virus to replicate after invading a cell, and τ 2 , which is the time it takes the immune system to recruit an effector cell after detecting an infected cell.The delay model is where σ = 1/τ 1 is the rate at which exposed cells become infective cells.We used 5 different models for population growth and infection propagation, each with its own set of hypotheses.Dependence on past values due to processes taking time was added to all of them using three different methods.For all three, their solutions approach the steady states at a slower rate than the one for the corresponding ODE model.While the fractional differential equations models incorporate dependence on the past in a different way, their simulations are similar if a value of α is chosen to make this happen.This value of α depends on the particular model and the value of the delay.

Results
This section presents the results for the models introduced in Section 2, including the numerical simulations.We show that the solutions using the different delay differential models are similar to each other and to the models with immature or exposed populations.We also show that even though the models given by FDEs do not have a dependence on the past based on common modeling assumptions, the additional parameters given by the order of the fractional derivatives enable the solutions to be similar to those obtained by DDEs or by introducing additional delaying populations.

Population Growth Models
We start with the Malthus ODE model (1).It is used to present the results in a simple form since it consists of only one population and one equation, and it is linear.The solution is The solution tends to 0 for b − µ < 0 and to infinity for b − µ > 0. The delayed Malthus model ( 2) has solutions of the form N(t) = exp(λt), with λ a solution of the characteristic equation λ = b exp(−λτ) − µ.This transcendental equation has no exact solutions, but λ > 0 if 1 ≥ exp(−λτ) ≥ µ/b.So, for b > µ, the solution also tends to infinity as t tends to infinity and to 0 otherwise.Some delay differential equations can be solved analytically using the method of steps, as described in [1,59].For example, for (2), taking and so on.The solutions become complicated very fast.For the fractional derivative model (10), an exact solution can also be calculated using fractional tools such as Mathematica [60].The solution is given by N is the Mittag-Leffler function, Γ is the Gamma function, and r = b − µ is the growth rate.Exact solutions of fractional differential equations, except for linear problems, are very few.In this example, from the definition of the Caputo derivative (9), if y(t) is of one sign, C 0 D α t y(t) has the same sign.So, if b − µ > 0, the solution of the fractional equation increases as t tends to infinity.
The Malthus model with an immature population is given in (5).The exact solution is a long formula, but the eigenvalues are given by which implies that the solution will tend to infinity as t tends to infinity for σ(4b + σ) > (2µ + σ) 2 .It can also be seen that by adding the two equations, the equation for the total population is the same as the equation for the ODE Malthus model ( 1).Numerical solutions were calculated for the four Malthus models.The ODEs and DDEs were solved using the julia package DifferentialEqns.jl[61,62], and the FDEs were solved using the package FdeSolver [63], which uses methods reviewed, for example, in [21].The same software was used for all the other models.Figure 2 shows the simulation of the different Malthus models (( 1), ( 2), ( 5), and ( 10)) for several values of the order α.The parameter values used are b = 0.16 (1/day), µ = 0.01 (1/day), τ = 14 (day), N(0) = 0.5, and σ = 1/τ (1/day), with Y[0] = 0 for (5).The values are not chosen for any specific population.The delay model and the model with an immature population have solutions that grow slower than the solution for the non-delay model.As can be seen, by changing the value of α, the fractional solution looks approximately like the delayed solution.As α decreases from 1, the solutions are more delayed.But there is no theoretical way to choose the α without comparing it to either empirical data or other calculated solutions.The fractional model involves memory (dependence on the past), but the memory does not come from modeling a population growth term in a common way (see ( 4)).So, the good agreement is due to the extra parameter in the FDE, α.Next, we studied the logistic model, which consists of one equation but is now nonlinear.The ODE logistic model ( 6) adds intra-species competition to the Malthus model.It has two steady solutions: N * = K, which is stable, and N * = 0, which is unstable, as can be seen by checking the sign of the derivative dN dt .A delay can be added to the death term (7).It has the same steady states and stability as the ODE model.A fractional derivative model obtained by replacing the time derivative with a fractional Caputo derivative model also has the same steady states.It also has the same stability, as can be checked by looking at the sign on the right-hand side.The steady solutions and their stability conditions for the two-population logistic models ( 8) have long expressions.The numerical simulations show similar behavior.The parameters used are b = 0.16 (1/day), r = 0.15 (1/day), µ = 0.01 (1/day), K = 1 (individual), τ = 14 (day), N(0) = 0.5, and σ = 1/τ (1/day), with Y(0) = 0.The graphs are shown in Figure 3. Again, it can be seen that as α decreases from 1, the FDE solution is more delayed and can be made similar to the delayed solution.The α that gives the best agreement with the delay model depends on the value of the delay.For this example, al pha is about 0.8.The model given by ( 6) is more realistic than the one given by (1) since the non-linear competition term keeps the solution bounded.Again, the fractional model produces similar solutions to the delay model provided that α is adequately chosen.

SIRS Models
The SIRS model is an epidemic model considering three populations: S(t), I(t), and R(t), satisfying three non-linear equations.The ODE model (11) has two equilibrium points: the disease-free equilibrium (DFE), S * = N, I * = 0, R * = 0, with N the total constant population, and the endemic one, S * = Nγ/ν, I * = (N − S * )ν/(γ + ν), R * = γI * /ν [38,64].The local asymptotic stability of the DFE can be determined by finding the basic reproduction number R 0 , which is usually calculated using the next-generation matrix method [53,54].For R 0 = β/γ < 1, the DFE equilibrium is locally asymptotically stable.The three delay models, Equations ( 12)-( 14), have the same equilibrium points with the same stability, as determined by linearizing about the DFE and analyzing the characteristic equation.The three models make biological sense.Their main difference is when the newly infected stops being counted as susceptible.For some epidemics, the first model may fit the data better than the others, but this may be due to factors not included in the models.For the FDE model, the equilibrium points are the same as for the ODE model.The stability can be analyzed by solving the linearized equation about the DFE, but the expressions for the solutions are not easy to analyze [11,65].For the SEIRS model ( 16), the DFE is the same, but the endemic one is different due to the extra compartment, E(t).The DFE is also locally asymptotically stable for R 0 = β/γ < 1, also calculated using the next-generation matrix method.Figures 4 and 5 show the simulations for SIRS models ( 11)-( 14) and ( 16), as well as (15), using various values of α.All the solutions for the integer-order models tend to the same endemic equilibrium point, but the delay and exposed models take more time.The solutions of the fractional models also tend to the same equilibrium.Again, for the fractional model, it can be seen that decreasing α slows the solution and that there is a value of α that gives a good fit with the delay model solutions, but the fit is not good for small time values.The best value for α is between 0.7 and 0.8.The parameter values used are N = 100 (individual), β = 1 (1/day), γ = 0.5 (1/day), ν = 0.4 (1/day), τ = 5 (day), and σ = 1/τ (1/day).These parameters are not based on an actual infection.The initial conditions are S(0) = 99, I(0) = 1, which were also used as the history for the delay models, and E(0) = 0 is additionally used for the SEIRS model.

Plant Virus with Vector Transmission Models
The plant virus model introduces the transmission through a vector instead of through direct contact.Now, the infection of the plant depends on the past since it takes time for the virus to spread inside the plant.But the infection of the vector does not.It is a more complicated model, and the objective is to see if the fractional models can still approximate the results of a DDE model by varying the value of α.Two different ways of introducing the delay were presented, (18) and (19).Only model ( 18) conserves the total plant population.It can be argued that the way of introducing the delay in (19) is acceptable for a model with a total population that is not constant, but it is clearly not a good choice for a constant population model.This can also be seen in Figures 6-9.The plots on the left of the figures correspond to the simulations using the models given by ( 17)-( 19) and (21).The plots on the right correspond to the simulation using (20) for different values of α.The parameters used in the simulations were taken from [43,50]: K = 100 (plant), β = 0.01 (1/(vector day)), β 1 = 0.01 (1/(plant day)), µ = 0.1 (1/day), d = 0.1 (1/day), m = 0.2 (1/day), and Λ = 10 (vector/day).These values are for illustration purposes only.All the models, including the fractional ones, have solutions that tend to the same endemic equilibrium, except for the model given by (19), which does not conserve the total plant population.For the FDE models, as expected, the solutions change slower than for the ODE model, and as α decreases from 1, they keep slowing.Assuming that the only dependence on the past is due to the delay in the plant infection, the FDE models can produce similar time trajectories.The value of α that gives the best agreement with the solutions from the delay model is around 0.9.

Within-Host Virus Propagation Models
Finally, we consider a more complicated model in which the populations can interact in several ways, and there are two terms of the equations for the susceptible and infected cells that have a delay.There is also a second delay.The values of the parameters used are Λ = 5 × 10 5 (cells/mL), µ x = 0.003 (1/day), β = 4 × 10 −10 (mL/(cells day)), µ y = 0.043 (1/day), µ v = 0.7 (1/day), B = 5.58, β xy = β/3 (mL/(cells day)), µ E = 0.5 (1/day), α 1 = 2.2 × 10 −7 (1/day), β E y = 0.6 × 10 −3 (mL/(cells day)), β yE = β Ey (mL/(cells day)), β Ev = 4 × 10 −10 (mL/(cells day)), β vE = β Ev (mL/(cells day)), s = 30 (cells/(mL day)), τ 1 = 1 (day), and τ 2 = 24 (day).These parameters are based on those in [55,56] for liver infection by the hepatitis B virus.The initial conditions are x(0) = Λ/µ x , y(0) = 0, v(0) = 0.33, and E(0) = s/µ E .Figures 10-13 show the simulations for the above-mentioned values of the parameters.The plots on the left of the figures correspond to the simulations using ( 22), (23), and (25).For the values of the parameters given above, the solutions tend to a chronic equilibrium state.The DFE is asymptotically unstable for the integer-order models, as was also shown in [66].The plots on the right correspond to the simulations using (24) for different values of the order of the derivative α.As expected, for the number of infected cells, for the fractional model, the rate of increase of the solutions increases as α increases toward 1.The optimal value of α is between 0.9 and 1.Since the solutions increased very slowly, the plots were truncated to clearly show the differences for a short time.Another simulation was conducted with β xy = 0, resulting in no transmission from cell to cell.In this case, the DFE is asymptotically stable, as illustrated in Figure 14.

Discussion
Many populations depend on their values in the past.The number of births depends on the population the maturation time ago.The number of newly infected depends on the number of contacts between an infected and a susceptible the infection time ago.Also, the immune system takes a certain time to react and start fighting an infection.Mathematical models for processes involving values in the past should incorporate this dependence and, therefore, produce solutions that are more realistic.Mathematical models are simplifications of reality.This can be due to keeping the model tractable or to the lack of information and data.ODE-based models assume that the population is homogeneously distributed and so avoid having spatial dependence.Usually, infection rates are assumed to be constant but they can vary with time, temperature, and the individual.There may be saturation effects not modeled by mass-action kinetics.Recovered individuals turning into susceptibles may retain some immunity and may infect others.Processes that take time should be modeled taking into account that time.As shown in Section 2, there are different ways of introducing delays into an ODE model, but they need to be introduced in a consistent manner.For example, model (19) does not conserve the total population and, therefore, should not be used.The maturation time is not a single time but has some variation.The same is true for the infection time.But usually, these variations are small, and using a single value is a good approximation.But if they are not, then a distributed delay model will be better.Also, there is the effect of population values on the history interval, [t 0 − τ, t 0 ], where t 0 is the initial time.But for epidemic or infection models, the values of the populations in those intervals can be assumed to correspond to the DFE, which is a reasonable assumption.There are many papers dealing with comparisons between a variety of epidemic models with delays and real data for diseases like dengue or COVID-19.
A few examples can be found in [67][68][69][70].Most show that the fit is good but do not compare their results with those of other methods.An exception is [70], where the results from a DDE model were compared with those from an ODE model.
A second way to introduce dependence on the past is to replace the time derivative with a fractional one.From the definition, for example, of the Caputo derivative at time t, it is a weighted average of the time derivatives at times from 0 to the current time t.As shown in Section 2, FDEs do not come from discrete-time conservation models by taking the limit as the time increment goes to zero.In models where the dependence on past times comes from the delay in a given process, the FDEs do not have a form similar to constant or distributed delay models.Furthermore, the effect of past times starts with the initial time of modeling, as opposed to models formulated using DDEs, where the influence extends over times t ∈ [t 0 − τ, t 0 ], with τ representing the delay.FDEs introduce memory, but there lacks a justification for the specific form of this memory.For epidemic models, it is true that susceptibles may remember previous infections, potentially leading to a different infection process.But this does not justify the form of the memory in the FDE.It is also true that many epidemic models based on FDEs fit data better than models based on ODEs, but this may be due to the additional parameter α.If there are susceptibles who previously had the infection and thus may react differently to a new infection due to their memory of the previous infection, an alternative is to add a new susceptible compartment.This can also be done for infectives and those who have recovered from previous infections.The main arguments for using FDEs in epidemic models are that they include memory and that they fit data better than models based on ODEs.This is true, but there is no justification that the memory that the FDEs provide has anything to do with any biological, physical, or chemical process with dependence on the past.Also, the value of the fractional order (or orders) is chosen to better fit the data.Mechanistic models should give a modeling reason or hypothesis on how to choose it.Statistical models, like linear regression, will find the values of the involved parameters to better fit a data set.But mechanistic models should do more than provide computational data that fit experimental data.Recently, there has been an awareness that fractional differential equations should also include delays to better model processes that have delays.This allows the fractional derivative to provide an additional memory effect.Some references supporting this notion are [71][72][73].
A third common way to introduce that processes take time, thus introducing dependence on the past, is to add new populations, immature populations for birth processes, or exposed populations or cells to epidemic or infection models.This approach has the advantage of being based on ODEs, and, for example, the next-generation matrix method can be used to determine if there will be an epidemic.Also, it does not require initial values for an interval prior to the initial value of t.But, just as a Malthus death model dN dt = −µN with µ > 0 can have individuals living for very long times, some exposed individuals may stay in the compartment for a long time.Of course, this can happen with a very small probability.
All the example models used have a dependence on past values based on the assumption that the dependence is due to the time it takes processes like infection or maturation to happen.Delay models, as well as exposed or immature models, explicitly include this dependence.Fractional differential equation models have a dependence on past values, but this dependence has a different form that is independent of the actual value of the delay or delays.But fractional differential equation models offer the flexibility of more parameters, specifically the fractional orders of the derivatives.By choosing these orders adequately, these models can produce similar results to models that explicitly include the delays.
In conclusion, the terms added to an ODE model to take into account dependence on past values need justification.The justification should be more than just saying there is memory or achieving a better fit.The processes in mechanistic mathematical models need to be based on realistic assumptions.For processes depending on past values due to a delay, DDEs are good choices since there is flexibility in choosing the form of the delay, whether constant or distributed.But models based on adding more populations and still relying on ODEs are simpler to analyze.When using fractional derivatives, modeling assumptions should also be clearly stated.
Future work will include models with stochastic delays, models with partial immunity, and comparisons of the predicted results with real data.

Figure 1 .
Figure 1.Flow diagram for the sequence of models used, from simplest to most complex.ODE is the ordinary differential equations model, DDE is the delay differential equations model, EXP is the model with immature or exposed populations, and FRAC is the fractional differential equations model.E indicates that the model has an exact solution, S indicates the local asymptotic stability of the steady states for the Malthus and logistic models, and of the disease-free equilibrium for the infection models.

Figure 2 .
Figure 2. Malthus models.(a) ODE, DDE, and immature population models.The legend "Both" refers to the sum of the mature and immature populations for model (5).(b) Fractional models for several values of α.

Figure 3 .
Figure 3. Logistic models.(a) ODE, DDE, and immature population models.The legend "Both" refers to the sum of the mature and the young or immature populations for model (8).(b) Fractional models for several values of α.

Figure 6 .
Figure 6.Plant vector virus models, susceptible plants.(a) ODE, DDE, and exposed population models.(b) Fractional models for several values of α.

Figure 10 .
Figure 10.Virus within-host models, susceptible cells.(a) ODE, DDE, and exposed population models.(b) Fractional models for several values of α.

Figure 11 .
Figure 11.Virus within-host models, infected cells.(a) ODE, DDE, and exposed population models.(b) Fractional models for several values of α.