Validity of Markovian modeling for transient memory-dependent epidemic dynamics

The initial transient phase of an emerging epidemic is of critical importance for data-driven model building, model-based prediction of the epidemic trend, and articulation of control/prevention strategies. In principle, quantitative models for real-world epidemics need to be memory-dependent or non-Markovian, but this presents difficulties for data collection, parameter estimation, computation and analyses. In contrast, the difficulties do not arise in the traditional Markovian models. To uncover the conditions under which Markovian and non-Markovian models are equivalent for transient epidemic dynamics is outstanding and of significant current interest. We develop a comprehensive computational and analytic framework to establish that the transient-state equivalence holds when the average generation time matches and average removal time, resulting in minimal Markovian estimation errors in the basic reproduction number, epidemic forecasting, and evaluation of control strategy. Strikingly, the errors depend on the generation-to-removal time ratio but not on the specific values and distributions of these times, and this universality will further facilitate estimation rectification. Overall, our study provides a general criterion for modeling memory-dependent processes using the Markovian frameworks.


INTRODUCTION
When an epidemic emerges, the initial transient phase of the disease spreading dynamics before a steady state is reached is of paramount importance, for two reasons [1][2][3][4][5][6][7][8][9][10][11][12][13][14].First, key indicators or parameters characterizing the underlying dynamical process and critical for prediction and articulation of control strategies, such as the generation time, the serial intervals and the basic reproduction number, are required to be estimated when the dynamics have not reached a steady state.Second, it is during the transient phase control and mitigation strategies can be effectively applied for preventing a large scale outbreak.Prediction and control depend, of course, on a quantitative model of the epidemic process, which can be constructed based on the key parameters estimated from data collected during the transient phase.In principle, since the dynamical processes underlying real-world epidemics are generally memorydependent in the sense that the state evolution depends on the history, a rigorous modeling framework needs to be of the non-Markovian type, but this presents great challenges in terms of data collection, parameter estimation, computation and analyses [15][16][17].The difficulties can be significantly alleviated by adopting the traditional memoryless, Markovian framework [18][19][20][21][22][23][24][25][26][27][28][29].An outstanding question is, under what conditions will an approximate equivalence between non-Markovian and Markovian dynamics hold during the transient phase of the epidemic?Additionally, another important issue remains, how are the errors of Markovian estimation determined?The purpose of this paper is to a comprehensive answer to these questions.
The COVID-19 pandemic has highlighted the need and importance of understanding disease spreading and transmission to accurately predict, control and manage future outbreaks through non-pharmacological interventions and vaccine allocation strategies [1][2][3][4][5][6][7][8][9][10][11].To accomplish these goals, accurate mathematical modeling of the disease spreading dynamics is key.In a general population, epidemic transmission occurs via some kind of point process, where individuals become infected at different points in time.It has been known that point processes in the real world are typically non-Markovian with a memory effect in which the distribution of the interevent times is not exponential [3,26,[30][31][32][33][34][35][36][37][38][39][40][41][42][43][44][45][46][47].For example, the interevent time distribution arising from the virus transmission with COVID-19 is not of the memoryless exponential type but typically exhibits memory-dependent features characterized by the Weibull distribution [3].Strictly speaking, from a modeling perspective, disease spreading should be described by a non-Markovian process.A non-Markovian approach takes into account historical memory of disease progression, mathematically resulting in a complex set of integro-differential equations in the form of convolution.
There are significant difficulties with non-Markovian modeling of memory-dependent disease spreading.The foremost is data availability.In particular, while standard epidemic spreading models are available, the model parameters need to be estimated through data.A non-Markovian model often requires detailed and granular data that can be difficult to get, especially during the early stage of the epidemic where accurate modeling is most needed [15].From a theoretical point of view, it is desired to obtain certain closed-form solutions for key quantities such as the onset and size of the epidemic outbreak, but this is generally impossible for non-Markovian models [18,48].Computationally, accommodating memory effects in principle makes the underlying dynamical system infinitely dimensional, practically requiring solving an unusually large number of dynamical variables through a large number of complex integro-differential equations [16,17].In contrast, in an idealized Markovian point process, events occur at a fixed rate, leading to an exponential distribution for the interevent time intervals and consequently a memoryless process.If the spreading dynamics were of the Markovian type, the aforementioned difficulties associated with non-Markovian dynamics no longer exist.In particular, a Markovian spreading process can be described by a small number of ordinary differential equations with a few parameters that can be estimated even from sparse data, and the numerical simulations can be carried out in a computationally extremely efficient manner [18][19][20][21][22][23][24][25][26][27][28][29].For these reasons, many recent studies of the COVID-19 pandemic assumed Markovian behaviors to avoid or "escape" from the difficulties associated with non-Markovian modeling [4][5][6][7][8][9][10][11].The issue is whether such a simplified approach can be justified.Addressing this issue requires a comprehensive understanding of the extent to which the Markovian approach represents a good approximation to model non-Markovian type of memory-dependent spreading dynamics, and specifically of the conditions under which the Markovian theory can produce accurate results that match those from the non-Markovian model.
There were previous studies of the so-called steady-state equivalence between Markovian and non-Markovian modeling for epidemic spreading.In particular, when the system has reached a steady state, such an equivalence can be established through a modified definition of the effective infection rate [21,24,25].From a realistic point of view, the equivalence limited only to steady states may not be critical as the transient phase of the spreading process before any steady state is reached is more relevant and significant.For example, when an epidemic occurs, it is of fundamental interest to estimate the key indicators such as the generation time (the time interval between the infections of the infector and infectee in a transmission chain), the serial interval (the time from illness onset in the primary case to illness onset in the secondary case), and the basic reproduction number (the average number of secondary transmissions from one infected person), but they are often needed to be estimated when the dynamics have not reached a steady state [3][4][5][6][7][8][9][10][11][12][13][14].It is the equivalence in the transient dynamics rather than the steady state that determines whether the transmission features in the early stages of a memory-dependent disease outbreak can be properly measured through Markovian modeling.Moreover, it is only during the transient phase control and mitigation strategies can be effective for preventing a large scale outbreak.Discovering when and how a non-Markovian process can be approximated by a Markovian process during the transient state is thus of paramount importance.To our knowledge, such a "transient-state equivalence", where the Markovian and non-Markovian transmission models produce similar behaviors over the entire transient transmission period, has not been established.In fact, the conditions under which the transient equivalence may hold are completely unknown at the present.
In this paper, we present results from a comprehensive study of how memory effects impact the Markovian estimations in terms of the errors that arise from the Markovian hypothesis.We consider both the steady-state and transient-state equivalences between non-Markovian and Markovian models.We first rigorously show that, in the steady state, a memory-dependent non-Markovian spreading process is always equivalent to certain Markovian (memoryless) ones.We then turn to the transient states and find that an approximate equivalence can still be achieved but only if the average generation time matches the average removal time in the memory-dependent non-Markovian spreading dynamics.Qualitatively, the equality of the two times gives rise to a memoryless correlation between the infection and removal processes, thereby minimizing the impact of any memory effects.We establish that the equality gives the condition under which Markovian theory accurately describes memory-dependent transmission.
One fundamental quantity underlying an epidemic process is the basic reproduction number R 0 .Our theoretical analysis indicates that, when the average generation and removal times are equal, the transient-state equivalence between memory-dependent and memoryless transmissions will minimize the error of Markovian approach in estimating R 0 and lead to its accurate epidemic forecasting and prevention evaluation.Another finding is that the generation-to-removal time ratio plays a decisive role in the accuracy of the Markovian approximation.Specifically, if the average generation time is smaller (greater) than the average removal time, the Markovian approximate will lead to an overestimation (underestimation) of R 0 and epidemic forecasting as well as the errors of the prevention evaluation, which can also be verified based on readily accessible clinical data of 4 types of real diseases.
Strikingly, the estimation accuracy is largely determined by the time ratio and hardly depends on the particular forms of time distributions or the specific values of the average generation and removal times.This property is of great practical significance, because it is in general challenging to obtain the detailed distributions of the generation and removal times in the early stages of the epidemic [49], but their average values can be reliably estimated even during the transient phase [50][51][52][53][54].Moreover, based on this property, we have developed a semi-empirical mathematical relationship that connects the errors in estimating R 0 with the generation-to-removal time ratio.This relationship holds practical value as it can be utilized to rectify errors in real-world scenarios.And the rectification of R 0 and epidemic forecasting can be accomplished through our web-based application [55].
Overall, our study establishes a general criterion for modeling memory-dependent processes within the context of Markovian frameworks.Once the condition for the existence of a transient-state equivalence between Markovian and non-Markovian dynamics is fulfilled, epidemic forecasting and prevention evaluation can be carried out using the Markovian model, again based solely on the data collected from the transient phase.

RESULTS
The overall structure of this work is depicted in Fig. 1.The section titled "Model" presents the Model building of the age-stratified Susceptible-Infected-Removed (SIR) spreading (Fig. 1a), highlighting the difference between the Markovian (memoryless) and non-Markovian (memory-dependent) dynamics (Fig. 1b).In the "Dynamical equivalence" section, we demonstrate the equivalence between Markovian (memoryless) and non-Markovian (memory-dependent) dynamics for steady state and transient dynamics, which will further lead to the accurate description of memorydependent dynamics by the Markovian theory (Fig. 1c).The section "Markovian approximation of memory-dependent spreading dynamics" analyzes the errors of the Markovian approach in estimating R 0 , epidemic forecasting, and prevention evaluation (Fig. 1d).

Model
We articulate an age-stratified SIR spreading dynamics model, in which the entire population is partitioned into various age groups with intricate age-specific contact rates among them.The distribution of population across different age groups is represented by an age distribution vector (p), with the age-dependent contact matrix (A) quantifying the transmission rates between different age groups.Both p and A can be constructed from empirical data [56,57]  Each individual belongs to one of the three states: susceptible (S), infected (I), or removed (R).When infected (i), a susceptible individual will switch into the I state (ii) and gain the ability to infect others (iii).An infected individual is removed (through recovery or death) with a probability (iv).(b) Non-Markovian versus Markovian process.The infection capacity of an infected individual is characterized by the infection time distribution ψ inf (τ ) and its removal can be described by the removal time distribution ψrem(τ ).For the non-Markovian process, the distributions can assume quite general forms, while the distributions are exponential for a Markovian process.(c) Equivalence between non-Markovian and Markovian processes: (i) steady-state equivalence holds under all conditions; (ii) transient-state equivalence only holds only when Tgen is equal to Trem.(d) Markovian estimation of memory-dependent process.(i) The initial phase of the Monte Carlo simulation is used to fit the parameters according to the Markovian theory.(ii) Important issues such as the estimation of R0, epidemic forecasting, and the evaluation of the vaccination strategies can be addressed by the theory.(iii) The remaining data generated by the Monte Carlo simulation is used to test the accuracy of the estimated R0, epidemic forecasting, and prevention evaluation.
convenience, to distinguish between the actual dynamical process and its theoretical treatment, throughout this paper we use the terms "memory-dependent" and "memoryless" to describe actual spreading processes and Monte Carlo simulations, while in various theoretical analyses, the corresponding terms are "non-Markovian" and "Markovian." The mechanism of disease transmission across different age groups and the recovery or death of infected individuals can be described by the SIR model, as illustrated in Fig. 1a, where the individuals are placed into three compartments: susceptible (S), infected (I), and removed (R).Susceptible individuals (S) have not contracted the disease and are at risk of being infected.Infected individuals (I) have contracted the disease and can infect others.Removed individuals (R) have recovered or died from the disease.There are two dynamical processes: (1) infection during which susceptible individuals become infected by others and transition to the I state so as to become capable of infecting others, as shown in Fig. 1a(i-iii); and (2) removal during which infected individuals recover or die from the disease transmission and transition to the R state, as shown in Fig. 1a(iv).The ability to infect others of an infected individual can be characterized by the infection time distribution, ψ inf (τ ), where τ denotes the time elapsed between the time the individual is infected and the current time, and the probability of the infection process occurring during the time interval [τ, τ + dτ ) is given by ψ inf (τ )dτ , as shown in Fig. 1b(i).Likewise, the removal process is described by the removal time distribution, ψ rem (τ ), where the probability of a removal occurring within the time interval [τ, τ + dτ ) is given by ψ rem (τ )dτ , as shown in Fig. 1b(ii).The time distributions of the infection and removal processes with memory effects are general, with the exponential distributions associated with the memoryless process being a special case of the memory-dependent process.
The generic memory-dependent SIR spreading dynamics can be described by a set of deterministic integro-differential equations: where s l (t), i l (t), and r l (t), respectively, denote the fractions of the susceptible, infected, and removed individuals in age group l.The term c l (t) = 1 − s l (t) = i l (t) + r l (t) represents the fraction of cumulative infections (including both infections and removals) with respect to the total population in age group l, while k is a parameter to adjust the overall contacts and n is the total number of age groups.The quantity ω inf (τ ) represents the hazard function of ψ inf (τ ), meaning the rate at which infection happens at τ , given that the infection has not occurred before τ .Ψ rem (τ ) is the survival function of ψ rem (τ ), meaning probability that the removal has not occurred by τ (see Method for detailed calculations for hazard and survival functions).When the infection and removal time distributions are known, Eqs.(1-3) provide an accurate description of generic (including memory-dependent and memoryless) SIR spreading Monte Carlo simulations in an age-stratified population-based system.As shown in Figs.2a-c, the theory is validated by the agreement between the numerical solutions of Eqs.(1-3) and the results from direct Monte Carlo simulations.(See Method for a detailed description of the Monte Carlo simulation procedure.)Eqs.(1-3) provide a general framework encompassing both non-Markovian and Markovian descriptions.If the infection and removal time distributions are exponential: ψ inf (τ ) = γe −γτ and ψ rem (τ ) = µe −µτ , Eqs. (1-3) can be reformulated into a Markovian theory and further simplified into a set of ordinary differential equations with the constant infection and removal rates γ and µ (A detailed derivation of these equations is presented in Supplementary Note 1) : While the non-Markovian and Markovian theories [Eqs.(1)(2)(3) and Eqs.(4)(5)(6), respectively] describe the memorydependent and memoryless Monte Carlo simulations, we focus on whether the Markovian theory can accurately capture memory-dependent dynamics and how memory effects influence its accuracy.For this purpose, we seek to establish the equivalence between Markovian and non-Markovian approaches for describing spreading dynamics.

Dynamical equivalence
Eqs. (1)(2)(3)(4)(5)(6) provide a base to study the steady-state equivalence and transient-state equivalence between non-Markovian and Markovian theories, where a steady state characterizes the long-term dynamics of disease spreading and a transient state is referred to as the short-term behavior prior to system's having reached the steady state.As illustrated in Fig. 1c, note that steady-state equivalence means that the two types of spreading dynamics attain identical steady states [21,24,25], whereas transient-state equivalence implies that two types of dynamics are consistent throughout the entire transmission period.Transient-state equivalence thus implies steady-state equivalence, but not vice versa.Since Eqs.(1-6) also provides a numerical framework for Monte Carlo simulations, the terms "steady-state equivalence" and "transient-state equivalence" not only describe the connection between non-Markovian and Markovian theories, but also illustrate the relationship between memory-dependent and memoryless processes.
Therefore, the equivalence between the two theories implies the equivalence between the two corresponding processes, and vice versa.

Steady-state equivalence
Eqs. (1)(2)(3) give the following transcendental equation for determining the steady state (see Supplementary None 2 for detailed derivation): where sl = lim t→+∞ s l (t) and rl = lim t→+∞ r l (t) denote the fractions of the susceptible and removed individuals in age group l at the steady state (note that sl = 1 − rl , because at steady state, no infection exists), while śl = s l (0) and ŕl = r l (0) represent the initial fractions of the susceptible and removed individuals in this age group.For non-Markovian dynamics, basic reproduction number R 0 can be determined by: For Markovian dynamics, R 0 is given by where Λ max is the maximum eigenvalue of the matrix kA • p, and • denotes a row-wise Hadamard product between a matrix and a vector (see Supplementary Note 2 for a detailed description).Since Eq. ( 7) applies to both non-Markovian and Markovian dynamics, an identical R 0 value in the two cases will result in equivalent steady states from the same initial conditions.Consequently, for a given non-Markovian spreading process, there exists an infinite number of Markovian models with the same steady state, as the R 0 value is only determined by the ratio of γ to µ, but not by either value.
As shown in Fig. 2d, memory-dependent and memoryless spreading dynamics that reach the same steady state with the identical R 0 value confirm the steady-state equivalence.Fig. 2e demonstrates that, even for R 0 ranging from 0.01 to 2, the equivalent memory-dependent and memoryless spreading dynamics still produce highly consistent steady states that can be calculated from Eq. ( 7), which share the same critical point of phase transition at R 0 = 1.

Transient-state equivalence
From the preceding section, we used the basic reproduction number R 0 , a fundamental metric quantifying the number of secondary infections generated by a single individual, to characterize the steady-state equivalence.Here, we propose to quantify the transient-state equivalence through the average generation time T gen that measures the "velocity" at which secondary infections occur.This time can be calculated as [3,58], is the generation time distribution.Effectively, T gen measures the average duration of disease transmission from an infected individual to the next generation of individuals.Likewise, the average infection time T inf and the average The solid brown, blue, and green curves represent the theoretical results of the susceptible, infected, and removed fractions, while the solid orange, red, and purple curves show the corresponding results of 100 independent Monte Carlo simulations.d: The red and blue curves, respectively, depict the removed fractions from the memory-dependent and memoryless Monte Carlo simulations of 100 independent realizations with steady-state equivalence.e: Red + and blue × markers, respectively, represent the steady-state removed fractions of memory-dependent and memoryless Monte Carlo simulations for different values of R0, where each marker is the result of averaging 100 independent simulations.The orange curve is the numerical calculations from Eq. ( 7), and the vertical dashed line denotes the critical point R0 = 1.f-g: For Tgen = Trem in the non-Markovian theory: the blue and green curves in f denote the susceptible and removed fractions, while the black × markers represent the inferred susceptible fractions calculated by substituting removed fractions in Eq. ( 11), which agrees with the susceptible curve calculated from Eq. ( 1).The red and blue curves in g denote the non-Markovian susceptible and removed fractions, while the orange and purple dashed curves are the corresponding curves of the Markovian transmission obtained from Eqs. (13)(14), which agree with the non-Markovian results.(The Euler-Lotka equation assumes exponential growth of a disease outbreak during the initial stage.As a result, the Markovian curves in g slightly deviate from the non-Markovian ones as the cumulative infections increase.)h-i: For Tgen = Trem in the non-Markovian theory, the inferred susceptible curve does not match the numerical result, and the susceptible and infected curves of the Markovian transmission obtained from Eqs. (13)(14) do not match the corresponding non-Markovian results.j: Five scenarios for the non-Markovian time-distribution setting (within each scenario, T inf and Trem are fixed): Weibull, T inf = 5, Trem = 7 (blue +); Weibull, T inf = 5, Trem = 5 (red ×); Weibull, T inf = 7, Trem = 7 (purple 2); log-normal, T inf = 5, Trem = 7 (green ); gamma, T inf = 5, Trem = 7 (orange 3).The value of Tgen is modified to adjust ln η for better visualization.
removal time T rem are defined as the mean values of the infection and removal time distributions: In calculating T gen , the individual's removal is taken into account, while T inf measures the average time of the first disease transmission of an infectious individual without factoring in removal.In the classical memoryless transmission with exponential distributions ψ inf (τ ) and ψ rem (τ ), the equality T gen = T rem holds.However, for memory-dependent spreading, the possible scenarios are: T gen = T rem , T gen < T rem , or T gen > T rem .Specifically, because T gen , T inf and T rem all represent the mean values of distributions, it is possible for T rem to be shorter than T gen or T inf in some situations.And our web-based application demonstrates the impact of parameters on the time distributions (infection, removal, and generation) as well as their average times [55].
For a non-Markovian spreading process, if the equality T gen = T rem holds, in the transient state we have which exhibits a memoryless transmission pattern similar to Markovian dynamics (see Supplementary Notes 3 for a detailed analysis).Intuitively, the equality T gen = T rem signifies that the infection and removal processes occur concurrently, which in turn leads to a memoryless relationship between the two processes, thereby minimizing the memory effects.Furthermore, to determine the corresponding Markovian parameters γ and µ of the Markovian transmission which is equivalent to the non-Markovian dynamics in the transient state, we need to utilize the Euler-Lotka equation [50,58,59]: where g denotes the growth rate of the non-Markovian dynamics and is another measure of how quickly the epidemic is spreading within a population.Therefore, we can calculate the values of the basic reproduction number, R 0 , and growth rate, g, in the non-Markovian dynamic by using Eqs.( 8), ( 10) and (12).Additionally, the Markovian form of ψ gen (τ ) according to Eq. ( 10) is µe −µτ , and the equivalent Markovian and non-Markovian dynamics in the transient state have the same values of R 0 and the equal values of g.By substituting ψ gen (τ ) = µe −µτ and the calculated R 0 and g into Eq.( 12), we can determine the value of µ.Furthermore, using Eq. ( 9), we can find the value of γ based on µ.Hence, the Markovian parameters γ and µ are determined as follows: And we provide visualizations that illustrate how the values of γ and µ are influenced by the distribution parameters in our web-based application [55].
As illustrated in Figs.2f-g, when the equality T gen = T rem holds for the non-Markovian dynamics, Eq. ( 11) holds, which can be seen by comparing the susceptible curve calculated from Eq. (1) to that inferred from Eq. ( 11), as shown in Fig. 2f.In this case, the Markovian spreading curves deduced from Eqs. (13)(14) closely align with the non-Markovian transient curves, as shown in Fig. 2g.However, as shown in Figs.2h-i, if the equality does not hold, the equivalence in transient states breaks down.It is important to note that the Euler-Lotka equation assumes an exponential growth of a disease outbreak and is only reasonable at the initial stage.Consequently, as the cumulative infections increase (Fig. 2g), the Markovian curves will exhibit slight deviations from the non-Markovian counterparts.Meanwhile, because the equivalent dynamics share the same R 0 , they will ultimately reach the same steady state, ensuring that deviations will diminish while they approach the steady state.
To evaluate, under different values of the generation-to-removal time ratio η ≡ T gen /T rem for non-Markovian dynamics we introduce a metric, ε, to quantify the difference from the corresponding Markovian results calculated from Eqs. (13)(14) (see Method for detailed definition of ε).Fig. 2j shows, for non-Markovian numerical calculations, five scenarios under various forms of time distributions ψ inf (τ ) and ψ rem (τ ) constrained by certain average infection and removal times: Weibull, T inf = 5, T rem = 7; Weibull, T inf = 7, T rem = 7; Weibull, T inf = 5, T rem = 5; log-normal, T inf = 5, T rem = 7; and gamma, T inf = 5, T rem = 7 (see Method for detailed definitions of Weibull, log-normal, and gamma distributions).The T gen value is adjusted to obtain different values of ln η.For ln η = 0, i.e., T gen = T rem , the "distance" ε between the transient states of non-Markovian and Markovian dynamics with parameters determined from Eqs. (13)(14) is minimal.Otherwise, ε increases as ln η deviates from zero.Remarkably, Fig. 2j shows that ε depends only on the ratio of T gen to T rem , not on the values of T gen , T inf , or T rem , and it is not affected by the specific form of the time distributions.
Furthermore, it is important to note that the condition T gen = T rem can guarantee transient-state equivalence between a non-Markovian dynamic and a Markovian one, but according to Eqs. (13)(14), it does not imply that the average generation and removal times of the non-Markovian dynamic must be equal to those of the equivalent Markovian one.For instance, if a non-Markovian dynamic satisfies the condition of transient-state equivalence and we keep its average generation and removal times fixed, altering the shape of the corresponding time distributions will change the transmission speed [58].This change, in turn, affects the infection and removal rates of the equivalent Markovian dynamic, leading to different average generation and removal times for the Markovian equivalent dynamic (see Supplementary 4 and 5 for a detailed analysis).

Markovian approximation of memory-dependent spreading dynamics
As illustrated in Fig. 1d, testing the applicability of Markovian theory for memory-dependent spreading dynamics requires three steps.The first step is fitting, where the memory-dependent Monte Carlo simulation data are divided into two parts: (a) a short initial period used as the training data for fitting the Markovian parameters in Eqs.(4)(5)(6), and (b) the remaining testing data for evaluating the performance of the Markovian model (see Methods for details of the fitting procedure).The second step is to use the fitted Markovian model for tasks such as estimating R 0 , predicting outbreaks and assessing the prevention effects of different vaccination strategies.The third step is testing, i.e., evaluating the accuracy of the Markovian model, e.g., by comparing the estimated and actual R 0 values, disease outbreaks and prevention effects.As real-world disease spreading is subject to environmental, social and political disturbing factors, for the fitting and testing steps, we conduct Monte Carlo simulations of stochastic memory-dependent disease outbreaks to generate the training and testing data.
Here, we first analyze the influence of η on the estimation of R 0 using the Markovian theory, and use the results to design two tasks to evaluate the applicability of the theory in epidemic forecasting and prevention evaluation of memory-dependent spreading.For comparison, we also generate the corresponding results from the non-Markovian theory in the two tasks.

Estimation of basic reproduction number
Estimating basic reproduction number R 0 is crucial for determining the ultimate prevalence of disease spreading and for assessing the effectiveness of various disease containment measures [50,58,59].When using the Markovian theory to fit the early-stage transmission of a memory-dependent process, a key parameter that can affect the estimation of R 0 is the ratio η.To develop an analysis, recall the basic principle for estimating R 0 : disease spreading dynamics can be viewed as a combination of two parallel processes: infection and removal.In particular, the infection process is the reproduction of the disease within each generation, where each infected individual generates an average of R 0 newly infected individuals in the subsequent generation after a mean time period T gen .In the removal process, infected individuals are removed from the spreading chain, where each generation takes an average time T rem to be removed.For a Markovian type of dynamics with constant γ and µ, the equality T gen = T rem holds.Consequently, during the Markovian fitting step, the average number of new infections upon the removal of a single infected individual is taken as the value of R 0 .For memory-dependent spreading, if the equality T gen = T rem holds, the memory-dependent spreading curves will possess an approximate memoryless feature so that R 0 can be still be estimated by counting the number of new infections at the time when the current generation of infections is removed, as shown in Fig. 3a.However, for T gen < T rem , more than one generation is produced while the current generation is removed, R 0 estimated by the Markovian theory will represent an overestimate, as shown in Fig. 3b.For T gen > T rem , less than one generation The relationship between Tgen and Trem influences the number of new infections when the current generation of infections is removed.For Tgen = Trem, the number of new infections is exactly R0.For Tgen < Trem, the number of new infections is greater than R0.For Tgen > Trem, the number of new infections is smaller than R0.d: Three distinct categories of disease spreading with the same value of R0: Tgen < Trem (blue curves), Tgen = Trem (red curves), and Tgen > Trem) (green curves), where the fractions of cumulative infected individuals (i.e., sum of infected and removed fractions) are calculated using 100 independent realizations.The predicted future evolution of the spreading dynamics by the Markovian theory with the fitted parameters are also shown: Tgen < Trem (blue + symbols), Tgen = Trem (red × symbols), and Tgen > Trem (green symbols), where the gray area marks the average cumulative infected fraction for selecting the training data (see Method for details).e: The relationship between ln (ln R0/ ln R0) (R0 represents the real basic reproduction number, while R0 denoted the estimated one) and ln η where the horizontal and vertical dotted lines show that the equality between Tgen and Trem results in an accurate estimation of R0 and the dashed line represents a linear fitting with the slope 1.55.Inset: the relation between ln R0 and η with the asymptotic behaviors: for η → 0, R0 → +∞ (ln R0 → +∞), and for η → +∞, R0 → 1 (ln R0 → 1).
is created during T rem , the Markovian theory will give an underestimate of R 0 , as shown in Fig. 3c.
Fig. 3d shows fitting curves (training data) from the early stage of memory-dependent spreading simulations that have identical R 0 values for T gen = T rem (red curves), T gen < T rem (blue curves), and T gen > T rem (green curves) from the Markovian theory.When the equality T gen = T rem holds, the Markovian theory with fitted parameters generates accurate the future evolution (red × symbols).For T gen < T rem , the outbreak in the initial stage is accelerated, resulting in an overestimation by the Markovian theory (blue + symbols).For T gen < T rem , the initial outbreaks are decelerated, leading to an underestimation by the Markovian approach (blue + symbols).
The above qualitative insights lead to a semi-empirical relationship between the Markovian-estimated basic reproduction number R0 and its actual value R 0 as: where a is a positive coefficient (see Method for a detailed derivation).The value of a is a crucial and constant parameter in Eq. ( 15), and it needs to be determined by fitting it to the data.Once this constant a is obtained, the actual value of R 0 can be derived by adjusting the estimated R0 based on Eq. ( 15), and more accurate steady state can be calculated by using Eq.(7).
Eq. ( 15) implies the relationship ln (ln R 0 / ln R0 ) = a ln η.We use the five scenarios specified in Fig. 2j for memorydependent Monte Carlo simulations.Fig. 3e shows the linear relationship between ln (ln R 0 / ln R0 ) and ln η, providing support for our qualitative analysis of the Markovian estimation.The estimation of R 0 also depends on the ratio η and is relatively insensitive to the particular forms of the time distributions or the specific values of T gen , T inf , or T rem .The results in the inset of Fig. 3e further confirm that the estimated R0 approaches one when T gen is much larger than T rem and tends to +∞ when T gen is much smaller than T rem .By fitting the available data, we have determined the value of a to be 1.55.After obtaining the value of a, we can develop our web-based application for rectifying R 0 and epidemic forecasting [55].

Epidemic forecasting
As suggested in Fig. 1d, we evaluate the efficacy of Markovian theory for epidemic forecasting.We use the initial period of Monte Carlo simulation data to fit parameters under both Markovian and non-Markovian hypotheses and then to predict future disease outbreaks.The remaining simulation data are leveraged to evaluate the accuracy of the Markovian and non-Markovian forecasting results.Regardless of the type of time distributions in the memorydependent Monte Carlo simulations (Weibull, log-normal, or gamma), the non-Markovian model fits the training data in a consistent manner, i.e., by selecting Weibull time distribution.
Figs. 4a-c show the evolution of the spreading dynamics from three types of memory-dependent Monte Carlo simulations with Weibull infection and removal distributions, where the shape parameters α inf and α rem are selected according to ln α inf = −0.3,ln α rem = 1.2 (Fig. 4a), ln α inf = 0.45, ln α rem = 0.45 (Fig. 4b), and ln α inf = 1.2, ln α rem = −0.3 (Fig. 4c), for T inf = 5 and T rem = 7.For the Weibull distributions, we have α inf < α rem , α inf = α rem and α inf > α rem , corresponding to T gen < T rem , T gen = T rem , and T gen > T rem , respectively.We compare the simulated cumulative infected fractions to those predicted by the Markovian and non-Markovian theories.In general, the non-Markovian theory provides more accurate predictions than the Markovian theory.For the specific parameter setting ln α inf = 0.45, ln α rem = 0.45 (i.e., T gen = T rem ), both theories yield a high accuracy.
The accuracy can be assessed through the forecasting error ε + that evaluates whether a theory overestimates or underestimates the steady-state cumulative infection, i.e., quantifying the extent of deviation between the results obtained from Markovian or non-Markovian theories and those derived from Monte Carlo simulations (see Method for detailed definition of ε + ).A plus value of ε + means overestimation while minus value indicates underestimation.We evaluate the accuracy measure ε + in the parameter plane of ln α inf and ln α rem , ranging from −0.3 to 1.2.Figs.4de show that the Markovian accuracy is sensitive to parameter changes: underestimated if α inf is greater than α rem (T gen > T rem ), overestimated when α inf is smaller than α rem (T gen < T rem ), and a high forecasting accuracy is achieved only for α inf = α rem (T gen = T rem ).In contrast, the non-Markovian theory yields highly accurate results in the whole parameter plane, with only a slight underestimation for α inf α rem (T gen T rem , this is primarily due to the increased difficulty in fitting simulation data, as the simulation parameters become increasingly unreasonable.).
Using the five scenarios specified in Fig. 2j for memory-dependent Monte Carlo simulations, we obtain the relationship between ε + and ln η, as shown in Fig. 4f.It can be seen that, in the Markovian framework, an overestimation arises for T gen < T rem , and an underestimation occurs for T gen > T rem .Only when T gen = T rem is an accurate estimate achieved.In general, the non-Markovian theory provides much more accurate forecasting than the Markovian theory, especially when T gen and T rem are not equal.The results further illustrate that the forms of time distributions or the specific values of T gen , T inf , or T rem have little impact on forecasting accuracy.
To establish the relevance of these results to real-world diseases, we obtain the ψ inf (τ ) and ψ rem (τ ) relations for four known infectious diseases, including COVID-19, SARS, H1N1 influenza, and smallpox, using the information in Refs.[3,[40][41][42][43][44][45][46][47].We then calculate the corresponding values of ε + and ln η based on the Markovian and non-Markovian approaches.As demonstrated in Fig. 4f, the positions of the four diseases in the (ln η, ε + ) plane are consistent with the results of our estimations.Because the data were from the reports of laboratory-confirmed cases incorporating the effects of the quarantine and distancing from susceptible individuals after the confirmation of the diagnosis, T gen of the four diseases are all smaller than the corresponding values of T rem , leading to some overestimation for the Markovian forecasting results.

Evaluation of vaccination strategies
In the development and application of a theory for disease spreading, assessing the effects of different vaccination strategies is an important task.Here we consider five prioritization strategies for vaccine distribution [6]: individuals under 20 years (denoted as m = 1), adults between 20 and 49 years (m = 2), adults above 20 years (m = 3), adults above 60 years (m = 4), and all age groups (m = 5), and implement these strategies in Monte Carlo simulations.Fig. 5a shows the results of epidemic evolution in comparison with those without any vaccination intervention (m = 0), where the shape parameters are chosen according to ln α inf = 0.45 and ln α rem = 0.45 (T gen = T rem ).Figs. 5b-c show the results from the Markovian and non-Markovian theories, respectively, with the corresponding fitted parameters for the vaccination strategies in comparison with those without vaccination (see Method for the detailed procedure of vaccination in the theoretical calculation).These results indicate that the Markovian and non-Markovian theories yield the correct epidemic evolution and future outbreaks under different vaccination scenarios.
To characterize the effectiveness of different vaccination strategies in blocking disease transmission, we introduce a vector, δ, whose m-th element quantifies the cumulative infected fraction with the m-th vaccination strategy in the steady state: δ m = cm , for m = 0, . . ., 5. Fig. 5d shows that the δ vectors from the Monte Carlo simulation, Markovian and non-Markovian theories from Figs. 5a-c, respectively.
We further introduce a metric, the so-called prevention evaluation error ε * , that gauges the ability of the Markovian and non-Markovian theories to estimate the total effectiveness of vaccination, i.e., measuring the disparity between the results calculated by the Markovian or non-Markovian theories and those obtained through Monte Carlo simulations considering various vaccination strategies (see Method for the detailed definition of ε * ).Figs. 5e-f show the average values of ε * of the two theories in the simulation parameter plane using 100 independent realizations, which are similar to those in Figs.4d-e, indicating that the error mainly comes from the R 0 estimation.In general, the Markovian theory performs well only in the diagonal area of the parameter plane where α inf = α rem , as shown in Fig. 5e, and the non-Markovian theory outperforms the Markovian counterpart in most cases, as shown in Fig. 5f.
Meanwhile, we assess the ability of both the Markovian and non-Markovian theories to detect the optimal vaccination strategy.We also define a quantity, optimization failure probability ε, to quantify the probability of a theory failing to identify the optimal strategy, i.e., that leads to the lowest cumulative infection among the five strategies (see Method for the detailed definition of ε).Figs.5g-h illustrate the results of ε for the two theories within the parameter plane (ln α inf , ln α rem ).While the non-Markovian theory still demonstrates superior performance, the Markovian approach proves capable of identifying the optimal strategy across a significantly larger parameter space compared to the Markovian results depicted in Figure 5e.
We obtain the relationships between ε * and ln η, as well as between ε and ln η, as shown in Fig. 5i-j with the same five time-distribution scenarios as in Fig. 2j.In all cases of Fig. 5i, ε * reaches a minimum for ln η = 0 and increases as ln η deviates from zero.The agreement of the results from the five scenarios further illustrates that the forms of time distributions or the specific values of T gen , T inf , or T rem play little role in the errors in vaccination evaluation.Fig. 5i also includes the values of ε * for the real-world infectious diseases COVID-19, SARS, H1N1 influenza, and smallpox, which are consistent with those from the non-Markovian and Markovian theories.Regarding the results depicted in Fig. 5j, it is observed that the non-Markovian theories consistently outperform the Markovian counterparts.On the other hand, within a wide range of ln η values around 0, the Markovian theories successfully identify the optimal vaccination strategy among various commonly employed ones.When the value of ln η significantly deviates from 0, Markovian theories become ineffective in determining the optimal strategy.(Note that on the left side of Fig. 5j, we only present the failures of Markovian theories to identify the optimal strategy in the Monte Carlo simulations with log-normal distribution.This is primarily due to the fact that the parameters associated with the Weibull and gamma distributions fall outside the acceptable range when we keep T inf and T rem fixed to modify ln η to a very low value.)Furthermore, we demonstrate that even when employing Markovian approaches, the optimal vaccination strategy can still be determined among the five strategies considered for the four distinct real diseases.
Conducting accurate evaluations in prevention serves as the sufficient condition of the successful identification of the optimal strategy.In comparison to the prevention evaluation errors ε * of Markovian theories, the optimization failure probability ε exhibits a wider range of ln η values that result in the lowest value.The lack of mathematical continuity among the five strategies is the primary reason for this.It indicates that there is no smooth transition or mathematical relationship connecting these strategies, resulting in the rank of the strategies not changing promptly when the value of ln η deviates from 0. Therefore, only significant errors from the Markovian theories can result in the failure to detect the optimal strategy.Based on this analysis, the extent to which ln η deviates from 0, leading to the failure of Markovian theories, as well as whether such failure will occur, depends on the selection of the tested strategies.

DISCUSSION
The COVID-19 pandemic has emphasized the importance of investigating disease transmission in human society through modeling.Empirical observations have consistently demonstrated strong memory effects in real-world transmission phenomena.The initial transient stage of an epidemic is critical for data collection, prediction, and articulation of control strategies, but an accurate non-Markovian model presents difficulties.In contrast, a Markovian model offers great advantages in parameter estimation, computation, and analyses.Uncovering the conditions under which Markovian modeling is suitable for transient epidemic dynamics is necessary.
We have developed a comprehensive mathematical framework for both Markovian and non-Markovian compartmentalized SIR disease transmissions in an age-stratified population, which allows us to identify two types of equivalence between Markovian and non-Markovian dynamics: in the steady state and transient phase of the epidemic.Our theoretical analysis reveals that, in the steady state, non-Markovian (memory-dependent) transmissions are always equivalent to the Markovian (memoryless) dynamics.However, transient-state equivalence is approximate and holds when the average generation and removal times match each other.In particular, when the average generation time is approximately equal to the average removal time, the disease transmission and removal of an infected individual exhibit a memoryless correlation, thereby minimizing the memory effects of the dynamical process.This results in highly accurate results from the Markovian theory that captures the characteristics of memory-dependent transmission based solely on the early epidemic curves.Our analysis also suggests that the Markovian accuracy is mainly determined by the value of generation-to-removal time ratio in disease transmission, where a larger-than-one (smaller-than-one) ratio can lead to underestimation (overestimation) of the basic reproduction number and epidemic forecasting, as well as the errors in the evaluation of control or prevention measures.The estimation accuracy primarily depends on this ratio, but is not significantly affected by the specific values and distribution forms of the various times associated with the epidemic.This property exhibits substantial practical importance, because the average generation and removal times can be readily assessed based on sparse data collected from the transient phase of the epidemic, but to estimate their distributions with only sparse data is infeasible.These results provide deeper quantitative insights into the influence of memory effects on epidemic transmissions, leading to a better understanding of the connection and interplay between Markovian and non-Markovian dynamics.
There were previous studies of the equivalence between Markovian and non-Markovian transmission in the SIS model [21,24,25].However, these studies addressed the steady-state equivalence rather than the transient-state equivalence.To our knowledge, our work is the first to investigate the transient-state equivalence of the SIR model.In addition, previous studies mainly examined the impact of the average generation time on the transmission dynamics, such as how the shape of the generation time distribution affects the estimation of R 0 [58] or the use of serial time distributions in estimating R 0 during an epidemic [50].There was a gap in the literature regarding how generation times affect the accuracy of different models.Our paper fills this gap by providing a criterion for using Markovian frameworks to model memory-dependent transmission based on the relationship between the average generation and removal times.
From an application perspective, our study suggests that the impact of the time distribution forms on Markovian estimation accuracy is minimal, making it easier to select models between Markovian and non-Markovian dynamics in the initial outbreak of an epidemic based only on the generation-to-removal time ratio.This insight is especially useful since detailed time distribution forms are often harder to detect than their corresponding mean values.In addition, we note that in previous studies, it was observed that in various scenarios, serial intervals, albeit with larger variances, are anticipated to possess a consistent mean value with the average generation time and are more straightforward to measure [3,[50][51][52][53][54].Given the practical difficulties in observing the generation time, our finding of minimal impact from the distribution forms suggests that the average serial interval can be utilized as a substitute of the average generation time to determine the applicability of the Markovian theories for modeling purposes without compromising accuracy, although numerous studies have indicated that replacing the generation time distribution with the serial interval distribution may affect the analysis of transmission dynamics [49,50,60].Meanwhile, based on the Eq. ( 15), if we determine the ratio of generation-to-removal time, the estimated R 0 obtained through the Markovian approach can be adjusted to approximate the true.And our web-based application showcases the demonstration of rectifying R 0 and epidemic forecasting [55].
Our study highlights the critical importance of accurately quantifying R 0 for achieving precise epidemic forecasting and prevention evaluation.A previous work [59] revealed that the value of R 0 depends on three key components: the duration of the infectious period (e.g., ψ rem (τ )), the probability of infection resulting from a single contact between an infected individual and a susceptible one (e.g., ψ inf (τ )), and the number of new susceptible individuals contacted per unit of time.However, given the practical limitations inherent in obtaining all three components, numerous methods have been developed for estimating R 0 .Although our work presents a specific approach, which fits the parameters of exponential or non-exponential time distribution by using the initial outbreak curves, it is not the only one available.For example, when contact patterns are unknown, R 0 can be estimated by fitting the growth rate g and the generation time distribution ψ gen (τ ), and then applying them in the Euler-Lotka equation [50,58,59].However, since the focus of our work is on epidemic forecasting and evaluation of prevention measure, R 0 can be directly calculated once ψ inf (τ ) and ψ rem (τ ) are fitted, without requiring the fitting of any additional quantity.The estimation of R 0 can also be achieved by using data in the steady state, such as the final size of an epidemic or equilibrium conditions [59].However, this method is not suitable for the transient phase where only early-stage curves are available.Utilizing the approach delineated in this paper is practically more appropriate for estimating R 0 .
While our study focused on transmission within the SIR framework, extension to SEIR or SIS models is feasible.While we emphasized the significance of the transient-state equivalence in disease transmission, transient dynamics are more relevant or even more crucial than the steady state in nonlinear dynamical systems [61].For example, in ecological systems, transient dynamics play a vital role in empirical observations and are therefore a key force driving natural evolution [62][63][64][65][66][67].In neural dynamics, transient changes in neural activity can mediate synaptic plasticity, a crucial mechanism for learning and memory [68][69][70].Therefore, the identification of suitable conditions for choosing between Markovian and non-Markovian dynamics may not be limited to transmission dynamics alone and may serve as a valuable reference for other fields as well.
Taken together, our study establishes an approximate equivalence between Markovian and non-Markovian dynamics in the transient state, assuming that time distributions follow Weibull forms (see Supplementary Note 3 for details).While the applicability of our findings to most synthetic and empirical distributions has been analyzed qualitatively, a quantitative analysis requires further studies.For extreme cases with non-Weibull distributions, the transmission should be evaluated using other specific methods.While we have provided a qualitative analysis of the mechanism underlying why time distribution forms have minimal impacts on the errors of Markovian estimations, a more rigorous theoretical analysis is needed and requires further exploration.In addition, due to the complexity of the nonlinear transmission, our study has produced a semi-empirical relationship to estimate the overestimation and underestimation of Markovian methods.Further research is required to develop a rigorous formula that can accurately predict these effects.

Monte Carlo simulation
In the simulation, we classify N individuals into n subgroups based on the age distribution p.The index of the subgroup to which an individual belongs is denoted by l (where 0 ≤ l ≤ n − 1), and the index of the individual within the subgroup is denoted by u (where 0 ≤ u ≤ p l N ).The state of the u-th individual in the l-th subgroup is represented by X lu , which includes the states S (susceptible), I (infected), W (recovered), and D (dead), where W and D both represent R (removed).For each individual, we also record the absolute time of infection and removal using two variables: T inf lu and T rem lu , respectively.The absolute time of the system is denoted by T , and we implement the total spreading simulation step by step using a finite time step ∆T as follows: i) Initialization: set T = 0, X lu for every individual is set to S. ii) Set infection seeds: choose a set of individuals as the infection seeds and the corresponding X lu are set to I, the corresponding T inf lu are set to 0, and T rem lu are set to a random value following the removal time distribution ψ rem (τ ).iii) Infection of one step: calculate the infection rate, ωinf lu (T ), of infected individual u in age group l during the current time step by The probability ωinf l (T ) of each susceptible individual in age group l being infected can be calculated by where I m is the index set of the infected individual in age group m.The number of the susceptible individuals being infected in a age group follows a binomial distribution B(s l (T )p l N, ωinf l (T )), where s l (T ) denotes the fraction of susceptible individuals in age group l at time T .Then generate a random number N l (T ) following this binomial distribution and set the N l (T ) individuals as I state.The corresponding T inf lu of the new infected individuals are set to the current T and T rem lu are set to T rem + T , where the random T rem follow the removal time distribution ψ rem (τ ).iv) Check if T rem lu of each infected individual is during the current time step.If this condition is satisfied, set their state to D with a probability of death and to W otherwise.Then let T = T + ∆T .v) Repeat the process iii) and iv), until no individual with I index exists.

Time distributions
In the Monte Carlo simulations, we employ three types of time distributions, i.e., Weibull, log-normal, and gamma, to describe the memory-dependent transmission process.
For Weibull time distribution, it follows: , where α and β denote the shape and scale parameters, respectively.The log-normal time distribution is defined as follows: ).
The gamma time distribution is expressed as follows: where Γ(•) denotes gamma function, while α and β represent the shape and scale parameters, respectively.
For each type of time distribution, denoted as ψ(τ ), it can be either ψ inf (τ ) or ψ rem (τ ).Similarly, the parameter α can take either α inf or α rem , and the parameter β can be either β inf or β rem .

Derivation of semi-empirical estimation of basic reproduction number
Intuitively, the period of T rem can accommodate T rem /T gen = 1/η time intervals of length T gen , corresponding to the result of 1/η generations of infections.This can lead to an exponential increase in the number of infections during T rem .This intuition suggests a relationship between the fitted basic reproduction number R0 and the actual R 0 , which can be expressed as an exponential function: , where f (•) is a monotonically increasing function that satisfies three conditions.First, f (1) = 1, indicating that R0 can be accurately estimated when T gen = T rem .Second, f (0) = 0, meaning that if T rem is an extremely small fraction of T gen , the transmission will take a long time to reach the steady state, causing the curve to be flat in the initial stage and potentially causing the Markovian fitting to produce the estimate R0 = 1.Third, f (+∞) = +∞, implying that if T rem is extremely large compared to T gen , the transmission will quickly reach the final prevalence, causing the Markovian fitting to give an extremely large estimate of R0 .Because the actual transmission process involves many complicated nonlinear relationships, identifying the specific form of the function f (•) is a challenging task.We thus assume where a is an unknown positive coefficient.This leads to Eq. (15).

Definition of errors
The difference ε between non-Markovian and the corresponding Markovian results calculated from Eqs. (13)(14) is defined as: ε = x,x * ∈{(s,s * ),(i,i * ),(r,r * )} x * − x 2,T x 2,T where the pairs (s, s * ), (i, i * ) and (r, r * ) correspond to the non-Markovian and Markovian susceptible, infected and removed curves, respectively.The 2-norm • 2,T on time duration T ensures that ε measures the "distance" between non-Markovian and the Markovian transient states.It is not appropriate to set T as the total transmission period because the cumulative infected fraction approaches the final value asymptotically, making it difficult to determine the exact time point of the steady state.To address this issue, we choose T as [0, tθ ], where tθ is the time when the cumulative infected fraction c(t) reaches the θ percentile point within the range that spans from its initial value (ć) to its final value (c).The value of θ in Fig. 2j is selected as 50 (see Supplementary Note 6 for more selection of θ and detailed analysis).
The forecasting error ε + that evaluates whether a theory overestimates or underestimates the steady-state cumulative infected fraction is defined as: where t denotes the time when the stochastic simulation reaches the steady state when no infection occurs in the population, c( t) and ĉ( t) are the cumulative infected fractions from theory and simulation, respectively.A positive value of ε + indicates overestimation, whereas a negative value indicates underestimation.The prevention evaluation error ε * , that gauges the ability of the Markovian and non-Markovian theories to estimate the total effectiveness of vaccination, is defined as: where δ is the result from Monte Carlo simulation and • 2 is the 2-norm of a vector.The optimization failure probability ε, which measures the probability that a theory fails to identify the optimal vaccination strategy, is defined as: ε = z l=1 ξ (l)  z , where ξ (l) = 0 if argmin σ(l) = argmin σ (l) 1 otherwise , σ(l) and σ (l) represent the vectors, σ and σ, for the l-th experiment, respectively, and z denotes the total number of experiments (in this paper, z is set to 100).Consequently, ε quantifies the fraction of experiments in which a theory fails to identify the optimal vaccination strategy, and serves as a measure of the probability of failure in optimizing the vaccination strategy.

Fitting method
Because the removal process is independent of the infection one, we divide the fitting method into two parts: removal parameter fitting and infection parameter fitting.Specifically, we use c * init and r * init to denote the cumulative infected and removed fractions in the initial stage of a Monte Carlo simulation.These two types of data are substituted into the Eq. ( 3) to fit the parameters of ψ rem (τ ).Likewise, we use c * l,init to denotes the cumulative infected fraction of age group l in the initial stage of a Monte Carlo simulation and s * l,init = 1 − c * l,init represents the corresponding susceptible fraction.After obtaining the removal parameters, these two types of data are put into Eq.(1) to fit the infection time distribution parameters.
In our study, we selected the curves of all states that occurred prior to the time point at which the cumulative infected fraction reached a specific percentile (e.g., 20%) situated between the initial and steady cumulative infected fractions, as the training data.Choosing a specific time period as the training data may not be appropriate, as it can result in an overabundance of data points for fitting due to some instances of fast transmission already having reached the steady state, while some instances of slow transmission may not have spread out yet, leading to too few data points.

Vaccination method
We assume that the individuals will build enough immune protection from the disease κ days after vaccination with the probability ρ.In Monte Carlo simulations, the susceptible individual who gets vaccinated and the corresponding time T vac are associated with the probability ρ, and the absolute time becomes T vac + κ.If this individual has not been infected, he/she will be set to a state called protected state, indicating that this individual is protected from the disease.
When a fraction of individuals in age group l get vaccinated, the detailed vaccination fraction υ l , susceptible time t vac and the fraction of susceptible individuals s l (t vac ) are recorded.When the absolute time reaches T vac + κ, the corresponding value of s l (t vac + κ) will be set as s l (t vac + κ)ρ υ l s l (tvac) .

FIG. 1 :
FIG.1: Overall structure of this work.(a) SIR Model.Each individual belongs to one of the three states: susceptible (S), infected (I), or removed (R).When infected (i), a susceptible individual will switch into the I state (ii) and gain the ability to infect others (iii).An infected individual is removed (through recovery or death) with a probability (iv).(b) Non-Markovian versus Markovian process.The infection capacity of an infected individual is characterized by the infection time distribution ψ inf (τ ) and its removal can be described by the removal time distribution ψrem(τ ).For the non-Markovian process, the distributions can assume quite general forms, while the distributions are exponential for a Markovian process.(c) Equivalence between non-Markovian and Markovian processes: (i) steady-state equivalence holds under all conditions; (ii) transient-state equivalence only holds only when Tgen is equal to Trem.(d) Markovian estimation of memory-dependent process.(i) The initial phase of the Monte Carlo simulation is used to fit the parameters according to the Markovian theory.(ii) Important issues such as the estimation of R0, epidemic forecasting, and the evaluation of the vaccination strategies can be addressed by the theory.(iii) The remaining data generated by the Monte Carlo simulation is used to test the accuracy of the estimated R0, epidemic forecasting, and prevention evaluation.

j 7 FIG. 2 :
FIG. 2:Steady-state and transient-state equivalence between Markovian and non-Markovian dynamics.a-c: The solid brown, blue, and green curves represent the theoretical results of the susceptible, infected, and removed fractions, while the solid orange, red, and purple curves show the corresponding results of 100 independent Monte Carlo simulations.d: The red and blue curves, respectively, depict the removed fractions from the memory-dependent and memoryless Monte Carlo simulations of 100 independent realizations with steady-state equivalence.e: Red + and blue × markers, respectively, represent the steady-state removed fractions of memory-dependent and memoryless Monte Carlo simulations for different values of R0, where each marker is the result of averaging 100 independent simulations.The orange curve is the numerical calculations from Eq. (7), and the vertical dashed line denotes the critical point R0 = 1.f-g: For Tgen = Trem in the non-Markovian theory: the blue and green curves in f denote the susceptible and removed fractions, while the black × markers represent the inferred susceptible fractions calculated by substituting removed fractions in Eq. (11), which agrees with the susceptible curve calculated from Eq. (1).The red and blue curves in g denote the non-Markovian susceptible and removed fractions, while the orange and purple dashed curves are the corresponding curves of the Markovian transmission obtained from Eqs.(13)(14), which agree with the non-Markovian results.(The Euler-Lotka equation assumes exponential growth of a disease outbreak during the initial stage.As a result, the Markovian curves in g slightly deviate from the non-Markovian ones as the cumulative infections increase.)h-i: For Tgen = Trem in the non-Markovian theory, the inferred susceptible curve does not match the numerical result, and the susceptible and infected curves of the Markovian transmission obtained from Eqs.(13)(14) do not match the corresponding non-Markovian results.j: Five scenarios for the non-Markovian time-distribution setting (within each scenario, T inf and Trem are fixed): Weibull, T inf = 5, Trem = 7 (blue +); Weibull, T inf = 5, Trem = 5 (red ×); Weibull, T inf = 7, Trem = 7 (purple 2); log-normal, T inf = 5, Trem = 7 (green ); gamma, T inf = 5, Trem = 7 (orange 3).The value of Tgen is modified to adjust ln η for better visualization.

d 7 FIG. 3 :
FIG.3: Estimation of R0. a-c: Mechanism of the R0 estimation.The red arrows represent the infection process of the next generation by the current generation, while the dashed arrows denote the removal of the current generation.The relationship between Tgen and Trem influences the number of new infections when the current generation of infections is removed.For Tgen = Trem, the number of new infections is exactly R0.For Tgen < Trem, the number of new infections is greater than R0.For Tgen > Trem, the number of new infections is smaller than R0.d: Three distinct categories of disease spreading with the same value of R0: Tgen < Trem (blue curves), Tgen = Trem (red curves), and Tgen > Trem) (green curves), where the fractions of cumulative infected individuals (i.e., sum of infected and removed fractions) are calculated using 100 independent realizations.The predicted future evolution of the spreading dynamics by the Markovian theory with the fitted parameters are also shown: Tgen < Trem (blue + symbols), Tgen = Trem (red × symbols), and Tgen > Trem (green symbols), where the gray area marks the average cumulative infected fraction for selecting the training data (see Method for details).e: The relationship between ln (ln R0/ ln R0) (R0 represents the real basic reproduction number, while R0 denoted the estimated one) and ln η where the horizontal and vertical dotted lines show that the equality between Tgen and Trem results in an accurate estimation of R0 and the dashed line represents a linear fitting with the slope 1.55.Inset: the relation between ln R0 and η with the asymptotic behaviors: for η → 0, R0 → +∞ (ln R0 → +∞), and for η → +∞, R0 → 1 (ln R0 → 1).

FIG. 4 :
FIG. 4: Epidemic Forecasting.a-c: Predicted evolution of the cumulative infected fraction (i.e., sum of infected and removed fractions) by the Markovian (orange dotted curves) and non-Markovian (red dashed curves) theories, in comparison with the Monte Carlo simulations with Weibull time distributions (blue solid curves), for three sets of simulation parameters, respectively.The results are the averages of 100 independent realizations with the standard deviations indicated by the shaded regions.The gray area marks the average cumulative infected fraction for training data selection.d-e: The forecasting errors ε + of Markovian and non-Markovian theories with respect to the memory-dependent Monte Carlo simulations in the parameter plane of ln α inf and ln αrem in the range [−0.3, 1.2].The green, blue and red squares mark the parameters of Monte Carlo simulations in a-c, respectively.f : The forecasting errors ε + from the Markovian and non-Markovian theories for different values of ln η under five scenarios of time distribution setting for Monte Carlo simulations.The corresponding estimations for a number of real-world diseases (COVID-19, SARS, H1N1 and Smallpox) are also included.

FIG. 5 :
FIG. 5: Evaluation of vaccination strategies.a-c: For simulation parameters chosen according to ln α inf = 0.45 and ln αrem = 0.45, the cumulative infected fraction (i.e., sum of infected and removed fractions) curves from Monte Carlo simulations and the corresponding Markovian and non-Markovian theories with fitting parameters for five vaccination strategies and the case of no vaccination.The average results are obtained from 100 independent realizations with the shaded regions representing the standard deviations.The gray area marks the average cumulative infected fraction for training data selection.d: Vector δ calculated from the results in Figs.5a-c.e-f : The prevention evaluation errors ε * of Markovian and non-Markovian theories for evaluating the effects of vaccination prevention in the parameter plane (ln α inf , ln αrem).The green squares mark the selected parameters in a-c.g-h: The optimization failure probabilities ε arising from the Markovian and non-Markovian within the parameter plane (ln α inf , ln αrem).The green squares mark the selected parameters in a-c.i: The prevention evaluation errors ε * from the Markovian and non-Markovian theories versus ln η under five scenarios of time distribution setting for Monte Carlo simulations.The estimated errors for four real diseases (COVID-19, SARS, H1N1 and Smallpox) are also shown.j: The optimization failure probabilities ε from the Markovian and non-Markovian theories against ln η in five different time distribution scenarios for Monte Carlo simulations.The optimization failure probabilities for four real diseases (COVID-19, SARS, H1N1 and Smallpox) are also presented.