Non-Markovian SIR epidemic spreading model of COVID-19.

We introduce non-Markovian SIR epidemic spreading model inspired by the characteristics of the COVID-19, by considering discrete- and continuous-time versions. The distributions of infection intensity and recovery period may take an arbitrary form. By taking corresponding choice of these functions, it is shown that the model reduces to the classical Markovian case. The epidemic threshold is analytically determined for arbitrary functions of infectivity and recovery and verified numerically. The relevance of the model is shown by modeling the first wave of the epidemic in Italy, Spain and the UK, in the spring, 2020.

Introduction

The ongoing pandemics of COVID-19, has claimed millions of human lives, caused stagnation of the global economy and excessive load on the healthcare systems throughout the world and changed the normal life. Mathematical models of epidemic spreading are important tools for predicting the effects that the pandemics can have on each segment of the society. They provide support for policy-makers to make adequate decisions in order to partially mitigate the consequences by planning various social distancing measures, preparation of healthcare facilities and appropriate adaptation of the economy.

The spectrum of mathematical models applied for the COVID-19 pandemic ranges from the simplest SIR to rather complex SIDARTHE [1], [2], [3], [4], [5], [6], [7], which are used for assessment of different aspects of the epidemics. One of the major features of these models is their Markovian nature, which considers transitions from one state to another to be independent on the past. As an example, when Markovian property is assumed to hold, an individual that has just become infected can proceed to recovered state with the same probability as another one which has been infected for longer period. This Markovian assumption, encapsulated in constant transition probabilities, or rates, makes the models easier to study analytically. The outcomes of these studies with Markovian approach offer some, and in certain instances satisfactory, assessment of the spreading dynamics. However, growing body of evidence, particularly for the COVID-19, suggests existence of latency period and certain infectivity patterns, with possibility for spreading the pathogen before onset of the symptoms, to which correspond functions that are rather distinct from the exponential distribution which the Markovian models rely on [8], [9]. Although adding one or more compartments for the Exposed, Asymptomatic, Presymptomatic, or Quarantined persons or considering various kinds of delay [10], [11], [12] address such observations to certain extent, they cannot systematically incorporate the observed distributions of the latency period and the healing process.

The non-Markovian setting is inherent in the pioneering works in the mathematical epidemiology by Ross [13], [14], Kermack and McKendrick [15], and in the related field of population dynamics by Böckh [16] and Lotka [17]. However, the more special and mathematically more tractable, Markovian approach has largely dominated in subsequent studies. In the recent time the non-Markovian framework has started to gain more attention in various settings. In one of the pioneering works on the subject [18], Gillespie algorithm is proposed as an adequate tool for numerical analysis of non-Markovian spreading models. The effects of the form of distribution of infection and curing (recovery) times on SIS epidemic model occurring on complex networks in continuous time has been analyzed in several studies [19], [20], [21], [22], [23], [24]. With the introduction of SI*V* model [25] it was suggested that non-Markovian spreading models have capacity to be extended to cover a wide variety of spreading sub-models and variants. Nontrivial distribution of infectious period in an integro-differential SIR model was considered in [26]. In a recent study, non-Markovian SIS model on complex networks, with arbitrary function for infectivity and recovery was proposed [27], in which control theory was successfully applied for determination of epidemic threshold. Another, novel key contributions in the theory of non-Markovian epidemic spreading models can be found in [28], [29]. In those works, with extensive theoretical work on models with integro-differential equations were obtained analytical results about the equilibria and the basic reproduction numbers. Our study adds determination of the epidemic threshold on base on the stability analysis for general distributions of infectivity and healing in a SIR model. By similar approach as in [27] we show how these functions determine the epidemic threshold. The relevance of the model, besides by numerical simulations, is verified by fitting to the observations of the first wave of the epidemic in Italy, Spain and the UK, in the spring, 2020. The predictions of the model are compared with those of the classical Markovian SIR model.

The paper is organized as follows. After providing initial setting of the model in Section 2, we introduce the discrete-time and continuous-time models in 3, 4, respectively, where we also derive the epidemic threshold relationships. The reduction to Markovian case of the model is presented in Section 5, while numerical simulations and discussions are given in Section 6. The paper concludes with Section 7.

Preliminaries

We consider SIR model that has three compartments: Susceptible - S, Infected - I and Recovered - R, with the usual transition S → I → R. Adequately, let the functions S(t), I(t) and R(t) denote the fractions of the population that are in the state S, I and R, at time t correspondingly, with S(t) + I(t) + R(t) = 1 being the conservation condition. The calculation of the fraction of infected individuals I(t), and the dynamics of S(t), I(t) and R(t) will be given below, separately for the discrete and continuous time. To capture the nontrivial dependence of the healing period and the different contagiousness of the infected individuals in different stages of the disease we introduce two functions. The first one, b(τ), captures the infectiousness intensity at which individuals that became infected before time τ are spreading the disease to the susceptible ones. Thus, by simply taking b(τ) = 0 for τ < T 0, one is able to introduce latency period of the infectiousness with length T 0. The second function is the healing one, g(τ), that denotes the probability that the healing takes period τ. To account for asymptomatic transmitters and existence of certain time window when presence of pathogen can be confirmed, one can introduce a reporting function r(τ). It is associated to the probability that the presence of the pathogen can be confirmed at moment τ after contraction with it. The asymptomatic cases are conveniently handled by normalizing the reporting function to value smaller than unity. In the literature, the two functions b(τ) and g(τ) are usually combined in single infectivity function. We pursue by considering discrete- and continuous-time models separately, and provide more details about these functions.

Discrete-time version

In this section we consider evolution in discrete time t and denote the fraction of individuals that have become infected within the continuous-time interval [t – 1,t] with I (t), where for simplicity the length of the interval was taken to be 1. This can be relevant for situations like those when newly infected cases are considered on daily basis. In such scenario, we have discrete-time healing function g(τ) and infection intensity one, b(τ), on which we put the constraint b(0) = 0. The probability that the individual will heal within the first τ time units is G(τ) = ∑ g(ν). We further assume finite duration T of the disease, what implies G(T) = 1 and for practical reasons introduce its complement , to denote the probability that individual has not healed yet within the first τ time units. The function g(τ) also has an interpretation as fraction of individuals that have contracted the disease within the same unit time interval, to become healed later within another unit time interval [τ – 1,τ]. Similar reasoning can be applied for the cumulative functions G(τ) and as well. On base on the classical SIR model, the proposed model of evolution of the compartments is given with the system

One can note that the infected individuals that have contracted the pathogen up to T periods before the current moment t, and which are not healed yet, can contribute to spreading of the disease, with appropriate intensity captured in the function b(τ). We note that in order to determine the fraction of all individuals that are in the infected compartment at given moment, I(t), one should sum over all that were infected in the past, but did not heal up to the given moment

To make the problem completely defined one has to specify the initial conditions for I (t). We assume that they are given for τ = T – 1, T – 2, …, 0. In general this model cannot be solved analytically and should be studied by application of numerical routines.

To get insight into the conditions when epidemic can emerge, one can determine the stability of the disease free state S ⁎ = 1, I ∗ = I ∗ = R ∗ = 0, that is an equilibrium point of the system. Its local stability is established by linearizing the dynamical equations (1) in its neighborhood. By making the linearization in vicinity of S ∗ = 1, I ∗ = R ∗ = 0, one can observe the dynamical evolution of the perturbations δS = S − S ∗, δI  = I  − I ∗, δR = R − R ∗. Under linearization, the perturbations are related with

Let us focus on the infected fraction and make Z-transform on the second equation in Eq. (3). To do so, multiply first both sides of that equation by z − and sum to obtain

By using the Z-transform of the fraction of the population that become infected at unit interval I (t), given as ℐ(z) = ∑ ∞ I (t)z −, the left hand side of Eq. (4) will become

Accordingly, the right-hand side of Eq. (4) can be rearranged as

By using substitution ν = t − τ, the last sum for τ ≤ −1 can be expressed aswhere we have introduced a function ℐ0(τ,  z) that corresponds to the initial conditions. Now, combining the relationships (5), (6), (7) one has

To shorten the notation, one can introduce the following two complex functions

The first one is simply the Z-transform ℰ(z) of the infectivity function , that is a combination of the infecting intensity and healing functions because . The second complex function ℰ0(z) is related to the initial conditions. Now, one has the following relationshipfrom where

From a result in theory of discrete linear time-invariant systems, a sequence (the impulse response of such system) is decaying if the poles of its Z-transform are within the unit circle [30]. Thus, when the poles of the function ℐ(z) of the complex function (11), or the roots of the polynomial z − ℰ(z) lie within the unit circle, the perturbation dies out at infinity. So, the epidemic threshold can be obtained by taking z = 1 in the denominator in Eq. (11), that results inwhich obviously depends on the functional forms of the healing and infection intensity functions.

We should finally note that any initial infection would not shift back the population to the disease-free state S = 1, I = R = 0, but to some endemic S ∗, I ∗ = 0, R ∗ = 1 − S ∗. However, if the conditions are not favoring epidemic both equilibria will be rather close S ∗ ≈ 1.

Continuous-time version

We will pursue similarly to the discrete-time approach, where the fractions of individuals within given compartment and the functions modeling the infectivity and healing are defined for continuous time t and we again assume finite healing period T. The fraction of infected individuals is conveniently modeled with the rate of infection, or the fraction of newly infected individuals I (t) within the infinitesimal interval (t – dt,t). The total fraction of infected persons is given with the integralwhich accounts for those that had become infected in the past and have not healed yet. Now, the dynamical evolution of the respective fractions is given with

One should note that in their original approach, the general version of the model by Kermack and McKendrick assumes dependence of the infectivity on the age of infection just as the last relationships (Eq. (14)) suggest [15], [31]. In order to determine whether the initial perturbation will grow into epidemics, one could focus on the second equation in the vicinity of the disease-free state S ∗ = 1, R ∗ = I ∗ = 0. Then, the perturbation of newly infected individuals will evolve aswhere it is assumed that in vicinity of the disease-free state S(t) ≈ 1. Now, make Laplace transform of the perturbation of the rate of infection, ℐ(s) = ∫0 ∞ δI (t)e − dt and use it in the last eq. (15). To do that, we will follow the same approach as in the discrete-time version. Multiply both sides with e − and integrate. The left hand side will result in the Laplace transform of δI (t), while the right hand one will be

The last integral can be expressed with

Now, one has

Similarly to the discrete-time case we can introduce the Laplace transform of the infectivity function and its initial conditions contribution

Finally, one obtainsfrom where the Laplace transform of the perturbation of the infection rate is

From the results of control theory, a continuous-time linear time-invariant system is stable if the poles of its transfer function, or Laplace transform of its impulse response have negative real part [30]. Thus, the perturbations δI (t) will decay if the poles of its Laplace transform ℐ(s) (Eq. (21)), or eigenvalues of the system (14) lie within negative half-plane Re{s} < 0. Then, the epidemic threshold can be obtained with s = 0 which leads tothat represents the relationship, which corresponds to the discrete-time case (12).

Markovian SIR model

In order to obtain the classical Markovian SIR model in discrete time, from the non-Markovian case (1), one should consider taking T → ∞, b(τ) = β and G(0) = 0, g(τ) = γ(1 − γ) for τ ≥ 1, where β and γ are constants. This further yields G(τ) = 1 – (1 − γ) and . First, one could observe that by using constant infectivity b(τ) = β in the first relationship of the model (1) and using Eq. (2) one will obtain the classical form for evolution of the susceptible population

Next, by implementing the condition G(0) = 0, and the relationship one can drop the first term in the sum in the recovered population in Eq. (1), and further obtainfrom where, for T → ∞, the recovered population evolves as

Finally, from the conservation relationship I(t) + S(t) + R(t) = 1, one can find that the infected fraction is given as

The relationships (23), (25), (26) represent the classical SIR model in discrete time.

As an example, in Fig. 1 we make a comparison between numerical solutions of the discrete classical SIR model given with Eqs. (23), (25), (26) and the non-Markovian form (1) that reduces to it for the infectiousness intensity function β(τ) = β and the healing one γ(τ) = γ(1 − γ), with T → ∞. The matching confirms that the classical model can be obtained as a special case of the more general non-Markovian model.

Fig. 1

Comparison between the discrete classical SIR model and the classical SIR - equivalent model obtained from the non-Markovian form, for β = 0.2, γ = 0.03. It is used rather large finite duration of the healing T = 150, as a proxy for T → ∞.

Similarly to the discrete-time version, to verify that the proposed continuous model is generalization of the classical, Markovian SIR model, one should consider two characteristics of the latter: 1. The infection rate is independent on the moment when the disease was contracted b(τ) = β; and 2. the duration of infectivity is infinite and exponentially distributed which implies that the healing function is g(τ) = γe −. We note that the respective cumulative distribution is G(τ) = 1 − e −, and accordingly . By using the functional form of the healing function, the total infectious population will be

Similarly, by using b(τ) = β, for the dynamics of the susceptible fraction one hasthat represents the corresponding relationship in the classical SIR model. Furthermore, by applying the functional form for the healing function, the dynamics of the recovered population will be as followsthat is the respective relationship in the classical SIR model. Finally, by using the conservation principle S(t) + I(t) + R(t) = 1, the total infectious fraction will evolve asthat is the remaining familiar relationship from the classical case. As a final note, we just mention that using respective forms for the infectivity and recovery functions for the Markovian case in the epidemic threshold relationships (12), (22), one will obtain the familiar threshold β t = γ.

Numerical experiments and discussion

Our numerical experiments with the proposed model were based on solution of the integro-differential equations for the continuous-time case. We have used the Euler method with step Δt = 0.01. For the functions for recovery and infection intensity were selected those that match the observations as evidenced from the literature. As a proxy for the infection intensity function we have used the incubation period distribution, that quantifies the period from contracting the pathogen to the onset of symptoms. It was observed that the infectivity starts nearly the onset of symptoms [32], so the moment of appearance of symptoms might be considered as start of the infectiousness. We did make this choice since there are available estimates of the incubation period distribution in the literature. As suggested in [8] the incubation period can be conveniently represented with Weibull probability density function W(τ;  α;  λ) = αλ(τλ) e −(, with parameters α = 2.04 and λ = 0.103. This function was further truncated to have support of 35 days and was normalized. We have also assumed that once becoming infectious, the infected person has constant capability of transferring the virus, and thus the distribution of the onset of infectiousness is exactly the infection intensity function. The daily recovering probabilities were modeled with log-normal probability density function , with parameters , σ 2 = ln (1 + σ 2/μ 2) chosen to match a mean value of μ  = 21 and standard deviation σ  = 6. The distribution is then normalized to 61 days, and time-shifted for 4 days in order to exclude immediate recovery. This results in the healing function g(τ) with mean recovery time of 25 ± 6 days, in the following fashion

This construct was based on the results from [33], [34], assuming that: 1. Onset of symptoms (on average) occurs after four days (the time shift); 2. It takes another 7–10 days from onset of symptoms to diagnosis confirmation and hospitalization; 3. Another 10–11 days, on average, are needed from hospitalization to recovery. The period of T = 65 days is considered in order to include even most extreme cases in which hospitalization exceeded 40 days.

Furthermore, we have chosen to decompose the infectivity function b(τ) = βB(τ) into scale parameter β and shape B(τ) that has the form of the above mentioned Weibull distribution with the appropriate truncation. The threshold value of the parameter β th was obtained from condition (22)

To confirm the value of the epidemic threshold we have varied the infectivity parameter β in vicinity of the critical value obtained from Eq. (32) and run the continuous-time model for total time equal to 5000. The final values of the susceptible and recovered fraction are plotted as function of the infectivity parameter in Fig. 2 . As one can see, once β is larger than its critical value, the epidemic emerges.

Fig. 2

Fractions of susceptible (red stars) and recovered (blue dots) individuals at the end of the epidemic as a function of the scaling of the infectivity function β given in terms of its threshold value βth. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

We have further tried to see what are the predictions of this model for the COVID-19 pandemic by using a value of the infectivity parameter β that nearly matches the growth patterns of the epidemic in the countries before countermeasures were applied. As was obtained in a detailed study [35], the doubling time in the first epidemic wave across different countries is approximately three days. For that reason, we have opted to make an experiment with the value β = 4.85β th that produces such growth. We have numerically verified that in the initial stage of the epidemic, the newly confirmed daily cases and the total number of infected individuals grow with the same rate, and have the same doubling time of about three days. Then, by running the model with β = 4.85β th for very long time, it was obtained that at the end less than 1 % of the population will remain susceptible! This result means that, if the doubling time is three days and free spreading of the virus is allowed, then nearly everyone, unless vaccinated (this refers to the original strain of the virus), would eventually contract the disease. This is a particular challenge of the model that should be addressed carefully. One approach is to make more appropriate choice of the healing and infectivity functions and the related parameters.

We have finally attempted to check how well the proposed model can explain the observed shape of the function of the reported cases. To do so, we have used the COVID-19 data from Our World in Data, for three European countries: Italy, Spain and the UK. Our focus was set on the first wave of the pandemic, since in its beginning no preventive measures were used and thus the model parameters could be considered as constant. We have opted first to make more detailed study of the epidemics in Italy, where the wave was the strongest. The other two counties were chosen to verify that the approach has general applicability.

There are three key dates for each country in this study: the start of continuous report of new cases, the lockdown start and the day of the peak of the reported cases. For each country these are summarized in Table 1 . We have used two different values of infectivity parameter β: one for the period before the lockdown start, while another one, smaller than the threshold β th, for the period that follows. The initial condition was set to I (0) = 10−7, that is one case in ten million inhabitants. We have chosen to apply detection of the infected individuals based on a function that has identical form as the infection intensity one, but which is delayed for certain number of days. This corresponds to situation that only those with symptoms are tested, and their appearance is delayed few days after the contraction of the virus. Also, there is certain delay that corresponds to the whole process from onset of symptoms, to visit to hospital to obtaining positive result. We note that this reporting, or testing function was normalized to 0.8 that corresponds to assuming existence of 20 % asymptomatic cases [36]. To reach good fit to the observations, we had to take the start of the simulation, as an assumed epidemic onset, to be D days before the beginning of the period when we compare the predicted daily cases with the actual data. Its exact value was tuned by fitting the logarithms of the daily detected cases from the simulation to the respective ones from the data. More precisely, we have looked for a shift D , that will result in minimal squared error of the following sum

Table 1

Key dates in 2020 for the first wave of the epidemic.

Country	Population	Start of study	Lockdown	Peak
Italy	60 million	February 21	March 10	March 21
Spain	47 million	February 25	March 14	March 25
UK	67 million	February 25	March 23	April 7

In the last relationship T is the duration of the period of comparison of the simulations with the observations, from the start of study (column three in Table 1) to the epidemic peak (column five in the same table). T is another free parameter that corresponds to the average number of days from contraction of the disease to reporting - reporting period. As can be seen, in the optimization procedure we have also varied the infectivity parameter β given in terms of its threshold value β th. Since the optimization procedure is computationally demanding we have chosen the value for β for the lockdown period as one that provides good fit for the whole first wave by visual inspection. As evidenced in Fig. 3 , the value β = 0.75β th is a good choice. We have used this value for the all three countries as constant. The values for the other parameters obtained by minimization of the squared error (Eq. (33)) for the three countries are summarized in Table 2 .

Fig. 3

Daily confirmed cases in the first epidemic wave in Italy in spring 2020 (in blue squares), compared to numerical simulations of the model. Confirmation function is delayed for five days, while the infectivity parameter in the lockdown period is: β = 0.5βth (green circles), β = 0.75βth (red stars), and β = βth (magenta crosses). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Table 2

Parameters quantifying the epidemics. Top three rows - non-Markovian model; the middle ones - Markovian with detection of non-healed cases only; bottom three rows - Markovian with detection of all cases. The asymptomatic cases are not reported for all three models.

Country	β	γ	De	Tr	ϵ
Italy	3.9	–	54	3	2.97
Spain	6.5	–	31	3	7.51
UK	4.7	–	33	1	8.93
Italy	0.64	0.41	50	14	2.30
Spain	0.75	0.42	33	14	6.66
UK	0.59	0.35	42	17	7.68
Italy	1.07	0.86	17	6	3.35
Spain	0.78	0.46	7	6	7.99
UK	0.52	0.29	6	5	9.23

For comparison, we have also considered the classical SIR model given by Eqs. (28), (29), (30), to see its performance for capturing the shape of the function of reported cases in the first wave. The number of new daily infections was obtained by integrating the first term of the right hand side of Eq. (30) asfor period of one day. The optimal parameters of this model were estimated identically as for the non-Markovian case, so that they correspond to minimal squared mismatch of the logarithms of the reported and predicted daily cases, as in the eq. (33). The epidemic was initialized in the same way with one in ten million individuals for all three countries. We have assumed again that only 80 % of the positive cases are detected and after the start of the lockdown the infectivity parameter was set as β = 0.95β th, since value 0.75 of the threshold one makes rather fast decline of the new cases. The detection of the new cases was considered to be delayed as in the non-Markovian model for certain number of days T . Since, in the Markovian approach immediate healing is possible, we have considered two scenarios. In the first one, we have assumed that all new cases, except the asymptomatic ones that account for 20 % are reported. In the second, that is more appropriate for the Markovian model we have considered that only those that have not healed yet are reported. Since the number of the individuals that have not recovered decays exponentially, the number of reported ones at day t was calculated as

In Fig. 4 are presented the theoretical predictions of the proposed non-Markovian and the Markovian models for the three countries, while in Table 2 are given the optimal parameters with the squared error ϵ as a measure of the prediction accuracy. As one can see, the Markovian model with delayed detection provides best match. However, as one can notice that this is achieved with very long delay of the reporting, that largely exceeds the observed incubation period. Another inconvenience is the fact that majority of the infected individuals would not be detected, because they would be healed when the reporting is delayed for so long time. The reports for COVID-19 do not support this result.

Fig. 4

Daily cases of the three countries Italy (top panel), Spain (middle), and the UK (bottom) in the first wave of the epidemic and the predictions of the non-Markovian and two versions of the classical Markov SIR model. The model parameters are given in Table 2. The meaning of the symbols is as follows: Blue squares - official data, red stars - non-Markovian model, green circles - Markovian model without considering healing before detection, magenta crosses - Markovian model with considering healing before detection. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Another peculiarity is the prediction of the peak of the reported cases, which is known to appear more than ten days after lockdown is introduced. For example, the peak of the theoretical curve for the non-Markovian model appears three days before the observed peak in Italy. The best-fit Markovian model misses the peak for one day only, but it features a sudden drop of the number of new infections. To remind, in this model the detection is delayed for two weeks after the contraction of the pathogen. The other Markovian model misses the peak for 11 days. Similar results are observed for the other two countries. Markovian model can predict the peak better, but with unreasonably long delay of detection.

Although providing natural framework for incorporation of observed distributions of the infectiousness of the infected individuals and the typical development of the disease, the proposed model has drawbacks as well. First, before using it, one needs to specify the functions modeling the infectiousness, healing and discovering the infected individuals. Their determination is a far from trivial task and needs careful analysis of epidemiological and medical data. As more complex one, the tuning of the model would need in general more data than the classical Markovian counterparts. Also, its full specification needs providing initial conditions that represent a high-dimensional vector, or an interval of values. How all these factors shape the outcome of the model, and how much is it robust to perturbations of any kind is unknown. We believe that their understanding could provide the epidemiologists with valuable information for better understanding of the possible outcomes of epidemics with pronounced non-Markovian nature.

Conclusions

The proposed general non-Markovian epidemic spreading model captures the typical patterns of the disease in person infected with SARS-CoV-2: delayed onset of symptoms and potential to infect the others and impossibility of immediate cure of those that will become sick. We have studied both discrete- and continuous-time versions and derived analytically the relationships for determination of the epidemic threshold. The model reduces to the classical SIR model with the corresponding choice of the functions of infection and healing. The theoretical analysis was supported by numerical confirmation of the epidemic threshold values. The good fit of the model to the real data shows its promising potential for application for modeling the spread of other infectious diseases. As compared to the classical SIR model this approach is able to reproduce the observations in more natural way. By introducing other appropriate functions one could possibly generalize this model to versions that include other compartments that correspond to hospitalized, quarantined, or deceased persons.

Although the epidemic threshold as key quantity was determined, we did not calculated the basic reproduction number R 0, that represents another important quantity. Furthermore, the relationship between the scaling of the infectivity function β/β th from one side and R 0 and the doubling time, from another should be explored as well. With this regard, we think that it is even more important to determine the herd immunity level needed to prevent the epidemic. Finally, analysis of epidemic spreading by nontrivial contact patterns, modeled with complex networks, and by incorporating the proposed approach could provide further insight in the evolution of the epidemics. These issues could provide better understanding of the non-Markovian setting in modeling the epidemic spreading.

CRediT authorship contribution statement

Lasko Basnarkov: Conceptualization, Methodology, Software, Writing – original draft, Writing – review & editing. Igor Tomovski: Conceptualization, Methodology, Writing – review & editing. Trifce Sandev: Conceptualization, Writing – review & editing. Ljupco Kocarev: Conceptualization, Supervision.

Declaration of competing interest

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Trifce Sandev reports financial support was provided by .