Simple Formulae, Deep Learning and Elaborate Modelling for the COVID-19 Pandemic

Abstract

Predictive modelling of infectious diseases is very important in planning public health policies, particularly during outbreaks. This work reviews the forecasting and mechanistic models published earlier. It is emphasized that researchers’ forecasting models exhibit, for large t, algebraic behavior, as opposed to the exponential behavior of the classical logistic-type models used usually in epidemics. Remarkably, a newly introduced mechanistic model also exhibits, for large t, algebraic behavior in contrast to the usual Susceptible-Exposed-Infectious-Removed (SEIR) models, which exhibit exponential behavior. The unexpected success of researchers’ simple forecasting models provides a strong support for the validity of this novel mechanistic model. It is also shown that the mathematical tools used for the analysis of the first wave may also be useful for the analysis of subsequent waves of the COVID-19 pandemic.

1. Introduction

The first period of the COVID-19 pandemic in Europe, often referred to as ‘the first wave’, provides an ideal situation for the retrospective study of the various mathematical formulations used for modelling the dynamics of this viral epidemic. Indeed, during this phase, many European countries employed similar measures for combating the spread of the epidemic, and moreover, the pandemic was dominated by a single viral strain. Such a study is not only useful for academic purposes, but also, importantly, could identify an accurate mathematical formalism that could be used in a variety of epidemiological phenomena.

In the framework introduced by Holmdahl and Buckee [1], epidemiological models are broadly divided into two categories: forecasting and mechanistic. The former models fit a specific formula to the data, and then attempt to predict the dynamics of the quantity under consideration. The main limitation of these models is that they usually remain valid for a short period of time, and specifically only if the epidemiological situation remains unchanged. For example, they can be used during a lockdown period, but will not make accurate predictions after the lockdown is lifted. In contrast to forecasting models, mechanistic models can make predictions even when the relevant circumstances change. The main limitation of the forecasting models is the difficulty of determining the parameters specifying them.

In what follows, the researchers first discuss specific forecasting and mechanistic models introduced by the collaborative efforts published in [2,3,4,5,6,7,8]. In particular, the researchers provide an explanation for the remarkable success of researchers’ simple forecasting models for predicting the dynamics of the first wave. The researchers then emphasize that the tools introduced for the analysis of the first wave can also be useful for analyzing subsequent waves of the COVID-19 pandemic. This is illustrated with a brief discussion of the second wave.

2. The Success of the Novel Formulae as Opposed to the Logistic Formula

A novel class of forecasting models was presented in [2,3]. Although, as correctly noted in [9], the forecasting models are ‘not well suited for long-term predictions,’ the researchers’ class of forecasting models provided highly accurate long-term predictions both for the number of deaths caused by COVID-19 and the number of reported infected cases. Specifically, by training researchers’ novel formulae (i.e., by fixing the constant parameters appearing in theses formulae) using data for the accumulative number of deaths in a given country until 1 May 2020, it was possible to make predictions for a period of more than 3.5 months. A typical example is shown in in Figure 1: for the COVID-19 epidemic in Italy, it is seen that there was no deviation between the curve depicting the number of deaths and the curve constructed via the model formula. Incidentally, although this formula was trained using dates during the lockdown period, its predictions remained accurate even after the easing of the lockdown conditions; possible explanations for this unexpected fact are discussed in [4].

In researchers’ previous work [2,3], the usual logistic formula used in epidemiology is replaced by a different formula, called the rational formula. The logistic and rational formulae are given, respectively, by

$$N\left(t\right)=\frac{N_f}{1+\beta e^{-kt}}$$

(1)

and

$$N\left(t\right)=\frac{N_f}{1+\beta {\left(1+dt\right)}^{-k} }.$$

(2)

In these formulae N(t) denotes the cumulative number of deaths (or reported infected cases) at time t; the parameters N_f, k, β, and d are constants that can be determined by ‘training’ the model using the given data.

In addition to the rational formula, researchers also used a simple variant of this formula, called birational [2,3]. The birational formula is similar to the rational formula, but the associated parameters are different before and after a fixed time point. Both the rational and birational formulae are based on a particular nonlinear ordinary differential equation, called the Riccati equation. This equation is specified by the constant parameter, N_f, denoting the final total number of deceased (or infected individuals), and by a time-dependent function, a(t). Remarkably, although this Riccati equation is a nonlinear equation that involves time-dependent coefficients, it can be solved in closed form. Its solution depends on N_f, a(t), and on a parameter related to the constant that arises in integrating this equation. In the particular case that a(t) is a constant, the explicit solution of the above Riccati equation becomes the classical logistic Formula (1). If a(t) = kd/(1 + dt), then it becomes Formula (2).

3. Deep Learning

In previous work [2], researchers also investigated deep learning algorithms, and in particular, the bidirectional long short-term memory (BiLSTM) network, which is a powerful generalization of recurrent neural networks (RNN) that can capture long-term dependencies while at the same time avoiding the problem of vanishing/exploding gradients. RNN can capture sequential information but cannot learn from long term dependencies. A common problem with RNN, particularly when many hidden layers are used, is the vanishing or the exploding of the gradient in the backpropagation algorithm occurring during the process of updating the weights. Hochreiter and Schmidhuber introduced in [10] the Long Short-Term Memory (LSTM) networks that can capture long term dependencies and at the same time avoid the problem of vanishing/exploding gradients. LSTM networks are a type of RNN, where a memory cell maintains its state over time; they use gates to decide whether to flow information in (keep) or out (forget) of the memory cell.

Graves and Schmidhuber introduced in [11] the BiLSTM networks; this development was motivated by the Bidirectional RNN networks, introduced earlier by Schuster and Paliwal in [12], where the training runs forwards and backwards using two separate RNN. Consequently, the main difference is the training sequence: in the LSTM, the sequence runs backwards, preserving information from the future, whereas in the BiLSTM training, the sequence runs backward and forward, preserving information from both the past and the future. The BiLSTM networks are well suited for time series prediction and can potentially completely capture the contextual information of the time series.

The rational and birational formulae are, in a sense, optimal, since such elaborate deep learning algorithms investigated in [2] could not improve the predictions obtained by using these simple analytic formulae.

What is the critical difference between the logistic and the rational formulae? Although, for large time, t, both the logistic and the rational formulae tend to a constant, N_f, namely, the plateau reached by the epidemic, this constant is approached with exponential and algebraic rates, respectively. Indeed, for large t, Equations (1) and (2) yield the following asymptotic expressions:

$$N\left(t\right) ~ N_f\left(1-\beta e^{-kt}\right), t\to \infty$$

(3)

and

$$N\left(t\right) ~ N_f\left(1-\frac{\beta }{d^k}{t}^{-k}\right), t\to \infty$$

(4)

It is interesting to note that given a set of data, it is always possible to choose the constant parameters that specify these two formulae, so that they can match the given data extremely well. However, when these formulae are used to predict the dynamics of the epidemic, it becomes clear that the curve obtained via the logistic formula immediately begins to deviate from the curve of the real data, whereas the rational formula yields excellent predictions.

Regarding the algebraic decay, it is worth noting the following observation: whereas physical phenomena are often accurately modelled by exponential processes, social phenomena usually require an algebraic formalism.

4. A Typical Susceptible-Exposed-Infected-Recovered (SEIR) Type Model

In summary, the analysis of the real data from several European countries during the first wave of the COVID-19 epidemic suggests that the long-time behavior of the dynamics of the underlying processes is characterized by an algebraic as opposed to an exponential decay. This raises the following important question: do the mechanistic equations used to model these dynamics exhibit large t algebraic decay? This important question is discussed below.

It is well known that the most widely used mechanistic models are the so-called susceptible-exposed-infected-recovered (SEIR) type models. A typical such model was analyzed in [5,6]. This model consists of six ordinary differential equations modelling the number of infected, asymptomatic, sick, recovered, hospitalized, and deceased. Numerous other investigators have analyzed similar models, which include the well know works of [13]. Different SEIR models use different number of ‘compartments,’ but the essential physiological hypothesis used for deriving these models is the same.

The large time asymptotics of the model of [5] are analyzed in [7], where it is shown that this SEIR model exhibits exponential decay. It is natural to expect that all other SEIR type models will exhibit a similar behavior, which is inconsistent with the summary at the beginning of this section.

By modelling in [7] a physiologically important mechanism, which was until recently ignored, an additional nonlinear effect was introduced in the six ordinary differential equations presented in [5]. Remarkably, for large t, the novel model introduced in [7] exhibits algebraic decay, consistent with the decay exhibited by the real data. It is noted that the new physiological mechanism used in [7] was first introduced in [8,14] in connection with a simpler model.

5. From the First to the Second Wave

RNA viruses, as exemplified by influenza, may acquire mutations that can have a significant effect on their biological behavior [15,16]. Viruses’ mutation rates are affected by various factors inherent to their replication such as the polymerase action, or external to the virus itself such as environmental effects. The precise contribution of these effects and how natural selection acts to optimize these rates remains unknown. The novel coronavirus, even before harboring genetic changes like the D614G polymorphism, was already fit enough to cause a pandemic.

The novel pandemic coronavirus SARS-CoV-2 acquired mutations early on [17], which altered some of its essential characteristics, affecting diagnosis [18], disease phenotype and severity [19,20], therapeutic interventions [21], and vaccine efficacy [22].

The exact public health impact of the early mutations remains debatable. By April 2020, strains carrying the D614G mutation [23,24] became prevalent in Europe, the USA, and Canada [23]. In particular, in a study in Houston, Texas, USA, the D614G variant strain gradually fully dominated the local epidemic [25]. The high prevalence of this strain was also confirmed in a UK study involving the analysis of 25,000 genomes [26]. In human infections, this strain was reported to be associated with a fitness advantage, as well as with greater transmissibility [23] and higher viral loads.

A difference in the growth rates between the higher levels of recorded infections and lower levels of respective number of deaths during the evolution of the pandemic, and especially during the second surge of the infection, is evident in the examination of the epidemic curves available from the European Centre for Disease Control [27,28]. A number of different factors may be associated with this ‘paradox,’ including increased testing, better health system capacity, earlier treatment, and increased levels of immunity in local populations. The possibility that the virus—due to structural differences—has evolved to be associated with less virulence, has been reported at the experimental level [29]. In this context, studies investigating the effect of earlier mutations including the D614G polymorphism showed that infections with D614G carrying strains were associated with higher viral loads in the upper respiratory tracts and a fitness advantage, leading to increased transmissibility [23,30,31,32,33,34]. Two large genomic analysis studies suggested otherwise [26,35]; they did not find convincing evidence with regards to infectivity or virulence [26,35]. The question of establishing whether this strain is indeed more infectious or less virulent remains crucial with regards to virus evolution to less virulent forms [36]. In the meantime, the world is faced with the emergence of new variants [37,38,39,40,41,42,43,44,45,46,47,48,49,50].

5.1. Predictions for the Number of Reported Infected Cases and the Number of the Deceased for the Second Wave in Portugal

In order to show that the techniques developed for the first wave can also be used for the analysis of subsequent waves, researchers applied the rational and logistic formulae to the cumulative number of deaths and reported infected cases for the second wave of Portugal. These formulae provided accurate predictions for the relevant dynamics (Figure 2). The data were fitted up to 27 December 2020, the last day before the first appearance of the UK lineage. Therefore, the researchers’ curves determine what the plateau for reported cases and deaths would have been without the emergence of the UK lineage. For the number of deaths, the logistic formula predicted a plateau on 21 February 2021 with 6002 deaths for the second wave and 7798 total deaths overall; the rational formula predicted a plateau on 2 March 2021 with 6432 deaths for the second wave and 8228 total deaths overall. The plateau is defined as the point when the rate of deaths for the given wave is 5% of the maximum rate. It should be noted that in order to get accurate predictions, the model should employ data at least up to a short period after the time when the rate of deaths reached a maximum. In Portugal, for the deaths of the second wave, this occurred on 6 December 2020.

For the number of reported infected cases for Portugal, the logistic formula predicted a plateau on 2 February 2021 with 364,526 cases for the second wave and 420,123 total cases overall; the rational formula predicted a plateau on 5 February 2021 with 369,040 cases for the second wave and 424,637 total cases overall. In Portugal, the rate of infections of the third wave reached a maximum on 21 November 2020.

5.2. Higher Infectivity and Lower Virulence

In order to compare the infectivity and virulence of the first two waves, additional techniques are needed. In this regard, first, researchers computed the rates dN/dt directly from the data, where N denotes the number of reported infected, in seven European countries (in order to avoid the noise of the real data, researchers computed these rates by fitting the data for the first and second waves to the rational and the logistic formulae, respectively). Second, researchers computed the function v(t) defined as the ratio of the cumulative number of deaths to the total number of reported infected cases that occurred three weeks prior to the date that the accumulated number of deaths was computed.

The computed rate of reported infected cases in the seven European countries analyzed was larger in the second wave, which supports the suggestion that the circulating strain of the virus of the second wave was more infectious. Figure 3a presents the relevant curves for Spain, which is prototypical of all cases examined. Similarly, the computed rate of deaths was much lower, seen in Figure 3b.

For the case of Portugal, the rate of deaths can also be computed indirectly by examining the logistic and rational models. By comparing the parameters of the fitted logistic formula for Portugal (Figure 2) between the first and the second waves, it is confirmed that the rate of deaths is lower in the second wave (the parameter k is 0.0778 for the first and 0.0564 for the second, whilst the parameter N_f is 1489 and 6076).

In order to examine whether death rates were slower during the evolution of the pandemic and especially in the last half of 2020, researchers computed v(t) using the following steps. First, researchers separated the reported infected cases and the deaths for each country in the two waves and calculated their plateau. Then, in each wave, researchers computed the function v(t) by dividing the cumulative number of deaths up to a given date, by the cumulative number of reported cases 3 weeks earlier. Researchers chose the period of 3 weeks because individual infections during the last 3 weeks are not expected to greatly affect the number of deaths [51].

The graphs of v(t) for Germany, Greece, Spain, UK, Italy, Czechia, France, and Belgium are shown in Figure 4. The values of v(t) after the situation gets stabilized are shown in

Table 1. Stabilization values for the function v(t) for the first and second waves, for various countries, for all population.

All Population	1st Wave		2nd Wave
	Date	v(t) (%)	Date	v(t) (%)
Germany	30/June/2020	4.86	26/November/2020	1.50
Greece	04/July/2020	6.18	02/January/2021	3.95
Spain	30/June/2020	11.72	25/December/2020	1.49
UK	30/June/2020	15.24	19/September/2020	2.55
Italy	30/June/2020	14.77	22/December/2020	2.54
Czechia	15/July/2020	3.33	31/December/2020	2.04
France	15/July/2020	18.62	18/December/2020	1.46
Belgium	14/July/2020	15.77	19/December/2020	1.78

. These results suggest that the pathogenicity of the second Autumn/Winter wave was smaller than the first for these countries. Various factors may affect this suggestion, including increased testing, better health system capacity, earlier treatment, and increased levels of immunity in local populations. Among these factors, perhaps the increased number of testing during the second wave is of particular significance. As shown below, this effect can be quantified.

Indeed, let I and R denote the number of infected and reported infected individuals, and let the subscripts 1 and 2 indicate first and second wave, respectively. In the analysis that follows, the researchers assume that I₁ = α₁R₁ and I₂ = α₂R₂. Differentiating these equations with respect to time and dividing the resulting expressions, the researchers find

$$\frac{\dot{I_2}}{\dot{I_1}}=\frac{{\alpha }_2}{{\alpha }_1}\frac{\dot{R_2}}{\dot{R_1}}.$$

(5)

Since the number of tests increased during the second wave, α₂ < α₁.

Thus, the researchers’ conclusions are valid provided that the growth of $\raisebox{1ex}{${\dot{R}}_2$}\!\left/ \!\raisebox{-1ex}{${\dot{R}}_1$}\right.$ is sufficiently large so it can compensate for the decrease in $\raisebox{1ex}{${\alpha }_2$}\!\left/ \!\raisebox{-1ex}{${\alpha }_1$}\right.$. Similarly, let ν^a, where the superscript ‘a’ denotes ‘actual,’ be defined by $v^a=\raisebox{1ex}{$d$}\!\left/ \!\raisebox{-1ex}{$I$}\right.$. Then, $v^a=\raisebox{1ex}{$v$}\!\left/ \!\raisebox{-1ex}{$a$}\right.$, and hence:

$$\frac{v_2^{\alpha }}{v_1^{\alpha }}=\frac{{\alpha }_1}{{\alpha }_2}\frac{v_2}{v_1}$$

(6)

The researchers’ conclusions remain valid, despite the increase in the number of tests during the second wave, provided that the decrease of $\raisebox{1ex}{$v_2$}\!\left/ \!\raisebox{-1ex}{$v_1$}\right.$ is sufficiently small to compensate for the growth of $\raisebox{1ex}{${\alpha }_1$}\!\left/ \!\raisebox{-1ex}{${\alpha }_2$}\right.$. The researchers have been able to find data for the ratio $\raisebox{1ex}{${\alpha }_2$}\!\left/ \!\raisebox{-1ex}{${\alpha }_1$}\right.$ only for Italy and UK; these numbers are quite small, namely 1.40 and 1.65, respectively [52,53]. Researchers speculate that similar rations are valid for the other counties considered in the researchers’ study. Such small ratios will yield slight modifications for the graphs in Figure 4 and the values in the Table 1, but the overall conclusions remain the same.

6. Conclusions

The researchers have reviewed different predictive models used in the collaborative efforts published in [2,3,4,5,6,7,8] and have emphasized that the success of the researchers’ forecasting models is due to the fact that, for large time, they capture the essential behavior of a recently introduced improved mechanistic models. Of course, mathematical modelling is only suggestive and cannot provide definitive proof of the different biological behaviors of the virus. In particular, the number of deaths and the rates of confirmed cases vary significantly across countries and there is significant dependency on the testing performed, as well as significant underreporting [54]. Moreover, the index of the probability of death per case of disease (Infection Fatality Rate) also varies from country to country. This is the result of a variety of factors, including, for example, differences in testing strategies and in case ascertainment. In addition, variance in population characteristics and comorbidities among countries and continents affects the observed morbidity and mortality (e.g., age is the most significant factor associated with increased mortality). Moreover, the percentage of cases involving specific high-risk sub-populations (such as those in nursing homes [55]) varies according to country. Furthermore, differences in timing of measures adopted, as well as differences in targeted control strategies across countries (including targeting the high-risk population like the middle-aged and over), affects death rates. Relating to this observation, the health system quality and capacity widely varies across nations and continents. The saturation of medical facilities may further affect mortality rates, especially in vulnerable high-risk populations.

Researchers have also shown that the mathematical modelling tools used in the first wave, as well as variations of these tools, can be useful for the analysis of subsequent waves. The researchers have illustrated this fact by analyzing the second wave. The researchers suggest that this wave is characterized by higher transmissibility and less virulence than the first wave. However, in these connections, it should be noted that several clinical and physiological factors may affect the disease course towards a better outcome. Among these are the increasing experience of doctors treating the disease and its complications, along with early diagnosis, small hospital volumes in early stages of the second wave that allowed for better care, and the use of new antivirals and immune modulation treatments. All these effects could potentially lead to milder clinical forms of the disease [56,57]. Thus, the value of the ratio proposed herein may tend to underestimate the virulence of new strains.

In conclusion, it should be noted that it is not possible to predict the clinical effect of new mutations. For example, although the second wave apparently involved a less virulent strain than the first, the third wave was more severe. On the other hand, the omicron variant is highly transmittable but also has very low virulence.