Identification of time delays in COVID-19 data

Introduction

The outbreak of COVID-19 in 2020 and into 2021 has led to a surge in interest in the mathematical modeling of the COVID-19 epidemic, as well as the mathematical modeling of epidemics generally. These models have taken many types. A full review of the relevant literature is beyond the current work. However, some general classes of models include data-driven and machine learning approaches (Barros et al. 2022; Bhouri et al. 2021; Jha, Cao, and Oden 2020; Linka et al. 2020; Viguerie et al. 2022a; Wang et al. 2020), models based on partial differential equation (PDE) systems (Albi et al. 2022; Bertaglia and Pareschi 2021a,b; Bertrand and Pirch 2021; Grave et al. 2021; Grave and Coutinho 2021; Viguerie et al. 2020, 2021), agent-based models (Zohdi 2020), and models based on ordinary differential equation (ODE) systems. This last category is by far the most common, with such articles numbering in the thousands. Some works of this type include e.g. (Calafiore, Novara, and Possieri 2020; Choi and Ki 2020; Ferguson et al. 2020; Gatto et al. 2020; Ivorra et al. 2020; Parolini et al. 2021; Piccolomini and Zama 2020; Remuzzi and Remuzzi 2020). We note that the categories listed above are not necessarily strict delimiters; indeed, many of the above cited works, as well myriad others, have characteristics of multiple of different model classes. One commonality that nearly all the presented models share is that they are, in some way, compartmental models in the SIR family owing to the seminal work of Kermack and Mackendrick (Ogilvy Kermack and McKendrick 1927) (see also: Breda et al. 2012; Murray 2007; Iannelli and Pugliese 2015).

One of the most important aspects of COVID-19 modeling, and epidemiological modeling in general deals with the presence of time delays in the data. Time delays certainly manifest themselves via the natural incubation period of the disease; this is reflected in the employed models commonly incorporating an exposed compartment, to model the period after exposure to the virus but before symptomatic infection (Bhouri et al. 2021; Gatto et al. 2020; Linka et al. 2020; Viguerie et al. 2020) also shown in Figure 1.

Figure 1:

Flow chart describing the evolution of the various compartments in a classical SEIRD model.

Such an approach is a reasonable one, though, there are several potential problems. First of all, the exposed compartment is not necessarily easily known from data, and hence must be estimated. This may be done, for example, through rule-of-thumb estimates based on the infected population, which result in high uncertainty (Albi et al. 2022; Grave and Coutinho 2021; Grave et al. 2021). Compounding this problem further are the well-known presence of asymptomatic patients, who may spread the disease while never showing appreciable symptoms, and, at least in the earlier stages of the pandemic, were almost certainly under-counted (Albi et al. 2022; Viguerie et al. 2020).

Another major problem with the exposed compartment in most SEIR-type formulations come from estimating the sojourn time of incubation period, or the length of time spent in the compartment. Under most common compartmental approaches, the sojourn time of the exposed period is implicitly assumed to follow an exponential distribution. In practice, this may not always be a realistic (Iannelli and Pugliese 2015) assumption.

In (Guglielmi, Iacomini, and Viguerie 2022; Viguerie et al. 2022a), the authors employed a delay-differential equation (DDE) model, in which the incubation period was modeled with a time-delay term, rather than as a separate exposed compartment. Such a modeling approach eliminated the need to estimate the exposed compartment, as well as the assumption of an exponentially-distributed sojourn time for the incubation period, instead assuming a Dirac-delta (constant) distribution. Other works incorporating DDE models into the analysis of COVID-19 include (Dell’Anna 2020; Devipriya, Dhamodharavadhani, and Selvi 2021; Kumar and Erturk 2020), while references for more general epidemiological models of this type may be found in (Brauer and Castillo-Chavez 2012; Buonomo, d’Onofrio, and Lacitignola 2008; Forde 2005; Iannelli and Pugliese 2015; Murray 2007; Takeuchi, Ma, and Beretta 2000).

Aside from the presence of incubation periods, delays in data may result from other sources: additional time from initial onset of symptoms to when a person decides to be tested, the availability of testing, the speed of test result processing, testing center capacity, and other factors may have a large influence on when a given case is properly identified in the data (Sarnaglia et al. 2021). The presence of such lags is important, particularly when attempting to model measured data, as such processes are intrinsic to the data (Bastos et al. 2019). Hence, when attempting to properly evaluate different intervention strategies, for instance, such delays must be considered in the modeling process. A model designed to fit measured parameters, while not taking into account the delays in measurement, may fail provide a wholly accurate description of the dynamics. In this sense, we may view such delays as no longer time from exposure to infection, but rather time from exposure to identification.

Quantifying both the magnitude and distribution of such time delays, considering the effects of both identification lag-time as well as the natural incubation period, is thus important for evaluating intervention efficacy and allocation of resources, and is the focus of the present work. In order to accomplish this task, we first introduce, as a forward problem, an SIRD (susceptible-infected-recovered-deceased) model, in which the time from infection to identification is not assumed to follow any particular distribution. At the continuous level, this is accomplished via a convolution with an unknown probability distribution function. We then approximate the continuous probability distribute function with a discrete approximation, reducing the convolution term to a weighted sum of constant time-delays. In this way, we obtain a linear problem in terms of the distribution weights, enabling us to employ an optimization procedure against measured data to identify the weights, thus approximating the distribution of the time delays present in the data. We note that the introduced method, while applied in the current work to COVID-19, is general in nature and such a technique may be employed in any situation in which a dynamical system and its associated measurements may exhibit time-delayed dynamics whose distribution is not known.

The article is outlined as follows. We first introduce the employed delay-differential equation model and discuss some of its characteristics at the continuous level. We then introduce its discrete analogue and describe the optimization process for the identification of the time delay distribution. We then validate the proposed technique on an interval of data for the COVID-19 outbreaks in Italy and Austria demonstrating its effectiveness for different cases under different modeling assumptions. Finally, we conclude with a summary of our main findings and possible directions for future research in this area.

Delay differential equation model

The delay differential model introduced in (Guglielmi, Iacomini, and Viguerie 2022) reads as follows:

Here, s(t), i(t), r(t), d(t) represent respectively the susceptible, infected, recovered and deceased compartments at time t. The parameter σ represents the time delay term which models the incubation period of the disease, as well as potential reporting delays. For a comprehensive description of the parameters and variables of the model, with the units of measure, we refer to Table 1. Moreover the living population is given by n = s + i + r + d.

We assume the following:

1. That the time scales considered are sufficiently short such that we do not need to consider new births or non-COVID-19 deaths;

2. That the population is sufficiently well-mixed in space, such that spatial dynamics (including the Allee effect on the transmission terms, see Viguerie et al. 2021) can be ignored;

3. That the time scales considered are sufficiently short such that we do not need to consider waning immunity;

4. That recovery and mortality rates follow an exponential distribution;

5. That the model depends on time-delays, arising from both natural (incubation period) and non-natural (measurement and reporting delays, and similar) causes, whose distribution is unknown.

where σ _min and σ _max denote the extrema of the time delays. Note that the presence of the convolution term:

where w(σ) ≥ 0 is a probability density function describing the distribution of the time-delay dependence present in the model. By definition,

Under the assumption (with δ(⋅) a Dirac-δ distribution):

for a certain constant σ ̂ ∈ σ min , σ max , one obtains the model introduced in (Guglielmi, Iacomini, and Viguerie 2022). In this sense, (5)–(8) can be considered a generalization of the model introduced in (Guglielmi, Iacomini, and Viguerie 2022).

The function w is unknown in general, and its identification in terms of the available measured data is the focus of the present work. To render the problem numerically tractable, we replace w(σ) in (5)–(8) with a piecewise-constant approximation, letting {σ _j } denote a partition of the interval, assumed uniform for simplicity: σ min , σ max and Δ = σ _j − σ _j−1:

where the χ _j are simple functions, defined such that:

From the mean value theorem, for each [σ _j−1, σ _j ), there exists a σ j ̃ ∈ [ σ j , σ j − 1 ) such that:

We then define the w ̃ j in (12) such that:

We then have:

where the final line follows from the simple approximation:

Hence, by defining

we obtain the approximation:

Note that:

To see this, observe from (14), (18):

From the preceding arguments, (16), replacing w(σ) in (5)–(8) with its piecewise approximation (12) gives the following system of ordinary delayed-differential equations:

It was shown in (Guglielmi, Iacomini, and Viguerie 2022) that, for k = 1 (and thus only one time-delay σ ̂ ), the model (22)–(25), exhibits stable behavior provided that:

which, for the more general case, easily extends to:

for each j.

The significance of replacing (5)–(8) with the above approximation (22)–(25) is that we now have a problem which is linear in w _j . This allows us to make use of an optimization procedure to identify the w _j , allowing us to estimate w(σ). Therefore, we are able to numerically estimate the dependence of the data on general time-delays, assuming only that w(σ) is sufficiently smooth over the relevant interval. In contrast, many common modeling approaches assume, for example, that such delays are exponentially distributed (as in most standard compartment models), or that the delays are constant (as in most basic delay-differential equation models). We remark while we consider only distributions ω assumed to not vary in time in the present work, the above framework can nonetheless accommodate time-dependent ω, provided that such time-dependence is sufficiently smooth.

Optimization procedure and methods

Denote the measured susceptible, infected, recovered, and deceased (from COVID-19) populations as s ̃ , i ̃ , r ̃ , d ̃ , respectively. We then define the following optimization problem: For a given β, ϕ _r , ϕ _d , σ _j , j = 1: k and w = [ w 1 , w 2 , … , w k ] T :

subject to

The novelty of the proposed method above lies in the approximation (12) of w(σ), allowing us to consider the nonlinear convolution term (9) as a weighted sum of constant-delay terms. Through this approximation, it is possible to obtain an optimization problem for the time delay that is linear in w _j . This is important, as attempting to optimize for w(σ) directly on (5)–(8) gives a highly nonlinear optimization problem that is likely untractable. This method of identifying a general distribution by approximating it as a weighted sum of constant time delays, then solving an optimization problem for the weights is, to the authors’ knowledge, novel. We note that this approach is highly general and is not restricted to models of COVID-19 or epidemiological models in general. Indeed, one may use an analogous approach to identify general sojourn times in compartmental models. Moreover, the accuracy and the detailed description of the delay depend on how k is chosen. This allows this method to be very applicable for other models as well.

We note that the model (22)–(25) depends heavily not only on the approximate delay-distribution w _j , but also on the contact rate β, recovery rate ϕ _r , and mortality rate ϕ _d . These parameters are also, in general, not known. In order to obtain more robust estimates for the w _j through the optimization procedure defined above, we resort to a Monte Carlo method. The idea of this class of algorithms is to repeat random sample from the distribution function of the random variable, and perform several runs of the simulation, one for each sample. The solution will be given by the mean of all the obtained result. For a detailed description of Monte-Carlo methods we refer to (Pareschi and Toscani 2013).

Here, the Monte Carlo method consists of considering the β, ϕ _r , ϕ _d as independent random variables sampled from a Gaussian distributions. The mean of the respective distributions is obtained via an empirical parameter fitting, with the standard deviation of each distribution taken to be 5 % the value of the mean. This value was chosen based on sensitivity analyses. We then run the optimization procedure (28a)–(28c) a total of 1,000 times, corresponding to 1,000 different parameter samplings. We then identify frequency with which different w _j appear in the optimization, as well as the mean value of each w _j over the sampled trial.

We note that the result optimization procedure (16a)–(16d) depends on the choice of input data. In general, we consider the infection and deceased data to be the most reliable such data and is emphasized in our calculations. As the easiest quantity to measure, the deceased compartment is assumed the most accurate. Further, the as the rate of new infections is the driving point behind policy measures, it is considered the second most accurate measure for our purposes. However, one may modify (16a) to include other compartments, or different subsets of the compartment set.

In terms of implementation, there are several important points to note. We define the optimization procedure using the MATLAB routine fmincon. In order to solve the delay differential equation system (22)–(25), we use the MATLAB function dde23. We note that dde23 uses adaptive grid points, and hence the temporal points of the solution obtained using dde23 do not correspond to those in our dataset in general. To circumvent this, we use deval to evaluate the dde23 solution at the desired temporal points.

Results

For our measured data, we consider the second wave of COVID-19 outbreak in Italy and Austria, as reported by the newspaper Il Sole 24 ORE (Coronavirus 2021) and AGES^[1] respectively. The choice of these two Countries is mainly due to the availability of data as current number of positive cases and daily updates on recovered people, which was the most difficult data to find.

We consider the dates ranging from September 11, 2020 to February 7, 2021. This date range was chosen to represent a 150-day span, in which different governmental policies were enacted, widespread testing was available, and before large-scale population vaccination. This choice was intentional, as we wish to avoid as many confounding factors as possible.

Italy

We consider the following distributions for the parameters

(29) β ∼ N . 1131 / n 0 , . 05 ( . 1131 / n 0 )

(30) ϕ d ∼ N 1 / 940 , . 05 ( 1 / 940 )

(31) ϕ r ∼ N ( 1 / 24 , . 05 ( 1 / 24 ) ) ,

where N(μ, ς) denotes a normal distribution with mean μ and standard distribution ς and n ₀ denotes the initial living population. For t > 73,^[2] where t ₀ = 0, we replace β by β/3 to model the introduced government restrictions, designed to curb the virus spread. The values of ϕ _d and ϕ _r remain unchanged throughout the considered interval. The mean values of these parameters were obtained via a preliminary parameter tuning. For the delays, we considered 12 possible values, from σ = 2 days to σ = 35 days, with 3 days of spacing between the different days. We also performed some experiments (not shown) using different discretizations over the same range (for instance, 10 possible values) to show the robustness of the optimization procedure, which did not appreciably change the conclusions drawn.

To examine the possible differences in delays among the different compartments, we perform the optimization four times: once considering the i and d compartments together in the objective function, and for the i, r and d compartments, considering each compartment individually. For each ensemble, we report the minimum, maximum, and mean L ² error over the simulation ensemble when compared to the measured data, as well as the mean value and frequency of the different weights w _j corresponding to the delays σ _j .

The aggregated results of the simulations in terms of the infected, recovered, and deceased compartments, when optimizing over the i and d compartments are shown in Figure 2. We see good qualitative agreement with the measured and simulated data across the three compartments of interest, with the major relevant trends being captured. In Table 2, we report the mean, minimum, and maximum of each error value computed over the simulation ensemble. We see particularly strong agreement in the susceptible and deceased compartments, reasonable agreement in the infected compartment (especially considering its rapid-changing trajectory), and somewhat worse performance in the recovered compartment. Given that the r compartment was not included in the optimization, this is unsurprising.

Figure 2:

Italy: Aggregated results of the infected, recovered, and deceased compartments over the range of simulations when optimizing over the i and d compartments. We see good agreement across all compartments with the measured data, and particularly strong agreement (less than 10 % in L ² norm) in the deceased compartment.

Table 2:

Table of relative errors in L ² norm against measured data for the different compartments over 1,000 simulations. The compartment towards which the optimization procedure is ran is shown in the columns.

		Optim. wrt i, d	Optim. wrt i	Optim. wrt r	Optim. wrt d
‖ s − s ̃ ‖ L 2 / ‖ s ̃ ‖ L 2	Mean	0.0051	0.0051	0.0025	0.0035
	Min	0.0022	0.0022	0.0016	0.0009
	Max	0.0104	0.0108	0.0041	0.0079
‖ i − i ̃ ‖ L 2 / ‖ i ̃ ‖ L 2	Mean	0.1101	0.1116	0.2347	0.1464
	Min	0.0798	0.0807	0.1580	0.0826
	Max	0.2757	0.2681	0.4672	0.4992
‖ r − r ̃ ‖ L 2 / ‖ r ̃ ‖ L 2	Mean	0.2602	0.2561	0.0586	0.1903
	Min	0.071	0.0650	0.0394	0.0440
	Max	0.5913	0.6054	0.1091	0.4769
‖ d − d ̃ ‖ L 2 / ‖ d ̃ ‖ L 2	Mean	0.0695	0.0684	0.0926	0.0430
	Min	0.0249	0.0266	0.0362	0.0170
	Max	0.2063	0.1659	0.1672	0.0864

Turning our attention to the influence of time delays, we report the mean and frequency of the different weights w _j in Figure 3. We see a heavy concentration of weights values, for both frequency and mean, in the range of σ = 8 to σ = 17 days. This is consistent with what has been reported in general media, as the incubation period of COVID is believed to have an average of 3–11 days (Byrne et al. 2020; Lauer et al. 2020; McAloon et al. 2020). When combined with natural delays coming from the testing and data reporting, a total time lag in the range of 8–17 days would appear consistent with this information.

Figure 3:

Italy: Left: the different values of the w _j by frequency of occurrence when optimizing over i and d. We see a strong concentration of weights in the range of σ = 8 to σ = 17 days. Right: the mean values of the weights over the simulation ensemble. These results mirror those of the frequency plot, showing a concentration between 8 and 17 days. For both frequency and mean, σ = 11 appears as the most represented delay.

We then repeat this experiment three more times: once each for the i, r, and d compartments considered individually in the optimization process. The error tables for the i, r and d compartments can be found in Table 2. The corresponding simulation outputs for each case are plotted in Figures 4, 6, and 8. The behavior is consistent with what one may intuitively expect; in each instance, the error behavior is minimized for the considered compartment, with the simulation variance greatly reduced for the isolated compartment. We note that one can observe from the relevant tables and figures that there is a much higher degree similarity between the i and d error behavior when compared to the r compartment, which appears to act more independently.

Figure 4:

Italy: Aggregated results of the infected, recovered, and deceased compartments over the range of simulations when optimizing over the i compartment. We see good agreement across all compartments with the measured data, and particularly strong agreement (less than 10 % in L ² norm) in the deceased compartment.

This is confirmed when looking at the frequency and mean distribution of the weights w _j in Figures 5, 7, and 9 for the i, r, and d compartments respectively. In the case of the i compartment, we see concentration between σ = 8 and σ = 17 days, similar to when i and d are considered jointly; however, for i alone, σ = 14 appears as the most dominant lag term, rather than σ = 11. The behavior appears similar but even more pronounced when considering d alone. In this instance σ = 14 is again the dominant delay term; however, the grouping is tightly clustered between σ = 11 and σ = 17, with terms outside this range showing very little influence.

Figure 5:

Italy: Left: the different values of the w _j by frequency of occurrence when optimizing over i. We see a strong concentration of weights in the range of σ = 8 to σ = 17 days. Right: the mean values of the weights over the simulation ensemble. These results mirror those of the frequency plot, showing a concentration between 8 and 17 days. For both frequency and mean, σ = 14 appears as the most represented delay.

Figure 6:

Italy: Aggregated results of the infected, recovered, and deceased compartments over the range of simulations when optimizing over the r compartment. We see good agreement in the r compartment; the agreement with the d compartment remains acceptable, though somewhat lesser than other optimization choices. The error in the infected compartment is notably higher. This different error behavior is reflected in the optimization results; when optimizing over r, the optimal time delays are different than when optimizing over i and d.

Figure 7:

Italy: Left: the different values of the w _j by frequency of occurrence when optimizing over r. We see a strong concentration of weights in the range of σ = 14 to σ = 20 days. Right: the mean values of the weights over the simulation ensemble. These results mirror those of the frequency plot, showing a concentration between 14 and 20 days. For both frequency and mean, σ = 17 appears as the most represented delay. These results are less similar to those obtained when one optimizes over i and d, and suggests different time-delays present in the r compartment as compared to other compartments.

Figure 8:

Italy: Aggregated results of the infected, recovered, and deceased compartments over the range of simulations when optimizing over the d compartment. We see good agreement in the d and r compartments, with somewhat less close agreement with the i compartment.

Figure 9:

Italy: Left: the different values of the w _j by frequency of occurrence when optimizing over d. We see a tight concentration of σ in the range of 11–17 days. Right: the mean values of the weights over the simulation ensemble. These results mirror those of the frequency plot, showing a concentration between 11 and 16 days. For both frequency and mean, σ = 14 appears as the most represented delay.

In contrast, when examining the case of isolated r in 7, we see a different distribution of delay terms. Indeed, in this instance, the dominant delay terms range from σ = 14 to σ = 20, and are very tightly concentrated within this range. This suggests a distinct and longer time scale when compared to the i and d, and may explain why the error behavior for the r compartment shows a distinctly different behavior in the plots and tables shown.

Austria

We consider the following distributions for the parameters

(32) β ∼ N . 6 / n 0 , . 05 ( . 6 / n 0 )

(33) ϕ d ∼ N 1 / 550 , . 05 ( 1 / 550 )

(34) ϕ r ∼ N ( 1 / 14 , . 05 ( 1 / 14 ) ) ,

where N(μ, ς) denotes a normal distribution with mean μ and standard distribution ς and n ₀ denotes the initial living population. Note that the parameters of the random variables are different from the ones related to Italy.

For t > 66, where t ₀ = 0 we replace β by β/4 to model the introduced government restrictions, namely the second hard lockdown in Austria.

As in the case of Italy, the optimization is performed four times: once considering the i and d compartments together in the objective function, and for the i, r and d compartments, considering each compartment individually. For each ensemble, we show, as before, the aggregate results of the infected, recovered and deceased compartments, the mean and the frequency of the weights w _j of the delay.

Optimizing with respect to i and d, Figure 10, leads to a good agreement across all compartments with the measured data with a particular strong agreement in the deceased compartment, as in Figure 2.

Figure 10:

Austria: Aggregated results of the infected, recovered, and deceased compartments over the range of simulations when optimizing over the i and d compartment. We see good agreement across all compartments with the measured data, and particularly strong agreement in the deceased compartment.

Concerning the weights, shown in Figure 11, both the frequency and the mean of σ are concentrated between 8 and 20 days, similar to Figure 3.

Figure 11:

Austria: Left: the different values of the w _j by frequency of occurrence when optimizing over i and d. We see a concentration of σ in the range of 8–20 days. Right: the mean values of the weights over the simulation ensemble. These results mirror those of the frequency plot, showing a concentration between 11 and 20 days. For both frequency and mean, σ = 14 appears as the most represented delay.

For the optimization with respect to the infected (Figures 12 and 13) and recovered (Figures 14, 15) compartments considered individually, we see a good agreement with the results obtained in the case of Italy, Figure 4, 5 and Figure 6, 7 respectively.

Figure 12:

Austria: Aggregated results of the infected, recovered, and deceased compartments over the range of simulations when optimizing over the i compartment. We see good agreement across all compartments with the measured data, and particularly strong agreement in the deceased compartment.

Figure 13:

Austria: Left: the different values of the w _j by frequency of occurrence when optimizing over i. We see more σ values in the range of 11–23 days. Right: the mean values of the weights over the simulation ensemble. These results mirror those of the frequency plot, showing a concentration between 11 and 20 days. For both frequency and mean, σ = 14 appears as the most represented delay.

Figure 14:

Austria: Aggregated results of the infected, recovered, and deceased compartments over the range of simulations when optimizing over the r compartment. We see good agreement in the r compartment; the agreement with the d compartment remains acceptable, though somewhat lower r than other optimization choices. The different error behavior is reflected in the optimization results; when optimizing over r, the optimal time delays are different than when optimizing over i and d.

Figure 15:

Austria: Left: the different values of the w _j by frequency of occurrence when optimizing over r. We see a tight concentration of σ in the range of 11–20 days. Right: the mean values of the weights over the simulation ensemble. These results mirror those of the frequency plot, showing a concentration between 11 and 16 days. For both frequency and mean, σ = 14 appears as the most represented delay These results are less similar to those obtained when one optimizes over i and d, and suggests different time-delays present in the r compartment as compared to other compartments.

Figure 16:

Austria: Aggregated results of the infected, recovered, and deceased compartments over the range of simulations when optimizing over the d compartment. We see good agreement in the d compartment, with somewhat less close agreement with the r and i compartments.

However, this can not be said in the case of the deceased compartments (Figures 16 and 17). Here the two cases are quite different, see Figures 9 and 17. This is possibly explained by differences in how the various data are reported in the different countries causing delay distributions to differ; however, we are not sure as to the exact cause.

Figure 17:

Austria: Left: the different values of the w _j by frequency of occurrence when optimizing over d. We see a concentration of higher weights for σ in the range of 14–23 days. Right: the mean values of the weights over the simulation ensemble. These results mirror those of the frequency plot, showing a concentration between 14 and 23 days, even though the presence of sigma values also in the first and the last time-bars is notable. According to the mean value, σ = 20 appears as the most represented delay.

Sensitivity of delay distribution to discretization

In order to ensure the validity of our method, it is important to ensure that the results are insensitive to the discretization of [σ _min, σ _max). We perform a brief sensitivity analysis, considering the data from Austria described in the previous subsection. We consider two additional simulations, optimized over the i and d compartments (corresponding to Figures 10 and 11 above). The setup is identical to the previous simulations, except we vary the resolution of the delay discretization:

1. Coarse-resolution case: we discretize the interval [2, 35) with Δ = 5 days.

2. Fine-resolution case: we discretize the interval [2, 35) with Δ = 1 day.

We remind the reader that in the results previously shown, Δ = 3 days. We show the results for the coarse-resolution case in Figures 18 and 19 and the fine-resolution case in Figures 20 and 21.

Figure 18:

Austria: Aggregated results of the infected, recovered, and deceased compartments over the range of simulations when optimizing over the i and d compartment using a coarser discretization of the delay term. The results are similar, qualitatively and quantitatively, to those obtained with a finer discretization.

Figure 19:

Austria: Left: the different values of the w _j by frequency of occurrence when optimizing over i and d when using a coarser discretization of the delay term. Right: average value of each weight over the simulation ensemble. The results are similar, qualitatively and quantitatively, to those obtained with a finer discretization.

Figure 20:

Austria: Aggregated results of the infected, recovered, and deceased compartments over the range of simulations when optimizing over the i and d compartment using a fine discretization of the delay term. The results are similar, qualitatively and quantitatively, to those obtained with a coarser discretization.

Figure 21:

Austria: Left: the different values of the w _j by frequency of occurrence when optimizing over i and d when using a fine discretization of the delay term. Right: average value of each weight over the simulation ensemble. The results are similar, qualitatively and quantitatively, to those obtained with a coarser discretization.

In both cases, we find that the solution of the system is similar to the one obtained in the preceding section. Similarly, the delay-distributions are also similar, both qualitatively and quantitatively, compared to the previous results. From our analysis, it appears that the resolution of the delay discretization does not result in significant changes to the qualitative or quantitative profile or the identified delay distribution.

We note that increasing the resolution of the delay discretization significantly increases the computational time necessary to complete the Monte Carlo analysis. On a MacBook Pro with Processor 2 GHz, Quad-Core, Intel Core i5, with Δ = 5 (and hence 7 wt), each individual simulation required an average of 1.178 s; the average simulation instead required 3.435 s for Δ = 3 (12 wt) and 48.994 s for Δ = 1 (34 wt). Based on this evidence, it appears that one should increase the resolution of the delay term only if absolutely necessary to understand the underlying distribution, as the necessary computational cost appears to scale quadratically with the number of weights. If a coarser resolution will suffice (as appears to be the case in the present example), the significantly shorter computational time make this a more attractive choice. The authors acknowledge, however, that the optimization algorithm employed in the present work are not necessarily the most efficient choice, and that a more carefully chosen optimization algorithm may provide better scaling behavior for increasing numbers of weights.

Conclusions and future directions

In the present work, we have attempted to identify time delays present in COVID-19 data by means of newly-introduced method. This method consists of discretizing k possible time delays σ _j , j = 1: k and expressing these terms as a weighted sum ∑ j = 1 k w j i ( t − σ j ) . This formulation then yields an optimization problem which is linear in the weights w _j . After a preliminary parameter fitting, we then solve an ensemble of optimization problems, finding the weights w _j which minimize the difference between simulated and measured data over different samplings of parameter values. We recover a distribution of delay terms for the different variables, allowing us to recover the inherent time-delays within the measured data.

We tested this approach on COVID-19 data in both Italy and Austria, and found that the time-delays differ among different data points, but generally speaking range between 8 and 17 days for i and d, with the most likely value between 11 and 14 days, and over a somewhat longer time scale (14–20 days) for the recovered compartment. This is consistent with our expectations given the incubation period of the disease and the natural time delays one may expect from the testing and data collection process. Our results confirm that proper evaluation of a given intervention for COVID-19 should allow for a sufficiently long time window to pass before conclusions are made. All of the identified distributions followed a bell-like distribution resembling a normal or lognormal distribution, and a sensitivity analysis suggests that this profile is independent of the discretization of the delay term.

There are several important directions for future work. Given the increased importance of waning vaccine efficacy over time, a similar analysis given vaccine data and data on breakthrough infections could be helpful in identifying similar such time-delay trends for waning immunity. While we have elected to avoid vaccines in the current work, preferring to establish a proof of concept on a more simplified problem, the approach offered is extremely general and may be applied to these important problems. Outside of COVID-19, time delays are an important part of other types of epidemiological modeling; similar methods could be used, for example, to improve back-calculation methods for estimating HIV incidence, in which incidence is estimated by reconstructing the distribution of the time-delays between infection and diagnosis (Hall et al. 2017; Song et al. 2017; Viguerie et al. 2022b; Xia et al. 2020). Further, given the generality of the proposed optimization procedure, it may have wide applications outside of epidemiology, and may be applied to any model where time delays are present, including applications in traffic modeling, pedestrian flow, remote sensing, and other fields. From the theoretical and methodological standpoint, an extension of the framework shown here for problems in which the delays are state-dependent is also an important direction of future research. Finally, we note that the general idea shown, in which general sojourn time-distributions are approximated by a weighted sum of delay-differential equations, may be applied to approximate systems featuring more complex delay-distributions.