Memory effects in disease modelling through kernel estimates with oscillatory time history

We design a linear chain trick algorithm for dynamical systems for which we have oscillatory time histories in the distributed time delay. We make use of this algorithmic framework to analyse memory effects in disease evolution in a population. The modelling is based on a susceptible-infected-recovered SIR—model and on a susceptible-exposed-infected-recovered SEIR—model through a kernel that dampens the activity based on the recent history of infectious individuals. This corresponds to adaptive behavior in the population or through governmental non-pharmaceutical interventions. We use the linear chain trick to show that such a model may be written in a Markovian way, and we analyze the stability of the system. We find that the adaptive behavior gives rise to either a stable equilibrium point or a stable limit cycle for a close to constant number of susceptibles, i.e. locally in time. We also show that the attack rate for this model is lower than it would be without the dampening, although the adaptive behavior disappears as time goes to infinity and the number of infected goes to zero.


Introduction
Memory effects are an important part of disease modeling (Roddam 2001;Sofonea et al. 2021;Liao et al. 2022;Espinoza et al. 2022;Ali et al. 2022;Spitzer et al. 2022;Dönges et al. 2022).Non-pharmaceutical interventions and individuals adapting their behavior in response to the news both fall under this category, and the challenge when implementing these in models is the non-local interaction in time, i.e., that the information about number of infected is delayed.
In this paper, we apply kernel methods from animal population dynamics to epidemiological models as the susceptible-infected-recovered SIR -model and the susceptible-exposed-infected-recovered SEIR -model.These models allow analytical calculation of equilibrium points and their respective stability properties.We follow a modeling tradition which is common in theoretical ecology based on a dynamical systems approach with a distributed time delay incorporated, see Murray (2002), Murray (2001), Roos (2014) and Cushing (2013) and the references therein.A notable feature is that the distributed time delay has the significant advantage of being Markovian.Especially the criterion for an outbreak, local stability in time (i.e., for approximately constant number of susceptibles), and the attack rate are of interest.It is known that using non-pharmaceutical interventions such as lockdowns allows control of the system (Bisiacco and Pillonetto 2021), but we here show that stability comes automatically from adaptive behavior.
First, we investigate the properties of an SEIR-model with added feedback on the activity based on the recent history of the number of infected.Second, we discuss the physiological of these memory effects.Technical details can be found in the appendices.where diag(S) denotes a diagonal matrix with the elements of S on the diagonal and the dot denotes differentiation with respect to time t, see for instance Arino et al in Arino et al. (2005).The parameter β is the disease transmission rate and γ is the recovery rate.The SEIR model in ( 1) is scaled to fractions of the total initial population N p such that

SEIR-model with memory effects and its basic properties
(2) where ν ∈ (0, 1) n is the fraction of the total population in group j.The population of susceptible at time t equals S N (t) = N p S(t), similarly for the exposed E N (t) = N p E(t), the infectious I N (t) = N p I (t), and for the recovered R N (t) = N p R(t).Subscript N in S N , E N , I N , and R N refers to the actual number of susceptible, exposed, infected and recovered.We will later look at special cases to simplify certain calculations.
The contact matrix β ∈ (0, ∞) n×n holds the rate of interactions between different groups.Such stratification could be ages (interactions between young and old), physical location (different cities or countries), or species (such as mosquitoes and humans) (Bastin 2012;Li 2018;Martcheva 2015).To study the effects of adaptive behaviour, we promote the contact matrix β to be a function of time and infection numbers in the following way: (4) That is, an integration kernel consisting of a linear series of the functions α k , where α k is proportional to a product of u (k−1) e −σ u with u = t − τ and an oscillating harmonic part.The first term is considered independent of time, and the second term is responsible for memory effects in the dynamics such as adaptive behavior.The 3-tensor kernel α may in general be very complicated as long as it satisfies (5) The kernel will typically decrease as u increases and otherwise the overall shape of β can be chosen to best represent the given data, e.g., whether it goes to 0 at the origin, which determines whether the feedback is immediate.In mathematical terms, we have a lot of freedom to choose a family of functions for the kernel.
We shall consider a kernel with an oscillating part for modelling e.g.seasonal variations or weekly variations due to shifts in behaviour between work and leisure time in the weekends.We use a kernel previously investigated by Ponosov et al. (2004) but now adding an oscillating term as follows The cc stands for complex conjugation of the preceding terms.We have here suppressed the matrix indices for readability, but all parameters can easily be given more indices if needed.We will also focus on a non-stratified model in this paper.
Note that we can consider a kernel with explicit time dependence α T (t, t − τ ) as long as it is of the form α T (t, t − τ ) = c(t) αT (t − τ ).This will allow c(t) to be pulled outside the integral and handled as a part of the rest of the differential equations.This is useful, either when looking at seasonal changes (Chowell et al. 2015;Viboud et al. 2016;Chowell et al. 2016) or sub-exponential growth (Brauer 2019).It makes sense to do this in conjunction with a time-dependent β 0 .To include oscillations in the kernel, depending on the delay time u=t-s, may seem less obvious.But the oscillations in the kernel depending on u could arise from people, who adjust their behavior today on their experience same time last year or same day last week.An example could be risky behavior due to gathering last weekend and in the current weekend people wish to behave less risky by staying home.From a mathematical point of view it is of interest that the linear chain trick can be extended to the case of an oscillating kernel as in Eq. (6).Furthermore, such oscillating kernels appear in physics and here we can mention the delayed Raman response in nonlinear optical fibers (Dudley and Taylor 2010).
The parameters σ ∈ (R + ) n×n are positive real numbers, and the positive integers k takes the values k = 1, 2, . . .N with N ∈ N + .Furthermore, c k , ω, ε k , μ k ∈ R n×n and n ∈ N.

Rewriting the integral kernel as a set of ODEs
In order to solve Eq. (1) numerically using Eq. ( 6) we apply the linear chain trick transforming the integro-differential equation into a set of ordinary differential equations.This transformation implies two advantages.First, the system of integro-differential equations in (1) can be solved numerically using ordinary differential equation solvers without invoking numerical methods for finding the integral parts.Secondly, stability analysis of equilibrium points for (1) can be conducted using methods from ordinary differential equations.We start by writing where we have introduced From the Eq. ( 4) we observe that we need to calculate integrals of the form where and the bar indicates the complex conjugate.
Applying the linear chain trick as presented in Ponosov et al. (2004) we can find differential equations for z ( j) k (t), j = 0, 1, 2, by differentiating the integrals in (10).The integrals in the delay differential Eq. (4) can thus be replaced by a set of differential equations for z ( j) k (t), using the particular form for α k in (6).Differentiating z For k = 1 we have from Eq. ( 8) that We can easily differentiate G k 0 in (8) with respect to t and use the definition of z We now continue by deriving a differential equation for z (1) From Eq. ( 8) using k = 1 we have 1 in ( 8) with respect to t and use the definition of z (1) 10) to obtain the differential equations governing z (1) (1) Differential equations for z (2) k (t) and accordingly we obtain Collecting the above, our aim is to solve the system of differential equations in ( 12), ( 14) and ( 15).Initial conditions are specified for S, I and R. The initial conditions for z k j , k = 1, 2, . . ., N , and j = 1, 2, 3, must also be specified.We have In these expressions we need to specify I (s) for s ∈ (−∞; 0], which is typically very difficult in realistic systems.Though this is of course also a problem if one wants to determine the initial conditions for I in an SEIR-model without adaptive behaviour.In both cases, fitting infection data to a polynomial is tenable (https://covid19.ssi.dk/analyser-og-prognoser/modelberegninger).In our simulations, we assume that I (s) < 0 for s ∈ (−∞; 0], simplifying this problem.

Numerics
For illustration, we restrict ourselves to the simpler form of the kernel That is, Eq. ( 8) for k = 2.For this specific kernel ( 17), the linear chain trick has the following form.We express it in the following real functions We will later set n = 1, but as some results are obtainable for general n, we keep it for now.The following differential equation for β in Eq. ( 4) is and similarly for the helper functions we get This gives us a set of differential equations that describes our system with the initial conditions β = β 0 and That is, we assume there has been no infection before t = 0. Numerical results illustrating the effect of the delay are shown in Fig. 1, using the parameter values in Table 1.We consider the scalar case of the dependent variables S, E, I and R corresponding to n = 1.

Equilibrium points and stability
The results for t → ∞ are of course still determined by depletion of susceptibles, but as long as the number of susceptibles is assumed to be roughly constant, the memory effects have interesting consequences, such as stability of the number of infections on short time scales.This explains why a contact number around 1 is observed more often in real-world systems than a traditional exponential model would suggest.We investigate these in detail with the example kernel in Eq. ( 17).
Leaving out the recovered state R through normalization and assuming approximately constant number of susceptibles S, the Jacobian is Table 1 Parameter values used in Fig. 1 The time unit is one day denoted d.These parameters are chosen to be illustrative rather than realistic Note that u, v, and β do not need indices, but β and β 0 do because we may want n > 1.

Disease free
Let us start with the trivial disease-free equilibrium point The Jacobian ( 22) reduces to Note that all blocks but S j (β 0 ) jk are proportional to the identity matrix, so we may diagonalize that block on its own and relate the eigenvalues x Sβ 0 of diag(S)β 0 to those of J DF .The largest eigenvalue (i.e. the one that potentially can be positive) is So the condition is x Sβ 0 < γ for the disease-free equilibrium point to be stable.This corresponds to a reproduction number above 1 in a normal SEIR-model.This is significant, because it shows that feedback of the kind (17) cannot change whether or not there will be an outbreak, only the severity of it.

Equilibrium point during outbreak
It turns out that there is also an equilibrium point during an outbreak The last condition on I is very difficult to solve in general, especially because the quantity n m= I * m is not invariant under diagonalization of β 0 .We therefore continue with a single group (i.e., n = 1) to see what properties can be divined in this case.This reduces the equilibrium point equations to We require I > 0, which means that the non-damped reproduction number has to be larger than 1, corresponding to an ongoing epidemic.(By non-damped we mean c 0 = 0 where there is no feedback.)It is clear that below this point, we transition to the disease-free equilibrium point.This assumes that the denominator in I * is positive.Note that negative denominator does not make physical sense, as we then get more infected when the contact rate is lowered.We also implicitly assume that I * S. Otherwise the change in susceptibles will play a role in the dynamics.
For the sake of stability analysis, we start by looking at a simplified case where ω = 0, and then treat the full version numerically.First we take the case ω = 0.As discussed above, this is a very physiologically relevant case.For simplicity, we also set = μ = 0 as these may otherwise simply be absorbed in c 0 .Here the equilibrium point for M = 1 is Note an important, but perhaps unsurprising feature here.As c 0 ∼ σ 2 for a normalized version of the kernel in ( 17), the stable level of infection I depends only on the integral of the kernel.So epidemic non-pharmaceutical interventions may be spread out over a period of time, or they may be strict and fast.It gives the same level of infection.
The variables u (0) , v (0) , u (1) , and v (1) drop out of the dynamics, and the relevant part of the Jacobian ( 22) therefore reduces to We will use the Routh-Hurwitz theorem to determine the stability of the system, specifically the formulation from Nordbø et al. (2007).The characteristic polynomial for the Jacobian is where These are all positive, because R N D = β 0 S/γ ≥ 1.Using the definitions from Nordbø et al. ( 2007) From the condition det(D 4 ) > 0 we find There is also a root of det(D 3 ) at this point, but as sign(a 1 ) = sign(a 3 ), there is still a Hopf bifurcation when this criterion is not satisfied (Nordbø et al. 2007).Note an interesting property of Eq. ( 34).The stability of the equilibrium point is independent of c 0 , and the interpretation seems to be the following: As long as the feedback is non-zero, it cannot change the stability of the equilibrium point, only the position.A sign change would give rise to an equilibrium point at negative I , which in turn would change the sign of the integral in Eq. ( 4).This is of course purely formal as I < 0 is unphysical.It is merely to explain the absence of c 0 in the stability condition.Numerics shows that the transition is between a stable equilibrium point and a stable limit cycle, see Fig. 2, left column.2 Phase space and stability for ω and β 0 using n = 1.We keep S constant to illustrate the local stability in time.We choose three set of parameters that are unstable, critical, and unstable for ω = 0, using Eq. ( 34), and then vary ω to analyse changes in stability.The parameters are γ = η = 1/3.5,c 0 = 5.5, σ = 0.5, S = 0.8, = μ = 0.5 for all the configurations and β 0 = 1.319, 1.519, 1.719 for the stable, critical, and unstable configurations respectively (top, middle, and bottom row respectively), corresponding to non-damped generation numbers of 3.70, 4.25, and 4.81 respectively.The initial conditions are chosen according to the equilibrium point from Eq. ( 28), but with I varied to show convergence of different paths.We have omitted runs that exit the interval E + I ∈ [0, 1] as they are unphysical.These typically display large oscillations.We also use max(β, 0) instead of β in the RHS of Eq. 21.We see that ω does indeed allow changes in stability, but mostly around the critical point

Numerical stability analysis of full model
We start by choosing parameters just above, below, and at the critical point in Eq. ( 34) for ω = 0.This shows the transition is between a stable point and a stable limit cycle.When increasing ω, the limit cycles increase radius to the point where the I -state exits the physical interval [0, 1].See Fig. 2.
We then investigate the effects of varying and μ for both the stable and unstable configurations.It turns out that the important transitions here happen at μ = 0 and = 0. See "Appendix A".

Physiological bounds on the parameters
In order to make realistic examples, let us briefly take a look at physiological scope of each parameter.Since we assume vanishing change in S in the section, we must first and foremost require This gives conditions on the feedback parameters through for the kernel (17 , which is the relation for no dampening.The red point are generated by uniformly sampling parameters in the intervals η, γ , so the kernel is always positive.We see that all points lie to the right of the black line, indicating a lower attack rate than one would expect from a non-damped system (c 0 = 0) for the general case and for = μ = ω = 0.These are necessary conditions, but unfortunately not sufficient ones if the trajectories of E and I go above 1 − S on their way to the equilibrium point.This may even be the case if I (0) < I * , see Fig. 2. Note that γ and η are typically equilibrium by the nature of the disease, though γ can be artificially lowered by testing the population and putting positive cases in isolation.In terms of physical size, both σ and ω have units of inverse time and it is unrealistic to have feedback that reacts faster than a day, unless the population have access to their own tests and report the results immediately (and accurately).We therefore consider σ, ω < 1 days −1 .
We also require β ≥ 0 at all times.A rough estimate can be obtained through Eqs. ( 4) and ( 17), where we for simplicity assume If the Eq. ( 21) are implemented as they stand, it is difficult to satisfy Eqs. ( 35) and ( 37) at the same time, because we only control c 0 and σ .However, as non-pharmaceutical interventions hardly are implemented to satisfy Eq. ( 37), we replace β by max(0, β) by hand on the RHS of Eq. ( 21).The numerics show that the trajectory still converges to either a stable limit cycle or a stable equilibrium point, see Fig. 2.

Attack rate
A different interesting question to ask is the effects of the feedback on the attack rate, i.e.R(∞).Note that at the end of the epidemic, the I -state is almost empty, and it is therefore tempting to think that the feedback does not affect the attack rate.This turns out not to be the case, as can be seen from the following.Start with the fraction of the S-and R-equations in the non-stratified model If we assume S(−∞) ≈ 1 and R(−∞) ≈ 0, the first terms simplify.Reapplying d R = γ I dt and shifting the τ -integral, we can rewrite the last term as where R N D,0 denotes the initial non-damped reproduction number.(As β(t = 0) = β 0 , we could also just call this the initial reproduction number.)While the second term is a complicated integral to handle in general, we can draw some conclusions from it.As long as α(τ ) ≥ 0, which is the case for our kernel when cos (arctan(μ/ )) + μ sin (arctan(μ/ )) ≥ −1, we can be certain that the integrand is non-negative too.This means that attack rate for a given initial reproduction number will be lower than one would expect for a standard SEIR-model.See Fig. 3 for a numerical check of this.A more general, but less transparent calculation for a stratified model (n > 1) can be found in "Appendix B".

Physiological interpretation of memory effects
While it is clear that non-pharmaceutical interventions as well as a population reading about high rates of infection in the news and therefore changing their behavior are obvious examples of memory effects, we would like to briefly discuss the physiological interpretation in more details.
If we start with ω = 0, we only have a dampening from the kernel in Eq. ( 17).So whenever the recent number of infectious has been high, the activity decreases.The behavioral part is not limited to individual choice.Governmental non-pharmaceutical interventions, such as local lockdowns or extra tests in areas with high rates of infection, of course restrict activity based on recent history of infection, but it may also come from increased focus in the media leading to more cautious behavior.It may also come naturally, i.e. without testing the population, from the acquired immunity.If the option of moving from S to R is included, a lingering but waning immunity would take the same kernel form, where infectiousness diminishes for some time.Note that self-isolation does not take this form as it is immediate in time, i.e. the individual isolates based on their own illness, not the previous illness of the population.
The interpretation of ω > 0 should be seen more exclusively as behavioral.It allows for a momentary increase in activity based on a previous wave of infection.Unless the disease weakens individuals and makes secondary infection more likely, this will not happen naturally.However, after a period of either lockdown or self-isolation, humans may be inclined to compensate socially, and thus a wave of infection may be followed by first a period of lower activity and then one of higher.In other words, as we are more social than logical beings, an oscillatory kernel may be very relevant for description of adaptive behavior.
Transmission of a disease may also vary over one week due to different social behaviour in workdays as opposed to weekends.In this case the period of the oscillating disease transmission is 7 days.

Conclusion
We have implemented a feedback mechanism in epidemic models and illustrated its properties.This is useful for modelling adaptive behaviour and non-pharmaceutical interventions, which are both central in epidemic control.
The feedback dampening cannot prevent an outbreak from happening, as the stability of the disease-free equilibrium point does not depend on the kernel.Instead, the feedback can create equilibrium points locally in time as well as decrease the attack rate, so the severity of the outbreak can be contained.
The existence of a stable equilibrium point on short timescales is noteworthy, because it explains how the effective generation number of many countries during the COVID19-pandemic stayed consistently close to 1, which is impossible in normal SIR-type models without fine-tuning.
We have also shown that the feedback affects the attack rate R(∞).This is important, because it shows that non-pharmaceutical interventions reduce the number of people that need to be infected on long time scales, instead of simply postponing the time of infection, which would otherwise be reasonable to assume.

Fig. 1
Fig. 1 Plots of SEIR-model with adaptive behaviour.Parameter values are taken from Table 1 and using n = 1, that is we consider S, E, I and R to be scalars.Top left: The classical variables S, E, I , and R. Bottom left: The auxiliary variables from Eq. (18) that pertain to the contact rate.Note that β stays positive for all times.Top right: Plot of the integral kernel from Eq. (17) used in the simulation.Bottom right: Comparison of the I -state for an SEIR-model with the same time parameters, but with and without adaptive behaviour.(The one without adaptive behaviour simply has c 0 = 0

Fig.
Fig.2Phase space and stability for ω and β 0 using n = 1.We keep S constant to illustrate the local stability in time.We choose three set of parameters that are unstable, critical, and unstable for ω = 0, using Eq.(34), and then vary ω to analyse changes in stability.The parameters are γ = η = 1/3.5,c 0 = 5.5, σ = 0.5, S = 0.8, = μ = 0.5 for all the configurations and β 0 = 1.319, 1.519, 1.719 for the stable, critical, and unstable configurations respectively (top, middle, and bottom row respectively), corresponding to non-damped generation numbers of 3.70, 4.25, and 4.81 respectively.The initial conditions are chosen according to the equilibrium point from Eq. (28), but with I varied to show convergence of different paths.We have omitted runs that exit the interval E + I ∈ [0, 1] as they are unphysical.These typically display large oscillations.We also use max(β, 0) instead of β in the RHS of Eq. 21.We see that ω does indeed allow changes in stability, but mostly around the critical point

Fig. 3
Fig. 3 Connection between the initial non-damped reproduction number R N D,0 and the attack rate R(∞)