The current COVID-19 pandemic highlights the relevance of reliable theoretical modelings able to provide correct predictions of epidemy propagation and extension. In the recent past, a great effort has been devoted to trying to give correct representation to such phenomena as infection latency before symptoms onset, the time required for healing, the role of asymptomatic individuals in epidemy propagation [1,2] and so forth. Several authors [3-8] have attempted to improve the likelihood of deterministic modelings trying to increase the number of model compartments or introducing a delay between the time of pathogen exposition and that of inclusion in the infected compartment [8-14]. Other works [15-17] have investigated the ability to predict epidemic dynamics of distributed delay models, which represent delays by introducing differential equations convolutions with some functions that model class transitions. Unfortunately, these modelings (delayed and time distributed) do not guarantee that time elapsed between exposition and illness is correctly distributed because they act on compartment populations rather than on fluxes among them. Despite expectations, these models do not control the time interval between entering and leaving a given class (compartment) and how this time lag is distributed when considering many class transitions. Moreover, the delay is usually introduced to model the incidence rate, while all other passages between different compartments only rarely are delayed. Nonetheless, also these class transitions need a realistic time distribution for the theoretical model to be able to provide realistic predictions of epidemic propagation. Evidently, this lack of likelihood can reduce the accuracy of quantitative predictions forecasted by these models when investigating epidemic evolutions. This point, which will be discussed in the following, has also been addressed in  although in the framework of stochastic modeling. Many authors have investigated the effects on the dynamical behavior of deterministic models when including in their equations a nonlinear incidence or recovery rate [19-23], proving the occurrence of bifurcations . The nonlinear terms originated trying to include into the model the effects of systemic phenomena such as the availability of hospital beds on the recovery rate, or restrictive measures dictated by governments. The mathematical formulation of these nonlinearities is however arbitrary, and no way exists to prove the physical foundation of different models of the same phenomenon. Din et al [25,26] have proposed the application of fractional order Atanganaa–Baleanu Caputo derivative to the theoretical modeling of the Hepatitis B epidemic, also differentiating acute infected from chronically infected individuals. An emerging approach is the application of data assimilation techniques to the standard or a modified version of the SIR model . This option leads to a modified Kalman filter, which would be specifically suited to solve the problem of parameter estimation stemming from epidemic data.
In this paper, we outline a new type of deterministic epidemic modeling that include five population compartments and introduces the concept of fluxes among them. Classes are considered time-invariant linear systems that connect the input to output fluxes, making it possible to account for the time evolution of each population class. The model guarantees a parametric temporal distribution of infection evolution that is dictated by the impulse response function associated with each class transition. Due to its mathematical structure, the new model prevents impossible phenomena (class transitions) such as some population fraction that persists into the infected class for a too short (e.g.: null), or too long (e.g.: years), time interval. The model can allow restored people to come back into the susceptible class, giving rise to successive epidemic peaks, and predicts the existence of “evanescent” epidemic waves circulating for a long time even for basic reproduction numbers less than a unit. Due to its mathematical structure, the model holds thirteen free parameters that can fit its behavior to different kinds of diseases.
The paper shows the outcomes of simulations executed to assess the performance of the new model in comparison with the SIR, and for the current COVID-19 pandemic.
The remainder of this paper is organized as follows. In Sect. 2.0 a self-contained description of the mathematical structure of the model is given, highlighting the main mathematical differences with previous epidemic modelings (e.g.: the SIR). Sect. 3.0 reports the software implementation of the model that has been adopted to perform simulations of epidemics shown in Sect. 4.0. Here three sets of simulations are shown that demonstrate the model's capability to foresee successive epidemic waves, its ability to account for some phenomena typical of the COVID-19 pandemic and its asymptotic behavior. Sect. 5.0 resumes the principal findings of this work.
Material and methods
The convolutional model of epidemics
- The new model assumes that the whole population can be divided into five classes, as specified below.
- The susceptible class s(t) contains those individuals who can be infected.
- Exposed individuals e(t) got in contact with the pathogen but do not have developed the related illness yet.
- Infected/sick class i(t) is constituted by persons who fell ill.
- Restored r(t); those people who recovered from infection of, or exposition to, the pathogen agent.
- Deceased class d(t). Individuals who died due to the infectious agent.
We assume the above functions of time are relative abundances (frequencies) of the corresponding class normalized to the whole population, with unitary sum s(t)+ e(t)+ i(t)+ r(t)+ d(t) = 1. Class transitions happen according to the scheme of Figure 1.
It is worth noting that a part of the exposed persons can turn into the restored class as being healthy carriers or quasi-asymptomatic, while a fraction of restored people can become newly susceptible.
Fluxes among classes in Figure 1 are input and output derivatives of the related class abundances, as exemplified in the relationship below for the flux
between the susceptible and exposed:
We assume that fluxes and input derivatives always are non-negative, while output derivatives are non-positive. With this convention, the algebraic sum of input and output derivatives equates to the total time derivative that can be written as a function of fluxes.
The system of equations (2) needs to be completed by additional equations that specify the rules for computing the six fluxes.
All fluxes, except
, are calculated according to a convolutional model that links the input signal (flux) of a class with its output, adding delay and spreading which are regulated by the Impulse Response Functions (IRFs) characteristic of each class transition. These IRFs are representative of specific biophysical phenomena, such as the incubation period of the disease (for e(t) to i(t) transition) or the healing time i(t) (to r(t)). IRFs are indicated in Figure 1 near the exit fluxes of each class.
Let us consider the Input/Output (I/O) relationship for a generic class with an input flux
and an exit flux
. The relative population abundance
entering the class between τ and τ + dτ is subsequently transferred into the output flux at time t according to the IRF
This equation implies that class transition is stationary (aka Time-Invariant, TI): delaying the input flux results in a delayed version of the output flux without affecting its profile. Assuming the epidemy broke out at time t = 0 (all fluxes are zero for negative times), we find the total output flux at the time t by integrating the input contributions over all the instants τ preceding t:
The right-hand side of Eq. (4) is the convolution between the input flux and the IRF for the considered class transition, meaning that class transitions behave like a TI Linear System (LS). The causality of this LS is guaranteed since admitted IRFs are null for negative values of their argument (input time τ in the future of the output time t). This feature allows us to arbitrarily extend the upper bound of the integral on the right-hand-side of Eq. (4). Conservation of the relative number of individuals throughout the class transition is guaranteed by the unitary area of function h(t). These conditions are recapped by the following equations.
We propose gaussian IRFs for our model, bell-shaped functions that are frequently adopted to fit epidemic data  and that avoid unrealistic wings for large t. However, the IRS's profile may be characteristic of the modeled disease and should be adapted accordingly. A typical instance of gaussian IRF is detailed in Eq. (6) where T is the average delay of class transition taking place and σ it's spreading.
Whenever the flux
is known, the iterative application of Eq. (4) to each class of Figure 1 leads to the determination of all other fluxes.
To ascertain the flux
, it is necessary to find the probability dP0 for a susceptible to come in contact with the pathogen in a given time interval dt. dP0 is given by the average frequency of encounters f0 between two individuals of that population multiplied by the conditional probability that one belongs to the susceptible class and the other to the exposed or infected classes. We assume that exposed and infected individuals can both transmit the infection, a characteristic that seems important for such diseases as COVID-19. This probability is outlined by the equation below:
A casual encounter between a susceptible and an individual of class i(t) or e(t) can happen in two distinct orders generating factor 2 on the right-hand side of Eq. (7), and the subsequent Eqs. (8,9). If the two individuals are indicated as A and B, the encounter between a susceptible and an infected takes place either as
, or as
, both with probability
. The probability that a randomly chosen individual belongs to the infected compartment is
, and the correction factor
is the total living population and expresses the condition that dead individuals cannot participate in encounters. Previous models [3,29] neglect this correction, underestimating the conditional probability of casual encounters. Only a fraction of these contacts, called transmissibility, can truly transmit the pathogen, and this factor should be lesser for exposed than infected people. Let b indicate the transmissibility for infected individuals, and ε . b, with
, be the transmissibility of exposed ones. Multiplying the two addenda of dP0 by the respective transmissibility, we obtain the people's amount
leaving the susceptible class within the time interval dt:
This relationship can be conveniently condensed as shown in Eq. (9):
Using this outcome, it is possible to write the equations that govern the time evolution of I/O fluxes.
Eqs. (2,10-11) form the new convolutional model described in this work. The four parameters Ae, Air, Aid, and Ar are area normalization constants of the related IRFs and are implicit functions of the delay and dispersion. u(t) is the Heaviside step function that makes IRFs null for negative values of their argument, α [1-α] is the fraction of exposed people who recover without developing illness [fall sick], and β [1-β] is the fraction of sick individuals who recover [die].
Remarkably, the output fluxes in the epidemic model discussed in this work are not represented by a decay process as happens in almost all deterministic epidemic modelings [3, 4, 15, 27, 29, 30]. It should be noted that the exponential decay is an improbable assumption for epidemic models because the output flux is proportional to the abundance of the considered class irrespective of when people entered it. Due to this essential feature, the convolutional model avoids impossible or implausible events such as the recovery or death of infected people immediately after their infection. In different wording, the model herein discussed can provide a realistic timing for the passage of population fractions among existing model compartments.
Let us note that while each class population is the probability that a randomly selected individual belongs to that class, the IRFs h(t-τ) are proportional to conditional probability functions; i.e.: the probability that an individual leaves the considered class at time t provided that he entered that class at time τ ≤ t. In different wording, IRFs behave like “time propagators” for fluxes among classes and, indirectly, for class probabilities(s(t), e(t), i(t), r(t), and d(t)).
Table 1 details the free parameters of Eqs. (10,11) necessary to fit the model to the epidemy of interest. The possibility that recovered persons turn back into the susceptible class can be precluded by taking γ = 0.
The two output fluxes of the class i(t) are regulated by different IRFs, hir (t) toward the restored class and hid (t) toward the deceased class. This means that this model can provide different timing of evolution for the fraction of infected people who die with respect to those who survive the disease.
The asymptotic behavior
Due to the IRFs‘ normalization, the overall death rate of the new model equates (1-α).(1-β). Whenever γ = 0, the asymptotic limit of susceptible and death classes (
are linked together with the death rate by the simple relationship shown below:
The convolutional model admits a quasi-stationary solution in which i(t) and e(t) are almost constant, with vanishing derivatives. Imposing this condition into Eqs. (2) leads to the equations below.
Substituting the second of these relationships into the first, we obtain:
Eq. (14) proves that in this regime the decrease of susceptibles is continuously absorbed by the restored and deceased classes, apart from the
flux transporting recovered persons back into the susceptible class. Simulations performed with our model show this asymptotic behavior for large t frequently.
Basic reproduction number
The estimation of the basic reproduction number  of the convolutive model is made complex by the circumstance that two classes, e(t) and i(t), can transmit the infection, giving rise to three outflows having different time constants. Using the same flux weighting shown in Eqs. (10-11) we write:
- The standard deviation of each IRF is far below the related delay; e. g.:
. Whenever the standard deviation is comparable to (or greater than) the delay, the integration domain only includes a small slice of the left flank of the IRF, hence the integral is dominated by contributions originating on the right flank of the IRF. In this case, exposed and infected individuals propagate the disease for an average time greater than the delay.
- γ = 0 and the overall death rate is small, thus no correction for the division by the
the factor is necessary.
- The epidemy is in its starting phase
A formal solution to epidemic equations
Whenever γ = 0 the fundamental equation of the convolutional model reduces to:
Introducing a new function of time p(t) defined as the primitive of
, and performing simple mathematical manipulations a formal solution of the model is obtained:
Expressing the other classes as a function of p(t) is rather complex and is neglected. Let us note that a similar solution can be obtained for almost any deterministic epidemic model, as for the SIR (Eqs. (18)) whose complete formal solution is shown by Eqs. (19):
Eqs. (17,19) is the quadrature of the related epidemic model but since the function p(t) remains unknown they do not permit us to deduce the epidemy evolution, thus they are formal solutions to the related problems. Nonetheless, they can be useful to interpret the outcome of simulations discussed in the next Sections.
Model implementation has been made by developing a custom software code utilizing the C/C++ programming language. The systems of coupled integrodifferential equations given by Eqs. (2,10,11) have been solved by adopting Euler’s method with a discrete time step ∆t far below the least standard deviation among the four IRFs of Eqs. (11). Discrete time evolution equations are recapped below, where the symbol * represents convolution, while the apex stands for time derivative.
Results and discussion
Simulations aim at demonstrating specific features of the convolutional model, comparing its performance to those achieved by the SIR model , and investigating some aspects of the COVID-19 pandemic. Unfortunately, available data for this epidemy show significant differences among different countries concerning sampling methods and classification of deaths, hence comparison between theoretical predictions and real data will be made at a qualitative level. All simulations have been computed starting the infection at t = 0 with a
exposed persons. The SIR model has also been solved utilizing Euler’s approach.
In this simulation, we have investigated an important issue of epidemics, i.e., the formation of oscillatory epidemic propagation modes that have been reported during the current COVID-19 epidemic frequently . We have selected the model-free parameters shown in Table 2 where time constants (T,σ) are expressed in days. The value γ = 0.8 has been chosen to investigate the model's ability to foresee successive epidemic waves that are summarized in Figure 2.
Let us note that the new model can reproduce oscillating temporal patterns of epidemics, which are forced by the high value imposed on the γ parameter. However, simulated oscillatory patterns always are damped, hence for a large time, the model reproduces the quasi-stationary regime dictated by Eqs. (13-14) notwithstanding the large γ value. The period of such oscillations seems to slightly increase with increasing time.
Simulation #2 aims at elucidating the principal difference between the standard SIR and the new epidemic convolutional model discussed in this work. The parameters of this simulation roughly represent the COVID-19 pandemic. The maximum incubation period is believed to not exceed 14 days [33-37] with an average between 5.1 days [35,36] and 8 days , so we have chosen
days. The fatality rate has been estimated between 2.2% , 3.8% , and 5.6%  (4% in our simulation), while the time between symptom onset and death ranged from about 2 weeks to 8 weeks in , was estimated to be 8 days in  and 11 days in 
days in this simulation). A recovery time of 2 weeks for patients with mild symptoms (3-6 weeks for patients in critical conditions) has been reported in , which we have condensed with
days. The remaining parameters of this simulation were: γ = 0 (
(2.2 in , between 2.8 and 3.9 in ).
The simulation outputs are summarized in Figures 3, 4, and 5. Figure 3 shows the time evolution of classes of both the convolutional and the SIR models. Ratios of corresponding classes of the two models are plotted versus time in Figure 4, while Figure 5 depicts the input (
) and output (
) fluxes of class i(t), together with the dynamic estimate of the Case Fatality Rate (dCFR), which was calculated as:
Let us note that fluxes (
) are proportional to the number of outcomes in the corresponding compartments within a short time interval, hence it is the dynamic (instantaneous) value of the CFR estimate defined as deaths/(deaths + recovered) that is adopted often [28,37]. We also note that the exposed category e(t) rarely is traced precisely, so the death rate is often estimated considering the outcomes of the infected class only disregarding the flux
Figures 3 and 4 demonstrate that the signals originated by the convolutional model are significantly delayed with respect to those calculated by the SIR model with the same R0. Due to its structure, the convolutional model foresees important trouble encountered when facing the Covis-19 epidemy: when the disease becomes apparent, i.e.: i(t) assumes non-negligible values and patient hospitalization begins to grow, the class e(t) already is significantly diffused reducing the possibility to restrain the epidemy. Moreover, the convolutional model exhibits a significant persistence of i(t) and e(t)signals after the epidemy peak that is not found with the SIR. This point is evident in Figure 4 where the ratio between the signals e(t) and SIR i(t) significantly grows for large t. In other words, the convolutional model does not exhibit a sharp epidemy stopping, the simulated epidemic wave slowly vanishes away after the peak, and a tiny fraction of the population remains infectious for a long time.
Figure 5 shows that the
flux toward the d(t) class (dashed line) resulting from simulations is larger in the initial phase of the simulated COVID-19 epidemy than near its termination when the predominant flux
(dotted line) is toward the r(t) class. In the epidemic terminating phase
is almost one order of magnitude greater than
, while at the beginning the simulated fluxes are comparable. The dCFR plot confirms this property that has been obtained with intrinsically constant mortality of the simulated epidemy. From a mathematical standpoint, the apparent reduction of infection gravity with increasing time is a direct consequence of the circumstance that the healing time
is on average greater than the dying time
. As a rule of thumb, the apparent dynamic recovery rate of class i(t) at time is proportional to
, while the dynamic death rate depends on the input flux
, a greater death rate with rarer recoveries will appear near the epidemy beginning with respect to its conclusion.
A similar phenomenon has been observed in China where “the overall CFR was higher in the early stages of the outbreak” , and its origin still is debated . Comparisons among the case-fatality rates in China, Italy, and other countries [40,41] showed a partially unexplained higher mortality of the COVID-19 epidemic in Italy. However, the data analyzed in  were sampled in the initial phase of the epidemic wave in Italy and near the top of the epidemy in China. Simulations performed with the convolutional model (Figure 5) suggest that the disparity of sampling time could imply a large difference in the apparent death rate, so the differences revealed by these works [28,33,40,41] might reflect the property
rather than possible differences in epidemy management or clinical response between the countries.
It is worth noting that recent virus epidemics have shown time-varying CFRs as reported in . The SARS epidemy in 2003 was found to have an increasing CFR, lesser in the initial phase of the epidemic wave with respect to its conclusion . We point out that the convolutive mathematical modeling discussed in this work can predict this characteristic by choosing
, i.e.: a recovering time on average shorter than the time to death.
In this Section, we show the aggregate outcomes of several simulations with many combinations of the model-free parameters, all with γ = 0. We have investigated the criterion for epidemy spontaneous stopping and the asymptotic value as a function of
For both the examined models (convolutional and SIR), it was possible to find a clear link between the asymptotic value
imposed on the simulations, as shown in Figure 6 and in Eq. (25).
Where C and
are positive constants that assume values close to one (best fit). Figure 6 clearly shows how carefully the relationship of Eq. (25) predicts the asymptotic limit
of the susceptible class with a correlation coefficient greater than 0.99.
Substituting into Eq. (25) the limit
from Eqs. (17) it is possible to find a simple relationship between function p(t) and R0, as shown in the next equations.
It is easy showing that the SIR model closely obeys the same equation:
Integrals in Eqs. (26,27) can be considered as a measure of the epidemy “strength”, thus these equations show how the impact of an epidemy depends on a few parameters such as R0, k and s(0)
It is believed that an epidemy does not break out if R0<1, and that it stops when its basic reproduction number falls below 1 . For the convolutional model, these properties seem untrue, and epidemy can survive for years although with vanishing rates. This circumstance is depicted in Figure 7, where class abundances are plotted versus time for an epidemic simulation characterized by R0= 0.946429,
. As can be seen, the convolutional model foresees a tiny-amplitude epidemic wave that propagates over years with almost vanishing e(t) and i(t) time derivatives, as dictated by Eqs. (13,14). This behavior should be ascribed to the delay and spreading of the output signals typical of the convolutive model: Delay allows e(t) and i(t) to grow even with a small R0 while spreading prevents their fast zeroing. We note that the role of R0 as epidemy stopping threshold seems untrue also for different epidemic models containing nonlinear terms .
The characteristics discussed above are relevant to developing strategies for epidemy mitigation or control.
We have shown the structure of a new epidemic model holding five population compartments, which introduces the concept of population fluxes among different classes. In the mathematical construction of the model, classes are considered time-invariant linear systems that connect the input to output fluxes by the convolution of the input flux with an IRF. The convolution with the normalized IRF introduces delay and spreading in the related class transition that can be thought of as a probabilistic temporal evolution of the epidemy. This mathematical structure guarantees that people leave a given class according to a lifelike temporal scheme (the IRF) starting from the time at which the same individuals entered the concerned class. The convolutional model does not adopt the unrealistic exponential decay, implicitly accepted by most deterministic models to compute the exit rate of a given compartment.
The model includes healthy carrier contribution to epidemy propagation, the option for recovered people to become newly susceptible, several free parameters that govern pathogen transmissibility, and time delay and spreading for any class transition admitted. The model can be useful to simulate the evolution of different types of epidemics and seems to be able to predict important features of the COVID-19 pandemic.
The performed analyses pointed out the following properties.
- The convolutional model foresees important trouble encountered when facing the Covis-19 epidemy: when class i(t) assumes non-negligible values and patient hospitalization begins to grow, the class e(t) already is significantly diffused, reducing the possibility to restrain the epidemy.
- Epidemic waves simulated with the convolutive modeling are significantly delayed in comparison with those calculated by the SIR model.
- The convolutional model can predict the formation of epidemic oscillatory patterns for simulations in which recovered people can become newly susceptible to losing their immunization.
- The new model admits a quasi-stationary regime in which the decrease of s(t) is continuously absorbed by the restored and deceased classes with vanishing derivatives of e(t) and i(t).
- The new model allows a tiny epidemic wave to propagate for a long time even when .
- The dynamic fatality rate of COVID-19 deduced from simulations performed with the convolutive model appears to be greater at epidemy beginning than in its stopping phase, a model property that can justify measured epidemiological data (behaviors) that recent investigations were unable to explain. According to the convolutional model, this epidemic property is the result of a long average recovering time . The convolutive model is even able to account for a CFR that apparently increases with increasing time, as in the case of SARS, a feature that would require the adoption of .
- We have defined a quadrature function p(t) whose asymptotic limit determines the epidemy strength (i.e.: its impact on the population) and which roughly matches the value R0 at epidemy beginning. We were able to find this quadrature function for both the convolutive and the SIR models.
Points 4 and 5 are less important per se, as their effects can be quite small from an epidemic point of view, but they prove the inexistence of a sharp epidemy stopping option. They indeed demonstrate that epidemics like COVID-19 can silently propagate over years even involving tiny fractions of the population whenever
, but being able to spring off every time the basic reproduction number rise.
The features outlined above should be carefully considered when planning strategies for epidemy mitigation or control.