Next Article in Journal
Spatiotemporal Variation of Urban Plant Diversity and above Ground Biomass in Haikou, China
Previous Article in Journal
Reducing Seed Shattering in Weedy Rice by Editing SH4 and qSH1 Genes: Implications in Environmental Biosafety and Weed Control through Transgene Mitigation
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Spectral Method in Epidemic Time Series: Application to COVID-19 Pandemic

1
Université Grenoble Alpes, AGEIS EA7407, F-38700 La Tronche, France
2
Université Bordeaux, IMB, UMR 5251, F-33400 Talence, France
3
CNRS, IMB, UMR 5251, F-33400 Talence, France
*
Author to whom correspondence should be addressed.
Biology 2022, 11(12), 1825; https://doi.org/10.3390/biology11121825
Submission received: 31 October 2022 / Revised: 6 December 2022 / Accepted: 8 December 2022 / Published: 14 December 2022
(This article belongs to the Section Theoretical Biology and Biomathematics)

Abstract

:

Simple Summary

This article aims to study the times series provided by data of the daily number of reported cases of COVID-19. During the COVID-19 pandemic, most people viewed the oscillations around the exponential growth at the beginning of an epidemic wave as the default in reporting the data. The residual is probably partly due to the reporting data process (random noise). Nevertheless, a significant remaining part of such oscillations could be connected to the infection dynamic at the level of a single average patient. Eventually, the central question we try to address here is: Is there some hidden information in the signal around the exponential tendency for COVID-19 data?

Abstract

Background: The age of infection plays an important role in assessing an individual’s daily level of contagiousness, quantified by the daily reproduction number. Then, we derive an autoregressive moving average model from a daily discrete-time epidemic model based on a difference equation involving the age of infection. Novelty: The article’s main idea is to use a part of the spectrum associated with this difference equation to describe the data and the model. Results: We present some results of the parameters’ identification of the model when all the eigenvalues are known. This method was applied to Japan’s third epidemic wave of COVID-19 fails to preserve the positivity of daily reproduction. This problem forced us to develop an original truncated spectral method applied to Japanese data. We start by considering ten days and extend our analysis to one month. Conclusion: We can identify the shape for a daily reproduction numbers curve throughout the contagion period using only a few eigenvalues to fit the data.

1. Introduction

Modeling an epidemic peak requires precise knowledge of the daily data corresponding to new cases. One of the aims of the paper is to extract the value of the average daily reproduction numbers. The daily reproduction numbers vary from individual to individual and from day to day during the period of contagiousness of an individual. These numbers depend on the age of infection, i.e., the number of days since the individual contracted the infectious disease.
From a discrete model of the evolution of new daily cases, we propose to evaluate the average number R 0 ( d ) of secondary infected individuals produced by a single infected individual on each day d since his infection. For this purpose, on the top of the dominant eigenvalue, we will estimate from the data other significant subdominant eigenvalues (complex), which explain the modulation of the growth and allow better adequacy of the model to the data.
For that purpose, we reconsider the discrete-time epidemic model with the age of infection presented in Demongeot et al. [1]. This model is a discrete-time version of the Volterra integral formulation of the Kermack–McKendrick model with age of infection [2]. The variation of the number of susceptible individuals S ( t ) is given each day t = t 0 , t 0 + 1 , , by
S ( t ) = S 0 d = t 0 t 1 N ( d ) , t t 0 ,
where S ( t ) is the number of susceptible individuals at time t, and N ( t ) is the daily number of new infected at time t. Throughout the paper, we use the following convention for the sum
d = k m = 0 , whenever m < k .
As a consequence, when t = t 0 , (1) gives
S ( t 0 ) = S 0 .
We assume for simplicity that the epidemic starts from a single cohort of infected at time t 0 , then the number of infectious individuals is given by
I ( t ) = Γ ( t t 0 ) I 0 + d = 1 t t 0 Γ ( d ) N ( t d ) ,
where I 0 is the number of infected individuals at time t 0 , and Γ ( d ) is the probability for an infected to be infectious after d day of infection. In particular, we have Γ ( 0 ) = 0 .
We assume that the number N ( t ) of new infected at time t is the product of the transmission rate τ ( t ) with the number S ( t ) of susceptible individuals and the number I ( t ) of infectious at time t. That is,
N ( t ) = τ ( t ) S ( t ) I ( t ) .
By replacing I ( t ) by the right hand side of (2) in (3), we obtain
N ( t ) = τ ( t ) S ( t ) Γ ( t t 0 ) I 0 + d = 1 t t 0 Γ ( d ) N ( t d ) .
Now assuming that τ ( t ) = τ 0 and S ( t ) = S 0 are constant (over a short period of time), then we define the daily reproduction numbers as
R 0 ( d ) = τ 0 S 0 Γ ( d ) , d 0 .
The quantity R 0 ( d ) is the average number of secondary infected produced by a single infected on the day d since infection (see [1] for more details). Therefore, the basic reproduction number is the following quantity
R 0 = d = 1 n R 0 ( d ) ,
where n is the maximal duration (in days) of the infection.
Moreover, when τ ( t ) = τ 0 and S ( t ) = S 0 are constant, Equation (4) becomes a linear discrete time Volterra integral equation
N ( t ) = R 0 ( t t 0 ) I 0 ( I ) + d = 1 t t 0 R 0 ( d ) N ( t d ) ( II ) , t t 0 ,
where (I) is the number of infected produced directly by the I 0 infected individuals already present on day t 0 , and (II) is the number of new infected individuals at time t produced by the new infected individuals since day t 0 .
If we consider the first terms of the discrete time Volterra Equation (6), we obtain
N ( t 0 ) = R 0 ( 0 ) I 0 , N ( t 0 + 1 ) = R 0 ( 1 ) I 0 + R 0 ( 1 ) N ( t 0 ) , N ( t 0 + 2 ) = R 0 ( 2 ) I 0 + R 0 ( 2 ) N ( t 0 ) + R 0 ( 1 ) N ( t 0 + 1 ) , N ( t 0 + 3 ) = R 0 ( 3 ) I 0 + R 0 ( 3 ) N ( t 0 ) + R 0 ( 2 ) N ( t 0 + 1 ) + R 0 ( 1 ) N ( t 0 + 2 ) ,
In practice, we can assume that R 0 ( 0 ) = 0 since infected individuals are not infectious immediately after being infected. Under this additional assumption, we obtain the system
N ( t 0 ) = 0 , N ( t 0 + 1 ) = R 0 ( 1 ) I 0 , N ( t 0 + 2 ) = R 0 ( 2 ) I 0 + R 0 ( 1 ) N ( t 0 + 1 ) , N ( t 0 + 3 ) = R 0 ( 3 ) I 0 + R 0 ( 2 ) N ( t 0 + 1 ) + R 0 ( 1 ) N ( t 0 + 2 ) ,
Therefore, (6) can be rewritten as a scalar delay difference equation
N ( t ) = R 0 ( 1 ) N ( t 1 ) + + R 0 ( t t 0 1 ) N ( t 0 + 1 ) + R 0 ( t t 0 ) I 0 , t t 0 .
Assume that the infectious period is n days. That is
R 0 ( a ) = 0 , a n + 1 .
Then by defining t 1 = t 0 + n + 1 , Equation (6) becomes
N ( t ) = d = 1 n R 0 ( d ) N ( t d ) , t t 1 ,
with the initial values
N ( t ) = N 0 ( t ) , t [ t 1 n , t 1 ] .
The goal of this article is to understand how to identify the daily reproduction numbers d 1 , , n R 0 ( d ) in (8) knowing t [ t 1 , t 2 ] N ( t ) on some finite time interval. This problem is particularly important to derive the average dynamic of infection at the level of a single patient.
One of the aims of this paper is to investigate the variations of the daily reproduction number d 1 , , n R 0 ( d ) during the period of contagiousness of infectious individuals.
This is not the case in influenza, as shown in simulated data [3] and in real infected animals, where we observe a U-shaped evolution of their viral load and symptoms as their body temperature during their contagiousness period. From there, it is possible to suspect a U-shaped variation in their ability to emit (aerosol transmission) the virus and, therefore, to contaminate it [4].
After the first asymptomatic period (without contagiousness), the daily reproduction number increases. After one to three days, this number decreases due to the action of the first defense of the innate immune system. But, the virus passes over this first immune defense, and the daily reproduction number increases again before the action of the second adaptive immune system. Then, after two to four days, the second adaptive immune response becomes fully effective. The combination of these biological mechanisms causes the daily reproduction numbers’ U- or M-shaped curve.
The literature about parameters identification for epidemic models with age of infection can be divided into two groups of articles depending on the assumptions made. The first group assumes that d Γ ( d ) is a given function and estimates the time dependent transmission rate t τ ( t ) . As a consequence, they obtain the instantaneous (daily or effective) reproduction number, which is
R 0 ( t ) = τ ( t ) S ( t ) d = 1 n Γ ( d ) .
We refer to [5,6,7,8,9,10,11,12] (and references therein) for more results about this subject.
The second group corresponds to the assumptions considered here. That is, we assume that t τ ( t ) = τ 0 and t S ( t ) = S 0 are constant functions (over a short period of time) and estimate the daily reproduction number. That is the case for the discrete time model in [13] and more recently for the continuous time model in [1]. The major default in [13] is that the estimated d R 0 ( d ) does not remain positive. We will have the same problem is Section 3.1 when we will use the full spectrum. In Section 3.2, to solve this problem, we introduce a method using the dominant and secondary eigenvalue only.
This article aims to investigate the shape of the distribution d R 0 ( d ) from the data of COVID-19. In Figure 1, we illustrate the notion of U- or M-shaped distribution.
The U and M shape distribution are well known in the context of influenza [3,4]. In Figure 2, we present some figures reflecting patients’ viral load for COVID-19.
Such U shape has not yet been systematically studied in COVID-19 data, but observations of the evolution of the viral load have been done in some patients and show this U shape. Figure 2 shows such a U-shaped evolution for the viral load in real patients [14].
The present work is directly connected to the original work of Peter Whittle in 1951 [15,16], who introduced the Auto Regressive Moving Average (ARMA) model, after the seminal paper on time series by N. Wiener [17],
N ( t ) = K ( 1 ) N ( t 1 ) + K ( 2 ) N ( t 2 ) + + K ( n ) N ( t n ) Auto regressive part + w ( t ) Moving average part ,
where N ( t ) is the size at time t of the population whose growth is forecasted, the kernel d K ( d ) has real values, n is the regression order, and here w ( t ) stands for a noise. Equation (10) has been extensively studied under the denomination of ARMA models by many authors [18,19,20,21,22,23,24].
Here, we propose a new approach based on the spectral properties of the population growth equation to capture information from data. Our goal is to estimate the shape of the daily reproduction numbers d R 0 ( d ) . Spectral methods are not new (see Priestley [20,25]), but it usually refers to Fourier transform with frequencies associated to various periods, corresponding to a fundamental period and its sub-multiples (harmonics). If we consider the auto regressive part only, the spectrum of the delay difference equation is determined by its characteristic equation
λ n = K ( 1 ) λ n 1 + K ( 2 ) λ n 2 + + K ( n 1 ) λ + K ( n ) .
The main idea in this article is to use these eigenvalues λ 1 , λ 2 , , λ n C (i.e., the solution of the characteristic equation) to identify the parameters K ( 1 ) , K ( 2 ) , , K ( n ) . The eigenvalues λ 1 , λ 2 , , λ n C are estimated by some separated method. In Section 2, we will see that when all the eigenvalues are non null and separated two by two, then we can compute the parameters K ( 1 ) , K ( 2 ) , , K ( n ) by using the eigenvalues only.
The idea of using eigenvalues in population dynamics goes back to Malthus [26], who, in 1798, first identified in a mixture of populations the one that would impose itself on the others, determined through the exponential growth of the largest exponent—this leading exponent having been called Malthusian parameter by Fisher [27]. The Malthusian growth seeming unrealistic, the saturation logistic term was introduced further by Lambert [28], and then extending the initial work by Euler [29], Lotka [30], Leslie [31], and Hahn [32] gave the current matrix form of the discrete population growth equations.
However, as far as we know, estimating the subdominant eigenvalues to characterize the system is new. So the key idea of this work is to use the dominant eigenvalue λ 1 and also the following pair of complex conjugated eigenvalues λ 2 , λ ¯ 2 as an estimator to reconstruct the kernel of the auto regressive part.
This work is motivated by the times series provided by data of the daily numbers of reported cases of COVID-19. During the COVID-19 pandemic, most people viewed the oscillations around the exponential growth at the beginning of an epidemic wave as the default in reporting the data. The residual is probably partly due to the reporting data process (random noise). Nevertheless, a significant remaining part of such oscillations could be connected to the infection dynamic at the level of a single average patient. Eventually, the central question we try to address here is: Is there some hidden information in the signal around the exponential tendency for COVID-19 data? So we consider the early stage of an epidemic phase, and we try to exploit the oscillations around the tendency in order to reconstruct the infection dynamic at the level of a single average patient.
We start by investigating the connection between a signal decomposed into a sum of damped or amplified oscillations and a renewal equation. The prototype example we have in mind is the following:
N ( t ) = A 1 e α 1 t + e α 2 t A 2 cos ω 2 t + B 2 sin ω 2 t + C , t t 1 n ,
where A 1 , A 2 , A 3 R , α 1 > 0 , α 2 R , and ω 2 > 0 .
In Figure 3, we illustrate a growing function with damped oscillations (i.e., α 2 < 0 ) and amplified oscillations (i.e., α 2 > 0 ). It is clear from Figure 3 that a periodic function can not represent such a signal, and extending such a signal by periodicity would be artificial. Indeed, the Fourier decomposition would only provide purely imaginary eigenvalues that would exclude a continuation of the exponential growth (i.e., eigenvalues with non-zero real parts). To apply wavelets theory (see, for example, in [33]), we need to extend the data for negative times by symmetry with respect to the initial time t = 0 , and we need a decreasing function ( α 1 < 0 and α 2 < 0 ).
Here, we are more interested in the model resulting from the data (i.e., R 0 ( d ) 0 , d = 1 , , n ) than in the fit to the data. The major problem with the Fourier method is that this method provides only eigenvalues with zero real parts (that is due to the periodicity required for this method). Such eigenvalues are well adapted to a periodic signal, but this is not suitable to describe, for example, an ever-growing function (as in Figure 3). Consequently, the Fourier method is not well adapted to derive non-negative daily reproduction numbers (i.e., R 0 ( d ) 0 , d = 1 , , n ).
Previous analogous approaches can be found in the seismic data modeling and statistical literature, like the Wiener–Levinson predictive deconvolution (Robinson [34], Peacock and Treitel [35], Robinson and Treitel [36]), which intends to estimate the minimum phase wavelet in the data, in particular in the case where the relatively weak sampling does not make it possible to affirm the Gaussian character of the errors (Walden and Hosken [37]). If the Gaussian character of the errors can be proven, another similar approach is that of the Geometric Brownian Motion (GBM) processes (Vinod et al. [38]) used, for example, in the analysis of financial data (Ritschel et al. [39]), which are based on the model of the solution of a stochastic differential equation, multiplied by a periodic component with a Gaussian noise.
The structure of this paper is as follows: Section 2 is devoted to the materials and methods. We recall some notions of matrices and spectra. We also present some phenomenological models that will be compared to the data. Section 3 contains the results. We fit the phenomenological models to the cumulative numbers of reported cases in Japan over 10 days and 30 days. We use the eigenvalues derived from the phenomenological model, and we identify the daily reproduction numbers by using: (1) all the spectrum (see Appendix B) and (2) part of the spectrum. The last section of the paper is devoted to the discussion and the conclusion. We present in the Appendices all the mathematical aspects of the paper (see Appendix A, Appendix B, Appendix C and Appendix D).

2. Materials and Methods

2.1. Identification of the Model

The Leslie matrix associated to the difference Equation (8) is
Biology 11 01825 i001
The characteristic equation of (11) is
λ n = d = 1 n R 0 ( d ) λ n d ,
for λ C , which is equivalent to (whenever λ 0 )
1 = d = 1 n R 0 ( d ) λ d .
The complex numbers satisfying the characteristic equation are called the eigenvalues of L.
In Appendix A and Appendix B, we discuss the identification problem of the daily reproduction numbers R 0 ( 1 ) , , R 0 ( n ) by using the eigenvalues of L. The main identification result of Appendix B corresponds to the formula (A3).
Definition 1.
We will say that L is a Markovian Leslie matrix if all the values d [ 1 , n ] R 0 ( d ) are non negative, and
d = 1 n R 0 ( d ) = 1 .

2.2. Phenomenological Model to Fit the Cumulative and the Daily Numbers of Reported Case Data

Due to Lemma A1 below, we propose the following phenomenological model to represent the data
CR ( t ) = CR 1 e λ 1 t + CR 2 e λ 2 t + CR 3 e λ 3 t + + CR m e λ m t ,
where CR 1 , , CR m C are non null, λ 1 = α 1 + i ω 1 , , λ m = α m + i ω m C are pairwise distinct, and m n .
Remark 1.
In the above formula, we allow the constant terms whenever λ n = 0 .
Assuming that the unit of time is one day, we have the following relationship between the cumulative number of cases CR ( t ) and the daily number of cases N ( t )
CR ( t ) = CR ( t 0 ) + t 0 t N σ d σ .
We deduce that the daily number of reported cases has the following form
N ( t ) = N 1 e λ 1 t + N 2 e λ 2 t + N 3 e λ 3 t + + N m e λ m t ,
where N 1 , , N m C are non null, and λ 1 , , λ m C are the same as in (14), and m n .
Since N ( t ) is obtained from CR ( t ) by computing the first derivative, we have the following relationship
N k = CR k λ k , k = 1 , , m .
Remark 2.
For the daily number of cases data t N ( t ) only a few eigenvalues will be tractable. For example, in Section 3.3, we will consider the following extension
N ( t ) = N 1 e λ 1 t + N 2 e λ 2 t + N 3 e λ 3 t + N 3 e λ 3 t + N 4 e λ 4 t + w ( t )
where w ( t ) will contain N 5 e λ 5 t + + N m e λ m t merged together with some random term.
Remark 3.
The identification of the eigenvalues λ 1 , , λ m as parameters of the phenomenological model is discussed in Section 3.3. So far, this problem for a finite time interval seems to be open.
We will first approach the data with the following phenomenological model.
Phenomenological model for the cumulative numbers of reported cases with λ > 0
We start with a first eigenvalue λ = e α > 0 , for some α R . The phenomenological model used to fit the cumulative numbers of reported cases has the following form
CR ( t ) = A e α ( t t 0 ) + C , for t [ t 0 , + ) ,
where A R , α R , and C R are real numbers.
For discrete times, it is equivalent to say that
CR ( n ) = A λ n + C , for n = 0 , 1 , 2 , .

By computing the first derivative of t CR ( t ) , we obtain a model for the daily number of cases of the following form
N ( t ) = A α e α ( t t 0 ) , for t [ t 0 , + ) .
Once the best fit of the above phenomenological model to the data is obtained, we can subtract this model to the data t CR Data ( t ) , then we obtain a first residual
Residual ( t ) = CR Data ( t ) CR ( t ) .
Next we will approach the residual with the following phenomenological model.
Phenomenological model for the cumulative numbers of reported cases with λ C
Assume that the eigenvalues are two conjugated complex numbers λ = e α ± i ω C , for some α R and ω 0 . The phenomenological model used to fit the cumulative numbers of reported cases has the following form
CR ( t ) = e α ( t t 0 ) A cos ω t t 0 + B sin ω t t 0 + C , for t [ t 0 , + ) ,
where α R , A R , B R , C R , and ω 0 are four real numbers.
For discrete times, it is equivalent to say that
CR ( n ) = A i B 2 λ n + A + i B 2 λ ¯ n + C , for n = 0 , 1 , 2 , .

By computing the first derivative of t CR ( t ) , we obtain a model for the daily number of cases of the following form
N ( t ) = e α ( t t 0 ) A ^ cos ω t t 0 + B ^ sin ω t t 0 , for t [ t 0 , + ) ,
where
A ^ = α A + ω B B ^ = ω A + α B A = α A ^ ω B ^ ω 2 + α 2 B = ω A ^ + α B ^ ω 2 + α 2 .
Remark 4.
When ω = 0 in (18), we obtain the previous model (15).

2.3. Cumulative and Daily Number of Reported Cases for COVID-19 in Japan

Here we use cumulative numbers of reported cases for COVID-19 in Japan taken from the WHO [40]. The data show a succession of epidemic waves (blue background color regions) followed by endemic periods (yellow background color regions). In Figure 4, black dots represent the data. The blue background color regions correspond to epidemic phases, and the yellow background color region to endemic phases. The region of interest to apply the method is between 19 October and 29 October 2020. This region is marked with light green vertical lines in the figure.

3. Results

3.1. Methods Applied to Ten Days Data

In this section, we will fit the phenomenological model (15) or (18) to the cumulative numbers of reported cases presented in the previous subsection. We consider a period of 10 days since the beginning of the third epidemic wave of COVID-19 in Japan. The period goes from 19 to 29 October 2020.
Step 1:
In Figure 5, we fit an exponential function (15) to the cumulative number of reported cases of COVID-19 in Japan between 19 and 29 October 2020.
In Figure 5, the best fit of model (15) is obtained for
A 1 = 2 . 8810 4 , C 1 = 6 . 4210 4 , and α 1 = 0.02 .
Hence,
λ 1 = exp ( α 1 ) = 1.02 .
Step 2:
Next, we consider the residual left after the previous fit,
Residual 1 ( t ) = CR ( t ) A 1 e α 1 t + C 1 .
In Figure 6, we fit the model (18) to the first residual function t Residual 1 ( t ) .
In Figure 6, the best fit of model (18) (i.e., minimizing the sum-of-squares error) is obtained for
A 2 = 138.16 , B 2 = 127.36 , C 2 = 11.88 , α 2 = 0.07 , and ω 2 = 0.95 .
The period associated to ω 2 is equal to P 2 = 2 π ω 2 = 6.609 days. This periodic phenomenon was observed in many countries (see for example [41]). Here,
λ 2 = exp ( α 2 + i ω 2 ) = 0.54 + 0.76 i ,
λ 3 = exp ( α 2 i ω 2 ) = 0.54 0.76 i .
By using
M = ( λ 1 1 λ 1 2 λ 1 3 λ 2 1 λ 2 2 λ 2 3 λ 3 1 λ 3 2 λ 3 3 ) ,
in (A3) below, with n = 3 , we obtain
( R 0 ( 1 ) R 0 ( 2 ) R 0 ( 3 ) ) = ( 2.09 1.96 0.87 ) .
Since
det ( M ) = 1.78 i ,
therefore, the components of M 1 are not too large, and the above result should not be too sensitive to the stochastic errors. The main problem in (22) is the second component 1.9625 , which is not making sense in this context.

3.2. Spectral Truncation Method Applied to Ten Days Data

In the previous subsection, the first two fits make perfect sense. However, adding more fits would be questionable because they become more and more random after a few steps. We could alternatively continue to fit the rest by using our phenomenological model, which would provide new eigenvalues.
The major problem in the previous section is that when we apply formula (A3) with all the eigenvalues, we obtain some R 0 ( 1 ) , , R 0 ( n ) with negative values. Instead here, we increase the dimension n of L, and we use only the eigenvalues λ 1 , λ 2 , λ 3 .

3.2.1. Re-Normalizing Procedure

Assume that λ 1 1 then by
N ¯ ( t ) = N ( t ) λ 1 t N ( t ) = λ 1 t N ¯ ( t )
where t N ( t ) is a solution of (8), we obtain the following normalized equation
λ 1 t N ¯ ( t ) = d = 1 n R 0 ( d ) λ 1 t d N ¯ ( t d ) , t t 1 ,
and by dividing the above equation by λ 1 t we obtain
N ¯ ( t ) = d = 1 n R ¯ 0 ( d ) N ¯ ( t d ) , t t 1 .
where
R 0 ( d ) = R ¯ 0 ( d ) λ 1 d , d = 1 , , n .
By using the procedure, we can always fix the dominant eigenvalue of L to 1 by imposing that L is Markovian (see Definition 1). Then we use the following re-normalizing procedure for the eigenvalues
λ 1 🟉 = λ 1 / λ 1 = 1 , λ 2 🟉 = λ 2 / λ 1 = 0.53 + 0.74 i , and λ 3 🟉 = 0.53 0.74 i .
In Figure 7, we fit these eigenvalues λ 2 🟉 and λ 3 🟉 with the spectrum of Markovian Leslie matrices L on a mesh. We observe that the fit improves when the dimension of L increases.
In Figure 8, we observe that, for n 3 , 5 , 6 , there is a unique set of eigenvalues λ 1 , λ 2 , λ 3 , , λ n of L (classified with decreasing real part) minimizing the distance | λ 2 🟉 λ 2 | and | λ 3 🟉 λ 3 | . This is no longer true for n = 7 .

3.2.2. Daily Basic Reproduction Numbers

In Figure 9, we plot the average distribution d R 0 ( d ) , standard deviation (blue region), and 95 % confidence interval.
In Figure 10, we plot the daily basic reproduction numbers R 0 ( d ) .
We can notice that following [42], the effective R 0 is between 1.06 and 1.14 on 19 October 2020, in Japan.

3.2.3. Applying the Model to Daily Number of Reported Cases

The model used to run the simulations is the following
N ( t ) = d = 1 6 R 0 ( d ) N ( t d ) , t t 0 + 6 ,
and according to the formula (17) and (20), with the initial condition
N ( t ) = A 1 ln λ 1 λ 1 t + e α 2 t [ A ^ 2 cos ( ω 2 t ) + B ^ 2 sin ( ω 2 t ) ] , t = t 0 , t 0 + 1 , , t 0 + 5 ,
with
A ^ 2 = α 2 A 2 + ω 2 B 2 and B ^ 2 = ω 2 A 2 + α 2 B 2 .
In (24)–(26) we use the parameter values estimated in Section 3.1.
In Figure 11, we plot the daily number of reported cases data from October 19 to November 19, 2020 (black dots) and from model (24) and (25) with the values of R 0 ( d ) obtained in Figure 10c (red dots).

3.3. Extension of the Spectral Truncation Method over One Month

In Figure 12, we apply respectively the AutoCorrelation Function (ACF) and Partial AutoCorrelation Function (PACF) to the daily number of cases for Japan from 19 October and 19 November 2020. It does not look like any standard cases. In the ACF, we observe the correlation is significant until 7 days, while in the PACF it is until 16 days.
  • Step 1: In Figure 13, we fit the model
    ϕ 1 ( t ) = A 1 e α 1 ( t t 0 ) + C 1 ,
    with the cumulative number of reported cases data between 19 October and 19 November 2020.
Figure 13. In this figure, we plot the cumulative number of reported cases data between 19 October and 19 November 2020 (black dots). We plot the best fit of the model (27) to the cumulative data (red curve).
Figure 13. In this figure, we plot the cumulative number of reported cases data between 19 October and 19 November 2020 (black dots). We plot the best fit of the model (27) to the cumulative data (red curve).
Biology 11 01825 g013
We obtain the following parameter values for the best fit
A 1 = 7 . 9310 3 , C 1 = 8 . 5510 4 , and α 1 = 0.05 .
  • Step 2: Next we define as before the first residual
    Residual 1 ( t ) = CR ( t ) A 1 e α 1 ( t t 0 ) + C 1 ,
    and we fit the Residual 1 ( t ) with the model
    ϕ 2 ( t ) = e α 2 ( t t 0 ) A 2 cos ω 2 t t 0 + B 2 sin ω 2 t t 0 + e α 3 ( t t 0 ) A 3 cos ω 3 t t 0 + B 3 sin ω 3 t t 0 + C 2 .
In Figure 14, we plot the cumulative number of reported cases data between 19 October and 19 November 2020.
The parameters of the phenomenological model ϕ 2 ( t ) obtained for the best fit are the following
A 2 = 55.21 , B 2 = 84.48 , A 3 = 391.57 , B 3 = 88.79 , C 2 = 7.68 ,
and
α 2 = 0.0501 , ω 2 = 0.91 , α 3 = 0.02 , ω 3 = 0.3 .
The periods associated to ω 2 and ω 3 are, respectively,
P 2 = 2 π ω 2 = 6.92 days , and P 3 = 2 π ω 3 = 21.24 days .
These periods are close multiples of 7 days.
Remark 5.
It is important to note that the period P 3 of 21 days is difficult to explain mechanically, but this value is the smallest value giving the best fit to the data. We tried to impose some upper bounds smaller than 21 days. In such a case, P 3 is always replaced by the upper bound. This is true for all constraints less that 21 days, and for each constraint larger than 22 days, we obtain P 3 = 21.24 days.
Remark 6.
It is important to note that α 1 = α 2 . That is because, during the fit, we impose that α 2 α 1 and α 3 α 1 . That is the condition coming from the Perron Frobenius theorem, in order to obtain
| λ 2 | | λ 1 | a n d | λ 3 | | λ 1 | .
This condition is coming from the fact that λ 1 must be the spectral radius of L and λ 2 , λ 3 belong to the circle centered at 0 and with the radius equal to the spectral radius of L (i.e., with a modulus less or equal to λ 1 ).
Eigenvalues associated to the model ϕ 1 ( t ) and ϕ 2 ( t ) : The first eigenvalue is
λ 1 = e α 1 = 1.05 .
The second pair of complex conjugated eigenvalues is
λ 2 = e α 2 cos ( ω 2 ) + i sin ( ω 2 ) = 0.65 + 0.83 i ,
and
λ 3 = e α 2 cos ( ω 2 ) i sin ( ω 2 ) = 0.65 0.83 i ,
and the modulus of λ 2 is
| λ 2 | = | λ 3 | = e α 2 = e α 1 = λ 1 = 1.05 .
The fourth eigenvalue is
λ 4 = e α 3 cos ( ω 3 ) + i sin ( ω 3 ) = 0.94 + 0.29 i ,
and the fifth eigenvalue is its conjugate
λ 5 = e α 3 cos ( ω 3 ) i sin ( ω 3 ) = 0.94 0.29 i ,
and the modulus of λ 4 is
| λ 4 | = | λ 5 | = e α 3 = 0.98 < 1.05 .
Using λ 2 and λ 4 as an estimator: Next we consider all the matrices L in which the component R 0 ( d ) is replaced by R ¯ 0 ( d ) , and we assume that
d = 1 n R ¯ 0 ( d ) = 1 .
The dominant eigenvalue of L is 1, and we look for matrices such that the second eigenvalue of L is close to
λ 2 🟉 = λ 2 / λ 1 ,
and the fourth eigenvalue of L is close to
λ 4 🟉 = λ 4 / λ 1 .
For realizing this approach, we minimize the
χ L = max d λ 2 🟉 , σ ( L ) , d λ 4 🟉 , σ ( L )
where
d λ 2 🟉 , σ ( L ) = min λ σ ( L ) | λ 2 🟉 λ | , and d λ 4 🟉 , σ ( L ) = min λ σ ( L ) | λ 4 🟉 λ | ,
where σ ( L ) is the set of all eigenvalues of L.
In Figure 15, we consider the d R ¯ 0 ( d ) such that the corresponding maximum satisfies
χ L ( R ¯ 0 ) inf R ^ 0 0 : R ^ 0 ( d ) = 1 χ L ( R ^ 0 ) + 10 2 .
We define
R 0 ( d ) = R ¯ 0 ( d ) λ 1 d , d = 1 , , n .
In Figure 16, we obtain a good description of the dynamic of infection at the individual level that confirms the one obtained over shorter periods. As expected, the average patient first loses its ability to transmit the pathogen, and after decreasing by day 1 to day 4, R 0 ( d ) increases between day 4 and day 7. Day 7 is a maximum. After day 7, R 0 ( d ) decays until day 9. Then a second peak arises, with a maximum on the day 14. We could explain this second peak by supposing that an important transmission of pathogen still exists from day 12 to day 16. We also obtain a third from day 19 to 23 with a maximum value on day 21.
In Figure 17, we plot the spectrum of the Leslie matrix L when d R ¯ 0 ( d ) corresponds to the average distribution (i.e., the red curve in Figure 15).
Recalling that, by definition, the basic reproduction number is
R 0 = d = 1 n R 0 ( d ) ,
we obtain the sum of the daily reproduction numbers (red curve in Figure 16)
R 0 = 2.13 .
In Figure 18, we plot a histogram for the values of the basic reproduction number obtained by summing the distributions d R 0 ( d ) from Figure 16.
Next, we consider
N ( t ) = d = 1 25 R 0 ( d ) N ( t d ) , t t 0 + 25 ,
and accordingly to the formula (17) and (20), with the initial condition for t = t 0 , t 0 + 1 , , t 0 + 25 , we have
N ( t ) = A 1 ln λ 1 λ 1 t + e α 2 t [ A ^ 2 cos ( ω 2 t ) + B ^ 2 sin ( ω 2 t ) ] + e α 3 t [ A ^ 3 cos ( ω 3 t ) + B ^ 3 sin ( ω 3 t ) ] ,
with
A ^ 2 = α 2 A 2 + ω 2 B 2 , B ^ 2 = ω 2 A 2 + α 2 B 2 , A ^ 3 = α 3 A 3 + ω 3 B 3 and B ^ 3 = ω 3 A 3 + α 3 B 3 .
In (24)–(26) we use the parameter values estimated in Section 3.1.
In Figure 19, we see the mean distribution d R 0 ( d ) permits to produce oscillations around the tendency for the daily number of cases. It is important to note that without the third peak in Figure 16 we do not obtain such a good correspondence between the model and the data.

4. Discussion

In this article, we start by investigating the connection between a signal decomposed into a sum of damped or amplified oscillations and a renewal equation. Namely, we connect the daily number of reported cases written as
N ( t ) = N 1 e α 1 t cos ω 1 t + i sin ω 1 t + + N n e α n t cos ω n t + i sin ω n t , t t 1 n ,
with the renewal equation
N ( t ) = d = 1 n R 0 ( d ) N ( t d ) , t t 1 .
In the context of epidemic time series, a spectral method usually refers to the Fourier decomposition of a periodic signal. In the present paper, the data are not periodic and are composed of an exponential function (Malthusian growth) perturbed with some damped oscillating functions. So we use complex numbers with non-null real parts. We refer to Cazelles et al. [33] for more results about time series.

4.1. Data over Ten Days

We can notice in Figure 9 and Figure 10 and Table 1 that the daily reproduction number as well as the instantaneous reproduction number are estimated. Concerning the instantaneous (or effective) reproduction number R e ( t ) [43,44] estimated by [42], which equals 1.1 at the 19th of October 2020, the best fit corresponds to n = 7 days (see (c) in Figure 9). This value of the duration of the contagiousness period is close to the values 6 or 7 days and are close to the values estimated from the virulence measured in [14,45,46]. In Figure 10, we always obtain a U -shaped distribution for the curve of daily reproduction numbers. This corresponds to the biphasic form of the virulence already observed in respiratory viruses, such as influenza, as recalled in the Introduction.
This temporal behavior of the contagiousness can correspond to the evolution of contagious symptoms like cough or spitting, which diminish during the innate immune response, followed by a comeback of the symptoms before the adaptive immune response (whenever the innate defense has been overcome by the virus). If the innate cellular immunity has been not sufficient for eliminating the virus, the viral load again increases, causing a reappearance of the symptoms before the adaptive immunity (cellular and humoral) occurs, which results in a transient decrease in contagiousness between the two immunologic phases. The medical recommendations are, in the case of U-shaped contagiousness, never to take a transient improvement for a permanent disappearance of the symptoms and to stay at home to avoid a bacterial secondary infection that is possibly fatal.
The estimation of the daily reproduction numbers in the COVID-19 outbreak constitutes an important issue. At the public health level, to publish only the sum of the daily reproduction numbers, that is, to say the basic reproduction number R 0 or the effective reproduction number Re, could suffice for controlling and managing the behavior of a whole population with mitigation or vaccination measures. At the individual level, it is important to know the existence of a minimum of the daily reproduction numbers, which generally corresponds to a temporary clinical improvement, after a partial success of the innate immune defense. This makes it possible to advise the patient to continue to respect his own isolation, prevention, and therapy choices (depending on his vaccination state) even if this transient clinical improvement has occurred. The present methodology allows also to estimate both the individual contagiousness duration in a dedicated age class and also its seasonal variations, which is crucial for optimizing the benefit–risk decisions of the public and individual health policies.

4.2. Data over One Month

Over one month, we obtain a daily reproduction number with three peaks. Each peak is centered respectively on 7 days, 14 days, and 21 days. These quantities coincide with the period of 7 days and 21 days obtained in Figure 14 in fitting the first residual when we subtract the exponential growth first fit to the cumulative data. As far as we understand the problem, that is the period of 21 days in the data, which induces the third peak. This third peak is very suspicious. Nevertheless, the data lead us to such a shape for the daily reproductive number. We also tried to run Figure 19 without the third peak, and we obtained a bad fit to the data, while with this third peak, the fit is good. One may also note that the 21-day period is insignificant for the ACF and the PACF in Figure 12.
Several possibilities exist to explain this strange shape for the daily reproduction number using the data over one month. One possible explanation is that the Japanese population should be subdivided into several groups having very different infection dynamics (at the level of a single patient). Here we have in mind the patient with a short infection period but high transmissibility (super spreaders) versus the patient with a long infection period with mild symptoms.
We suspect that such a shape for the daily reproduction number could be attributed to the time since infection to report a case. The daily number of reported cases would be obtained from N ( t ) , and the daily number of new infected cases by using the following model
D ( t ) = f d = 1 q K ( d ) N ( t d ) ,
where the integer q 1 is the maximum number of days needed to report a case, f [ 0 , 1 ] is the fraction reported, and K ( d ) 0 is the probability to report a case after d days. Therefore, we must have
d = 1 q K ( d ) = 1 .

4.3. Perspectives and Conclusions

In the present paper, we only consider the Japanese data in the exponential phase of the third epidemic wave.
The case of Japan seems emblematic to us, as it corresponds to a wave of well-identified new cases following a clearly characterized endemic phase. The exponential growth phenomenon being transitory, this explains the relatively limited duration of the sampling, which corresponds to a period in days during which the epidemiological parameters (such as the transmission rate) can be considered as constant. It is in such circumstances where the Gaussian nature of the errors is difficult to prove, due to the small sampling, such that similar methods based on wavelets have been proposed (Walden and Hosken [37]).
The method of the present paper should be applied to several countries for each epidemic wave to obtain a more systematic study. For the moment, over one month, we obtained a shape for the daily reproduction number that follows the data very well. However, we are suspicious about the third peak. We suspect that the default of our analysis is coming from the model itself. Such a question has been recently studied by Ioannidis and his collaborators in [47], and we believe that we are facing such modeling difficulties.

Author Contributions

Conceptualization, J.D. and P.M.; methodology, P.M.; software, P.M.; writing—original draft preparation, J.D and P.M.; writing—review and editing, J.D and P.M.; All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

There is no subject involved in the present study.

Data Availability Statement

No data were produced for this study.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. Non Identifiability Result

From Formula (13), we deduce that the characteristic (12) has exactly one positive solution. By the Perron–Frobenius theorem applied to the Leslie matrix L defined by (11), we know that (by considering the norm of linear operator) the spectral radius of L
r ( L ) : = lim n + L n L ( R ) 1 / n > 0 ,
is the unique positive solution of (12). Moreover, all the remaining eigenvalues have a modulus smaller or equal to r ( L ) . We refer to ([48], Chapter 4), for more results about this subject.
Non identifiability result: Let λ 🟉 > 0 and N 🟉 0 . Then
N ( t ) = N 🟉 λ 🟉 t t 1 , t t 1 ,
is a known solution of (8) if and only if λ 🟉 is a solution of the characteristic equation.
Assume that d [ 1 , n ] R 🟉 ( d ) 0 is given, and satisfies
d = 1 n R 🟉 ( d ) > 0 .
Then if we define
R 0 ( a ) = R 🟉 ( a ) d = 1 n R 🟉 ( d ) λ 🟉 d , a = 1 , , n ,
we deduce that the equation (12) is satisfied for λ = λ 🟉 , and N ( t ) = N 🟉 λ 🟉 t t 1 is a solution of (8). We conclude that a single function N ( t ) = N 🟉 λ 🟉 t t 1 is not enough to identify R 0 ( 1 ) , R 0 ( 2 ) , R 0 ( 3 ) , , R 0 ( n ) .

Appendix B. Identifiability Result

Assumption A1.
Assume that λ 1 , , λ n C are nonzero complex numbers, and are separated two by two. That is,
λ i 0 , i = 1 , , n .
and
λ i λ j , w h e n e v e r i j .
Remark A1.
Since the coefficients of the characteristic Equation (12) are all real, we could also impose that the conjugate of each eigenvalue belongs to the spectrum. That is
λ ¯ i λ 1 , , λ n , i = 1 , , n .
However, that is not necessary in this subsection.
Remark A2.
When all the eigenvalues are real, the above assumption will be satisfied if and only if λ 1 , , λ n R are nonzero real numbers which are pairwise distinct. Up to a permutation, that is
λ i 0 , i = 1 , , n ,
and
λ 1 < λ 2 < < λ n .
Lemma A1.
Let Assumption A1 be satisfied. Assume that each λ i C satisfies the characteristic Equation (12). Then the Leslie matrix L defined by (11) is diagonalizable (and invertible); moreover, for each U 1 , U 2 , , U n C ,
U ( t ) = U 1 λ 1 t + U 2 λ 2 t + + U n λ n t , t t 1 n ,
is a solution of (8). That is to say,
U ( t ) = d = 1 n R 0 ( d ) U ( t d ) , t t 1 .
Identification of the components U i from the values of t N ( t ) : Assume that the values of N ( t ) are given for t = t 1 , , t 1 + n 1 . We claim that we can compute U 1 , U 2 , U 3 , , U n C . Indeed,
N ( t 1 ) = U 1 λ 1 t 1 + U 2 λ 2 t 1 + + U n λ n t 1 , N ( t 1 + 1 ) = U 1 λ 1 t 1 + 1 + U 2 λ 2 t 1 + 1 + + U n λ n t 1 + 1 , N ( t 1 + n 1 ) = U 1 λ 1 t 1 + n 1 + U 2 λ 2 t 1 + n 1 + + U n λ n t 1 + n 1 ,
can be rewritten as the system
( N ( t 1 ) N ( t 1 + 1 ) N ( t 1 + n 1 ) ) = ( λ 1 t 1 λ 2 t 1 λ n t 1 λ 1 t 1 + 1 λ 2 t 1 + 1 λ n t 1 + 1 λ 1 t 1 + n 1 λ 2 t 1 + n 1 λ n t 1 + n 1 ) ( U 1 U 2 U n ) .
The determinant of the above Vandermonde-like matrix
det ( λ 1 t 1 λ 2 t 1 λ n t 1 λ 1 t 1 + 1 λ 2 t 1 + 1 λ n t 1 + 1 λ 1 t 1 + n 1 λ 2 t 1 + n 1 λ n t 1 + n 1 )
= λ 1 t 1 λ 2 t 2 λ n t n det ( 1 1 1 λ 1 1 λ 2 1 λ n 1 λ 1 ( n 1 ) λ 2 ( n 1 ) λ n ( n 1 ) ) ,
therefore,
det ( λ 1 t 1 λ 2 t 1 λ n t 1 λ 1 t 1 + 1 λ 2 t 1 + 1 λ n t 1 + 1 λ 1 t 1 + n 1 λ 2 t 1 + n 1 λ n t 1 + n 1 ) = λ 1 1 λ 2 1 λ n 1 1 i < j n λ i 1 λ j 1 .
Therefore, under Assumption A1, this determinant is non null, and we obtain the following result.
Proposition A1.
Let Assumption A1 be satisfied. Then we can compute the components U 1 , , U n in function of the given elements of the trajectory N ( t 1 ) , , N ( t 1 + n 1 ) by solving the linear system (A1), and
( U 1 U 2 U n ) = ( λ 1 t 1 λ 2 t 1 λ n t 1 λ 1 t 1 + 1 λ 2 t 1 + 1 λ n t 1 + 1 λ 1 t 1 + n 1 λ 2 t 1 + n 1 λ n t 1 + n 1 ) 1 ( N ( t 1 ) N ( t 1 + 1 ) N ( t 1 + n 1 ) ) .
Identification of the component R 0 ( d ) from the λ i : By assuming that each λ i is a solution of the characteristic Equation (12), we obtain
1 = R 0 ( 1 ) λ 1 1 + R 0 ( 2 ) λ 1 2 + + R 0 ( n ) λ 1 n , 1 = R 0 ( 1 ) λ 2 1 + R 0 ( 2 ) λ 2 2 + + R 0 ( n ) λ 2 n , 1 = R 0 ( 1 ) λ n 1 + R 0 ( 2 ) λ n 2 + + R 0 ( n ) λ n n ,
which rewrites in the matrix form as
( 1 1 1 ) = ( λ 1 1 λ 1 2 λ 1 n λ 2 1 λ 2 2 λ 2 n λ n 1 λ n 2 λ n n ) ( R 0 ( 1 ) R 0 ( 2 ) R 0 ( n ) ) .
Under Assumption A1 the Vandermonde-like matrix
( λ 1 1 λ 1 2 λ 1 n λ 2 1 λ 2 2 λ 2 n λ n 1 λ n 2 λ n n )
is invertible, because
det ( λ 1 1 λ 1 2 λ 1 n λ 2 1 λ 2 2 λ 2 n λ n 1 λ n 2 λ n n ) = λ 1 1 λ 2 1 λ n 1 det ( 1 λ 1 1 λ 1 ( n 1 ) 1 λ 2 1 λ 2 ( n 1 ) 1 λ n 1 λ n ( n 1 ) )
hence
det ( λ 1 1 λ 1 2 λ 1 n λ 2 1 λ 2 2 λ 2 n λ n 1 λ n 2 λ n n ) = λ 1 1 λ 2 1 λ n 1 1 i < j n λ i 1 λ j 1 0 .
Therefore, we can compute the component of the map d [ 1 , n ] R 0 ( d ) by solving a linear system involving the eigenvalues of the characteristic equation.
Theorem A1.
Let Assumption A1 be satisfied. Then the following properties are equivalent
(i)
The set λ 1 , , λ n is the spectrum of the Leslie matrix L defined in (11).
(ii)
Each element of λ 1 , , λ n satisfies (A2).
(iii)
The elements λ 1 , , λ n satisfy
( λ 1 1 λ 1 2 λ 1 n λ 2 1 λ 2 2 λ 2 n λ n 1 λ n 2 λ n n ) 1 ( 1 1 1 ) = ( R 0 ( 1 ) R 0 ( 2 ) R 0 ( n ) ) .
In Figure A1, we plot all the spectrum’s location for Markovian Leslie matrices on a mesh. We can observe the changes of location of the spectrum depending of the dimension n. It seems that the spectrum is fielding more and more the unit circle in C when the dimension increases. We refer to Kirkland [49] for more results going in that direction.
Figure A1. We plot all the spectrum’s locations for Markovian Leslie matrices on a mesh whenever n = 3 in (a), n = 4 in (b), n = 5 in (c), and n = 7 in (d). Here the dominant eigenvalue is always 1, and we can see the corresponding isolated blue dot. The blue region corresponds to the spectrum of Markovian Leslie matrices whenever R 0 ( n ) = 0 . The red region corresponds to the spectrum of Markovian Leslie matrices whenever R 0 ( n ) > 0 .
Figure A1. We plot all the spectrum’s locations for Markovian Leslie matrices on a mesh whenever n = 3 in (a), n = 4 in (b), n = 5 in (c), and n = 7 in (d). Here the dominant eigenvalue is always 1, and we can see the corresponding isolated blue dot. The blue region corresponds to the spectrum of Markovian Leslie matrices whenever R 0 ( n ) = 0 . The red region corresponds to the spectrum of Markovian Leslie matrices whenever R 0 ( n ) > 0 .
Biology 11 01825 g0a1
Continuous dependency of the component R 0 ( d ) with respect to the λ i : Define the set Ω C n of all the elements Λ = λ 1 🟉 , , λ n 🟉 C n satisfying Assumption A1. For each Λ = λ 1 🟉 , , λ n 🟉 Ω , we define
M ( Λ ) = ( λ 1 1 λ 1 2 λ 1 n λ 2 1 λ 2 2 λ 2 n λ n 1 λ n 2 λ n n ) , Λ = λ 1 , , λ n Ω .
Theorem A2.
Consider a sequence Λ m = λ 1 m , , λ n m m 0 Ω , and a point Λ 🟉 = λ 1 🟉 , , λ n 🟉 Ω (i.e., all satisfying Assumption A1). Assume that
lim m + Λ m = Λ 🟉 ,
then
lim m + R 0 m = R 0 🟉 ,
where
R 0 m = M ( Λ m ) 1 ( 1 1 1 ) , m N , and R 0 🟉 = M ( Λ 🟉 ) 1 ( 1 1 1 ) .
Proof. 
We have
( 1 1 1 ) = M ( Λ m ) R 0 m , n N , and ( 1 1 1 ) = M ( Λ 🟉 ) R 0 🟉 .
Subtracting the two above quantities, we obtain
0 = M ( Λ m ) R 0 m M ( Λ 🟉 ) R 0 🟉 ,
which is also equivalent to
0 = M ( Λ m ) R 0 m M ( Λ 🟉 ) R 0 🟉 R 0 m M ( Λ 🟉 ) R 0 m ,
hence,
R 0 🟉 R 0 m = M ( Λ 🟉 ) 1 M ( Λ m ) M ( Λ 🟉 ) R 0 m .
Setting
L m = M ( Λ 🟉 ) 1 M ( Λ m ) M ( Λ 🟉 )
and so
lim m + L m = 0 M n ( C ) .
Now since
R 0 🟉 R 0 m = L m R 0 🟉 L m R 0 🟉 R 0 m ,
we deduce that
R 0 🟉 R 0 m L m L ( C n ) R 0 🟉 + L m L ( C n ) R 0 🟉 R 0 m ,
hence, for all m 1 large enough (i.e., satisfying L m L ( C n ) < 1 )
R 0 🟉 R 0 m L m L ( C n ) 1 L m L ( C n ) R 0 🟉 ,
and the proof is completed. □

Appendix C. Identification of the Phenomenological Model

Here we assume that the daily number of reported cases has the following form
N ( t ) = N 1 e λ 1 t + N 2 e λ 2 t + N 3 e λ 3 t + + N m e λ m t ,
where N 1 , , N n C are non null, and λ 1 , , λ n C are pairwise distinct.
If we assume to know t N ( t ) for all positive integer values t = 0 , 1 , 2 , , then we can compute the discrete Laplace transform
L N λ = t = 0 e λ t N ( t ) ,
which is well defined for all λ C such that
Re λ > max i = 1 , , n Re λ i .
By using (A5), we obtain
L N λ = p = 1 m N p 1 e λ p λ ,
whenever Re λ > max i = 1 , , n Re λ i .
Let k 1 , , m be an integer such that
Re λ k = max i = 1 , , n Re λ i ,
we obtain
lim λ λ k Re λ > Re λ k | L N λ | = + .
The Laplace transform could be used to identify the unknown parameters λ 1 , , λ m in (A5). Then by combining this idea with linear regression of t e λ k t , we could identify the parameters N k , then step by step compute all the parameters of N ( t ) in (A5).
In practice, we only know t N ( t ) on a finite time interval t = 0 , 1 , 2 , , L . In that case, we can define the truncated Laplace transform as
L N λ = t = 0 L e λ t N ( t )
and we have by (27)
L N λ = p = 1 m N p 1 e λ p λ L + 1 1 e λ p λ .
The Laplace transform does not permit to detect the eigenvalues λ k (we tested without success some examples with values of complex numbers coming from the present article). Identification of the eigenvalues λ k , whenever t N ( t ) is known only on a finite time interval, seems to be an open intriguing question.

Appendix D. About Residual 2 (t) in Section 3.3

In Figure A2, we observe that average of Residual 2 ( t ) = Residual 1 ( t ) ϕ 2 ( t ) is close to 0, but its histogram does not have the shape of a normal distribution. So, there might be some residual information in Residual 2 ( t ) .
Figure A2. In this figure, we plot Residual 2 ( t ) .
Figure A2. In this figure, we plot Residual 2 ( t ) .
Biology 11 01825 g0a2

References

  1. Demongeot, J.; Griette, Q.; Maday, Y.; Magal, P. Kermack-McKendrick model with age of infection starting from a single or multiple cohorts of infected patients. arXiv 2022, arXiv:2205.15634. [Google Scholar]
  2. Kermack, W.O.; McKendrick, A.G. Contributions to the mathematical theory of epidemics: II. Proc. R. Soc. Lond. Ser. B 1932, 138, 55–83. [Google Scholar]
  3. Chao, D.L.; Halloran, M.E.; Obenchain, V.J.; Longini, I.M., Jr. FluTE, a publicly available stochastic influenza epidemic simulation model. PLoS Comput. Biol. 2010, 6, e1000656. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  4. Itoh, Y.; Shichinohe, S.; Nakayama, M.; Igarashi, M.; Ishii, A.; Ishigaki, H.; Ishida, H.; Kitagawa, N.; Sasamura, T.; Shiohara, M.; et al. Emergence of H7N9 Influenza A Virus Resistant to Neuraminidase Inhibitors in Nonhuman Primates. Antimicrob. Agents Chemother. 2015, 59, 4962–4973. [Google Scholar] [PubMed] [Green Version]
  5. Alvarez, L.; Colom, M.; Morel, J.D.; Morel, J.M. Computing the daily reproduction number of COVID-19 by inverting the renewal equation using a variational technique. Proc. Natl. Acad. Sci. USA 2021, 118, e2105112118. [Google Scholar] [PubMed]
  6. Alvarez, L.; Morel, J.-D.; Morel, J.-M. Modeling COVID-19 incidence by the renewal equation after removal of administrative bias and noise. Biology 2022, 11, 540. [Google Scholar]
  7. Demongeot, J.; Griette, Q.; Magal, P. SI epidemic model applied to COVID-19 data in mainland China. R. Soc. Open Sci. 2020, 7, 201878. [Google Scholar] [CrossRef]
  8. Griette, Q.; Demongeot, J.; Magal, P. What can we learn from COVID-19 data by using epidemic models with unidentified infectious cases? Math. Biosci. Eng. 2021, 19, 537–594. [Google Scholar] [CrossRef]
  9. Griette, Q.; Demongeot, J.; Magal, P. A robust phenomenological approach to investigate COVID-19 data for France. Math. Appl. Sci. Eng. 2021, 2, 149–218. [Google Scholar]
  10. Nishiura, H. Time variations in the transmissibility of pandemic influenza in Prussia, Germany, from 1918–19. Theor. Biol. Med. Model. 2007, 4, 20. [Google Scholar] [CrossRef] [Green Version]
  11. Nishiura, H.; Chowell, G. The Effective Reproduction Number as a Prelude to Statistical Estimation of Time-Dependent Epidemic Trends. In Mathematical and Statistical Estimation Approaches; Epidemiology, G., Chowell, J.M., Hyman, L.M., Bettencourt, A., Castillo-Chavez, C., Eds.; Springer: Dordrecht, The Netherlands, 2009; pp. 103–121. [Google Scholar]
  12. Bakhta, A.; Boiveau, T.; Maday, Y.; Mula, O. Epidemiological forecasting with model reduction of compartmental models. application to the COVID-19 pandemic. Biology 2020, 10, 22. [Google Scholar] [CrossRef] [PubMed]
  13. Waku, J.; Oshinubi, K.; Demongeot, J. Maximal reproduction number estimation and identification of transmission rate from the first inflection point of new infectious cases waves: COVID-19 outbreak example. Math. Comput. Simul. 2022, 198, 47–64. [Google Scholar] [CrossRef] [PubMed]
  14. Pan, Y.; Zhang, D.; Yang, P.; Poon, L.L.M.; Wang, Q. Viral load of SARS-CoV-2 in clinical samples. Lancet Infect. Dis. 2020, 20, 411–412. [Google Scholar] [CrossRef] [PubMed]
  15. Whittle, P. Hypothesis Testing in Time Series Analysis. Almquist Wicksell 1951.
  16. Whittle, P. Prediction and Regulation; English Universities Press: London, UK, 1963. [Google Scholar]
  17. Wiener, N. Extrapolation, Interpolation, and Smoothing of Stationary Time Series; MIT Press: Cambridge, MA, USA, 1949. [Google Scholar]
  18. Chan, K.S.; Tong, H. A note on certain integral equations associated with non-linear time series analysis. Probab. Th. Rel. Fields 1986, 73, 153–158. [Google Scholar]
  19. Lim, K.S.; Tong, H. A statistical approach to difference-delay equation modelling in ecology—Two case studies. J. Time Ser. Anal. 1983, 4, 239–267. [Google Scholar] [CrossRef]
  20. Priestley, M.B. Spectral Analysis and Time Series; Academic Press: Cambridge, MA, USA, 1981. [Google Scholar]
  21. Ramsay, J.O. Monotone Regression Splines in Action. Stat. Sci. 1988, 3, 425–441. [Google Scholar] [CrossRef]
  22. Ramsay, J.; Hooker, G. Dynamic Data Analysis: Modeling Data with Differential Equations; Springer: New York, NY, USA, 2017. [Google Scholar]
  23. Tong, H. Non-Linear Time Series: A Dynamical System Approach; Oxford University Press: Oxford, UK, 1990. [Google Scholar]
  24. Tuan, P.D. The estimation of parameters for autoregressive moving average models. J. Time Ser. Anal. 1984, 5, 53–68. [Google Scholar]
  25. Priestley, M.B. Evolutionary spectra and non-stationary processes. J. R. Stat. Soc. Ser. 1965, 27, 204–229. [Google Scholar] [CrossRef]
  26. Malthus, T.R. An Essay on the Principle of Population as It Affects the Future Improvement of Society, with Remarks on the Speculations of Mr. Godwin, M. Condorcet, and Other Writers; J. Johnson: London, UK, 1798. [Google Scholar]
  27. Fisher, R.A. The Wave of Advance of Advantageous Genes. Ann. Eugen. 1937, 7, 353–369. [Google Scholar] [CrossRef] [Green Version]
  28. Lambert, J.H. Beytrage Zum Gebrauche Der Mathematik Und Deren Anwendung; Verlage des Buchladens der Realschule: Berlin, Germany, 1765; p. 72. [Google Scholar]
  29. Euler, L. Recherches générales sur la mortalité et la multiplication du genre humain. MéMoires L’AcadéMie Des Sci. Berl. 1767, 16, 144–164. [Google Scholar]
  30. Lotka, A.J. Relation between birth rates and death rates. Science 1907, 26, 121–130. [Google Scholar] [CrossRef]
  31. Leslie, P.H. On the use of matrices in certain population mathematics. Biometrika 1945, 33, 183–212. [Google Scholar] [CrossRef] [PubMed]
  32. Hahn, G.M. Mammalian cell populations. Math. Biosci. 1970, 6, 295–315. [Google Scholar] [CrossRef]
  33. Cazelles, B.; Chavez, M.; Magny, G.C.D.; Guégan, J.F.; Hales, S. Time-dependent spectral analysis of epidemiological time-series with wavelets. J. R. Soc. Interface 2007, 4, 625–636. [Google Scholar] [CrossRef] [Green Version]
  34. Robinson, E.A. Predictive deconvolution of time series with application to seismic exploration. Geophysics 1967, 32, 418–484. [Google Scholar] [CrossRef]
  35. Peacock, K.L.; Treitel, S. Predictive deconvolution: Theory and practice. Geophysics 1969, 34, 155–169. [Google Scholar] [CrossRef]
  36. Robinson, E.A.; Treitel, S. Geophysical Signal Analysis; Prentice-Hill, Inc.: Englewood Cliffs, NJ, USA, 1980. [Google Scholar]
  37. Walden, A.T.; Hosken, J.W.J. The nature of non-Gaussianity of primary reflection coefficients and its significance for deconvolution. Geophys. Prosp. 1986, 34, 1038–1066. [Google Scholar] [CrossRef]
  38. Vinod, D.; Cherstvy, A.G.; Wang, W.; Metzler, R.; Sokolov, I.M. Nonergodicity of reset geometric Brownian motion. Phys. Rev. E 2022, 105, L012106. [Google Scholar] [CrossRef]
  39. Ritschel, S.; Cherstvy, A.G.; Metzler, R. Universality of delay-time averages for financial time series: Analytical results, computer simulations, and analysis of historical stock-market prices. J. Phys. Complex. 2021, 2, 045003. [Google Scholar] [CrossRef]
  40. Data from WHO. Available online: https://COVID19.who.int/WHO-COVID-19-global-data.csv (accessed on 20 July 2022).
  41. Demongeot, J.; Oshinubi, K.; Rachdi, M.; Seligmann, H.; Thuderoz, F.; Waku, J. Estimation of Daily Reproduction Numbers during the COVID-19 Outbreak. Computation 2021, 9, 109. [Google Scholar] [CrossRef]
  42. Powered by the Institute of Global Health, Faculty of Medicine, University of Geneva and the Swiss Data Science Center, ETH Zürich-EPFL. Available online: https://renkulab.shinyapps.io/COVID-19-Epidemic-Forecasting/_w_850fb011/?tab=jhu_pred&country=Japan (accessed on 20 July 2022).
  43. Cori, A.; Ferguson, N.M.; Fraser, C.; Cauchemez, S. A new framework and software to estimate time-varying reproduction numbers during epidemics. Am. J. Epidemiol. 2013, 178, 1505–1512. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  44. Scire, J.; Nadeau, S.; Vaughan, T.; Brupbacher, G.; Fuchs, S.; Sommer, J.; Koch, K.N.; Misteli, R.; Mundorff, L.; Götz, T.; et al. Reproductive number of the COVID-19 epidemic in Switzerland with a focus on the Cantons of Basel-Stadt and Basel-Landschaft. Swiss Med. Wkly. 2020, 150, w20271. [Google Scholar] [PubMed]
  45. Kawasuji, H.; Takegoshi, Y.; Kaneda, M.; Ueno, A.; Miyajima, Y.; Kawago, K.; Fukui, Y.; Yoshida, Y.; Kimura, M.; Yamada, H.; et al. Transmissibility of COVID-19 depends on the viral load around onset in adult and symptomatic patients. PLoS ONE 2020, 15, e0243597. [Google Scholar] [CrossRef]
  46. Kim, S.E.; Jeong, H.S.; Yu, Y.; Shin, S.U.; Kim, S.; Oh, T.H.; Kim, U.J.; Kang, S.J.; Jang, H.C.; Jung, S.I.; et al. Viral kinetics of SARS-CoV-2 in asymptomatic carriers and presymptomatic patients. Int. J. Infect. Dis. 2020, 95, 441–443. [Google Scholar] [CrossRef]
  47. Ioannidis, J.P.; Cripps, S.; Tanner, M.A. Forecasting for COVID-19 has failed. Int. J. Forecast. 2022, 38, 423–438. [Google Scholar] [CrossRef]
  48. Ducrot, A.; Griette, Q.; Liu, Z.; Magal, P. Differential Equations and Population Dynamics I: Introductory Approaches; Springer Nature: Berlin, Germany, 2022. [Google Scholar]
  49. Kirkland, S. On the spectrum of a Leslie matrix with a near-periodic fecundity pattern. Linear Algebra Its Appl. 1993, 178, 261–279. [Google Scholar] [CrossRef]
Figure 1. In this figure, we illustrate the notion of U shape distribution in (a) and M shape distribution in (b). Recall that R 0 ( d ) represents the ability of patients to transmit the pathogen after d days since they were infected. The U shape or M shape distribution means that patients can transmit the pathogen since the beginning of their infection. Then they become less infectious in the middle of the infected period. Finally, they become infectious again at the end of the infected period. The only difference between U and M shape distribution is to include days 0 and 8 and R 0 ( 0 ) = R 0 ( 8 ) = 0 in the plot.
Figure 1. In this figure, we illustrate the notion of U shape distribution in (a) and M shape distribution in (b). Recall that R 0 ( d ) represents the ability of patients to transmit the pathogen after d days since they were infected. The U shape or M shape distribution means that patients can transmit the pathogen since the beginning of their infection. Then they become less infectious in the middle of the infected period. Finally, they become infectious again at the end of the infected period. The only difference between U and M shape distribution is to include days 0 and 8 and R 0 ( 0 ) = R 0 ( 8 ) = 0 in the plot.
Biology 11 01825 g001
Figure 2. Viral load in COVID-19 real patients [14]. In (a), the red curve corresponds to the throat swab and the blue curve corresponds to the sputum. In (b), the curves correspond to several patients (A), (B), and (C).
Figure 2. Viral load in COVID-19 real patients [14]. In (a), the red curve corresponds to the throat swab and the blue curve corresponds to the sputum. In (b), the curves correspond to several patients (A), (B), and (C).
Biology 11 01825 g002
Figure 3. We plot an exponentially growing function with (a) damped oscillations and (b) amplified oscillations.
Figure 3. We plot an exponentially growing function with (a) damped oscillations and (b) amplified oscillations.
Biology 11 01825 g003
Figure 4. In this figure, we plot the daily number of reported cases for COVID-19 in Japan.
Figure 4. In this figure, we plot the daily number of reported cases for COVID-19 in Japan.
Biology 11 01825 g004aBiology 11 01825 g004b
Figure 5. In this figure, the black dots correspond to the cumulative numbers of reported cases of COVID-19 in Japan between 19 October and 29 October 2020 (black dots). The red curve corresponds to the best fit of model (15) to the cumulative numbers of reported cases.
Figure 5. In this figure, the black dots correspond to the cumulative numbers of reported cases of COVID-19 in Japan between 19 October and 29 October 2020 (black dots). The red curve corresponds to the best fit of model (15) to the cumulative numbers of reported cases.
Biology 11 01825 g005
Figure 6. In this figure, the black dots correspond to the function t Residual 1 ( t ) from 19 October and 29 October 2020 (black dots). The red curve corresponds to the best fit of model (18) to Residual 1 ( t ) .
Figure 6. In this figure, the black dots correspond to the function t Residual 1 ( t ) from 19 October and 29 October 2020 (black dots). The red curve corresponds to the best fit of model (18) to Residual 1 ( t ) .
Biology 11 01825 g006
Figure 7. We plot the spectrum of the Markovian Leslie matrices L (red dots) when n = 3 , 5 , 6 , 7 , (respectively in (ad)) giving the best match to the secondary eigenvalues λ 2 🟉 and λ 3 🟉 (green dots). We observe that the best fit of the two secondary eigenvalues remain far away from λ 2 🟉 and λ 3 🟉 for n = 3 , then get closer for n = 5 , and are very close for n = 6 and n = 7 .
Figure 7. We plot the spectrum of the Markovian Leslie matrices L (red dots) when n = 3 , 5 , 6 , 7 , (respectively in (ad)) giving the best match to the secondary eigenvalues λ 2 🟉 and λ 3 🟉 (green dots). We observe that the best fit of the two secondary eigenvalues remain far away from λ 2 🟉 and λ 3 🟉 for n = 3 , then get closer for n = 5 , and are very close for n = 6 and n = 7 .
Biology 11 01825 g007
Figure 8. We plot the spectrum of the Leslie matrix L (red dots) when n = 3 , 5 , 6 , 7 , (respectively in (ad)) giving the best match to the secondary eigenvalues λ 2 🟉 and λ 3 🟉 (green dots). The red dots correspond to the spectrum of L for all the possible matrices L, having their second pair of eigenvalues close to the minimal distance to λ 2 🟉 and λ 3 🟉 .
Figure 8. We plot the spectrum of the Leslie matrix L (red dots) when n = 3 , 5 , 6 , 7 , (respectively in (ad)) giving the best match to the secondary eigenvalues λ 2 🟉 and λ 3 🟉 (green dots). The red dots correspond to the spectrum of L for all the possible matrices L, having their second pair of eigenvalues close to the minimal distance to λ 2 🟉 and λ 3 🟉 .
Biology 11 01825 g008
Figure 9. In this figure, we use the distributions d R 0 ( d ) minimizing the distance | λ 2 🟉 λ 2 | and | λ 3 🟉 λ 3 | whenever n = 7 . In (a), we plot the average distribution d R 0 ( d ) (red curve), standard deviation (blue region), and 95 % confidence interval (light blue region). In (b), we plot the 24 distributions d R 0 ( d ) . In (c), we give a histogram with the multiple values of R 0 = d = 1 n R 0 ( d ) . We observe that some of the d R 0 ( d ) are similar to the case n = 6 , with a maximum on day d = 6 , but on average the maximum value is on day 7.
Figure 9. In this figure, we use the distributions d R 0 ( d ) minimizing the distance | λ 2 🟉 λ 2 | and | λ 3 🟉 λ 3 | whenever n = 7 . In (a), we plot the average distribution d R 0 ( d ) (red curve), standard deviation (blue region), and 95 % confidence interval (light blue region). In (b), we plot the 24 distributions d R 0 ( d ) . In (c), we give a histogram with the multiple values of R 0 = d = 1 n R 0 ( d ) . We observe that some of the d R 0 ( d ) are similar to the case n = 6 , with a maximum on day d = 6 , but on average the maximum value is on day 7.
Biology 11 01825 g009
Figure 10. We plot the daily basic reproduction numbers R 0 ( d ) obtained for n = 3 in (a), n = 5 in (b), n = 6 in (c), and n = 7 in (d). The distribution for n = 7 corresponds to the red curve in Figure 9.
Figure 10. We plot the daily basic reproduction numbers R 0 ( d ) obtained for n = 3 in (a), n = 5 in (b), n = 6 in (c), and n = 7 in (d). The distribution for n = 7 corresponds to the red curve in Figure 9.
Biology 11 01825 g010
Figure 11. In this figure, we plot the daily number of reported cases data from October 19 and November 19, 2020 (black dots) and from model (24) and (25) with the values of R 0 ( d ) obtained in Figure 10c (red dots).
Figure 11. In this figure, we plot the daily number of reported cases data from October 19 and November 19, 2020 (black dots) and from model (24) and (25) with the values of R 0 ( d ) obtained in Figure 10c (red dots).
Biology 11 01825 g011
Figure 12. Autocorrelation Function (ACF) (left hand side) and Partial Autocorrelation Function (PACF) (right hand side) applied to the daily number of cases for Japan between 19 October and 19 November 2020.
Figure 12. Autocorrelation Function (ACF) (left hand side) and Partial Autocorrelation Function (PACF) (right hand side) applied to the daily number of cases for Japan between 19 October and 19 November 2020.
Biology 11 01825 g012
Figure 14. In this figure, we plot the cumulative number of reported cases data between 19 October and 19 November 2020 (black dots). We plot the best fit of the model (30) to the cumulative data (red curve).
Figure 14. In this figure, we plot the cumulative number of reported cases data between 19 October and 19 November 2020 (black dots). We plot the best fit of the model (30) to the cumulative data (red curve).
Biology 11 01825 g014
Figure 15. In this figure, we consider the case n = 25 . We plot the distributions of daily basic reproduction numbers d R ¯ 0 ( d ) corresponding to the distributions having some secondary eigenvalues and fourth eigenvalues at a distance less than 10 2 to the best match. The red curve is the average distribution d R ¯ 0 ( d ) . The blue region corresponds to the standard deviation around the mean distribution. The light blue region corresponds to the 95 % confidence interval.
Figure 15. In this figure, we consider the case n = 25 . We plot the distributions of daily basic reproduction numbers d R ¯ 0 ( d ) corresponding to the distributions having some secondary eigenvalues and fourth eigenvalues at a distance less than 10 2 to the best match. The red curve is the average distribution d R ¯ 0 ( d ) . The blue region corresponds to the standard deviation around the mean distribution. The light blue region corresponds to the 95 % confidence interval.
Biology 11 01825 g015
Figure 16. In this figure, we consider the case n = 25 . We plot the distributions of daily basic reproduction numbers d R 0 ( d ) = R ¯ 0 ( d ) λ 1 d , where R ¯ 0 ( d ) is the red curve in Figure 15.
Figure 16. In this figure, we consider the case n = 25 . We plot the distributions of daily basic reproduction numbers d R 0 ( d ) = R ¯ 0 ( d ) λ 1 d , where R ¯ 0 ( d ) is the red curve in Figure 15.
Biology 11 01825 g016
Figure 17. In this figure, we consider the case n = 25 . We plot the spectrum of the Leslie matrix L (red dots) when d R ¯ 0 ( d ) corresponds to the average distribution (i.e., the red curve in Figure 15).
Figure 17. In this figure, we consider the case n = 25 . We plot the spectrum of the Leslie matrix L (red dots) when d R ¯ 0 ( d ) corresponds to the average distribution (i.e., the red curve in Figure 15).
Biology 11 01825 g017
Figure 18. In this figure, we consider the case n = 25 , and we plot a histogram for the values of the basic reproduction number obtained by summing the distributions d R 0 ( d ) from Figure 16.
Figure 18. In this figure, we consider the case n = 25 , and we plot a histogram for the values of the basic reproduction number obtained by summing the distributions d R 0 ( d ) from Figure 16.
Biology 11 01825 g018
Figure 19. In this figure, we plot the daily number of reported cases data between 19 October and 19 November 2020 (black dots). The red curve corresponds to ϕ 1 + ϕ 2 , and the green dots correspond (34) and (35) whenever R 0 ( d ) comes from the average distribution (i.e., the red curve in Figure 15). We observe a very good match between the green dots and the red curve (the phenomenological model).
Figure 19. In this figure, we plot the daily number of reported cases data between 19 October and 19 November 2020 (black dots). The red curve corresponds to ϕ 1 + ϕ 2 , and the green dots correspond (34) and (35) whenever R 0 ( d ) comes from the average distribution (i.e., the red curve in Figure 15). We observe a very good match between the green dots and the red curve (the phenomenological model).
Biology 11 01825 g019
Table 1. The above reproduction numbers are obtained by using the formula R 0 = d = 1 n R 0 ( d ) .
Table 1. The above reproduction numbers are obtained by using the formula R 0 = d = 1 n R 0 ( d ) .
n3567
R 0 1.021.041.061.07
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Demongeot, J.; Magal, P. Spectral Method in Epidemic Time Series: Application to COVID-19 Pandemic. Biology 2022, 11, 1825. https://doi.org/10.3390/biology11121825

AMA Style

Demongeot J, Magal P. Spectral Method in Epidemic Time Series: Application to COVID-19 Pandemic. Biology. 2022; 11(12):1825. https://doi.org/10.3390/biology11121825

Chicago/Turabian Style

Demongeot, Jacques, and Pierre Magal. 2022. "Spectral Method in Epidemic Time Series: Application to COVID-19 Pandemic" Biology 11, no. 12: 1825. https://doi.org/10.3390/biology11121825

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop