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

We introduce an extended SEIR infectious disease model with data assimilation for the study of the spread of COVID-19. In this framework, undetected asymptomatic and pre-symptomatic cases are taken into account, and the impact of their uncertain proportion is fully investigated. The standard SEIR model does not consider these populations, while their role in the propagation of the disease is acknowledged. An ensemble Kalman filter is implemented to assimilate reliable observations of three compartments in the model. The system tracks the evolution of the effective reproduction number and estimates the unobservable subpopulations. The analysis is carried out for three main prefectures of Japan and for the entire country of Japan. For these four communities, our estimated effective reproduction numbers are more stable than the corresponding ones estimated by a different method (Toyokeizai). We also perform sensitivity tests for different values of some uncertain medical parameters, like the relative infectivity of symptomatic/asymptomatic cases. The regional analysis results suggest the decreasing efficiency of the states of emergency.

Introduction

The outbreak of coronavirus 2019 (COVID-19) had a huge impact on human society, and its devastating effects are still present more than 18 months after its official acknowledgment. One of the specific characteristics of this disease increases the difficulties for any scientific research based on infectious disease models, namely the existence of asymptomatic cases. Their existence, but also their potential infectiousness, blur the full picture of the disease spreading inside the society [1]. The simultaneous existence of mild symptomatic cases, which often remain undocumented, also generates an additional challenge for any epidemiological investigations. For example, one study discussing the potential fraction of undocumented cases at the beginning of the outbreak ends up with a ratio of 86% undocumented infections [2].

There exists now a huge literature related to COVID-19, but only a small number of papers are using techniques of data assimilation, as we shall do. For that reason, we present only a few references closely related to our investigations. To get a better understanding of the spread of COVID-19, some studies start with a specially designed infectious disease model with appropriate structures. Those elaborated models may accurately describe the process in reality but they also bring more uncertainties and unknown parameters. For example in [3], a variables model is adopted with each variable representing a subpopulation. The structure related to asymptomatic, pre-symptomatic, detected, and undetected cases are all properly described by the structure. However, the inference is done with only observations of compartments among all variables. Another example is [4], where an extended SEIR model with age-classes is used to estimate the posterior distribution of parameters and make short period predictions. The distribution of mild, severe and fatal cases for each age class is introduced as pre-defined parameters. Also, a matrix which describes the relative transmissions between each pair of age groups is predefined.

Despite these elaborated models, most of the investigations are still performed on the common infectious disease models, namely the SIR and SEIR models. Such models are easier to implement and contain less parameters, but they also sacrifice some details of the process. For example the SEIR model is used to study the dynamical behavior of COVID-19 outbreaks at the regional level in [5], and the SIR model was used to track the effective reproduction number for countries in [6]. A slightly more elaborated model is also introduced in [7], namely the SITR model with standing for under treatment, for inferring infection numbers and parameters values.

For dealing with real data, one of the recent tools employed for the study of infectious diseases is data assimilation. The techniques developed for data assimilation have the ability to optimally meld dynamical systems with noisy and imperfect observations. They also provide forecasts and estimations of parameters and variables [8]. As examples of investigations studying the details of the implementation of these techniques to epidemic models, let us mention [9], [10]. The previously cited work [3] also uses statistical data assimilation for identifying the measurements required for getting accurate states and parameters estimations. Finally, in [11] an ensemble Kalman filter is adopted to forecast the COVID-19 pandemic in Saudi Arabia with an extended SEIR model including vaccination.

Let us now present our research. First, we construct an extended SEIR model which takes into account one of the main specificities of COVID-19: the simultaneous existence of an asymptomatic/pre-symptomatic population and of a symptomatic population. These two infected populations have different characteristics which are encoded with different parameters. The model is also constructed such that it can be fed with the data of only three compartments, namely the Hospitalized (or treated) population, the Recovered population, and the Deceased population, naturally identified by the letters , , and , respectively. Note that the population or are coming from the population , and these three populations are considered as fully recorded by health services. Except for the first couple of weeks at the beginning of the epidemic, these data are also considered as the most reliable ones. With this model and these data, our main aim is to provide information about the effective reproduction number and about the population of asymptomatic/pre-symptomatic or undetected symptomatic, with uncertainty ranges. For these investigations, we shall also mainly concentrate on Tokyo, but provide additional analysis with two additional prefectures (Osaka and Kanagawa) and with the entire country of Japan. Let us mention that targeting a specific community gives the opportunity to rely on local parameters and to use information provided by local medical research or health services. We also emphasize that for data assimilation, we use the ensemble transform Kalman filter [12] which is efficient compared with the standard extended Kalman filter.

As a result of our investigations we get an effective reproduction number for the four communities (Tokyo, Osaka, Kanagawa, Japan), as well as a time dependent estimate for the number of asymptomatic/pre-symptomatic in these four populations. A comparison between our estimated effective reproduction number and the same statistic provided by Toyokeizai [13] shows the reliability of our strategy. As an asset of our approach, we test the sensitivity of our model and its outcomes to the values of uncertain parameters borrowed from the literature. These experiment’s results show a robust performance of the strategy used in our investigations. Another clear result from our computation of the effective reproduction number is the decay of the effectiveness of the states of emergency as their number increases. This effect is clearly visible independently for the three prefectures and for the entire country. Note that a similar effect is also visible in the use of public transportation, see for example [14]. Other outcomes of our investigation are gathered in Sections 4, 5.

Let us finally describe the structure of this paper. In Section 2, we recall the standard SIR and SEIR models, and introduce the extended SEIR model. We also provide information for the computation of the effective reproduction number, and discuss the constant parameter values and the observations. In Section 3, we explain the technical details of ensemble transform Kalman filter. Readers familiar with data assimilation and Kalman filters can skip this section without any loss for the application to infectious diseases. In Section 4, background information about the experiments, and subsequently technical information, are provided. The main results of our investigations are also presented in this section. The discussion and the comparisons between the different communities are provided in Section 5. We finally prepare the ground for future projects.

Extended SEIR model

In this section, we recall a few information on the SIR and SEIR models, and introduce the extended SEIR model. We also discuss medical parameters, and provide some explanations about the observations.

The SIR and SEIR models

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 is divided into three mutually exclusive sub-populations called compartments: the susceptible hosts , the infectious hosts , and the recovered hosts . For any given time , , , and describe the number of individuals in each respective compartment, and they satisfy . Individuals in compartment can be infected by individuals in compartment through direct contact. Once a successful contact happens, the infected individual leaves compartment and becomes a member of compartment . The number of newly infected individuals per unit time is given by , where the parameter is called the transmission coefficient. Similarly, once a patient recovers, this person leaves compartment and goes to compartment immediately. The transfer between and is controlled by the recovery rate . The transfer rate, namely the number of transfer individuals per unit time, is given by .

Fig. 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.

Fig. 1

Transfer diagram for the SIR model.

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

Fig. 2 shows the process of the SEIR model, and the following differential system describes its dynamic:

Fig. 2

Transfer diagram for the SEIR model.

The extended SEIR model

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 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 to the health service [15, 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 the reported data are not very accurate. In such a situation one needs to use proper strategies to analyze the observation data.

In order to meet the special features of COVID-19, we separate the compartment of the SEIR model into two compartments and . These compartments correspond to asymptomatic/pre-symptomatic and to symptomatic states, respectively. Both and can infect . As shown in Fig. 3, the transmission processes related to and to have transmission coefficients and respectively. Once an individual in gets infected, this person becomes a member of , and then moves to when it becomes infectious. In , part of these individuals (the asymptomatic) will never generate any symptoms. They will thus spend some time in , and then recover. The corresponding compartment is denoted by . In contrast, the other part of individuals in will develop symptoms, and then move to . In , 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 corresponds to those recovered individuals who have not been registered. The identified persons in will then move to the compartment , 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 . Finally, the patients in will recover, and thus move to the compartment , or will decease and end up in the compartment . As for the SIR or SEIR models, we assume a permanent immunity, which means that there is no flow from , , or to . Also, the total number of individual is constant, namely at all time one has

Fig. 3

Transfer diagram for the extended SEIR model. Compartment of the SEIR model is divided into two compartments (asymptomatic/pre-symptomatic) and (symptomatic) .

Based on Fig. 3 and on the explanations provided above, the differential system of the extended SEIR model reads: where and are some medical parameters. Note that some of them are time dependent.

The reproduction number

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

To study 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 while in the extended SEIR model it means that . 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 is that it corresponds to a threshold parameter, namely if , the DFE is locally asymptotically stable, while if , the DFE is unstable and an outbreak is possible.

The precise expression for 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 [16]. In this reference, the precise expression for is then provided for some classes of models, and the extended SEIR model meets the assumptions of the staged progression model, as presented in [16, Sec. 4.3]. For the extended SEIR model, the expression for then reads: where and are the initial transmission coefficients at time .

To understand the above expression, observe that corresponds to the fraction of individuals of going to the compartment , while and define the average times an infected individual spends in compartments and 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 and .

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

The medical parameters

It clearly appears in Fig. 3 and in the corresponding system (3) 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 [19, Sec. 1.4] and to [20, 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 or the death rate , but others are difficult to estimate, as for example the transmission coefficients and . 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 [4], [21]. For and , the percentage parts should sum up to . For the parameters and , some information can be found in the surveys [15], [22] and the health services website [23]. For parameter and , 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.

Table 1

Medical parameters.

Parameter	Estimation	Source	Remark
ɛ	(3 days)^−1	[27]	incubation period, not contagious
δ	83% × (2 days)^−1(95% CI 80% to 86%)	[24], [27]	proportion of pre-symptomatic ×(duration of pre-symptomatic)^−1
τ_H	78%×(8.3 days)^−1 (∼2020/5/31) 78%×(5.2 days)^−1 (2020/6/1 ∼)	[15], [28]	proportion of detected symptomatic ×(days of symptoms onset)^−1
γ_a	17% × (9 days)^−1(95% CI 14% to 20%)	[24], [26]	proportion of asymptomatic ×(duration of asymptomatic)^−1
γ_s	22% × (7 days)^−1	[15], [22]	proportion of undetected symptomatic ×(duration of symptomatic)^−1
γ_H	computed by observations	[13]	discussed in Section 4
γ_D	computed by observations	[13]	discussed in Section 4

Let us stress that the value assigned to some of these parameters has evolved during our investigations. 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 [24], [25], [26] 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 ) has been shorten a little bit, and the last 2 days of this incubation period have been moved to .

One very delicate question is the relation between and , 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 [24] 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 [25].

The observations

As introduced in Section 1, one of the aims of this study is to estimate the unknown parameters and , and the unobservable compartments and . Because of relation (5), it is clear that for the parameters only 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 Fig. 3, we highlight the three compartments with observations in yellow, namely , and . 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 (approximate mean of the population of 2020 and 2021). Subsequently, other prefectures in Japan are also considered. As already mentioned, the compartment corresponds to the identified individuals either hospitalized or treated at home or in a hotel. On the other hand, and 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.

Ensemble transform Kalman filter

In this section, we introduce the main tool employed in the study, namely the ensemble Kalman filter. The analysis of parameter values and unobservable compartments is performed based on it.

The method Ensemble Transform Kalman Filter (ETKF) used in this paper has been introduced and developed in [12], [29]. It is based on the Ensemble Kalman Filter [30], [31], [32], but it provides analyses in a more efficient way.

Let us consider a discrete time state space model given for any time by where is an -dimensional state vector and an -dimensional observational vector, is the operator which defines the dynamics of the state, is an operator corresponding to the observation model, and provides the observation error. All these vectors or operators are explicitly dependent. The vector represents the state of the dynamical system, while is called the noisy observation, both at time . The observation error follows a normal distribution , where is an observation error covariance matrix. In this framework, the general question is: given a list of noisy, unbiased observations , how can one find the best estimates for ?

When the operators and are linear, they can be represented by matrices, denoted also by and . In such a situation, and with the assumption of a Gaussian initial distribution, the linear model leads to the propagation of Gaussian distributions for the states. More precisely, assume that at some time , there is a prior estimate of the state with its covariance estimate (also called the background error covariance). Then the Kalman filter [33], [34] provides an optimal solution to minimize the variance of its uncertainty. For this reason, it is known as the minimum variance estimator. Formally, the prior estimate is updated by using the information of the observation at time as follows [35, Eq. (1), (2), (8) & (10)]: where is called the Kalman gain, the analysis, and its estimated covariance. The superscript denotes the transpose of a matrix. With the analysis, one can use the forecast model (6) to evolve the analysis and get the prior estimate and the covariance at the next time step.

When the operators and are nonlinear, the Kalman filter is not applicable in its original form, but a simple extension with linear approximation is known to be effective. The straightforward nonlinear extension of the Kalman filter is known as the extended Kalman filter (EKF). In this case, and in Eqs. (8) to (11) are replaced by their linear approximations respectively. For low dimensional models, the cost of EKF is low. However, when the model complexity and the state dimension increases, one may consider the ensemble Kalman filter. In this case, an ensemble of model states is used to estimate the covariance matrices. Namely, let be an ensemble of estimates of the model states at time . By letting this ensemble evolve according to the model (6), one gets the forecast or background ensemble at time , where The mean of the ensemble can be used as the most probable estimate of the state, with its covariance estimated by the sample covariance of the ensemble: where is an matrix with the th column .

The next step is to update the forecast ensemble mean and its covariance by using the information of the observation . The updated ensemble mean is denoted by and its covariance by . To get them, one first needs to determine an analysis ensemble . Similar to the background ensemble, the mean and covariance of the analysis ensemble are given by In our application the operator is linear, so that is given by a matrix. In this case and for ETKF, the computations can be written as follows [35, Eq. (15): where the matrix is called the analysis error covariance in ensemble space. For (16), the background perturbation is transformed into the analysis perturbation with the weight . By adding to each column of , one gets the analysis ensemble that satisfies (12). Note that the ETKF is a deterministic filter since no randomly perturbed observations are used in the computation [36], and it is also a square-root filter because its takes the power of the matrix  [37].

To avoid a variance underestimation, caused for example by the limited ensemble size, an artificial inflation of the ensemble spread is usually applied. In [35], several ways to perform the variance inflation are discussed. In our study, we apply multiplicative inflation and additive inflation. For multiplicative inflation, the background error covariance is multiplied by a tunable factor before the analysis. The multiplicative inflation can be thought as a procedure to increase the influence of the current observations on the analysis. For additive inflation, since the extended SEIR model is a conserved closed system, it is then applied by adding a random vector to the background ensemble. The random vector can be sampled from the normal distribution with mean 0 and some covariance matrix. This covariance matrix is assumed to be proportional to the background error covariance by a tunable parameter  [35].

For the current study with a low dimensional model, we could use EKF or EnKF. However, by considering the possibility of increasing the model complexity in future, we opted for ETKF with 50 ensemble members.

Experiments

In this section, we explain the design of the experiments and illustrate the outcomes by using observations from Tokyo Metropolis. These observations (, , and ) are obtained from the website [13]. Discussions and comparisons with other prefectures are provided in the subsequent section. As mentioned in Section 2.5, health officials started providing data from February 2020 but they included some uncertainties caused by delay or by policy changes. Our experiments start with the observations from March , 2020, when the record of data became more systematic.

The experiments, using data from Tokyo

For the state space model mentioned in Section 3, the dynamical system (6) is described by the differential system (3). To estimate the parameter , we use the augmented state by adding one more equation to the system (3), assuming persistence for . Namely, 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 has a correlation with the observations, will be updated together with the other states. In system (3), the unit of time represents one day. Thus, the one day forecast from day is obtained by integrating the system (3) on with the initial values equal to the analysis on day . To avoid negative values in the analysis step, all the compartments are transformed to scale with base . The lower case will be used for these new variables, as for example , and the corresponding equation becomes . The DA analysis update directly applies to the -transformed values.

Now the forecast value of is a -dimensional vector Note that the compartment 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 introduced in (7) is defined by the projection which sends the forecast vector (17) on , 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 . Different values of were tested, until values comparable to the data of the compartments on day zero (March , 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 . Note that each member is a -dimensional vector containing the initial conditions for the compartments and for the parameter .

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

Fig. 4

Data assimilation flow-chart.

As mentioned in Section 2.4, parameters and are estimated by assimilating the observations. Consider equation in model (3). Since the time is discrete, we can rewrite this relation as which leads to Similar computation can be generated for also, namely, The daily values of both parameters are shown in Fig. 5. For the implementation of these parameters for the computation of the -day forecast, we have used a slightly smoothed version obtained by a -day convolution with the symmetric weights . The effect is to decrease the amplitude of the weekly oscillations.

Fig. 5

Daily estimation of the parameters and .

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 , . Thus, the -day forecast is obtained with the parameters of the previous day, each of them perturbed by a normal distribution , where corresponds to the value of this parameter. These perturbations are independent and randomly generated for each -day forecast.

The last initial setting for data assimilation is the observation error covariance matrix . We assume that the covariance of the errors between different observations is zero, namely, the matrix is diagonal. Since we observe three states, is a 3 × 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 for any time , representing the observation error variance: Namely, for one has as 68% CI and as 95% CI. Considering that all the observations have been transformed into scale, that is, the observations in original scale are distributed in as 68% CI and as 95% CI.

Before showing the analysis result, a scatter plot of the values of and of for the different ensemble members on a certain day is provided in Fig. 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 Fig. 6, a positive correlation is observed between and . Thus we can use the observation of for the estimation of .

Fig. 6

Scatter plot of and on August 2, 2020.

In Fig. 7 we finally provide the analysis value for from March 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 is computed by using together with the information contained in (4) and in (5). The red curve shows the provided by Toyokeizai.net for reference. The gray 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 as indicated in Table 1. Additional discussions about this graph will be provided in the subsequent section.

Fig. 7

Analysis value for , and comparison with the value of Toyokeizai.net.

The analysis results for the compartments , , , , , with confidence intervals are also presented in Fig. 8. For compartments , , , the red dots represent the observations of Tokyo. For all pictures, the -axis is in scale.

Fig. 8

The 6 compartments for Tokyo.

Technical discussion, using data from Tokyo

In Fig. 7, one first notices that the similarity between the analysis and the reference provided by Toyokeizai.net. The Toyokeizai on day 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 [38]. Note that this formula is based on the daily confirmed cases, which experiences quick changes everyday. On the other hand, the analysis and its confidence interval is computed by formula (4) which involve the analysis ensemble of and . The average over the ensemble member provides a smoother evolution of , which is certainly closer to a true statistic’s evolution. Note that a third approach for the computation of is available in [39]: it involves an agent-based model together with a particle filter.

In Fig. 8, we show the analysis mean of compartments , , with the 68% CI (dark orange) and 95% CI (light orange), and the same results for compartments , and with the daily observations represented by red dots. One notes that the analyses of , , 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 and 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 with and 1.0, we run the experiments independently (always with the data from Tokyo) and the analysis mean of are shown in Fig. 9. The patterns of each are similar, especially after May , 2020. However, at the beginning of the outbreak and during the first state of emergency, all the analysis value of experience quick changes but with different speed. In Fig. 10 we also provide the analysis mean of the different compartments, since they also depend on the ratio between and . The same conclusion is obtained: except for the first three months, the possible ratio do not generate any noticeable difference.

Fig. 9

Comparison of analysis results with different ratio between and .

Fig. 10

Analysis mean of compartments with different ratio .

In Section 4.1, the error covariance matrix is defined as a diagonal matrix where is the identity 3 × 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 equals to with , , , , and . For these investigations the ratio between and is fixed as 0.58. In Fig. 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 Fig. 11(b), and one notices that the spread increases with the , but for larger than , the increase is less clear. We also show the analysis mean of compartments , , , , and in Fig. 12(a). Again, the analysis results of these compartments with different observation error covariance are quite close to each other. In Fig. 12(b), the spread is computed by the standard deviation of ensemble members of the transformed compartments. The relations of spread are quite clear for the three compartments with observations, and there are some overlaps for spread of , and , but in general the spread increases with .

Fig. 11

Analysis mean and spread of with different observation covariance.

Fig. 12

Analysis mean (a) and spread (b) of the compartments with different observation covariance.

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 cannot 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 and . Compared to the previous settings, in these two scenarios more people become asymptomatic and recover without showing any symptoms, as illustrated in Fig. 13(a). On the other hand, the three curves for are close to each others, but the one corresponding to only 50% of symptomatic is slightly higher than the other two, see Fig. 13(b). In order to understand this, let us recall that the relation between the transmission coefficients and 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, will also be bigger.

Fig. 13

Analysis mean of state (a) and (b) with different ratios between symptomatic and asymptomatic.

Discussion, comparisons, and future projects

Let us start by a few easy observations about the results obtained in the previous section.

The analysis values of quickly decrease during the first state of emergency, and this parameter is lower than around the middle of May 2020. This indicates that the disease has started to die out. During the second state of emergency, analysis values of had a much slower decay, but nevertheless in early February 2021, successfully drops below . However it starts to increase from March while the state of emergency had not been released. In the third state of emergency, the analysis had the slowest decay, compared to the previous two states of emergency. During the first several weeks of the fourth state of emergency, the effective reproduction number continued increasing, even at a faster speed. It is only after the end of the first week of August that a decay finally took place.

For , and , one finds that the three values increase at the beginning of the first state of emergency, then they slow down and change the direction to decrease. Similar behaviors can be observed in the second state of emergency, and the values start to increase around the middle of March 2021 while the starts to increase at the beginning of March. There is clearly a delay between the changes of and their effect on the three compartments , and .

Let us now provide a comparison with two different regions from Japan, namely Osaka and Kanagawa, and also provide the analysis for Japan as a whole. The experiments for Osaka, Kanagawa, and Japan start from March 2020, March 2020, March 2020 respectively. The observed compartments of these regions are the same as those for Tokyo, namely , and . The way to determine the initial setting is the same as that for Tokyo and we use the same setting for the parameters which have been presented in Table 1. The daily values for the parameter and for each region are also computed by using (18), (19). The ratio between and is 0.58, and the observation error covariance matrix is assumed to be .

In Fig. 14 we provide the analysis values of with 68% CI and 95% CI for each region. The additional red curves are again given by Toyokeizai.net. The gray shaded regions mark the respective periods of state of emergency. Though each region declared some states of emergency at different periods of time, they all have the common feature that the recent declarations are less efficient than the first ones. In Figs. 15, 16, and 17, we provide the analysis results for compartments , , , , and for each region.

Fig. 14

Analysis value of for Osaka, Kanagawa and Japan.

Fig. 15

Analysis results of Osaka.

Fig. 16

Analysis results of Kanagawa.

Fig. 17

Analysis results of Japan.

In Fig. 18, the three analysis curves for Tokyo, Osaka and Kanagawa are drawn simultaneously. The analysis for Japan is represented by red dots. One observes that the general trends for each region are similar, but some delays are also visible. For example in early April 2020, in Tokyo reached the first peak then started to decline, while in Kanagawa reached the first peak in the middle of April, and then start to decline. Subsequently, these two regions are quite correlated, but Osaka has a slightly independent behavior, both in the summer 2020 and in April 2021.

Fig. 18

for each region.

Based on these observations, one can suspect that the movement of populations might be related to the evolution of for regions which are close to each other. The study of the effect of movements can be carried out by constructing a proper model with multi-communities connected to each others. Some parameters may be tuned, for example for simulating movement restrictions. Such simulations could be useful for public health officials or for local governments to understand and estimate different disease-control strategies. We plan to work in this direction in the future.