Model framework

This is a computational study. Based on the existing knowledge of the clinical progression of COVID-19 infection, we developed a new transmission model by using ordinary differential equations and then applied this model to simulate disease evolution. Humans are classified into different categories according to their infection status: susceptible (S), latent (E), asymptomatic infection (I _b), clinical infection (I _c), vaccinated (V), recovered (R) and dead (D). The sum of these classes (except D) equals the total population size, i.e. N = S + E + I _b + I _c + Ic + R + V.

Newly infected individuals are generated on account of exposure to infectious individuals. A susceptible population (S) can be infected at rate λ ₁ through contact with infectious individuals (I _b and I _c). It then transforms into a latent state (E) and becomes infectious after an incubation period 1/η. Parts of them (accounting for β proportion) become clinically infected, and the others become asymptomatic. People can acquire immunity through vaccination and then move from susceptible state to vaccinated state. Known that vaccine only reduces the risk of infection, those vaccinated population (V) will be infected at a lower rate λ ₂ after exposure to infectious population. Asymptomatic infected individuals (I _b) are not easily detected and treated, but they can recover after period 1/δ. Clinically infected individuals (I _c) recover after an average treatment period 1/γ. The corresponding flow chart is shown in Figure 1, which is formulated by a differential system in the model (1).

(1)$$\left\{\matrix{\displaystyle{{dS( t ) } \over {dt}} = A-\lambda_1f\displaystyle{{S( t ) } \over {N( t ) }}( {I_c( t ) + \alpha I_b( t ) } ) -( {\vartheta + d} ) S( t ) , \hfill \cr \displaystyle \hskip-1pc {{dE( t ) } \over {dt}} = \lambda_1f\displaystyle{{S( t ) } \over {N( t ) }}( {I_c( t ) + \alpha I_b( t ) } ) + \lambda_2f\displaystyle{{V( t ) } \over {N( t ) }}( {I_c( t )} \cr \quad\quad\quad\quad{+\; \alpha I_b( t ) } ) -( {\eta + d} ) E( t ) , \;\hfill \cr \displaystyle{{dI_b( t ) } \over {dt}} = ( {1-\beta } ) \eta E( t ) -( {\delta + d} ) I_b( t ) , \;\hfill \cr \displaystyle{{dI_c( t ) } \over {dt}} = \beta I_b( t ) -( {\gamma + \mu + d} ) I_c( t ) , \;\hfill \cr \displaystyle{{dR( t ) } \over {dt}} = \delta I_b( t ) + \gamma I_c( t ) -dR( t ) , \;\hfill \cr \displaystyle{{dV( t ) } \over {dt}} = \vartheta S( t ) -\lambda_2f\displaystyle{\vartheta \over {N( t ) }}( {I_c( t ) + \alpha I_b( t ) } ) -dV( t ) . \hfill} \right.$$

Fig. 1. Flow diagram of COVID-19 transmission.

The combined effect of vaccine $\vartheta = \theta v$ is the product of the vaccination rate v and the vaccine potency θ. Here our model equalised the effectiveness of different vaccines in the world and regarded them as a whole.

The adaptive behaviours are considered in the following ways. During the COVID-19 pandemic, people usually change their behaviour and take preventative steps to reduce infection risk, such as social distancing and wearing masks. Individuals' self-protection awareness would intensify as the infection cases increase, and then they would adopt stricter measures against infection. We called such performances as adaptive behaviours. They are a collective designation of protection actions which are subject to people's self-protection awareness. Given that adaptive behaviours can reduce infection risk, they are quantified as a coefficient of infectivity, which is a decreasing function of case number, written as f(I,  R). According to the expression given in reference [Reference Eksin17], we set f = (1 − (I + ϕR)/N)^k. Parameter k measures the magnitude of adaptive behaviour, in which k = 0 indicates no adaptive behaviour, and a larger k indicates that people are more sensitive and cautious to the infection situation. Parameter ϕ ∈ [0,  1] represents the sensitivity of humans to the number of historical cases, where a larger ϕ indicates a stronger long-term awareness of protection in the population.

Parameter setting and data analysis

The total population of the world is set to be N = 7.898 billion. Birth constant A was estimated by using model dN/dt = A − dN to fit the world population during 2012 and 2021. The vaccination rate v is determined by the proportion of the average vaccinated population. The vaccine efficacy is first set to be θ = 60% for validating the model, and then we changed its value and run the model for observing the outcomes of different vaccine efficacies. The numbers of vaccinated people and reported cases per day were extracted from the website (see Supplementary Material) [18]. The time span is from 22 January 2020 to 18 July 2022, with a total of T  =  909 days. Considering the existence of changes in population behaviour and virus variation during epidemic evolution, infection rates (λ ₁,  λ ₂) and adaptive parameters (ϕ,  k) vary during three periods with different main strains. Model parameters are shown in Tables 1 and 2.

Table 1. Parameter introduction

Table 2. Parameter changes before and after vaccination

We applied the proposed model to fit the global COVID-19 infection data from 22 January 2020 to 18 July 2022. To take into account the different infectivity from different viral strains, we divided the time period into three segments according to the changes of main virus strains: initial strain, Delta and Omicron, and then we estimated their respective infectivity. MCMC algorithm was used to estimate unknown parameters (including k, λ ₁,  λ ₂,  ϕ and the initial values of E), which was achieved by using the FME package in R. Considering the protective effect of vaccine, we assumed that vaccination halves the infection rate, that is λ ₂ = 0.5λ ₁. The initial values of model variables were set equal to the reported number in initial time. We examined the model performance quantified by minimum fit error criterion. The error is determined by ${\rm Error} = \sum\nolimits_{t = 1}^n {{( {{\hat{C}}_t-C_t} ) }^2/n}$, where C _t ($\hat{C}_t$) is the numbers of reported cases (modelling cases) on day t.