The classic susceptible-exposed-infected-recovered (SEIR) model has been applied to case analysis of a variety of infectious disease outbreaks. Among them, susceptible groups, exposed groups, infected groups and recovered groups transform into each other. However, due to the high contagiousness and huge number of infections, the model seems not to work well when applied to the COVID-19 epidemic. Thus, some improvements are made in our model with reference to the work of Li et al. (15).

The infectious disease cases can be divided into two categories: confirmed cases, which are infection cases that are included in official statistics after infection; and unconfirmed cases, which are infection cases that were not included in official statistics because of the mild symptoms. To model an epidemic, the infection rate is one of the most important parameters. It is the number of people that one patient can infect in one day and represents the speed of transmission. Another variable with a similar meaning is the reproductive number R, which represents the total number of people a patient may infect. In our modeling, we have R = αβD + (1 − α)μβD, where α is the proportion of confirmed cases over all cases, β is the infection rate of a confirmed case in one day, μ is the reduction factor of the infection rate of an unconfirmed case compared to that of a confirmed case, and D is the average duration of an infection. When focusing on the early stages of an infectious disease outbreak, most studies using infectious disease models (15–17) assumed that the infection rate of the disease is a fixed value, which is reasonable during the free transmission stage. However, this assumption became problematic once local governments began to respond to the epidemic. It can be expected that after the free spread of infectious diseases in the early stages of the outbreak, the infection rate should continue to decline at a certain rate as the government's control measures advance and the public's awareness of prevention and control increases. Therefore, in studying the outbreak of the COVID-19 epidemic, we divide the spread of the epidemic into two stages: the first stage is the free spread stage, and the second stage is when the government implemented control measures. The infection rate as a function of time kept unchanged in the first stage, which is β₀, but gradually decreased in the second stage. It would be reasonable and convenient to assume that the infection rate of the epidemic continued to decrease at a fixed rate, which implies that the reproduction of the epidemic also decreased at the same rate. Thus, we set up the model M as follows:

where S_i, E_i, Iir, Iiμ, and N_i are respectively the susceptible population, exposed population, confirmed infected population, unconfirmed infected population, and total population in City i, β₀ is the basic transmission rate in the free transmission stage, τ is the attenuation rate of the daily infection rate, so the infection rate of a confirmed case at time t is β0τt, and Z is the average latency period.

The approximate solutions of the above simultaneous nonlinear equations can be obtained using the fourth-order Runge-Kutta method. It is noted that the values on the right-hand sides of these equations in Equation (1) are all determined by random samples from Poisson distributions with appropriate parameters to ensure robustness in the actual calculation.

To evaluate the performance of the proposed model in epidemic modeling, we examine the accuracy of model parameter estimates by setting up synthetic epidemics with different parameter values. We simulate two sets of bursts with different parameter values, one with larger parameter values τ = 0.8, α = 0.5, μ = 0.6, Z = 3.0, D = 3.0, denoted by Ξ₁, and the other with smaller parameter values τ = 0.6, α = 0.3, μ = 0.5, Z = 2.5, D = 2.5, denoted by Ξ₂. The initial number of confirmed cases is generated based on a discrete uniform distribution between 0 and 1000.

As mentioned in the previous section, before the first iteration of the algorithm, we need to select an initial range for each parameter. In light of Li et al. (15), we let the upper and lower bounds of S be S_max = S_min = s, and the lower bounds of E, I^r, and I^μ is Emin=Iminr=Iminμ=0, where s is the total population of a city. However, since our inferences involve the decay phase of the epidemic, it is no longer appropriate to set upper bounds on E, I^r, I^μ to fixed values as in (15). In order to be closer to the actual situation, based on the approximate ratios of E, I^r, and I^μ given in Li et al. (15), we let Emax=4×o0,obs, Imaxr=o0,obs, and Imaxμ=3×o0,obs using the most recent number of new daily confirmed cases o_0,obs. By taking into account the impact of reporting delays, the initial number of confirmed cases o_0,obs is set to the total number of confirmed cases in the two weeks closest to the first day of inference.

We now consider the global parameters. In our simulation studies, we find that the estimates of global parameters are sensitive to the initial ranges given beforehand. If we adopt the method in Li et al. (15), we need to set a relatively large initial range for each global parameter, which would result in that if Ξ₁ is the true set of global parameters, X_[g] can be well estimated, and if Ξ₂ is the true set of global parameters, only the estimates of α and D have small biases. After appropriately narrowing the initial range of each global parameter, we obtain an estimate of X_[g] with acceptable biases when its true value is Ξ₁ or Ξ₂.

The simulation results, shown in Table 1 and Figure 1, demonstrate that in the model we build, the algorithm can estimate the parameters well given an appropriate initial range. This fact allows us to apply it to real data in the framework of the model.

In this subsection, we apply the methodology in Section 2 to the 2020 outbreak in Chinese cities. Because many cities had too few infections then, which would lead to incorrect results, we only consider 136 cities with at least 20 cases there.

Since the Chinese government issued a national-level response on January 23, 2020, and implemented direct measures such as the lockdown of Wuhan, this date divides the 2020 epidemic into two time periods: the period before January 23, 2020, when the epidemic was in the stage of free transmission, and the period after January 23, 2020, the attenuation stage, where the spread of the epidemic was gradually restricted. Our analysis is mainly on the second stage, because there are already many studies on the first stage.

In the attenuation stage, the base value of the transmission rate is 1.12, which is the infection rate we have estimated in the first stage using the SEIR framework [consistent with Li et al. (15)], and the actual infection rate decreased at a fixed rate of τ over time, which, as explained in Section 2.1, represents the percentage reduction in the transmission rate per day. It is obvious that the smaller the τ, the faster the transfer rate will drop. Table 2 displays the global parameter estimates for fifteen cities.

By Table 2, it can be seen that the infectious attenuation rates of different cities are mostly in the range of [0.7, 0.8], which means that their reproductive numbers R_es drop below 1 after four days and falls to a low level (0.4) a week later. We believe this fact shows that the epidemic has basically been brought under control.

After obtaining the global parameter estimates, we generate the daily number of new confirmed cases by putting the global parameter estimates into our SEIR framework, where the state variables are generated as before. We repeat it 100 times and compute the average simulated daily number of new confirmed cases. Figure 2 displays the actual daily number of new confirmed cases and the fitted curve of the average simulated daily number of new confirmed cases. For comparison with our results, we also apply the modified SEIR model in Li et al. (15) to the epidemic in the corresponding city, and compute the average simulated daily number of new confirmed cases in 100 replicates. The results are shown in the Figure 2, from which it can be seen that the average simulated daily number of new confirmed cases of the modified SEIR model in Li et al. (15) is only in good agreement with the actual number of cases at the early stage of the outbreak, while the average simulated daily number of new confirmed cases from our varying coefficient SEIR model closely matches the number of daily new confirmed COVID-19 cases observed each day throughout the outbreak.

To show the impact of reporting delays, in Figure 3, we display fitted curves of the average simulated daily number of new confirmed cases by the varying coefficient SEIR model with and without considering reporting delays for Beijing and Shanghai. It can be seen that if the reporting delay is considered, the fitted curves are closer to the real data for both cities. Ignoring reporting delays will cause the estimated number of new confirmed cases to peak earlier than the peak of the true number of new confirmed cases, resulting in a non-negligible estimation bias.

After obtaining the estimate of τ, a common interest is to find out the factors associated with it. Generally speaking, people believe that cities with better economic development conditions are better able to control the epidemic. Strong local governments can effectively improve the efficiency of epidemic control. On the other hand, the attenuation rate may be affected by the location of the city, and some other factors. For demonstration purposes, we use the real data example from the previous subsection.

We consider the following factors: the total local GDP in 2019 and the number of permanent urban residents, which represent the impact of urban development on the attenuation rate; the GDP growth rate in 2019, which represents the capacity of local governments. In addition, we also measure the impact of geographical location by the city's latitude and longitude and the city's distance from Wuhan. Because the size of a local outbreak may affect the difficulty of controlling the outbreak, we include the total number of local cases as a factor.

Since the effect of the attenuation rate τ on the infection rate and the number of daily new confirmed cases is obviously nonlinear, it does not seem appropriate to assume that the attenuation rate is linear with local factors. To find out which local factors have non-negligible correlation with τ, we perform the Spearman correlation test. Considering that the attenuation rate τ is positive, we use logτ in the Spearman correlation test. The Spearman correlation test results between each local factor and logτ based on 136 cities are shown in Table 3, in which the local GDP in 2019, urban latitude and longitude, and urban population in 2019 were all respectively correlated with logτ, except for the significance level of urban longitude, which is <0.1, the significance levels are all <0.05. The remaining issue is deciding whether to include urban longitude as a local factor correlated with logτ. Note that the normality assumption for both logτ and urban longitude is acceptable at the 0.01 significance level after removing a small number of outliers (<5) by Shapiro-Wilk's test. Hence, we test the Pearson correlation between logτ and urban longitude and find that the Person correlation coefficient is 0.1998 and the p-value is 0.0197, indicating a non-negligible correlation between logτ and urban longitude. Therefore, we have included urban longitude in the analysis below as well.

It can be seen from our analysis that in those 136 Chinese cities, the attenuation rate has a tendency to increase as the regional longitude or latitude increases, which may be affected by the regional climate due to the corresponding changes in longitude and latitude. Furthermore, cities with better economic development tend to exhibit lower attenuation rates, consistent with general thinking that developed economy often indicates strong government execution and perfect public health system. More populous cities have lower attenuation rates, which, contrary to popular belief, may be due to the high correlation between urban population and the local GDP. The Spearman correlation analysis shows that the correlation between the urban population and the local GDP is as high as 0.79, indicating that the correlation between the urban population and the epidemic attenuation rate is affected by the local GDP, and hence urban population can be ignored when modeling the relationship between attenuation rates and other factors. Therefore, we fit a log-linear model of the attenuation rate based only on longitude, latitude, and local GDP. The results are shown in Table 4.

From Table 4, it is obvious that the local GDP and the longitude and latitude of a city are significant in the log-linear model. The parameter coefficients are –0.0347, 0.00284, and 0.00231, respectively, that is, for every 10 degrees increase in latitude or longitude, the regional attenuation rate increases by 1.0288 times or 1.0234 times. This may be caused by the climate difference between coastal and inland regions in China, similar conclusions can be found in the work of Wu et al. (25) and Srivastava (26). For every 1 trillion yuan increase in the region's total GDP, the attenuation rate would be reduced to 0.9659 times its original value, indicating that the local development level might have an impact on the speed of epidemic containment, which is consistent with Bambra et al. (27).

Although the p-value of less than 0.01 for the F-test indicates that the model is significant, the fitting performance of our log-linear model is still not good enough, with an R² of 0.216, suggesting that the explanatory variables we included in the log-linear model only explain part of the attenuation rate, and the rest may be local subjective factors such as the enthusiasm of the government and residents to control the epidemic. For city i, if we set the attenuation rate inferred by the varying coefficient SEIR model to τ1i and the attenuation rate estimated by the log-linear model to τ2i, then we can evaluate this city's ability to control the epidemic after eliminating objective factors by τ1i-τ2i. That is to say, for controlling the epidemic, City i is considered to have a strong ability to do so if τ1i-τ2i is low, otherwise, City i is considered to have a poor ability to achieve it.

From the real data analysis in Section 3.2 above, we can see that the varying coefficient SEIR model performs well in modeling China's first epidemic in 2020. We are also interested in seeing how the varying coefficient SEIR model performs on other epidemic data. Therefore, in this subsection, we apply the varying coefficient SEIR model to COVID-19 data from other regions in Asia. Subsequently, we also use a log-linear model to fit the attenuation rate. We then compare the attenuation rate estimate from the varying coefficient SEIR model to the prediction from the log-linear model.

An outbreak of COVID-19 infections occurred in Seoul, the capital of South Korea, in early August 2020, and the local government announced on August 16 an elevated quarantine response level to deal with the crisis. Local measures were as harsh as those in China: intercity transport was closed, public gatherings were banned, and schools were required to hold online classes. Using the varying coefficient SEIR model, the attenuation rate of the epidemic in Seoul after August 16 is estimated to be 0.7089, and the log-linear model predicts 0.7038. These two numbers are very close.

In early April 2020, a large-scale epidemic occurred in Tokyo, the capital of Japan, and the Japanese government declared a state of emergency on April 7, 2020. Based on the daily case data after April 7, 2020, the attenuation rate was 0.6869, compared to the log-linear model's estimate of 0.6309, which means that Tokyo was actually declining more slowly than it should be, most likely because Japan's control measures are much less stringent than those in China and South Korea. Instead of locking down the city, the government only restricted restaurant hours and public events, and urged people to stay indoors. We display the daily numbers of new confirmed cases and the fitted curves of the average daily numbers of new confirmed cases simulated by the varying coefficient SEIR model for outbreaks in Tokyo and Seoul in Figure 4.

India does not release city-specific epidemic data, so we have to use nationwide data to make a rough estimate. The second outbreak in India, which began in April 2021, was a major outbreak in the country. Many local governments announced strengthened prevention and control measures from May 6, 2020. Therefore, we regard this date as the dividing line of the epidemic. By using the approach in Section 2, India's attenuation rate after May 6, 2020, is estimated to be 0.6919. In terms of the GDPs, longitudes and latitudes of several major cities in India, India's attenuation rate would have been around 0.65, which implies that the prevention and control measures implemented in India were not effective enough. This result is consistent with the conclusions of other articles (8, 9) analyzing the effectiveness of control measures in various countries. We believe that during the 2020 epidemic, the control measures taken by China and South Korea were much more effective compared to other countries.