A reductive analysis of a compartmental model for COVID-19: data assimilation and forecasting for the United Kingdom

I. INTRODUCTION

The highly contagious SARS-CoV-2 has infected more than six million people worldwide since its first detection in China on December 31 [1]. The novel coronavirus is the fourth wave in the class of coronaviruses. In less than two months, the virus has spread all over the world, posing serious threats to healthcare systems and economies. The alarming speed of transmission, the virulence of the disease, and the unprecedented high proportion of fatalities even in countries with high healthcare indices have raised questions about what kind of interventions are appropriate for a given setting. The wide variability in infected numbers and fatalities in different countries and settings has also brought into sharp focus a debate about the underlying causes of the variability. In the absence of any treatment for the disease and non-availability of vaccines in the near future, policy makers have resorted to standard epidemiological interventions, such as social distancing, isolation, contact tracing, and quarantining, and more recently a complete lock-down.

At a basic level, the purpose of all non-pharmacological interventions is to control disease transmission by limiting the proportion of the population exposed to the virus as much as possible. Furthermore, inherent in the process of implementation of these interventions are delays at each stage. The delay time-scales are specific to the particular intervention.

The importance of mathematical models describing the spreading dynamics of infectious diseases has been recognized since early days [2]. In particular, the fact that timely models that include realistic features have often been helpful in decision making on healthcare issues is well recognized [2–4]. In the short period since the emergence of the coronavirus, there have been several mathematical models [5–17], to name a few. Some of these models attempt to evaluate the contribution from different transmission routes, such as contact tracing and isolation [10, 11, 18], travel restrictions [16, 17] social distancing [19, 20], lock-down measures [16, 21, 22], and a combination of these interventions [12, 23–26]. These models broadly fall into three categories, deterministic, stochastic, and simulations. Several new mathematical techniques used in different disciplines have been employed to gain insights, which would not be possible with the traditional approaches in the field. These include the human mobility model [16], differential evolution[19], heuristic optimization technique [19], stochastic agent-based discrete time simulation [27], supply chain risk simulations [28], etc.

One class of epidemiological models attempts to describe the transmission dynamics by partitioning the population into smaller subsets based on the disease status such as the susceptible, exposed, infected, quarantined, recovered, etc [5–9]. Most models of this kind ignore age-dependent infection and fatality rates, and the heterogeneous spatial distribution of the population. In a sense, these models describe the evolution of the mean response of each type of population. Despite these limitations, these models have the ability to include several realistic features.

In the compartment type of models, the disease status of individuals changes with the development of the disease, i.e., transitions occur between two compartments either due to interaction of the infected with the susceptible or due to interventional actions. These models include delay time-scales inherent in the dynamics of transmission, for instance, the period spent in quarantine and the time required for tracing individuals exposed to the infected. These models have the ability to include several realistic features, such as the response of the population to interventional measures. However, generally, inclusion of more and more realistic features requires a larger number of partitions. Then, the number of differential equations increases and so does the number of parameters, making calibration of the parameters difficult [5–9, 13].

Motivated by the complexity of such models, we have devised a compartment-based model having the susceptible, infected, quarantined, traced, recovered, and the deceased populations. The susceptible and the infected form the core populations in the sense that it is through these two populations that inward/outward transitions occur with other populations. However, even this simple model contains several parameters, leaving numerical solutions as the only option for further analysis. Considering the fact that the model equations form a coupled set of nonlinear differential equations, we adopt mathematical methods drawn from nonlinear dynamical systems most appropriate for analyzing such equations [29]. We hypothesize ‘accessible population for infection’ to be a small fraction of the total population. The validity of this hypothesis can be appreciated by noting that the purpose of interventions is to minimize the exposure of the population to virus transmission, thereby limiting the spread of infection. We further assume that the order of magnitude of the accessible population is similar to that of the infected population. This assumption is made more quantitative. This, together with the structure of the model equations, allows us to decouple them into two equations. These two equations further reduce to a logistic type of equation for the total infected population with well defined parameters namely, the ‘testing rate’ and ‘contact transmission rate’ parameters [30, 31]. The equation can be solved analytically, thereby allowing for a clear interpretation of the parameters controlling the growth and inhibiting factors. The validity of the ‘accessible population’ hypothesis and the efficacy of the reduced logistic model is demonstrated by the ease of fitting the cumulative number of infections and daily new cases for the United Kingdom (UK). The procedure further allows us to forecast the progression of the disease. Using this information and calibrating the relative importance of the various transition rates (equivalently the associated parameters), we optimize the parameter values specific to the UK. This procedure allows us to delineate the various time scales contributing to the various regions in the time development of the disease. Using this, we numerically solve the full model equations. The calculated total infected population fits very well with the available data for the UK [32]. (UK does not publish data on the recovered and the active populations.) The model exhibits a turning point in the active infected population around May 15. However, since the rate of slowing down beyond the turning point is poor, the projected end time of the epidemic would be around late July with the predicted saturation level of the total number of infections ~ 3.52 × 10⁵ assuming lock-down conditions continue.

II. THE MODEL

The total population N is partitioned into the susceptible S, active infected I, quarantined Q, those traced T after being exposed to the infected, recovered R, and the deceased D. The respective populations are denoted by N_s, N_i, N_q, N_tr, N_r and N_d.

Testing is one of the standard protocols used for identifying the infected. If α_s is the rate of testing per day per million and p_s is the probability of testing positive, then α_sp_sN_s is the transition rate from S to I. Infected individuals coming into contact with the susceptible class can transmit the virus. If β_i is the transmission rate per contact, p_i is the probability of transmission of the disease and F(d_i) is a distance-dependent interaction, then, p_iβ_iF(d_i)N_iN_s is the transition rate from S to I. Considering the fact that one of the primary routes of transmission is through airborne aerosols generated by the infected, larger separation is known to reduce the risk of transmission [19, 20, 33]. This distance-dependence of F(d_i) is generally expressed as . However, in the present context where we will be dealing with a lock-down situation for most part of the progression of the disease, we set F(d_i) = 1.

During testing, some individuals would always exhibit mild or ambiguous symptoms. These are identified as pre-symptomatic. If the probability of finding the pre-symptomatic is p_q, then, α_sp_qN_s transition out of S to Q. Subsequently, when tested again, say after a quarantine duration [5, 34, 35], some of them may either test positive with a probability q₁ or negative with a probability (1 − q₁). If positive, the transition out of Q (to I) is q₁λ_qN_q. Here, 1/λ_q is the quarantine duration, usually of the order of the incubation period [34, 35]. Similarly, if tested negative, the transition rate out of Q into S is (1 − q₁)λ_qN_q. The total loss rate to N_q is λ_qN_q.

Tracing those exposed to the infected, and testing to find if they are infected, are important steps in controlling the spread of the disease. Inherent in tracing such individuals are delays in tracing. Such delays cause increased transmission of the disease. If p_t is the probability of tracing such individuals, then, α_tp_tN_s is the transition rate from S to T. Subsequently, individuals testing positive will move to the infected compartment I with a probability q₂ and the rest with a probability (1 − q₂) move to S. The total transition out of T is equal to λ_tN_tr, where 1/λ_t is the time taken to trace the individuals. (There is also another possibility, namely, some individuals may show mild symptoms. Then, there would be a transition into Q. For the sake of simplicity, we have ignored this route.) Finally, the outward transitions from I are the recovery and death rates respectively, γ_r N_i and κ_dN_i.

Collecting these terms, we have the following set of coupled nonlinear ordinary differential equations

(Here, we have suppressed F(d_i) factor since it has been set equal to unity. Definitions of the notations are given in Table I). Note that the total infected population is given by N_t = N_i + N_r + N_d.

To begin with, we highlight a few features of the model equations. Our model, much as other compartment-type models, has several parameters. However, several of these are directly measurable and therefore can be obtained from the literature. A few others are related to testing protocols and again can be obtained from the literature or from relevant open sources [32]. For instance, α_sp_s, α_sp_q and α_tp_t are directly related to testing rates and therefore, these are known for a given situation. A few other rate parameters such as λ_q, λ_t, λ_r and κ_d, are inversely related to measurable time-scales, such as the duration of quarantine τ_q, time required for tracing τ_t, time for recovery starting from illness τ_r and the time from illness to death τ_d, respectively [35, 36].

The present model includes two delay loops defined by Eqs. (3) and (4). These delays are natural to the implementation of the protocols. For instance, once quarantined, subsequent tests are conducted after quarantine duration to identify if quarantined individuals test positive or negative. Similarly, delays in tracing individuals are common. A more transparent way to describe these delay loops is through the integral representation of Eqs. (3) and (4), which forms the definitions of the two populations N_q and N_tr, respectively. For instance,

When the kernel K(t) is modeled using an exponential form with a single time scale 1/λ_q, i.e., , one can easily verify that differentiating Eq. (7) (using the Leibniz rule) leads to Eqs. (3). The convoluted nature of the integral physically implies that those quarantined earlier will leave the quarantine sooner than those quarantined later.

Equations (1-6) constitute a set of coupled nonlinear differential equations. A standard procedure for further analysis of such equations is through numerical integration. While this is a necessary step, here we adopt a reductive analysis of the model equations by exploiting the fact that there are two main populations, namely, the susceptible (Eqs. 1) and the infected (Eq. 2). Furthermore, Eqs. (5,6) are essentially decoupled from the rest (transitions to R and D are from I). These two features suggest that Eqs. (1,2) can be decoupled from the rest of the equations. We refer to the decoupled equations as the reduced model equations. Since the two equations can be further reduced to a logistic-type equation (referred to as the reduced logistic equation), it can be analytically solved. As we shall see, analysis of this equation provides insights that prove to be useful for the analysis of the full model Eqs. (1-6). (We shall often refer to Eqs. (1-6) as “full model equations” to avoid confusion.)

III. CONCEPT OF ACCESSIBLE POPULATION: THE REDUCED MODEL

We now introduce the concept of ‘accessible population for transmission of the disease’. To appreciate this concept, consider the spreading dynamics of a contagious disease in the absence of any interventions. Then, in principle, the entire population is exposed to the disease, and it may spread to the entire population (barring the possibility of the population acquiring herd immunity). In this case, the entire population is the accessible population. However, since no Government would like to see the entire population infected, interventional measures are enforced precisely to mitigate the risk of transmission and limit the population exposed to the disease to a minimum. In this case, the accessible population is expected to be a small fraction of the total population. These two limiting cases can be accommodated by hypothesizing the accessible population to be a fraction of the total population, where the fraction is determined by whether interventional measures are imposed or not.

Consider dropping all terms except α_sp_sN_s and p_iβ_iN_iN_s in Eqs. (1,2). Then, these two equations get decoupled from the rest of the equations. Further, because all other inward/ outward transitions are removed, the character of the compartment I changes from the active infected to the cumulative infected I_t with N_t denoting the corresponding population. Then, we have

Noting that we have N_t + N_s = constant. Without loss of generality, we set N_t + N_s = N_s(0), the total population. Then, we get a single equation governing the cumulative infected population, given by

Equation (11) has the well known form of logistic equation, extensively studied in the context of population dynamics [37], with a notable difference, namely, the parameters a, b and c have a well defined interpretation as discussed above. We refer to Eq. (11) as the reduced logistic equation. (For brevity we often refer to α_sp_s and α_iβ_i as testing and contact transmission rates, respectively.)

We begin with a few observations on the relative magnitudes of the model parameters in the absence and the presence of interventions. Consider a situation when there are no constraints. Then, one should expect that the testing rate (α_sp_s) to be low (compared to the lockdown period) due to absence of any guidelines from policy makers. Similarly, since infected individuals carry on with their routine activity, the number of contact transmissions is high and hence, the contact transmission rate (p_iβ_i) is expected to be high (compared to when interventions are in place). Then, the total accessible population denoted by N_a(0) is the entire population of the region or the country, i.e., N_a(0) = N_s(0). In contrast, when interventions are in place, testing rates are high to ensure identification of the infected, therefore, α_sp_s is high. In this situation, since the mobility of individuals is restricted, the number of contacts is severely limited, i.e., p_iβ_i will be small. Therefore, the accessible population N_a(0) is expected to be small compared to the total population N_s(0). These two limiting cases can be written as . These qualitative statements about the accessible population will be made quantitative by carrying out a detailed analysis of Eq. (11).

Consider the initial growth of Eq. (11) by dropping the quadratic term. Then, we have

The solution is given by where N_t(0) is the initial number of infections. As can be seen, the growth rate given by b ≈ p_iβ_iN_s (0) depends on N_s(0), the total population. Therefore, the growth rate can be high. In addition, the pre-factor for the exponential growth term (in Eq. 16) depends not only on N_t(0) but also on c/b = α_sp_s/p_iβ_i. Thus, the initial growth depends on relative magnitudes of N_t(0) and a_sp_s/p_iβ_i.

It is straightforward to obtain the solution of Eq. (11). (See Appendix for details.) Here, it is adequate to consider the solution in terms of the parameters a, b, and c, given by

We now examine two limiting cases. For short times, N_t tends to (since the denominator is dominated by b/a = N_s(0)), consistent with the short time solution given by Eq. (16). For long times however, N_t tends to b/a = N_s (0), the total population.

The self-limiting nature of Eq. (17), a characteristic feature of logistic equations, is evident from the fact that N_t tends to N_s(0). In other words, the entire population becomes accessible for transmission of the disease. Clearly, the situation can only represent the growth of infection in the absence of any kind of interventions.

On the other hand, the effect of all interventions is to limit the contact transmission rate, thereby limiting the proportion of the exposed population to the disease to a small fraction. It is this that we call the accessible population. In other words, the accessible population N_a(0) is of the same order as the infected population. This can be written as , where is a small fraction. It must be noted here that the value of the fraction depends on the nature of interventions in force.

However, within the scope of the reduced logistic model, the evolution of N_t is independent of the values of the parameters α_sp_s and p_iβ_i during the absence or presence of interventions. As a consequence, the asymptotic value of the cumulative infected population is always N_t = N_s (0), the entire population. Therefore, demonstrating the accessible population is a small fraction of the total population is outside the scope of Eq. (11) and the full model Eqs. (1- 6). An independent way of demonstrating is desirable.

A. Quantitative estimate of the accessible population

Since the factor is not well determined, there is a necessity to get a better estimate of this parameter or the accessible population N_a (0). Assuming that the accessible population is of the order of N_t, we assume that . This is equivalent to using N_a(0) in place of N_s(0) in Eq. (17). Then we numerically evaluate the dependence of N_t on the parameters on N_a (0), α_sp_s, and p_iβ_i. Given the fact that the disease evolves, we expect that the accessible population N_a(0) also evolves with time and in the early stages of evolution, N_a (0) will be small, even in the absence of interventions.

Consider the dependence of N_t on N_a (0), keeping α_sp_s and p_iβ_i fixed. We find that even for relatively large values N_a(0), N_t grows exponentially; for intermediate values, a near saturation value is reached in a relatively short duration of 10-15 days; and for small values, the saturation value is not reached even after 100 days. These features are illustrated in Fig. 1(a) in plots (i-iii) for N_a (0) = 8 × 10⁵, 2.8 × 10⁵ and N_a (0) = 1.45 × 10⁵ respectively, keeping p_iβ_i = 3.3913 × 10⁻⁷ and α_sp_s = 1.0 × 10⁻⁴. We have also examined the influence of p_iβ_i,

• Download figure
• Open in new tab

FIG. 1.

(Color online) (a) Plots of N_t for decreasing values of N_a (0) as a function of time: (i) N_a (0) = 8 × 10⁵, (ii) 2.8 × 10⁵, and (iii) 1.45 × 10⁵ respectively, keeping p_iβ_i = 3.3913 × 10⁻⁷ and α_sp_s = 1.0 × 10⁻⁴ fixed. (ii) shows decreasing N_a(0) by a factor of 5.51 leads to slow increase in N_t. (iv) Plot of N_t for p_iβ_i = 4.752 × 10⁻⁷, keeping Na(0) = 2.8 × 10⁵ and α_sp_s = 1.0 × 10⁻⁴. Smaller values of p_iβ_i take longer time for N_t to grow as is clear (see (ii) and (iv)). (v) Plot of N_t for α_sP_s = 1.1 × 10⁻³, keeping N_a(0) = 2.8 × 10⁵ and p_iβ_i = 3.3913 × 10⁻⁷ fixed. Increase in α_sp_s leads to a faster initial growth seen in (v) and (iii). Also shown is the cumulative number of infections (●) for the UK. (b) Figure shows the two-phase evolution of the disease. The inset shows the good fit using Eq. (17) with the cumulative infected cases for the UK (●) prior to March 23, 2020. Parameter used are N_a(0) = 3.27 × 10⁵, N_t(0) = 13, α_sp_s = 1.0 × 10⁻⁵, p_sβ_i = 6.4622 × 10⁻⁶. Post lock-down period: The four curves (i-iv) correspond to the four iterations of N_a(0) values. (i) N_a(0) = 5.0 × 10⁵,α_sp_s =0.0, (ii) N_a(0) =4.0 × 10⁵, α_sp_s = 1.8 × 10⁻³, (iii) N_a(0) = 3.0 × 10⁵, and (iv) N_a(0) = 2.8 × 10⁵, with α_sp_s = 4.83 × 10⁻³ for latter two. The initial value of N_t(0) = 5687 on March 23, 2020. (c) Plots of daily new cases for the UK and the calculated daily new cases using the reduced logistic model predicted. keeping N_a(0) = 2.8 × 10⁵ and α_sp_s = 1.1 × 10⁻³. Smaller values of p_iβ_i, N_i take longer time for the infection (N_t) to grow. This feature can be seen from the curves (iv) and (iii) for p_iβ_i = 4.7522 × 10⁻⁷ and 3.3913 × 10⁻⁷ respectively. We have also examined the growth dependence of N_t on a_sp_s, keeping the other two parameters fixed. The dependence of N_t on this parameter is similar to that on p_iβ_i. The curve (v) taken together with (ii) shows that increasing α_sp_s also leads to faster initial growth of N_t. In the same plot, we have also plotted the total number of infected cases (●) for the UK.

A careful scrutiny of the total coronavirus cases (●) in the UK shows that it is similar- both in magnitude and shape -to the plot of N_t corresponding to N_a(0) = 0.28 × 10⁶ marked (ii) shown in Fig. (1a). This similarity suggests two important points. First, noting that the UK is under lock-down, one expects that the accessible population is a small fraction of the total population, and therefore we see that the order of magnitude of the accessible population N_a(0) used is comparable to that of the infected population N_t shown in curve (ii). The figure also shows that as much as all populations dynamically evolve during the development of the pandemic, N_a(0) also keeps evolving with time. Second, the similarity in shape of the UK data (●) with the sigmoidal shape of the logistic solution raises a question whether the similarity is accidental. If not, can this be used to fit the UK data?

IV. THE FULL MODEL

One of the challenges of compartmental models is the difficulty associated in making accurate predictions, mainly attributable to the uncertainties in obtaining proper estimates of the parameters [13–15, 35, 36]. For the same reason, forecasting is even more challenging. Often, several factors may also contribute to the same parameter, making it difficult for proper interpretation. In our model, however, several parameters in Eqs. (16) are related to measurable quantities. For instance, the parameters α_sp_s, a_sp_q and a_tp_t respectively represent rates of testing positive, rates identified as pre-symptomatic, and tracing rate of those exposed to the infected. Similarly, parameters λ_q, λ_t, λ_r and κ_d are inversely related to quarantine duration τ_q = 1/λ_q and time required for tracing τ_t = 1/λ_t, time from illness to recovery τ_r = 1/γ_r, and time from illness to death T_d = 1/κ_d. Though these quantities are country/region-specific, their values have been estimated in the literature [5, 6, 8, 35, 36, 38–42]. Some values are also available in the public domain [32, 43]. One parameter that is hard to estimate is the contact transmission rate β_i, which is already estimated in the context of the reduced logistic equation.

A. Calibration of relative strengths of the parameters: insights into disease evolution number of days

However, the dynamical evolution of a nonlinear coupled set of equations such as Eqs. (1- 6) is necessarily complex. Therefore, in the absence of appropriate values relevant for the country/region, a systematic method of finding optimized values of parameters that fit the data under consideration requires calibration of all parameters in the model. This will also help us to delineate the different time-scales participating in Eqs. (1- 6). For instance, several of these parameters represent the growth or decay rates of these populations. Our model contains eight time-scales which are inversely related to the transition rates. These are β_i, α_sp_s, α_sp_q, α_tp_t, λ_q, λ_t, λ_r and κ_d. These time scales control how these populations evolve with time. From the structure of Eqs. (1-6), it is clear that N_i, N_q and N_tr exhibit peaks as a function of time (days) whereas N_r and N_d grow monotonically. However, at what point of time do the peaks appear in these populations with the progression of the disease cannot be easily determined since these are coupled nonlinear differential equations where the evolution of any population depends on the evolution of all other populations [29]. More importantly, if one is interested in fitting the model predicted growth of populations, which convey the disease status (such as the total infected, active, recovered and dead), estimating the relative proportion of the populations as the disease evolves is necessary. Further, the total infected population commonly used to convey the disease status in daily briefings has contributions from all populations. Therefore, delineating and determining at what points of time each of these populations contribute to the total populations would provide required insight into further analysis.

Following the recently developed method in the area of plasticity [44–46], we investigate the influence of the parameters to identify the relative importance of the transition rates. Since it is a multi-parameter space, we vary each parameter, keeping all other parameters fixed at the reference set of values listed in Table II. The results are illustrated using plots of the active infected population N_i. The dotted curve shown in Fig. 2 is the reference curve corresponding to the reference set of parameters given in Table II. As in the case of the reduced logistic model, the growth of N_i depends sensitively on the contact transmission parameter p_iβ_i. (Note that in our model, this is the only parameter that directly contributes to growth of infections.) A 50% increase in the parameter induces a substantial increase (38%) in the peak height with the position shifting towards earlier time (by 25 days) as is clear from the curve (i). Noting that the position of the peak, i.e., the turning point of N_i is indicative of slowing down of the rate of infection, increasing peak height of N_i suggests slowing down of infections occurs at higher values of infections, where as a shift towards shorter times implies that slowing down occurs earlier. A similar effect but of lesser magnitude (20%) is seen when testing rate α_sp_s is increased by a factor 8 as is clear from (ii).

View this table:

• View inline
• View popup
• Download powerpoint

Table II.

Post lock-down period: Select set of parameter values serving as a reference set. q₁ = q₂ = 0.08 and p_i = 0.1.

• Download figure
• Open in new tab

FIG. 2.

(Color online) Calibration of parameters for identifying the relative importance of transitions contributing to N_i by varying one parameter, keeping all others parameters fixed at reference values listed in Table II. The dotted curve is the reference plot for the active infected population N_i corresponding to the values in the Table II. (i) Plot of N_i for a 50% increase in p_iβ_i showing a substantial increase in the peak height (38%) and a shift by 25 days. (ii) Similar effects are seen when α_sp_s is increased by a factor 8. (iii) The plot shows a substantial decrease in the peak height with a marginal shift in its position when α_sp_q is increased by a factor of 2. (iv) The plot shows a decrease in the peak height as is increased by 60%. (v) Similar effect is seen when the recovery rate γ_r is increased by 60%.

In contrast, an increase in the quarantining rate α_sp_q by a factor 2 decreases the peak height comparable in magnitude to that induced by β_i with a shift in the peak position away from the origin by 20 days (see iii). A similar effect, but of lesser magnitude, is seen when α_tp_t is increased (not shown). The reduction in the peak height of N_i is understandable because the total inward transition is α_sp_s + α_sp_q + α_tp_t, and therefore increasing one of these changes the relative weights. Physically, the decrease in the peak height (commonly referred to as ‘flattening the curve’) accompanied by a shift in the peak position for longer times implies that as more individuals are quarantined, the disease control is facilitated.

We have also investigated the dependence of the recovery (γ_r) and death rate (κ_d) parameters on N_i. An increase in γ_r by 60% decreases the peak height substantially as is clear from (iv). Similarly, increase in the death rate κ_d leads to a decrease in the peak heights of N_i as is clear from (v). Clearly, increased recovery rate, a desirable feature, also leads to a decrease in the number of active infections. This feature is also intuitively obvious. On the other hand, increased death rate, though not desirable, also leads to decrease in N_i. We have also investigated the influence of other parameters and find that N_i is relatively insensitive to these parameters. Noting that any change in the parameter values relative to those corresponding to the reference curve changes the peak position and height, we conclude that the parameters listed in Table II are close to the optimized values.

The above analysis on the relative importance of the transition rates (equivalently the corresponding parameters), and the accompanying discussions, can now be used to delineate the contributions from different time scales participating in the transmission dynamics of the virus. As discussed above, the direct transitions from S to I, namely, p_iβ_i, to a lesser extent α_sp_s, control the initial growth of N_i. The same also to the initial growth of the total infected population N_t. Similarly, the delayed routes, namely, the quarantine S and tracing T contribute to the mid region of the evolution of N_i and to therefore to the total population N_t also. Recall that the turning point of N_i is controlled by the balance between all inward transitions (from S to I, and Q and T to I) and outward transitions (I to R and D compartments). Then, the peak position of N_i can be identified with the inflection point of N_t. Therefore, the time development beyond the point of inflection of N_t is controlled by a balance of all inward and outward transitions. The insights from the above analysis identify three distinct stages in the evolutionary period of the disease, namely, the initial growth period, the mid developmental period and the final approach to saturation. As we shall see, this identification will be helpful in obtaining a good fit to the UK data.

B. Data assimilation and forecast

Having demonstrated that the two direct transition rates p_iβ_i and α_sp_s are the dominant contributions to the growth of N_i and having assessed the relative importance of other transitions, we now consider the solution of the full model Eqs. (1- 6) with a view to obtaining the best possible fit with the COVID-19 United Kingdom data. Attempt will also be made to forecast the future progression of the disease.

Recall that the spread of coronavirus in the UK falls into two phases of development. During the first phase prior to March 23, 2020, there were no constraints and the disease transmission was free. After the lock-down date, the transmission is restricted. Therefore, the model parameters and the initial conditions relevant for the two phases are different. As in the reduced model, we assume that the dynamics of the disease transmission is limited by the accessible population N_a(0) and not by the total population N_s(0), i.e., . Note that both N_i and therefore N_a(0) depend on whether there are any interventions or not, and also on the type of interventions.

Consider the period between January 31 and March 23 corresponding to the initial phase. Here, we note that prior to the lock-down date, there would be no quarantining and tracing procedures and therefore, it is adequate to solve Eqs. (1,2,5,6), More specifically, we ignore the delayed inward transitions into I. Furthermore, in the first few days of the development of the disease, we may assume that the total number of the infected cases N_t is equal to the active infections N_i. Finally, since, we plan to fit the model solution with the UK data [32], publicly available coronavirus data for the total number infected, active infected, recovered and the dead are useful in further optimizing the parameters. Unfortunately however, only the total numbers of the infected and the dead are made available in the UK.

Now we are in a position to solve the relevant equations for the first phase. As discussed in Section III B, we use February 27, 2020 as the starting day for the first phase evolution of Eqs. (1,2,5,6). The local growth rate of 0.25849/day (on the starting day) over 8 days obtained from the log-linear plot of the cumulative infected cases for the UK is equated with the model growth rate given by p_iβ_iN_a(0) to fix β_i = 6.4622 × 10⁻⁶ by using the initial value of N_a(0) = 4.0 × 10⁵. Further, using the initial conditions for N_t(0) = N_i(0) = 13, N_r(0) = 0 and N_d(0) = 0, we solve Eqs. (1,2,5,6) from February 27 to March 23 by choosing a value for α_sp_s that gives the best fit to the data for the period. (Here, α_sp_s = 3.6 × 10⁻⁶ and the values of other relevant parameters are those listed in Table II.) The model-predicted total infected population N_t (continuous curve) along with the data points (●) is shown in the inset of Fig. 3(a). Clearly, the match is seen to be very good. Also shown is a plot of active infections N_i (dotted curve). Equations (1,2,5,6) also provide the values of N_i, N_r and N_d on March 23, 2020. These are N_i = 5407, N_r = 400, N_d = 285.

• Download figure
• Open in new tab

FIG. 3.

(a) Inset: Plot of the total infected population (continuous curve), along with the UK data (●) from February 27, 2020 till March 31. Also shown is the active infected N_i (dotted curve). Post Lock-down period: Plots of the total infected population (curve marked i) along with the corresponding data for the UK (●) from March 23, 2020. Curve marked (ii) shows the active infected population. (b) Plots of the reported daily new cases for the UK along with the model predicted daily new cases. (c) Plots of the infected (i), quarantined (ii), traced (iii), recovered (iv) and deceased (v), starting from March 23, 2020. The values of the parameters used are given in Table II.

Now we consider the solution of Eqs. (1- 6) with a view to obtaining the best fit for the UK data for the period starting from March 23, 2020. Here again, we first find the growth rate from the data and equate it with the model growth rate. Here, we note that the effect of lock-down is expected to manifest after some time. For this reason, we use a 17-point slope in the log-linear plot. Equating 0.13/day with p_iβ_iN_a(0) and using N_g(0) = 3.75 × 10⁵ we get β_i = 3.4693 × 10⁻⁶. The initial values used for evolving Eqs. (1- 6) are N_i(0) = 5407, N_q(0) = 0, N_tr(0) = 0, N_r(0) = 0, N_d(0) = 0. (The reason for using zero initial conditions for N_q(0), N_tr(0), N_r(0) and N_d(0) is that the values of these populations would not be recorded during the first phase. However, using the values obtained from the first phase for N_r and N_d makes little difference. Note that N_i(0) = 5407 is smaller than the total number of infected cases. Again, using N_i(0) = 5687 does not alter the results.) The parameters values used are those listed in Table II. Figure 3(a) shows plots of the calculated total infected population N_t and the total infected cases in the UK (●). Clearly, the fit is very good. The asymptotic saturation value of the total number of infections N_t turns out to be 3.52 × 10⁵. Also shown is the active infected N_i labeled (ii). Further, the plot of the model predicted active infected population N_i (ii) exhibits a turning point (peak) around May 15. Subsequent decrease in N_i is seen to be slow. At this rate of slowing-down, the model predicts that a near saturation value of 3.52 × 10⁵.

Now consider the estimation of the end time of the epidemic. A natural candidate for estimating this is the daily new infections. Considering the close fit of the model predicted cumulative infected population Nt with the UK data (see Fig. 3a), the model calculated daily new infections can be compared with the UK data for the daily new cases. This is shown in Fig. 3(b). As can be seen the general profile of the model predicted daily new infected cases matches well with the published data. Recall that the end time of the epidemic is defined as the time at which no new cases are reported. This definition is clearly impractical for forecasting since the state of saturation is always reached asymptotically. For this reason, we use an arbitrarily small value, say 300 new cases as the end date of the epidemic. Then, the estimated end time is late July.

To the best of our knowledge, we are not aware of any model that fits the COVID-19 data for any country as has been done here, particularly over such long periods with the ability to forecast the future progression of the disease. However, there have been some efforts to fit data for the initial periods [8, 21, 22, 24, 26, 38, 47]. In view of the good fit for the UK data, the model can be used to fit the COVID-19 data for other countries and also to forecast the progression of the disease.

V. SUMMARY, DISCUSSION AND CONCLUSIONS

Recent literature has focused on abstracting the effect of various types of interventions through epidemiological models to make projections of how the disease progresses under different conditions. Recall that one limitation particularly applicable to the deterministic compartmental models is the difficulty in getting proper estimates of the parameters, particularly when the number of compartments is large. In this respect, simpler models with fewer compartments have an advantage. However, several factors may contribute to a single parameter. This is also a model-dependent feature. Therefore the ability of such parameters to represent the mitigating efficacy of interventions appears limited (see below). Furthermore, the number of parameters in such models is not necessarily small, making numerical solutions often the only choice. Therefore, any method - whether mathematical or conceptual - which simplifies analysis and easy interpretation is welcome.

Motivated by this, we hypothesize accessible population for transmission of the disease that can be the total or a small fraction of the total population depending on whether the transmission dynamics evolves in the absence or presence of interventions. Indeed, the effect of lock-down is evident in all counties where the disease has been controlled or nearly eliminated. At the mathematical level, we introduce a decoupling scheme to aid mathematical analysis that also helps easy interpretation. The model equations have been devised in such a way that the susceptible and active infected populations form the main populations. The decoupling is affected by dropping all inward and outward transitions excepting the direct transitions (p_iβ_iN_a(0) and α_sp_sN_a(0)). Because, all outward transitions from I are ignored under this decoupling, the active infected population N_i takes the role of the cumulative infected population N_t. The simplicity of the reduced logistic equation (Eq. 11) allows easy identification of the growth and inhibiting factors in terms of the dominant growth factors (direct inwards transitions or parameters). Surprisingly, this simple equation provides a good fit to the reported cumulative number of infections for the United Kingdom, as is clear in Fig. 1(b).

The full model Eqs. (1- 6) contain several parameters whose range has been estimated in a number of studies[8, 21, 22, 24, 38, 40–42, 47]. However, when it comes to explaining or capturing the growth characteristics for a specific country, optimized parameters suitable for the situation are required. Following [44–46], we have determined the relative importance of the various transition rates (equivalently the associated parameters) subject to the constraint that the parameter values provide the best fit for the given data. In this work, we have made use of publicly available data on the total infected and daily new infected cases for the United Kingdom [32].

Figure 3(a) shows the fit obtained for the period till March 23, 2020 (shown in the inset) and for the period beyond. Clearly, the fit is seen to be very good for both the period till the lock-down date and the period thereafter. Comparing Fig. 3(a) with Fig. 1(b) for the reduced logistic map, we see that while the fits in both cases are comparable, the projections of the future progressions are significantly different. The saturation value predicted by the full model (shown in Fig. 3(a) is close to 3.52 × 10⁵, whereas that predicted by the reduced logistic equation in Fig. 1(b) is ~ 2.9 × 10⁵. Conventionally, the end time of epidemic is defined as the day on which no new infections are reported. However, since the approach to the end point is generally slow, a better variable to predict the end point of the epidemic is by comparing the model computed daily new case with the reported data on the daily new cases for the UK. Then, the end time of the epidemic predicted by the full model turns out be late July (see Fig. 3b). In contrast, the end time for the epidemic predicted by the reduced logistic model is late June. Clearly, the results obtained from the full model emphasize the limitations of the reduced model. A natural question is: what are the underlying causes?

The fact that the reduced logistic model provides a good fit also means that the major contributing factors for the growth of infection are included in Eq. (2). To see this, consider Eqs. (1-6). The growth of N_i(t) or more appropriately the daily new infected cases has two types of inward transitions, namely, direct and delayed. However, the dominant direct transition from S to I given by p_iβ_i controls the growth rate of N_i. The other direct transition α_sp_s into I also contributes to a lesser extent. Now, consider the delayed inward transitions to I coming from Q and T. These transitions are smaller in magnitude and contribute to sub-exponential growth of N_i in time. More importantly, the turning point in N_i or is due to a competition between the growth factors (all inward transitions) and the outward transitions (recovery and fatality terms). Further, since the time evolution beyond the turning point of is controlled by a balance between all inward and all outward transitions, the approach towards the state of no infections or the saturation value of N_t is slow in our case. These features are clear from Fig. 3(a,b). Note that the fit, till June 1, is just two weeks beyond the turn point of . and it has a long way to evolve to the end point of the epidemic.

These arguments clarify two features of the data fit obtained using the reduced logistic equation. Because the total number of infected cases N_t has the dominant growth contributions, the good fit is not surprising. On the other hand, growth dynamics beyond the turning point in the daily new infected cases (i.e., ) is controlled by a balance between growth factors (all inward transitions) and inhibiting factors (the rate of recovery and dead). However, these competing time scales are absent in the logistic equation. This clarifies why the projected saturation value of N_t and the end time of the epidemic is not well captured.

A few comments above model are in order. As stated in the introduction, the standard compartmental models ignore both age dependence and the heterogeneous spatial distribution of the population, and therefore cannot abstract the aspects that depend fundamentally on these two features. For the same reason, the populations in these models represent only the mean response of each of these populations. From this point of view, the good fit obtained by the full model, to lesser extent by the reduced logistic model, may come as a surprise. However, the success story of the mean field approaches has been established in physics literature [29, 44–46] and the limitations of the approach has also been well established. To this extent, the success of the present model can be attributed to the way the model equations are structured and the underlying nonlinear dynamical methods used for analysis.

Here, we mention that the parameter values used for fitting have been obtained using an optimization procedure subject to the constraint that the optimized values should fit the data for the cumulative infected cases for the UK. Although, each parameter is varied within the range of values estimated in the literature (and open sources), since the optimization has been carried out subject to only one constraint, the values may not be unique. From this point of view, more number of constraints such as the data for the active infected, recovered, quarantines etc., would be helpful in removing the non-uniqueness of the optimized values, at least partially.

As stated in the text, the contact transmission rate parameter β_i is one of the crucial parameter in the model because this parameter largely controls the time development of the disease. In our approach, this parameter has been estimated by using the local initial slope in the log-linear plot of the UK data. The slope itself, however, depends on whether the disease development occurs in the absence or presence of interventions, which in turn depends on the accessible population N_a(0). Since the model growth rate depends on p_iβ_iN_a(0), the value of β_i depends inversely on N_g(0) corresponding to the absence or presence of interventions. However, the value reported by Tang et al. [6] is two orders smaller than that estimated in our paper. However, in both cases, the value of β_i is inversely proportional to the relevant population used for modeling. In this context, we mention that there is only independent estimate both under free evolution [5]. Such an independent estimates are desirable.

One point that needs some discussion is about the values of the recovery γ_r and the death rates obtained from the optimization procedure used that fit the UK data very well (see Table II). These rates are inversely related to the time duration between detection of illness till recovery and death respectively. While these values are within the published range of values [40–42], they appear to be on the lower side of the mean. However, in principle, the parameter values depend on the structure of model equations. In our model, the recovery and fatality are outward transitions from a single compartment, namely, I to R and D. This feature is clearly because our model equations were structured to have only two core populations. However, the recoveries can occur if other kinds of compartments (populations) are included [5]. Thus, the best fit for lower values may be the result of simplicity of the model.

In conclusion, the simple compartmental model not only provides a good fit to the United Kingdom coronavirus data but also makes concrete long term predictions for the future. We believe that these results have been made possible due to the reductive approach adopted here.