Analysis of COVID-19 in Japan with Extended SEIR model and ensemble Kalman filter

The SIR model is a deterministic epidemic model expressed by a system of differential equations. Its construction is based on Kermack–McKendrick theory, and it describes the transmission process of an infectious disease. A given population of size $N$ is divided into three mutually exclusive sub-populations call compartments: the susceptible hosts $S$, the infectious hosts $I$, and the recovered hosts $R$. For any given time $t$, $S(t)$, $I(t)$, and $R(t)$ describe the number of individuals in each respective compartment, and they satisfy $S(t)+I(t)+R(t)=N$. Individuals in compartment $S$ can be infected by individuals in compartment $I$ through direct contact. Once a successful contact happens, the infected individual leaves compartment $S$ and becomes a member of compartment $I$. The number of newly infected individuals per unit time is given by $\beta I(t)S(t)/N$, where the parameter $\beta$ is called the transmission coefficient. Similarly, once a patient recovers, this person leaves compartment $I$ and goes to compartment $R$ immediately. The transfer between $I$ and $R$ is controlled by the recovery rate $\gamma$. The transfer rate, namely the number of transfer individuals per unit time, is given by $\gamma I(t)$.

Figure 2.1 shows the process of the SIR model, and the following differential system describes its dynamic:

Note that the system (1) assumes a permanent immunity once recovered. In other words, no individual can be infected a second time. The model also assumes a constant total population: there is no inflow to the system, or outflow from the system.

The SEIR model is very similar to the previous model, but with an additional compartment $E$ between the compartments $S$ and $I$. The individuals in $E$ are exposed, namely they have been contaminated, but they are not infectious yet. The time spent in $E$ corresponds to the incubation period, before becoming infectious. The transfer rate from $E$ to $I$ is given by $\varepsilon E(t)$.

Figure 2.1 shows the process of the SEIR model, and the following differential system describes its dynamic:

The extended SEIR model has been developed based on some specific features of the COVID-19 epidemic, as described now. In early 2020, health services already noticed that some infected individuals, who did not show any symptom, could spread the disease. These persons correspond either to asymptomatic hosts or to pre-symptomatic hosts. In the first cohort, people will never show any symptom, while the second cohort corresponds to individuals who will show some symptoms subsequently. Clearly, asymptomatic hosts and pre-symptomatic hosts are difficult to be identified by health services, even though they play an important role in the spread of the disease. Symptomatic individuals are also infectious, but they can be more easily identified precisely because of their symptoms. Thus, if symptomatic individuals correspond to the compartment $I$ of the above SIR or SEIR models, then there is no compartment left for the asymptomatic or the pre-symptomatic individuals.

The existence of COVID-19 infectious spreaders who do not show (yet or at all) any symptom already brings a lot of uncertainties to health services. In addition, whenever the symptoms are relatively mild, some symptomatic individuals do hesitate to report the health service [31, p. 11]. As a consequence, daily new cases are under-reported. Combining this effect with some delayed information, some inaccurate tests (especially in the first phase of the epidemic), and other reasons that we are not aware of, it turns out that that the reported data are not very accurate. In such a situation one needs to use proper strategies to analyse the observation data.

In order to meet the special features of COVID-19, we separate the compartment $I$ of the SEIR model into two compartments $I_{a}$ and $I_{s}$. These compartments correspond to asymptomatic ​/ ​pre-symptomatic and to symptomatic states, respectively. Both $I_{a}$ and $I_{s}$ can infect $S$. As shown in Figure 2.2, the transmission processes related to $I_{a}$ and to $I_{s}$ have transmission coefficients $\beta_{a}$ and $\beta_{s}$ respectively. Once an individual in $S$ gets infected, this person becomes a member of $E$, and then moves to $I_{a}$ when it becomes infectious. In $I_{a}$, part of these individuals (the asymptomatic) will never generate any symptoms. They will thus spend some time in $I_{a}$, and then recover. The corresponding compartment is denoted by $R_{a}$. In contrast, the other part of individuals in $I_{a}$ will develop symptoms, and then move to $I_{s}$. In $I_{s}$, a majority of individuals will be identified by health services, but as mentioned above a minority will recover without being identified by any health services. Compartment $R_{s}$ corresponds to those recovered individuals who have not been registered. The identified persons in $I_{s}$ will then move to the compartment $H$, which corresponds to hospitalized or treated patients. The ones staying at home but under medical control, or the ones isolated in a hotel, are all included in the compartment $H$. Finally, the patients in $H$ will recover, and thus move to the compartment $R$, or will decease and end up in the compartment $D$. As for the SIR or SEIR models, we assume a permanent immunity, which means that there is no flow from $R_{a}$, $R_{s}$, or $R$ to $S$. Also, the total number $N$ of individual is constant, namely at all time one has

Based on the Figure 2.2 and on the explanations provided above, the differential system of the extended SEIR model reads:

where $\beta_{a},\beta_{s},\varepsilon,\delta,\gamma_{a},\gamma_{s},\gamma_{D}$ and $\gamma_{H}$ are some medical parameters. Note that some of them are time dependent.

The basic reproduction number, denoted by $R_{0}$, can be interpreted as the average number of secondary cases generated by on primary case in a susceptible population. It is commonly admitted that $R_{0}=1$ is a threshold value, as explained below. We also refer to [8, 16] for more explanations and more precise statements.

To study $R_{0}$ for different models, a general definition of the basic reproduction number is introduced based on the notion of disease free equilibrium (DFE). In a DFE the population remains in the absence of disease. For example, in the SIR or SEIR models, it means that $I=0$ while in the extended SEIR model it means that $I_{a}=I_{s}=0$. In this context, the basic reproduction number can be defined as the number of new infections produced by a typical infectious individual in a population at a DFE. The main feature of $R_{0}$ is that it corresponds to a threshold parameter, namely if $R_{0}<1$, the DFE is locally asymptotically stable, while if $R_{0}>1$, the DFE is unstable and an outbreak is possible.

The precise expression for $R_{0}$ is clearly model dependent, but numerous examples are available in the literature. For example, let us consider a general compartmental disease transmission model. Such models are represented by a system of ordinary differential equations, and under natural assumptions it can be shown that these models have a DFE, see [8]. In this reference, the precise expression for $R_{0}$ is then provided for some classes of models, and the extended SEIR model meets the assumptions of the staged progression model, as presented in [8, Sec. 4.3]. For the extended SEIR model, the expression for $R_{0}$ then reads:

where $\beta_{a}$ and $\beta_{s}$ are the initial transmission coefficients at time $0$.

To understand the above expression, observe that $\delta/(\delta+\gamma_{a})$ corresponds to the fraction of individuals of $I_{a}$ going to the compartment $I_{s}$, while $1/(\delta+\gamma_{a})$ and $1/(\gamma_{s}+\tau_{H})$ define the average times an infected individual spends in compartments $I_{a}$ and $I_{s}$ respectively. Thus, each term on the R.H.S. of (4) represents the number of new infections generated by an infectious individual during the time spent in the compartments $I_{a}$ and $I_{s}$.

In contrast, the definition of the effective reproduction number $R_{t}$ at time $t$ is much more delicate. The various challenges and possible errors have recently been discussed in [15], to which we refer for a thorough discussion. In our setting, we shall keep the interpretation of $R_{t}$ as the average number of secondary cases generated by one primary case. This approach is possible because the transmission coefficients at time $t$ are available in our approach, and therefore one can compute $R_{t}$ with (4) and the corresponding $\beta_{a}$ and $\beta_{s}$ at time $t$. Additional information on $R_{t}$ will be provided in Section 5.

It clearly appears in Figure 2.2 and in the corresponding system (2.2) that several parameters control the flows between the compartments. The values of these parameters may result in very different behaviors of the model. We refer for example to [22, Sec. 1.4] and to [4, Chap. 2] for general discussions of model behaviors and the role of parameters. In our setting, some parameters are easy to evaluate, as for example the recovery rate $\gamma_{H}$ or the death rate $\gamma_{D}$, but others are difficult to estimate, as for example the transmission coefficients $\beta_{a}$ and $\beta_{s}$. Let us also stress that some parameters depend on location. Since our experiment is based on data from Japan, see Section 4, the medical parameters are chosen accordingly. For that reason, we use research or survey results specific to Japan, or even more precisely to specific prefectures in Japan.

In Table 1 we list the estimated values of some parameters for the extended SEIR model, and provide the sources of information. Several parameters in the table have the form of the product of a percentage and the inverse of a time length. A similar setting for the construction of the parameters can be found for example in [9, 21]. For $\delta$ and $\gamma_{a}$, the percentage parts should sum up to $1$. For the parameters $\tau_{H}$ and $\gamma_{s}$, some information can be found in the surveys [31, 37] and the health services website [38]. For parameter $\gamma_{H}$ and $\gamma_{D}$, instead of using constant value estimated by health services, we shall use the observation data to estimate their values on a daily basis. The details will be discussed in Section 4.

Let us stress that the value assigned to some of these parameters has evolved during the last 12 months. For example, the ratio of asymptomatic individuals was thought to be very high at the beginning of the epidemic (up to $80\%$), but some recent research [3, 5, 32] have revised this ratio to $17\%$ to $20\%$. Our knowledge about the length of infectious periods has also evolved, and the possible values are spread over a rather long interval. In Table 1 we list some mean values, but in the simulations additional uncertainties are implemented. Since pre-symptomatic patients become infected before the appearance of symptoms, the incubation period (encoded in $\varepsilon$) has been shorten a little bit, and the last 2 days of this incubation period have been moved to $I_{a}$.

One very delicate question is the relation between $\beta_{a}$ and $\beta_{s}$, namely between the transmission coefficient for asymptomatic ​/ ​pre-symptomatic and the transmission coefficient for symptomatic. For our investigations, we shall rely on the result of the systematic review [3] which asserts that the relative risk of asymptomatic transmission is $42\%$ lower than that for symptomatic transmission. As a consequence, we shall fix

This factor $0.58$ is slightly smaller but of a comparable scale compared to earlier investigations, see for example [5].

As introduced in Section 1, one of the aims of this study is to estimate the unknown parameters $\beta_{a}$ and $\beta_{s}$, and the unobservable compartments $I_{a}$ and $I_{s}$. Because of relation (5), it is clear that for the parameters only $\beta_{s}$ has to be evaluated. The method that we are going to introduce in Section 3 requires observations of some compartments in the extended SEIR model. In Figure 2.2, we highlight the three compartments with observations in yellow, namely $H$, $R$ and $D$. The data corresponding to these compartments may not be perfect, and for the analysis based on these data we shall take some uncertainties into consideration. More precisely, some random noise will be added to the observations.

For our experiment, we shall firstly and mainly use the data about Tokyo, with a total population $N=13^{\prime}955^{\prime}000$ (approximate mean of the population of 2020 and 2021). Subsequently, other prefectures in Japan are also considered. As already mentioned, the compartment $H$ corresponds to the identified individuals either hospitalized or treated at home or in a hotel. On the other hand, $R$ and $D$ describe the accumulated number of recovered and deceased individuals. The information about these three compartments are very easy to find for Tokyo, but also for any region in Japan. One can check for example the website of Ministry of Health, Labour and Welfare, prefecture’s websites or some mass communication companies’ websites to get more details. The information is provided on a daily basis. Note however, that these records were not very accurate at the beginning of the outbreak. This was caused by the delay of reports, but also by some changes in the policy for the collect of information. Our analysis will take this source of uncertainties into consideration.

For the state space model mentioned in Section 3, the dynamical system (6) is described by the differential system (2.2). To estimate the parameter $\beta_{s}$, we use the augmented state by adding one more equation $\frac{d\beta_{s}}{dt}=0$ to the system (2.2), assuming persistence for $\beta_{s}$. Namely, $\beta_{s}$ will stay at the same value during the time integration process, and will be updated during the analysis step of the data assimilation. In other words, if $\beta_{s}$ has a correlation with the observations, $\beta_{s}$ will be updated together with the other states. In system (2.2), the unit of time $t=1$ represents one day. Thus, the one day forecast from day $n$ is obtained by integrating the system (2.2) on $[n,n+1]$ with the initial values equal to the analysis on day $n$. To avoid negative values in the analysis step, all the compartments are transformed to $\log$ scale with base ${\mathrm{e}}$. The lower case will be used for these new variables, as for example $s(t):=\log S(t)$, and the corresponding equation becomes $\frac{ds(t)}{dt}=\frac{1}{S(t)}\frac{dS(t)}{dt}$. The DA analysis update directly applies to the $\log$-transformed values.

Now the forecast value of $x_{t}$ is a $9$-dimensional vector

Note that the compartment $s(t)$ is not explicitly included since it can be deduced from the conservation relation (2). Before performing the DA analysis update, the above vector has to be projected to the observation space by the observation system (7). In our setting, the operator $H_{t}$ introduced in (7) is defined by the projection which sends the forecast vector (19) on $\big{(}h(t),r(t),d(t)\big{)}^{T}$, namely on the three compartments with observations.

To initiate the experiments, initial values for all compartments and parameters have to be provided. When these initial values are far away from the observations, the system goes through an unreasonable transition period. To avoid or reduce the unreasonable transition, preliminary experiments were performed for determining suitable initial values. Namely, we run the system for a few days, starting with a few individuals in some compartments and a presumed value for $\beta_{s}$. Different values of $\beta_{s}$ were tested, until values comparable to the data of the compartments $H,D,R$ on day zero (March $6^{\mathrm{th}}$, 2020) were obtained. The corresponding values for the different compartments were then chosen as initial conditions. However, independent normal distributed random errors are also added to create 50 ensemble members $\{x^{a(i)}_{0},i=1,\dots,50\}$. Note that each member $x^{a(i)}_{0}$ is a $9$-dimensional vector containing the initial conditions for the compartments and for the parameter $\beta_{s}$.

The ensemble on day zero are then integrated by the model (2.2) for one day and the ETKF will update the one day forecast by using the observations. The procedure of the experiments is shown in Figure 4.1.

As mentioned in Section 2.4, parameters $\gamma_{H}$ and $\gamma_{D}$ are estimated by assimilating the observations. Consider equation $\frac{dR}{dt}=\gamma_{H}H$ in model (2.2). Since the time is discrete, we can rewrite this relation as $R(t+1)-R(t)=\gamma_{H}(t)H(t)$ which leads to

Similar computation can be generated for $\gamma_{D}$ also, namely,

The daily values of both parameters are shown in Figure 5. For the implementation of these parameters for the computation of the $1$-day forecast, we have used a slightly smoothed version obtained by a $7$-day convolution with the symmetric weights $\frac{1}{64}\big{(}1,6,15,20,15,6,1)$. The effect is to decrease the amplitude of the weekly oscillations.

For the integration process, we have also included some uncertainties to all parameters: to the ones presented in Table 1, but also to the values of the parameters $\gamma_{H}$, $\gamma_{D}$. Thus, the $1$-day forecast is obtained with the parameters of the previous day, each of them perturbed by a normal distribution $N\big{(}0,(M/10)^{2}\big{)}$, where $M$ corresponds to the value of this parameter. These perturbations are independent and randomly generated for each $1$-day forecast.

The last initial setting for data assimilation is the observation error covariance matrix $\textbf{R}_{t}$. We assume that the covariance of the errors between different observations is zero, namely, the matrix $\textbf{R}_{t}$ is diagonal. Since we observe three states, $\textbf{R}_{t}$ is a $3\times 3$ diagonal matrix. Clearly, the bigger the variance, the weaker contribution the corresponding observation will make to the update.

For simplicity, we assume that the observation error covariance is independent of time. The diagonal elements are chosen as $\big{(}\log(1.3)\big{)}^{2}$ for any time $t$, representing the observation error variance: Namely, for $x\in\{h,r,d\}$ one has $\big{[}x-\log(1.3),x+\log(1.3)\big{]}$ as $68\%$ CI and $\big{[}x-2\log(1.3),x+2\log(1.3)\big{]}$ as $95\%$ CI. Considering that all the observations have been transformed into $\log$ scale, that is, the observations in original scale are distributed in $\big{[}X/(1.3),(1.3)X\big{]}$ as $68\%$ CI and $\big{[}X/(1.69),(1.69)X\big{]}$ as $95\%$ CI.

Before showing the analysis result, a scatter plot of the values of $h$ and of $\log(\beta_{s})$ for the different ensemble members on a certain day is provided in Figure 6. Indeed, for the parameter estimation, one assumption is the linear relation between the ensemble of parameters and the ensemble of compartments with observations. As shown in Figure 6, a positive correlation is observed between $h$ and $\log(\beta_{s})$. Thus we can use the observation of $H$ for the estimation of $\beta_{s}$.

In Figure 7 we finally provide the analysis value for $R_{t}$ from March $6^{\mathrm{th}}$ 2020 to the middle of October 2021. The black curve represents the mean value of the ensemble. The two shaded regions represent the $68\%$ CI (dark orange) and $95\%$ CI (light orange). The analysis value of $R_{t}$ is computed by using $\beta_{s}$ together with the information contained in (4) and in (5). The red curve shows the $R_{t}$ provided by Toyokeizai.net for reference. The grey shaded regions represent the periods of state of emergency in Tokyo. Note that the jump at the end of May 2020 is due to an abrupt change in the value of $\tau_{H}$ as indicated in Table 1. Additional discussions about this graph will be provided in the subsequent section.

The analysis results for the compartments $E$, $I_{a}$, $I_{s}$, $H$, $R$, $D$ with confidence intervals are also presented in Figure 8. For compartments $H$, $R$, $D$, the red dots represent the observations of Tokyo. For all pictures, the $y$-axis is in $\log_{10}$ scale.

In Figure 7, one first notices that the similarity between the analysis $R_{t}$ and the reference $R_{t}$ provided by Toyokeizai.net. The Toyokeizai $R_{t}$ on day $t$ is computed by the following formula:

where the average generation time is assumed to be 5 days and the reporting interval is assumed to be 7 days. The formula is a simplified version of a maximum likelihood estimation for the effective reproduction number provided by H. Nishiura [28]. Note that this formula is based on the daily confirmed cases, which experiences quick changes everyday. On the other hand, the analysis $R_{t}$ and its confidence interval is computed by formula (4) which involve the analysis ensemble of $\beta_{s}$ and $\beta_{a}$. The average over the ensemble member provides a smoother evolution of $R_{t}$, which is certainly closer to a true statistic’s evolution. Note that a third approach for the computation of $R_{t}$ is available in [35]: it involves an agent-based model together with a particle filter.

In Figure 8, we show the analysis mean of compartments $E$, $I_{a}$, $I_{s}$ with the $68\%$ CI (dark orange) and $95\%$ CI (light orange), and the same results for compartments $H$, $R$ and $D$ with the daily observations represented by red dots. One notes that the analyses of $H$, $R$, $D$ fit the observations appropriately except at the beginning of the observation period. This may be related to the initial guess and the model imperfections.

As mentioned in Section 2.4, the ratio between $\beta_{a}$ and $\beta_{s}$ is fixed at $0.58$ based on the information provided by the literature. Since this ratio is quite uncertain, we performed similar investigations with different ratios varying from $0.1$ to $1.0$. For $k=\frac{\beta_{a}}{\beta_{s}}$ with $k=0.1,0.3,0.58,0.8$ and $1.0$, we run the experiments independently (always with the data from Tokyo) and the analysis mean of $R_{t}$ are shown in Figure 9. The patterns of each $R_{t}$ are similar, especially after May $31^{\mathrm{st}}$, 2020. However, at the beginning of the outbreak and during the first state of emergency, all the analysis value of $R_{t}$ experience quick changes but with different speed. In Figure 10 we also provide the analysis mean of the different compartments, since they also depend on the ratio between $\beta_{a}$ and $\beta_{s}$. The same conclusion is obtained: except for the first three months, the possible ratio do not generate any noticeable difference.

In Section 4.1, the error covariance matrix $\textbf{R}_{t}$ is defined as a diagonal matrix $\big{(}(\log(1.3)\big{)}^{2}I$ where $I$ is the identity $3\times 3$ matrix. We also run a sensitivity test for different choices of the covariance matrix, and checked if it has an effect on the analysis results. We choose $\textbf{R}_{t}$ equals to ${\rm sd}^{2}I$ with ${\rm sd}=\log(1.1)$, $\log(1.3)$, $\log(1.5)$, $\log(2.0)$, and $\log(2.5)$. For these investigations the ratio between $\beta_{a}$ and $\beta_{s}$ is fixed as $0.58$. In Figure 11(a), the analysis means are provided, and they are clearly close to each other except at the beginning stage of the outbreak. The different spreads are shown in Figure 11(b), and one notices that the spread increases with the ${\rm sd}$, but for ${\rm sd}$ larger than $\log(1.5)$, the increase is less clear. We also show the analysis mean of compartments $E$, $I_{a}$, $I_{s}$, $H$, $R$ and $D$ in Figure 12(a). Again, the analysis results of these compartments with different observation error covariance are quite close to each other. In Figure 12(b), the spread is computed by the standard deviation of ensemble members of the $\log-$transformed compartments. The relations of spread are quite clear for the three compartments with observations, and there are some overlaps for spread of $e$, $i_{a}$ and $i_{s}$, but in general the spread increases with ${\rm sd}$.

As shown in Table 1, the proportion of symptomatic and asymptomatic used in our investigation is $83\%$ and $17\%$, respectively. These values were determined based on sources mentioned in Section 2.4. However, since asymptomatic cases are very difficult to detect, and since we can not be fully confident in the above ratio, a sensitivity test is necessary. To do this, we choose two different combinations which are $70\%$ and $30\%$, and $50\%$ and $50\%$, and keep the observation error covariance matrix at $\big{(}(\log(1.3)\big{)}^{2}I$ and $k=0.58$. Compared to the previous settings, in these two scenarios more people become asymptomatic and recover without showing any symptoms, as illustrated in Figure 13(a). On the other hand, the three curves for $R_{t}$ are close to each others, but the one corresponding to only $50\%$ of symptomatic is slightly higher than the other two, see Figure 13(b). In order to understand this, let us recall that the relation between the transmission coefficients $\beta_{a}$ and $\beta_{s}$ is set to be $0.58$, it means that the asymptomatic cases are less infectious than the symptomatic and symptomatic cases. Thus, when more people do not show any symptoms, the transmission coefficient has to be bigger to create enough cases to fit the observations. As a consequence, $R_{t}$ will also be bigger.