Here is an overview of the transmission mechanism of COVID-19 in the population (Figure 1):

It should be noted that the birth rate of newborns and the natural mortality rate of the population were ignored in this study.

The first subject is the susceptible population S(t). Because the birth rate of newborns was ignored, there is no source of the susceptible population, and the output includes contacts of susceptible people and asymptomatic undiagnosed people, asymptomatic diagnosed people, and symptomatic infected people. The infection rates are r₁, r₂, r₃, respectively.

Followed by the asymptomatic untested population A(t), whose sources are contacts among susceptible people and asymptomatic undiagnosed people, asymptomatic diagnosed people, and symptomatic infected people. Then, the input includes r₁S(t)A(t), r₂S(t)A_D(t), r₃S(t)I(t), and the outputs are the asymptomatic undiagnosed people will be diagnosed as asymptomatic people with probability a, be diagnosed as a symptomatically infected population b, and with probability c to develop into a cured population.

The source of the asymptomatic diagnosed population A_D(t) is the asymptomatic undiagnosed population that will be diagnosed as asymptomatic diagnosed population with probability a. The output includes the development of symptomatic infection population with probability d and the cure with probability g.

The sources of symptomatic infected population I(t) include the asymptomatic undiagnosed people with probability b to develop into symptomatic infected population and asymptomatic diagnosed people with probability d to develop into symptomatic infected population. The output includes part of the symptomatically infected people are cured with probability f₁, and part of them died with probability f₂.

The source of the recovered population R(t) includes three parts: first, the asymptomatic undiagnosed population is cured with probability c, the second is the asymptomatic diagnosed population is cured with probability g, and the third is the symptomatic infected population is cured with probability f₁.

Without considering the natural death of the population, the source of the death population D(t) is the development of the symptomatic infection population with probability f₂.

Therefore, this study proposed a mathematical model with Susceptible (S), Asymptomatic and undiagnosed (A), Asymptomatic and diagnosed (A_D), Symptomatic and infected (I), recovered (R) and death (D) groups.

In the above equation (1), S(t), A(t), A_D(t), I(t), R(t), D(t) represent the susceptible population, asymptomatic undiagnosed population, asymptomatic diagnosed population, symptomatic infected population, recovered population, and death population in the United States at time t, respectively.

We infer the model parameters based on the baidu data in United States from February 22, 2020 (day 1) to January 10, 2021 (day 324). The data from February 22, 2020 to December 1, 2020, which were used for the training set, are provided in Table 1 and were turned into fractions over the whole United States population (~3 ×10⁸) for simulations.

The estimated parameter values are based on the gathered data which are the number of currently infected individuals I(t), the number of diagnosed individuals who recovered R(t) and the number of death D(t) because of SARS-CoV-2 virus.

The nonlinear least square (NLS) method is regarded as the most basic way to estimate unknown parameters for ordinary differential equation model and facilitate to implement algorithm (24). Therefore, the NLS method was adopted to find the parameters that locally minimize the sum of the squares of the errors. During the course of the simulation, the parameters have been updated based on the successive measures at different stages.

Firstly, let P = (r₁, r₂, r₃, a, b, c, d, g, f₁, f₂), and I(t, P), R(t, P) and D(t, P) are the numerical solution of I(t), R(t) and D(t). Then, collect the data of currently infected individuals I₀, I₁, I₂, …, I_n−1, the number of diagnosed individuals who recovered R₀, R₁, R₂, …, R_n−1 and the number of death D₀, D₁, D₂, …, D_n−1. The initial values of six kinds of individuals are S₀, A₀, A_D0, I₀, R₀, D₀.

We assumed that I(t_i, P), R(t_i, P) and D(t_i, P) are the numerical solution at t_i and the parameter vector is P = (r₁, r₂, r₃, a, b, c, d, g, f₁, f₂). Next, we need to find the best parameter vector P~=(r~1,r~2,r~3,ã,b~,c~,d~,g~,f~1,f~2) to minimize the following equation:

Furthermore, based on previous literature, give the initial value of parameter vector P0=(r10,r20,r30,a0,b0,c0,d0,g0,f10,f20) and set the bound. Use the fmincon function to estimate the approximate range of each parameter and the estimated parameters were regarded as new initial values. Lastly, the lsqnonlin function was employed to achieve the best fitting effect and get the optimal parameters.

Similar to our previous research (25), extended Fourier Amplitude Sensitivity Test (EFAST) was adopted to search the sensitive parameters. The first order sensitivity index (S_i) and full order sensitivity index (S_Ti) are calculated by the following equations:

Based on reference (26), the sensitivity value (S_Ti) larger than 0.5 was defined as a sensitive parameter; otherwise, it was defined as a non-sensitive parameter.

The basic reproduction number reflects the size of the virus transmission capacity. The larger the basic reproduction number, the stronger the virus transmission ability; the smaller the basic reproduction number, the weaker the transmission ability. Therefore, the study on basic reproduction number is very necessary.

Based on equation (1), the equilibrium (S*,A*,AD*,I*,R*,D*)=(S*,0,0,0,R*,D*) was calculated firstly, and the Jacobian matrix around the equilibrium is:

The corresponding Jacobian determinant is:

Let

where

Therefore, we get G(λ) is as follows:

The system is positive and ||G(λ)||∞=G(0)=N(0)D(0), where ||G(λ)||_∞ is the H_∞ norm of G(λ).

Based on Hurwitz criterion (27) and conference (16), all roots in the left-hand plane if, and only, if S^*<G(0). Therefore, the basic reproduction number is:

The attack rate denotes the proportion of all infected individuals to the total population, and the infected individuals include the asymptomatic and undetected, the asymptomatic and detected, and the symptomatic and infected. The method for calculation of the attack rate is:

Since the development trend of the epidemic varies with the different control strategies adopted by countries, the parameters of model were estimated by piecewise fitting method in this study. Based on Baidu data and the above NLS algorithm, the model parameters obtained by piecewise fitting are as follows:

The parameters from days 1 to 50 are:

r₁ = 0.049477, r₂ = 1.3392e−05, r₃ = 0.27981, a = 0.23686, b = 0.71424, c = 0.0233352, d = 0.79452, g = 0.0040538, f₁ = 0.0061802, f₂ = 0.009, and the basic reproduction number is R₀ = 6.0064.

The parameters from day 51 to day 100 are:

r₁ = 0.049477, r₂ = 1.3392e−05, r₃ = 0.27981, a = 0.13686, b = 0.009424, c = 0.007352, d = 0.007452, g = 0.0018538, f₁ = 0.0061802, f₂ = 0.003, and the basic reproduction number is R₀ = 6.2908.

The parameters from days 101 to 200 are:

r₁ = 0.049477, r₂ = 1.3392e−03, r₃ = 0.0131, a = 0.13686, b = 0.008424, c = 0.0053352, d = 0.003052, g = 0.0018538, f₁ = 0.0061802, f₂ = 0.0047, and the basic reproduction number is R₀ = 1.8003.

The parameters from days 201 to 250 are:

r₁ = 0.049477, r₂ = 1.3392e−03, r₃ = 0.0131, a = 0.13686, b = 0.004424, c = 0.0053352, d = 0.003052, g = 0.003538, f₁ = 0.006, f₂ = 0.003, and the basic reproduction number is R₀ = 1.4888.

The parameters from day 251 to day 324 are:

r₁ = 0.049477, r₂ = 1.3392e−03, r₃ = 0.0281, a = 0.13686, b = 0.00424, c = 0.0123352, d = 0.0138, g = 0.00538, f₁ = 0.0062, f₂ = 0.0032, and the basic reproduction number is R₀ = 3.2698.

Based on the parameter estimation of the above five stages, the fitting results are obtained, as shown in Figure 2 (from February 22, 2020 to December 1, 2020, 284 days in total). The black curve represents the simulation result of the mathematical model, and the red curves are the collected data of the current infected, recovered and death, respectively. The correlation between the number of current infected simulated by the mathematical model and the collected was 98.85%, the correlation of the recovered was 99.84%, and the correlation of the death was 99.54%.

In addition, we have also calculated the basic reproduction number in 5 stages. The basic reproduction number in the first stage is R₀ = 6.0064, and the basic reproduction number in the second stage is R₀ = 6.2908, indicating that the transmission capacity of the epidemic is very strong in the first two stages. If it is not controlled, it will quickly spread throughout the United States. In the third stage, the basic reproduction number is R₀ = 1.4888. Compared with the first two stages, the transmission capacity of COVID-19 has been significantly reduced. In the fourth stage, the basic reproduction number further drops to R₀ = 1.4888, indicating that with joint efforts of the US government and the people, the spread of the epidemic has further declined. In the fifth stage, the basic reproduction number has risen again to R₀ = 3.2698, indicating that the transmission capacity of COVID-19 has increased, and that it is necessary to take measures to prevent and control it.

Next, in order to verify the accuracy of the model established in this research, we continued to collect data from December 2, 2020 to January 10, 2021 (40 days in total) (marked with blue curves). The data are provided in Table 2.

The red curves are from February 22, 2020 to December 1, 2020, a total of 284 days, and the black curve is the simulation result of the model, as shown in Figure 3. We use the following equation to calculate the accuracy of model prediction:

where x_i is the real value, y_i is the forecasting value, and n is the number of all data that needs to be forecast.

The results show that the prediction accuracy rate of the current infected population I(t) is 98.33%, the prediction accuracy rate of the recovered population R(t) is 98.2%, and the prediction accuracy rate of the death population D(t) is 98.96%. This indicates that the model established in this study has high accuracy, and that the model can accurately predict the trend of the epidemic in the United States, and provides a guarantee for the following applied research.

Therefore, without taking into account the vaccine and the recovery of the survivors, we used the above model to make a long-term prediction of the epidemic in the United States (as shown in Figure 4). The prediction results showed that the susceptible population continued to decline, reaching the lowest level around the 1,560th day (4.45% of the total population) and continued until the end. The asymptomatic undiagnosed population reached a peak on the 518th day (1.53% of the total population), and then gradually decreased, dropping to 0% on the 1367th day. The asymptomatic diagnosed population reached the peak on the 570th day (10.32% of the total number), and gradually decreased, and the dropped to 0% on the 1,440th day. The infected population reached a peak around the 679th day (17.92% of the total number), gradually decreased, and then dropped to 0% around the 1,852th day. The number of recovered people was increasing, reaching a peak around the 1,533th day (91.87% of the total number), and was stable at this value. The number of deaths also increased gradually, reaching a peak around the 1,619th day (3.16% of the total number), and remained stable.

In addition, the attack rate was also studied and shown in Figure 5. The numerical simulation result indicated that the attacked individuals increased from the beginning and reached a peak around the 624th day, and that the maximum attack rate was 0.281, which is different from the asymptomatic and undetected, the asymptomatic and detected, and the symptomatic and infected.

Next, parameter sensitivity analysis was performed to screen the factors that significantly affected the development of the epidemic. In order to overcome the influence of coupling between parameters on the results, we adopted a robust, highly calculated, and low-sample required global sensitivity analysis method, the Extended Fourier Amplitude Sensitivity Test (EFAST) method.

Based on the method of global sensitivity analysis, we obtained the first-order and full-order sensitivity indexes of 10 parameters in the mathematical model (Table 3 and Figures 6A,B). Based on reference, the sensitivity value (S_Ti) higher than 0.5 was defined as a sensitive parameter; otherwise, it was defined as a non-sensitive parameter (26).

Therefore, we theoretically analyzed the parameters that affect the epidemic. The sensitivity parameters are r₁ = 0.8317 (the infection rate of asymptomatic undiagnosed population to susceptible population), r₃ = 0.6427 (the infection rate of symptomatic infection population to susceptible population), d = 0.6163 (probability of asymptomatically diagnosed population turning into symptomatically infected population), g = 0.521 (recovered rate of asymptomatic diagnosed population), and f₁ = 0.8247 (recovered rate of symptomatically infected population), non-sensitive parameters are r₂ = 0.2353 (the infection rate of the asymptomatic diagnosed population to the susceptible population), a = 0.1664 (the probability of the asymptomatic undiagnosed population becoming asymptomatic diagnosed population), b = 0.1731 (the probability of asymptomatic undiagnosed population becoming symptomatically infected), c = 0.1757 (the recovered rate of asymptomatic undiagnosed population), and f₂ = 0.169 (the death rate of symptomatic infection).

In the following, we used numerical simulation method to verify the correctness of the above global sensitivity analysis results. That is, the gradient change of parameters is used to study the impact of the output results of the susceptible population, the asymptomatic undiagnosed population, the asymptomatic diagnosed population, the symptomatic infected population, the recovered population, and the death population. For example, the initial value of parameter r₁ is 0.05, and it increases to 0.1 in steps of 0.01 to study the dynamics trend of each type of population.

Numerical simulation results showed that when parameter r₁ increases from 0.05 to 0.1 in steps of 0.01, the number of susceptible populations is greatly reduced, and that the number of asymptomatic undiagnosed population, asymptomatic diagnosed population, symptomatic infected population, recovered population, and death population all increases significantly with increase in r₁ (as shown in Supplementary Figure 1). This is consistent with the sensitivity of r₁, calculated by our theory, of 0.8317 (Table 3). Therefore, the parameter r₁ is a sensitive parameter that affects the epidemic.

When parameter r₂ increases from 0.0013 to 0.0018 in steps of 0.0001, the number of susceptible people, asymptomatic undiagnosed population, asymptomatic diagnosed population, symptomatic infected population, recovered population, and death population has no significant changes (Supplementary Figure 2). This is consistent with the sensitivity of parameter r₂ of 0.2353 (Table 3), indicating that parameter r₂ is an insensitive parameter, and that the change in parameter r₂ has no significant effect on the epidemic.

When parameter r₃ increased from 0.016 to 0.021 in steps of 0.001, the number of susceptible populations significantly reduced, and the number of asymptomatic undiagnosed population, asymptomatic diagnosed population, symptomatic infected population, recovered population, and death population all increased significantly with the increase in r₃ (Supplementary Figure 3). This is consistent with the sensitivity of r₃ of.6427 (Table 3). Therefore, parameter r₃ is a sensitive parameter that affects the epidemic. The sensitivity of parameter r₃ (0.6427) is lower than that of parameter r₃ (0.8317). The numerical simulation results also verify that the change in parameter r₁ has a greater impact on the epidemic than parameter r₃, which preliminarily proves that the parameter sensitivity obtained by our theoretical calculation is correct.

When parameter a increased from 0.13 to 0.18 in steps of 0.01, except for a certain degree of change in the asymptomatic undiagnosed population, the number of susceptible people, asymptomatic diagnosed population, symptomatic infected population, recovered population, and death population did not change significantly (Supplementary Figure 4), which is consistent with the sensitivity of parameter a of 0.1664 (Table 3), indicating that the parameter is a non-sensitive parameter with no significant impact on the epidemic.

When parameter b increased from 0.004 to 0.014 in steps of 0.002, the number of susceptible people, asymptomatic undiagnosed population, asymptomatic diagnosed population, symptomatic infected population, recovered population, and death population had no significant changes (Supplementary Figure 5), indicating that parameter b is a non-sensitive parameter, which is consistent with the sensitivity of parameter b of 1731 (Table 3).

Similarly, when parameter c increased from 0.005 to 0.01 in steps of 0.001, there was no significant change in the 6 groups of people (as shown in Supplementary Figure 6), indicating that parameter c is a non-sensitive parameter, which corresponds to the sensitivity of parameter c of 0.1757 (Table 3).

When parameter d increased from 0.013 to 0.018 in steps of 0.001, there were significant changes in the asymptomatic undiagnosed population, symptomatic infected population, and death population (Supplementary Figure 7), indicating that parameter d is a sensitive parameter, which is consistent with the sensitivity of parameter d of 0.6163 (Table 3).

When parameter g increased from 0.004 to 0.009 in steps of 0.001, in addition to the cured population, the number of susceptible populations, asymptomatic undiagnosed populations, asymptomatic diagnosed populations, symptomatic infected populations, and death populations all changed significantly (Supplementary Figure 8). It means that parameter g is a sensitive parameter that affects the epidemic, which is consistent with the sensitivity of parameter g of 0.521 (Table 3).

When parameter f₁ increased from 0.006 to 0.011 in steps of 0.001, the number of susceptible populations, asymptomatic undiagnosed population, asymptomatic diagnosed population, symptomatic infected population, and death population all changed significantly, except for the recovered. Besides, the amplitude of change is greater than the parameters d and g (Supplementary Figure 9). This is because the sensitivity of parameter f₁ is 0.8247, which is significantly greater than the sensitivity of parameters d and g (0.6163 and 0.521) (Table 3).

When parameter f₂ increased from 0.0002 to 0.0007 in steps of.0001, except for the obvious changes in death population, there was no significant change in the susceptible population, asymptomatic undiagnosed population, asymptomatic diagnosed population, symptomatic infected population, and recovered population (Supplementary Figure 10), indicating that parameter f₂ is a non-sensitive parameter, which is consistent with the sensitivity of parameter f₂ (0.169) (Table 3).

In summary, we have verified the sensitivity indexes of the 10 parameters obtained by theoretical calculations through numerical simulation. The results showed that gradient change in the sensitive parameters can significantly affect the development of the epidemic, and that change in the insensitive parameters has little impact on the epidemic (only a small impact on a certain group of people or no significant impact at all). Thus, the accuracy of the parameter sensitivity analysis was proved. In addition, the higher the sensitivity of the parameter, the greater the impact of parameter changes on the epidemic.

In addition, the influence of parameter changes on attack rate was also studied. The simulation results (Figure 7A) showed that when increase parameter r1 from 0.05 to 0.1 in steps of 0.01, the attack rate increases significantly. which was shown in Figure 7A, where the color and black arrow have the same meaning as above. Similarly, increase r₃ from 0.016 to 0.021 in steps of 0.001 (Figure 7C), increase g from 0.004 to 0.009 in steps of 0.001 (Figure 7H), and increase f₁ from 0.006 to 0.011 in steps of 0.001 (Figure 7I), the attack rates are all affected significantly. What the difference is that with the increase in r₁ and r₃, the attack rate also increases, and when parameters g and f₁ increase, the attack rate decreases significantly. This is consistent with our sensitive analysis.

Besides, if r₂ increases from 0.0013 to 0.0018 in steps of 0.0001, the attack rate has no significant difference (Figure 7B). Similarly, increase a from 0.13 to 0.14, 0.15, 0.16, 0.17 and 0.18 (Figure 7D), increase parameter b from 0.004 to 0.006, 0.008, 0.01, 0.012, and 0.014 (Figure 7E), increase parameter c from 0.005 to 0.006, 0.007, 0.008, 0.009, and 0.01 (Figure 7F), there was no significant change in their attack rates. The non-sensitive parameters are consistent with our theoretical analysis.

However, when d increases from 0.013 to 0.018 in steps of 0.001, the change is not particularly dramatic. This is because with the increase in d, the asymptomatic undiagnosed population and symptomatic infected population increased, and the asymptomatic diagnosed population decreased.