Introduction
Since the pioneer study of [Reference Kermack and McKendrick1] in epidemiology, SIRD (Susceptible S-Infected I-Recovered R-Death D) compartmental models are often used for modelling of infectious diseases, and are now used for the recent Covid-19 pandemic in China [Reference Massad2–Reference Zhou4]. SIRD models can be used at the patch scale or over a large area [Reference Sattenspiel and Dietz5–Reference Vaidya10], and can be coupled to probabilistic approaches [Reference Kadowaki11, Reference Aswi12] and various processes such as population mobility [Reference Bertuzzo13–Reference Massad, Coutinho and Wilder-Smith16]. When used at the scale of a patch, or of a country when the latter is considered as one homogeneous patch, SIRD models can be viewed similar to lumped hydrological models representing the links between several water cycle compartments [Reference Beven17] with rain as input (equivalent to the input S) and river flow, infiltration and evaporation as outputs (equivalent to the outputs R and D). Different hydrological modelling approaches were developed on the basis of an interaction between two categories of ‘functions’ called ‘production function’ and ‘transfer function’ [Reference Beven17, Reference Moussa, Chahinian and Bocquillon18]. The production function separates the rain into runoff, infiltration and evaporation, and therefore distributes the mass balance between different compartments. Based on mass conservation, the transfer function is considered as a filter which transforms an input signal into an output signal, and calculates the time distribution of the production function outputs. In epidemiology, we can define by analogy a ‘production function’ to represent the relationship S–I which calculates the part of S contributing to I, and a ‘transfer function’ to represent the relationships I–R–D.
For the Covid-19 pandemic, the daily number of infected cases I(t), recoveries R(t) and deaths D(t) are freely available on websites of governmental institutions [19–Reference Dong, Du and Gardner21]. From graphics, we observe that the relationships I–R and I–D have similar shapes to those obtained from hydrologic transfer functions. First, mass conservation is verified because the total number of I(t) is equal to the total numbers of R(t) and D(t). Second, the output signals R(t) and D(t) can be derived from the input signal I(t) using a mathematical formulation identical to the transfer function describing the physical advection-dispersion processes: a lag time translation (for the advection process) and an attenuation of the peak (for the dispersion process). Although on the basis of analogical properties, transfer functions have been successfully used in other disciplines such as astronomy, geophysics, soil mechanics, meteorology, oceanography, traffic simulation or biological flows [Reference Singh22, Reference Cheviron and Moussa23], to our knowledge, no applications of the transfer functions were conducted for modelling I–R–D relationships in compartmental models in epidemiology.
This study aims to adapt transfer functions used in hydrology for modelling I–R–D relationships in epidemiology. I–R–D models are generally described by partial differential equations [Reference Feng, Xu and Zhao24–Reference Krylova and Earn26] where parsimonious first-order kinetics remains to be largely used because of their simplicity, robustness and their low number of parameters [Reference Murray27, Reference Gallos and Fefferman28]. In hydrology, under the hypotheses that the system is linear and time invariant, simplified versions of the transfer function were developed like the ‘unit hydrograph’, that is equivalent to a convolution with a kernel function representing the response to a Dirac input [Reference Singh, Mishra and Jain29]. Furthermore, a linear reservoir is a one-parameter kernel function (or unit hydrograph) with an exponential decrease equivalent to first-order kinetics used in epidemiology (dR/dt = γI or dD/dt = δI with γ and δ as parameters). However, one-parameter unit hydrograph is not adapted to represent processes such as flood routing in channels where at least two-parameter unit hydrographs are needed to describe both advection and dispersion processes.
In epidemiology, first-order kinetics remains to be the simplest model [Reference Massad2, Reference Massad, Coutinho and Wilder-Smith16], and a question arises whether one-parameter first-order kinetics is sufficient to model at the country scale the I–R–D relationships for Covid-19, and whether simulations could be improved by using other mathematical formulations such as the two-parameter kernel functions. Various unit hydrographs have been developed [Reference Singh, Mishra and Jain29], but one is of interest for possible applications in epidemiology: the Hayami model [Reference Hayami30] used as a unit hydrograph resolution of the diffusive wave equation with parameters having a physical interpretation [Reference Moussa31–Reference Moussa and Majdalani34]. For applications in epidemiology, questions arise on the adaptability of the Hayami model to simulate the relationships I–R–D, the interpretation of model parameters and how could these parameters be used to compare different case studies.
This study attempts to make the Hayami solution of the diffusive wave equation available for modelling the I–R–D relationships in epidemiology. Applications were implemented on the actual 2020 Covid-19 pandemic. First, we present the I–R–D model, the Hayami model, its adaptation for compartmental models in epidemiology, and an application example to simulate R and D (given I) in the case of China where the pandemic is ending. Then, applications were implemented on data from 24 countries all round the world, with different levels of evolution of the pandemic, contamination rates and recoveries strategies. Finally, we show how the model parameters can be used to compare recoveries and deaths in different countries, and discuss the usefulness and the limitations of the model. The Supplementary material presents the theory of the unit hydrograph, and shows a comparison between the observed and calculated R(t), D(t) an A(t) for 24 countries.
The I–R–D model
The I–R–D model developed herein is based on the unit hydrograph theory in hydrology adapted for applications in epidemiology. Inputs are the observed number of daily infected cases I o(t) and the observed mortality ratio μ, and outputs are the calculated daily recovered cases R c(t) and the calculated daily death cases D c(t). The performances of the I–R–D model are obtained by comparing the observed daily recovered cases R o(t) to the calculated R c(t), and the observed daily death cases D o(t) to the calculated D c(t).
The distribution I o(t) is separated into two series: IR(t) representing the part of I o(t) that would recover, and ID(t) the part of I o(t) that would die (Fig. 1). We suppose a simple proportionality relationship between I o(t) and each of IR(t) and ID(t) according to:
Fig. 1. I–R–D model structure. The input I o(t) is the observed number of daily infected cases divided into IR(t) and ID(t) proportional to the mortality ratio μ. t is the time expressed in days. The I–R model calculates the daily number of recoveries R c(t), with a performance NSER. The I–D model calculates the daily number of deaths D c(t), with a performance NSED. The active cases A c(t) are calculated with a performance NSEA.
The I–R–D model is divided into two models: I–R and I–D. On the basis of the unit hydrograph theory [Reference Beven17, Reference Singh, Mishra and Jain29] (see details in the Supplementary material) we have
where the symbol ‘*’ represents the convolution relation, and uR(t) and uD(t) [T −1] are mathematical kernel functions with an integral equal to 1. For a Dirac (unit) input at t = 0, the output has the same mathematical equation as uR(t) (or uD(t)). In hydrology, the functions uR(t) and uD(t) are called unit hydrographs [Reference Singh, Mishra and Jain29]. Section ‘Adapting the Hayami kernel function to epidemiology’ shows how these functions can be calculated.
The model also enables us to compare the observed cumulated number of recoveries R to(t) to the calculated R tc(t), the observed cumulated number of deaths D to(t) to the calculated D tc(t), and the observed daily active cases A o(t) to the calculated A c(t) with the below equation:
The performances of the I–R model, the I–D model and the active cases are calculated using the Nash–Sutcliffe efficiency (NSE) criteria [Reference Nash and Sutcliffe35] traditionally used in hydrology, with the below equations, respectively:
with Rm the mean value of R o(t), Dm the mean value of D o(t), Am the mean value of A o(t) and n the number of time steps. The criteria NSE is inferior to 1 with an optimum corresponding to NSE = 1. NSE = 0 corresponds to using the mean value of observations as a benchmark predictor and is regularly used as a benchmark to compare models [Reference Schaefli and Gupta36, Reference Moussa37]. NSE < 0 indicates that the model is a worse predictor than the mean of observations. However the NSE is an indicator of performance which is sensitive to peak values. Increasingly an alternative metric, the Kling–Gupta efficiency (KGE) [Reference Gupta38] is used instead with the below equations, respectively
where rR is the linear correlation between R o(t) and R c(t), σR ,o is the standard deviation of R o(t), σR ,c is the standard deviation of R c(t), μR ,o is the mean of R o(t) and μR ,c is the mean of R c(t):
where rD is the linear correlation between D o(t) and D c(t), σD ,o is the standard deviation of D o(t), σD ,c is the standard deviation of D c(t), μD ,o is the mean of D o(t) and μD ,c is the mean of D c(t):
where rA is the linear correlation between A o(t) and A c(t), σA ,o is the standard deviation of A o(t), σA ,c is the standard deviation of A c(t), μA ,o is the mean of A o(t) and μA ,c is the mean of A c(t).
Knoben et al. [Reference Knoben, Freer and Woods39] analysed the relationships between the NSE and the KGE, and concluded that NSE and KGE values cannot be directly compared because their relationship is non-unique. They also showed that using the mean observations as a benchmark does not result in KGE = 0 but KGE = −0.41. Thus KGE > −0.41 indicate that a model improves upon the mean observation benchmark.
In the following, we note the input i(t) for IR(t) or ID(t) and the output o(t) for R c(t) or D c(t). Although i(t) and o(t) are continuous functions in time in hydrological models, herein i(t) and o(t) are discrete functions representing the distribution of the number of cases per day, with t an integer representing the day number. In hydrology, different approaches can be used to calculate u(t) [Reference Singh, Mishra and Jain29].
The exponential model, similar to the solution of the first-order kinetics used in compartmental models in epidemiology, is a unit hydrograph considered as a reference model as below:
where k [T] is a parameter representing linear reservoir retention.
Adapting the Hayami kernel function to epidemiology
Saint-Venant [Reference Saint-Venant40] formulated the system of partial differential equations (continuity and momentum) to describe one-dimensional, gradually-varied, unsteady flow in rivers. The solution of the Saint-Venant equations has given rise to a number of numerical methods because no analytical solution is available. In practical applications of flood routing in natural channels, the acceleration terms in the Saint-Venant equations can be neglected, and the system is reduced to one parabolic equation, the diffusive wave equation [Reference Moussa31, Reference Moussa and Bocquillon32]:
where x [L] is the length along the channel, t [T] is the time and the celerity C(Q) [LT −1] and the diffusivity D(Q) [L 2T −1] are functions of the discharge Q [L 3T −1]. Let I(t) and O(t) be, respectively, the upstream inflow and the downstream outflow. In the particular case of a semi-infinite channel, no physical downstream boundary condition exists, and C(Q) and D(Q) constant, the diffusive wave equation has as solution the analytical Hayami (1951) equation [Reference Hayami30]:
where u(t) is the Hayami kernel function:
with two parameters θ [T] the lag time and z [dimensionless] a shape parameter. When z tends to zero, the Hayami equation has comparable shape as the exponential model equation (7). For hydraulic applications, equation (10) is generally written as a function of two physically based parameters, changing the two-parameter (θ, z) into (C, D), with C = L/θ (L being the length) and D = L 2/(4θz). The use of the two-parameter (C, D) is more adapted for hydraulic application because both parameters have physical significance. The use of the two-parameter (θ, z) is more adapted for hydrologic applications on catchments because θ represents the lag time between i(t) and o(t), and can be thus estimated. However, z is empirical and needs to be calibrated. For applications in epidemiology, the parameter θ represents the lag time between either IR(t) and R o(t), or ID(t) and D o(t), and hence can be easily interpreted, whereas the parameter z has no significance for the relationships I–R and I–D. This is why equation (10) needs to be rewritten according to two parameters that can be more easily interpreted in epidemiology. For that, i(t) and o(t) can be considered as distributions characterised respectively as below (Fig. 2):
• GI and GO: the centres of gravity.
• TI and TO: the abscissae of GI and GO which represent the means of i(t) and o(t).
• sI and sO: the standard deviations of i(t) and o(t), and hence sI 2 and sO 2 are the variances.
Fig. 2. Example of the distributions of the input i(t) and the output o(t). GI and GO are, respectively, the centres of gravity of i(t) and o(t). TI and TO are the abscissae of GI and GO representing the means of i(t) and o(t). sI and sO are the standard deviations of i(t) and o(t). θ is the time delay between the two centres of gravity.
The values of TI and TO depend on the time origin (t = 0) whereas sI and sO do not depend on time origin. Following Moussa [Reference Moussa31, Reference Moussa and Bocquillon32], we have
The parameter θ represents the time delay between the two centres of gravity GI and GO. Let τ be
The new parameter τ has a time dimension, and is proportional to the square root of the difference of the variance of o(t) and i(t). Combining equations (11) and (12), we have the relationships:
Substituting equation (13) in (10) we obtain the expression of the Hayami unit hydrograph function of (θ, τ) instead of (θ, z):
The Hayami u(t) given by equation (14) is applicable for the continuous function data. But, the Covid-19 pandemic problem dealing with I–R–D and active variables are observed in discrete time interval of 1 day. This necessitates the conversion of the continuous function u(t) to discrete 1-day time interval. To obtain at the discrete u(t) of the Hayami model at Δt time interval (with Δt = 1 day herein), we use the discrete u(t) of the Hayami model using a similar approach as the method described by [Reference Todini and Bossi41] (see details in the Supplementary material). Another solution consists of subdividing the daily time step into finer numerical time steps under the hypothesis of uniform distribution of data. The finer the time step, the closer one gets to a continuous function of time (see the sensitivity analysis for Δt = 10 min, 1 h, 3 h and 1 day in Section ‘Sensitivity analysis of the I–R–D Hayami model to the time step of calculation’). The I–R–D model proposed herein was developed in Matlab® and integrated in the MHYDAS hydrological model (MHYDAS-IRD) [Reference Moussa, Voltz and Andrieux42].
The two parameters (θ, τ) are related to the means and variances of the two distributions i(t) and o(t) by equations (10) and (11). The advantage of equation (14) using (θ, τ) instead of (θ, z) is that both parameters (θ, τ) can be estimated from data analysis without any calibration. However, the calibration of (θ, τ) can be also undertaken in order to improve the model performances.
Figure 3a shows the shape of the output distribution for different values of θ. A small θ value gives a sharp output distribution with a small time delay between the input and the output, whereas a large θ value gives a damped output distribution with a large time delay between the input and the output. Figure 3b shows the output distribution for different values of τ. A small τ value gives a damped and more asymmetric output distribution with a small time delay between the input and the output, whereas a large τ value gives a sharp and more symmetric output distribution with a large time delay between the input and the output. Figures 3c and d show the cumulative distribution function (CDF) of Figures 3a and b: they represent the time evolution of the total number of o(t) when i(t) is a Dirac delta function.
Fig. 3. Examples of the Hayami unit hydrograph u(t). (a) For τ = 10 days and different values of θ. (b) For θ = 20 days and different values of τ. (c) CDF of (a) and (d) CDF of (b).
In the applications, we compare the reference one-parameter exponential model (equation 7) with kR the parameter of uR(t) and kD the parameter of uD(t), to the two-parameter Hayami model (equation 14) with (θR, τR) the parameters of uR(t) and (θD, τD) the parameters of uD(t). The mortality ratio μ is supposed to be known from the observed data. For both models, the parameters of uR(t) are calculated separately from those of uD(t). The performances are NSER and KGER for I–R, NSED and KGED for I–D and NSEA and KGEA for the active cases (equations (5) and (6)). The parameters are calibrated using an iterative automatic trial-and-error method maximising the NSE criteria. This method can be used for both exponential and Hayami models, for datasets (I o(t), R o(t) and D o(t)) that are partial (pandemic is evolving) or complete (pandemic is ending).
Applications on the Covid-19 pandemic
Data
Applications were implemented on the actual Covid-19 pandemic which occurred in Wuhan, Hubei Province, China, in December 2019. Data are available from Worldometer [20] Covid-19 data used by Johns Hopkins CSSE [19], governmental institutions and many others [Reference Dong, Du and Gardner21]. Daily data for I o(t), R o(t) and D o(t) are available from 1 January until 19 May 2020 for 24 countries: Australia (denoted AU), Austria (AT), Belgium (BE), China (CN), Cuba (CU), Czechia (CZ), Denmark (DK), France (FR), Germany (DE), Iceland (IS), Iran (IR), Italy (IT), Japan (JP), Malaysia (MY), New Zealand (NZ), Romania (RO), Slovakia (SK), South Korea (KR), Spain (ES), Switzerland (CH), Thailand (TH), Turkey (TU), USA (US) and the whole world (World). The latter case enables us to study an overall average trend in the whole world.
Figure 4 shows the characteristics of the data for all 24 countries, where for each country It is the total number of infected cases, Rt the total number of recoveries and Dt the total number of deaths. Figure 4a shows the large range of variability of It with 2 × 103 (IS, NZ, SK) < It < 5 × 106 (World). We define also the index (Rt + Dt)/It representing the evolution of the pandemic at the date of data availability 19 May 2020: the index is close to 1 when the pandemic is ending (e.g. CN, IS, NZ) and lowest than 0.5 (e.g. BE, US, World) when the pandemic is still actively evolving. Figure 4b shows the large range of variability of the mortality ratio μ = Dt/It with 0.005 (IS) < μ < 0.16 (BE). Probably different data collection strategies have been carried out among the different countries. However, the I–R–D model was applied separately to each country, and therefore the calibrated sets of parameters can be considered as descriptors of the I–R and I–D relationships for each country.
Fig. 4. Data characteristics for the 24 countries, where for each country It is the total number of infected cases, Rt the total number of recoveries and Dt the total number of deaths. Data are available from 1 January to 19 May 2020. (a) Pandemic evolution index (Rt + Dt)/It function of It and (b) mortality ratio: μ = Dt/It.
We compare the exponential and the Hayami I–R–D models to simulate the relationships I–R, I–D and the active cases. As we observe noisy signals for I o(t), R o(t) and D o(t), simulations can be improved by smoothing these signals. For both the exponential and the Hayami models, we also compare four smoothing strategies: without smoothing, 3-days moving average, 5-days moving average, and 7-days moving average. A sensitivity analysis was also undertaken to study the impact of the time step of calculation on the calibrated parameters and the model performances, by subdividing the daily time step into finer numerical time steps (Δt = 10 min, 1 h, 3 h and 1 day). For the I–R model, we note NSER, NSER 3, NSER 5 and NSER 7, the values of NSE respectively for the four smoothing strategies. For the I–D model, we note NSED, NSED 3, NSED 5 and NSED 7. For the active cases, we note NSEA, NSEA 3, NSEA 5 and NSEA 7. For the I–R model, we note KGER, KGER 3, KGER 5 and KGER 7, the values of KGE respectively for the four smoothing strategies. For the I–D model, we note KGED, KGED 3, KGED 5 and KGED 7. For the active cases, we note KGEA, KGEA 3, KGEA 5 and KGEA 7.
First, we present an example of application and the sensitivity analysis on China where the pandemic is ending. Then, we show the results for all 24 countries, for both the exponential and Hayami models, and for the three smoothing strategies. The Supplementary material shows the comparison between the observed {R o(t), D o(t), A o(t)} and the calculated {R c(t), D c(t), A c(t)} for all countries.
Application on China
First, we study in detail the case of China where the pandemic was ending by the mid of April 2020. The total number of infected I is 82 123 cases, recoveries R is 78 042 cases and deaths D is 3324 cases and the mortality ratio is μ = 0.041. Figures 5 and 6 show the observed {R o(t), D o(t), A o(t)} and the calculated {R c(t), D c(t), A c(t)} for respectively the exponential and the Hayami models after calibration of the parameters using data with 3-days moving average. For the I–R model, the calibrated parameters of uR(t) are for the exponential model kR = 26.5 days with fair performance NSER 3 = 0.648 (Fig. 5a), and for the Hayami model θR = 21.8 days and τR = 10.3 days (very close to those calculated from equation (10) θR = 20.3 days and τR = 9.5) with excellent performance NSER 3 = 0.99 (Fig. 6a). The calibrated values of the lag time kr and θR are comparable to cure rate for infections (17–20 days) obtained in [Reference Zhou4]. For the I–D model, the calibrated parameters of uD(t) are for the exponential model kD = 9.3 days with very good performance NSED 3 = 0.93 (Fig. 5b), and for the Hayami model θD = 10.1 days and τD = 11.1 days with excellent performance NSED 3 = 0.95 (Fig. 6b). When comparing the observed and calculated active cases A(t), we observe that the exponential model gives good results NSEA 3 = 0.82 (Fig. 5c) whereas the Hayami model gives much better results NSEA 3 = 0.95 (Fig. 6c). Finally, when comparing the cumulated observed and calculated number of recoveries, and the cumulated observed and calculated number of deaths, we observe also fair performance of the exponential model (Fig. 5d) and excellent performance of the Hayami model (Fig. 6d).
Fig. 5. I–R–D exponential model application on Covid-19 in China (with 3-days smoothing average of data) using the calibrated parameters. (a) Comparison of the observed and calculated recoveries; (b) comparison of the observed and calculated deaths; (c) comparison of the observed and calculated active cases; (d) comparison of the cumulated observed infected cases, the observed and calculated recoveries, and the observed and calculated deaths. Data are available from 1 January 2020 to 19 May 2020.
Fig. 6. I–R–D Hayami model application on Covid-19 in China (with 3-days smoothing average of data) using the calibrated parameters. (a) Comparison of the observed and calculated recoveries; (b) comparison of the observed and calculated deaths; (c) comparison of the observed and calculated active cases; (d) comparison of the cumulated observed infected cases, the observed and calculated recoveries, and the observed and calculated deaths. Data are available from 1 January 2020 to 19 May 2020.
These results show that the one-parameter exponential model is sufficient and give comparable results to the Hayami model only for modelling the relationship I–D. However, for modelling the relationship I–R and the active cases, the Hayami model highly improves the performances in comparison with the exponential model, presenting also the advantage that its parameters can be easily estimated from the analysis of the data distributions of I(t) and R(t). Figure 7 shows the uR(t) and uD(t) obtained with the calibrated parameters for China. For uR(t) (Fig. 7a and the corresponding CDF in Fig. 7c), we observe different behaviours for the two models. The exponential model cannot reproduce a unit hydrograph that rises to a maximum and then falls down, which explains the fair performance of the exponential model. However, the Hayami model succeeds because it enables us to take into account a lag time translation (comparable to the advection process) and an attenuation of the peak (comparable to the dispersion process). For uD(t) (Fig. 7b and the corresponding CDF in Fig. 7d), we observe very comparable results for both models.
Fig. 7. Comparison of the unit hydrographs calibrated for the exponential and Hayami models for the Covid-19 in China: (a) uR(t); (b) uD(t); (c) CDF of (a) and (d) CDF of (b).
Sensitivity analysis of the I–R–D Hayami model to the time step of calculation
This section presents a sensitivity analysis of the Hayami I–R–D model on China. Compared to applications in hydrology where the discharge is a continuous function of time, the major difference for epidemiological applications of the Hayami I–R–D model is that in epidemiology the data are a discrete time function representing the daily number of cases. However, all unit hydrograph approaches, including the Hayami model, can be easily adapted for applications on discrete or continuous data as discussed in Section ‘Adapting the Hayami kernel function to epidemiology’ [Reference Todini and Bossi41, Reference O'Connor43]. Around the reference data time step Δt = 1 day, we conduct a sensitivity analysis by subdividing the daily time step into finer time steps (Δt = 10 min, 1 h and 3 h) under the assumption of a uniform distribution of the number of cases. The finer the time step Δt, the closer one gets to a continuous function of time. For the data used, we also compare four smoothing strategies: without smoothing, 3-days moving average, 5-days moving average and 7-days moving average.
Table 1 shows the values of the calibrated parameters of the I–R model (θR and τR) and the corresponding criteria functions (NSER and KGER), the calibrated parameters of the I–D model (θD and τD) and the corresponding criteria functions (NSED and KGED), and the criteria functions corresponding to the actives cases (NSEA and KGEA). For the fine time steps (10 min, 1 h and 3 h), the calibrated parameters and the corresponding performance criteria remain very close in comparison with those obtained with the reference time step of 1 day. This is due to the fact that the use of a fine time step only allows tending to continuous time functions (i.e. Δt = 10 min) but using the same data as for Δt = 1 day. However, we observe an improvement in model performances and a slight change in the values of the set parameters when the original data are smoothed (21.5 < θR < 21.8 days and 10.3 < τR < 10.5 days). The improvement of model performances comparing the data without smoothing to 7-days moving average smoothed data are as follows: NSER = 0.966 and 0.992, KGER = 0.603 and 0.619, NSED = 0.915 and 0.978, KGER = 0.966 and 0.987, NSEA = 0.997 and 0.998, KGEA = 0.739 and 0.740. As the model results are more sensitive to smoothing data then to the use of finer time steps, in the following we limit the applications to compare the performance of the models using the daily time step for the four smoothing strategies (without smoothing, 3-days, 5-days moving and 7-days moving average).
Table 1. I–R–D Hayami model application on Covid-19 in China for different time steps analysis by subdividing the daily time step into finer time steps (Δt = 10 min, 1 h, 3 h and 1 day) and four smoothing strategies (without smoothing, 3-days moving average, 5-days moving average and 7-days moving average): the calibrated parameters of the I–R model (θ R and τR) and the corresponding criteria functions (NSER and KGER), the calibrated parameters of the I–D model (θ D and τD) and the corresponding criteria functions (NSED and KGED), and the criteria functions corresponding to the actives cases (NSEA and KGEA).
Application on 24 countries
This section aims to show how the methodology can be applied, and how the parameters could be used to analyse the variability among different countries. For a part of the 24 countries, the pandemic is still evolving, and the series I o(t), R o(t) and D o(t) are still incomplete. Hence, the parameters of both the exponential and the Hayami models calibrated on the available data, may slightly change when the pandemic ends. For each country, we obtain generally similar values of the calibrated parameters for the four smoothing strategies (without smoothing, 3-days moving average, 5-days moving average and 7-days moving average) but the performances (NSE and KGE) may vary drastically with the smoothing strategy.
Figure 8 shows a comparison of the NSE (Figs 8a–c) and KGE (Figs 8d–f) of the exponential and the Hayami models for the four smoothing strategies for all studied countries. For the I–R model, the 7-days moving average improves the performances NSE (Fig. 8a) of both the exponential and the Hayami models in comparison with the two remaining smoothing strategies. Moreover, the Hayami model performs much better than the exponential model: for the 5-days moving average, NSER 5 > 0.90 for four countries with the exponential model, and for 20 countries with the Hayami model (Fig. 9a). For the I–D model, Figure 8b shows that for both the exponential and the Hayami models, the 3-days, 5-days and 7-days moving average give similar results but better than without smoothing. The Hayami model gives comparable and slightly better results than the exponential model: for the 5-days moving average, NSED 5 > 0.90 for 12 countries with the exponential model, and for 16 countries with the Hayami model (Fig. 9b). For the active cases, Figure 8c shows that the three smoothing strategies give similar results for both the exponential and the Hayami models, because the active cases is a cumulated function, and consequently is less sensitive to noisy data. However, the Hayami model gives better results than the exponential model: for the 5-days moving average, NSEA 5 > 0.95 for nine countries with the exponential model, and for all 24 countries with the Hayami model (Fig. 9c).
Fig. 8. For the I–R model (a and d), the I–D model (b and e) and the active cases (c and f), comparison of the exponential model (denoted Exp) and the Hayami model (Hay) for four different smoothing strategies: without smoothing, 3-days moving average (3d), 5-days moving average (5d), and 7-days moving average (7d). The values of the NSE (a, b and c) and KGE (d, e and f) performance criteria are classified by ascending order for the 24 countries.
Fig. 9. Comparison of the performances NSE (a, b and c) and KGE (d, e and f) of the exponential (denoted Exp) and the Hayami models (smoothing data with 5-days moving average) for : the I–R model (a and d), the I–D model (b and e), and the active cases (c and f).
Similar results are obtained with the KGE. The Hayami model performs much better than the exponential model: for the 5-days moving average, KGER 5 > 0.90 for two countries with the exponential model, and for 16 countries with the Hayami model (Fig. 9d); KGED 5 > 0.90 for six countries with the exponential model, and for 12 countries with the Hayami model (Fig. 9e); KGEA 5 > 0.90 for seven countries with the exponential model, and for 14 countries with the Hayami model (Fig. 9e).
Finally, Figure 10 shows the Hayami calibrated parameters for the 24 countries for both the I–R model (θR and τR) and the I–D model (θD and τD). The parameters range are: 11 days (IR) < θR < 89 days (US), 1.5 days (IR) < τR < 99 days (US), 4 days (DK, MY) < θD < 30 days (KR), 1 day (NZ) < τD < 27 days (KR). High values of θR (or θD) correspond to a long time interval between infection and recoveries (or deaths). High values of τR (or τD) correspond to a high difference of the variance between I o(t) and R o(t) (or I o(t) and D o(t)) as given in equation (9) and shown in Figure 2. Some particular cases can be examined such as Iran which has the shortest recovery period of θR = 11 days with small τR = 1 day in comparison with results of other countries. A small value of θR can be due either to a late date of detection of infected cases and consequently a shorter recovery period, or to the recovery criterion with, for example a single negative test instead of two for recovered patients, which shortens the contamination period. The small value of τR is an indicator that both I(t) and R(t) have the same dynamic with low diffusion. Conversely, France and the US data are characterised by a large τR (respectively, 82 and 99). This is mainly due to the fact that the R(t) curves were rising and incomplete by mid-May 2020. Consequently, the calibration of the I–R–D model on the NSE criterion favouring the high values of the observation can induce a poor estimation of the parameters. This is not the case where the pandemic is over (e.g. Austria, China, Germany, Iceland, New Zealand, etc.), where small values of τR (between 6 and 10 days) are calibrated.
Fig. 10. Hayami calibrated parameters for the 24 countries for: (a) I–R model (θR and τR) and (b) I–D model (θD and τD).
The four parameters (θR, τR, θD and τD) are descriptors of the relationships between the distributions I(t), R(t) and D(t). Consequently, they depend on various factors such as the measurement strategy of each country, the health policy, the population density, the presence and date of closure of the local airport, etc. Given the heterogeneities of countries, and the time evolution of health policies, it could be hard to explain as well. However, these parameters remain useful indicators on data acquisition, recoveries strategies and pandemic evolution. They are also useful to compare and classify countries and regions.
Finally, Figure 11 shows the overall very good performances of the I–R–D Hayami model (with 5-days moving average data smoothing) for the 24 countries, comparing the observed R o(t) and the calculated R c(t) recoveries, and the observed A o(t) and the calculated A c(t) active cases. Calibrated parameters and model performances will probably change when additional data will be available. However, these applications should be seen as first tests that show that the Hayami solution of the diffusive wave equation, and more generally the unit hydrograph theory, can be easily adapted to compartmental I–R–D models in epidemiology. In hydrology, the input i(t) and the output o(t) signals are continuous functions, whereas I o(t), R o(t) and D o(t) functions in epidemiology are discrete functions. Despite discontinuities in data, and the different methods among countries to define the observed I o(t), R o(t) and D o(t), the I–R–D model provides simulations mostly qualified as excellent when slightly smoothing the noisy data. The Hayami model is easy-to-use and parsimonious with only two parameters for each of uR(t) and uD(t). One main advantage of the Hayami model is that both parameters (θ, t) can be estimated using equation (11) from the analysis of I o(t), R o(t) and D o(t) when the pandemic ends, avoiding the calibration procedure. The parameters, either calibrated or estimated, are useful simple describers to compare the temporal evolution of recoveries and deaths in different countries.
Fig. 11. Comparison between the observed R o(t) and the calculated R c(t) recoveries, and between the observed A o(t) and the calculated A c(t) active cases, for 24 countries using the Hayami I–R–D model (smoothing data with 5-days moving average).
Conclusion
On the basis of the analogy between SIRD and compartmental models in hydrology, this study makes mathematical formulations developed in hydrology available for modelling in epidemiology. We adapt the ‘transfer functions’ generally used in hydrological modelling to compartmental I–R–D models in epidemiology in order to simulate the relationships between the number of infectious I(t), the number of recovered R(t), the number of death cases D(t) and the number of active cases A(t). Simplified approaches of the transfer functions such as the unit hydrograph are easy-to-use and parsimonious with a low number of parameters. We compare the one-parameter exponential model usually used in SIRD epidemiologic model to the two-parameter physically based Hayami model solution of the diffusive wave equation. Applications were implemented on the recent Covid-19 pandemic.
The application on 24 countries enables us to compare the performances of the two models. The exponential model gives very good performances for modelling the relationship I–D, but fair performances for modelling I–R and the number of active cases. For I–R, the Hayami model improves significantly the performances with excellent performances for all variables. The Hayami model presents also the advantage that its parameters can be easily estimated from the analysis of the data distributions of I(t), R(t) and D(t).
The Hayami model is parsimonious with only two parameters which are useful simple describers to compare the temporal evolution of recoveries and deaths in different countries with different contamination rates and strategies for recoveries.
These first results illustrate the interest of adapting mathematical formulations developed in a physical discipline like hydrology for applications in epidemiology. This allows epidemiology to benefit from the numerous advances in hydrology, and provides epidemiological modellers simple and easy-to-use parsimonious tools that have been evaluated in the literature, and could possibly make a modest contribution to the complex modelling exercise in epidemiology.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0950268821001011
Acknowledgements
The authors thank Claude Bocquillon at University of Montpellier, Jean-Claude Davidian at Montpellier SupAgro, Claude Casellas at University of Montpellier and Jean-Paul Lhomme at IRD for fruitful discussions and for their comments on this study. The authors thank the Editor Eduardo Massad and the reviewers for their constructive comments.
Financial support
This is a purely theoretical work.
Conflict of interest
None.
Ethical standards
Not applicable (this is a purely theoretical work with no human subject involved).
Data availability statement
Datasets for this research are available on the Worldometer's website https://www.worldometers.info/coronavirus/.
Open access
