An algebraic formula, deep learning and a novel SEIR-type model for the COVID-19 pandemic

The most extensively used mathematical models in epidemiology are the susceptible-exposed-infectious-recovered (SEIR) type models with constant coefficients. For the first wave of the COVID-19 epidemic, such models predict that at large times equilibrium is reached exponentially. However, epidemiological data from Europe suggest that this approach is algebraic. Indeed, accurate long-term predictions have been obtained via a forecasting model only if it uses an algebraic as opposed to the standard exponential formula. In this work, by allowing those parameters of the SEIR model that reflect behavioural aspects (e.g. spatial distancing) to vary nonlinearly with the extent of the epidemic, we construct a model which exhibits asymptoticly algebraic behaviour. Interestingly, the emerging power law is consistent with the typical dynamics observed in various social settings. In addition, using reliable epidemiological data, we solve in a numerically robust way the inverse problem of determining all model parameters characterizing our novel model. Finally, using deep learning, we demonstrate that the algebraic forecasting model used earlier is optimal.


Introduction
In the framework of Holmdahl & Buckee [1], epidemiological models are broadly divided into two categories: forecasting and mechanistic. The former fit early-time data with a specific empirical formula, which is then used to predict the time evolution. The limitation of the forecasting models is that they usually remain valid for only a specific, short, period of time, during which the epidemiological situation is unchanged. For example, if a forecasting model is valid during part of a lockdown period, this model is not expected to make accurate predictions after the lockdown is lifted. As noted in Zhou et al. [2], forecasting models are 'not well suited for long-term predictions'. Mechanistic models based on differential equations circumvent this limitation and can be used to make predictions even when the relevant circumstances change. Their limitation is that it is difficult to solve the associated inverse problem, namely, to determine the relatively large number of parameters involved in these models from the knowledge of epidemiological data.
The first period of the COVID-19 pandemic, known as 'the first wave', was characterized by the fact that a single viral strain dominated the pandemic. This setting provides an ideal situation for the retrospective study of the various mathematical formulations (both forecasting and mechanistic) used in viral epidemics models. Based on this first wave, we previously introduced [3,4] a novel class of forecasting models, which provided accurate long-term predictions for the number of deaths in several European countries caused by the epidemic (figure 1). This was accomplished by using a simple algebraic formula as opposed to the standard exponential formula (logistic model). It is important to emphasize that the logistic formula, as well as similar formulae (see for example [5]) provide accurate fitting for the data of a fixed period. In particular, all these formulae can be used to model the data up to 1 May 2020. However, after fixing the parameters of these formulae using the data of a fixed period, these models failed to provide accurate predictions. On the other hand, after 'training' the algebraic formula of for the same period (i.e. after fixing the relevant parameters using data until 1 May 2020), this formula did provide accurate predictions for a period extending longer than 3.5 months, namely until 1 September 2020.
We show typical comparisons in figure 1, where we plot data for deaths during the first wave in Italy, as well as the predictions associated with three models, the logistic and the two novel models of [3,4]. The model routinely used in epidemiology is the logistic defined in (1.1) below; the rational and the birational formulae are given, respectively, by equations (1.2) and (1.3) below: where X is a fixed parameter X usually taken to be T. Here, D(t) denotes the cumulative number of deaths at time t, with the constant parameters D F , β, γ and δ determined by 'training' (fitting) the model with the data. Similar considerations are valid for the parameters defining the birational formula (1.3), where the time T corresponds to the inflection point of D(t) (details are given in [3,4]). Equation (1.1) shows that, for large t, the expression D(t) approaches its final value D F exponentially. In the case of the rational and birational models, the approach is algebraic. In this sense, it is better to refer to the rational formula as algebraic [3,4]. All three expressions in (1.1-1.3) are particular solutions of a specific differential equation (ODE) of the Riccati type, uniquely specified by the constant D F and by a function a(t). The formulae in (1.2) and (1.3) correspond to the case where a(t) is a rational function of t, whereas the logistic model corresponds to the case that a(t) is a constant. The logistic, rational and birational curves are shown along with real data in figure 1 above.
As shown in figure 1, the birational formula closely matches the actual data throughout the time interval considered, in contrast with the other two curves (rational and logistic). It is interesting that even though 'training' occurred during the lockdown period, its predictions remained accurate after the easing of the lockdown conditions. Possible explanations for this unexpected agreement are discussed in [6].
The fact that the correct asymptotic decay is algebraic rather than exponential suggests a social behaviour origin. Indeed, while most, although not all, physical and biological phenomena exhibit exponential decay in the distribution of key variables, or in the behaviour of correlation functions, this is not the case for social phenomena, where one usually observes algebraic dependences and power laws, perhaps because of selfsimilarity or fractal geometry notions [7]. Indeed, examples of power laws in social settings are ubiquitous, from economics to the distribution of income or of the populations of cities. While such behaviour may also be observed in natural phenomena, e.g. in long-range correlations, or near percolation thresholds in phase transitions [8], the key point regarding epidemiological models is the following: because epidemiological quantities are also affected by social behaviour, e.g. the relaxation of precautions as the epidemic wanes, it is likely that the apparent success of the birational equation is due to such phenomena. Using the mechanistic susceptible-exposed-infectious-recovered (SEIR) model and previously reported results in [9], we will demonstrate this assertion quantitatively in the following section.
The most extensively used mechanistic models in epidemiology are variations of the classical SIR model. This is based on the assumption of a constant total population, consisting of three subpopulations, namely, susceptible, infected and recovered, and three implicit conditions, that the control volume where the epidemic occurs is fixed, there is no influx of additional populations, and spatial gradients do not exist. Following chemical reaction engineering process analogies, Ramaswamy et al. [10] generalized the SIR model to a spatio-temporal model which relaxed all these restrictions; this led to a set of partial differential equations in terms of population densities. This generalization predicted the onset of spatio-temporal waves, as well as other characteristics inherently absent in the SIR model. In addition, it led to the emergence of a single dimensionless number that characterizes the process, the so-called reproduction number R 0 , quantitatively expressed as the ratio of the kinetics of infection to those of recovery [10]. Importantly, R 0 incorporates a combination of physiological and biological parameters, and also reflects measures of social dynamics (e.g. via facial coverings, spatial distancing, areal density, etc.).
If R 0 is constant, SIR predicts that any infection wave wanes asymptotically following an exponential behaviour. In previous work [9], we examined the question of what modifications of the model would lead to an algebraic rather than an exponential asymptotic behaviour, as appears to be the actual case in the first COVID-19 wave. Given that the only relevant parameter in SIR models is R 0 , we showed that for an algebraic behaviour to be obtained, it is necessary that beyond a certain time, R 0 must become an increasing function of the infected fraction denoted by i (ratio of the number of infected individuals to the total). The particular dependence was found in [9] to be where the parameter a is a positive constant, and the exponent n was found to satisfy n > 2. Under these conditions, the asymptotic decay is algebraic and follows a power law in time with an exponent equal to − n/ (n − 1). Equation (1.4) shows that in such cases, as the infection wanes, the reproduction number R 0 is not constant, but rather increases in an algebraic manner with respect to i. The conclusion that the emergence of a power law in time requires that R 0 must increase as the infection wanes, is consistent with the social tendency to relax restrictions (e.g. decrease spatial distancing, increase areal densities, etc.) as the epidemic wanes. Motivated by the fact that the actual data are best matched using an algebraic dependence and that Luhar et al. [9] provide a way to incorporate such behaviour in a simple SIR model, we will extend here this methodology to a more elaborate SEIR-type model, that also includes hospitalized and deceased patients data, neither of which are included in the simpler SIR model of [9].
royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 10: 230858 3 The paper is organized as follows: in the Formulation and results section, we present a modification of the SEIR model that accurately gives rise to an algebraic asymptotic behaviour. We then fit the experimental results on the epidemiological data for the deceased of the first wave in Portugal to the algebraic and birational formulae (1.2) and (1.3). These formulae together with the new model are then used to solve the inverse problem and to determine the various parameters of our full mechanistic novel model. Note that, in a certain sense, the mathematical expressions for our forecasting models are optimal, since an elaborate deep-learning algorithm does not seem to yield better results. In solving the inverse problem, all parameters of the new SEIR-type model are determined from the knowledge of the deceased and hospitalized data. It is shown that the numerical solution matches well the analytical asymptotic results. Furthermore, as elaborated in the Discussion section, this solution provides a confirmation of the accuracy of the solution of the inverse problem.

Formulation and results
The SEIR model previously discussed in several publications [11,12] contains the following populations: and D(t), which denote susceptible, infected, asymptomatic, sick, hospitalized, recovered and deceased individuals, respectively. In this notation, I(t) refers to individuals infected, but not yet infectious, since it is generally assumed that infectiousness requires the elapse of a certain incubation period (typically 5 days) [12]. Then, that individual will either become sick or asymptomatic, both of which are infectious. Importantly, this means that a susceptible individual joins the population I(t) only because of encounters with individuals from populations A(t) or S(t) (but not from I(t)). The various entities are interdependent, as indicated in figure 2: if the daily rates at which an infected person becomes either sick or asymptomatic are s or a, respectively, then each day sI(t) or aI(t) individuals leave population I(t) and enter S(t) and A(t), respectively. Likewise, if asymptomatic individuals recover at a daily rate r 1 , each day r 1 A(t) individuals leave the asymptomatic population and enter the recovered population, R(t). Sick individuals either recover at a rate r 2 or they become hospitalized at a rate h. The hospitalized patients also have two possible outcomes: either they recover at a rate r 3 , or they are deceased at a rate d. Finally, susceptible individuals convert to infected, but not yet infectious, populations, because of their encounter with either asymptomatic or sick individuals (we assume that due to special precautions, hospitalized populations cannot infect) at corresponding rates c 1 A and c 2 S, where c 1 and c 2 denote the respective transmission rates (per person in the respective pools). It is important to note that from all these coefficients, it is only c 1 and c 2 that can be largely affected by social dynamics, the other being largely independent.

The basic formulation
The previous description leads to the following seven ordinary differential equations: royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 10: 230858 where superscript (n) indicates nth time derivative. Total conservation of individuals leads to where T is the sum of all individuals in the control volume prior to the pandemic. Eliminating Σ from equations (2.6) and (2.7) leads to We will next proceed with the assumption that the coefficients a, r 1 , s, h, r 2 , r 3 and d are constants. Then, for arbitrary, and not necessarily constant, values of c 1 and c 2 , equations (2.1-2.7) can be reduced to a system of two ODEs. The final result reads as follows: Populations H, S and I can be expressed in terms of D, via the following equations: Also, equations (2.1) and (2.2) can be written in the form ð2:12Þ Integrating this equation and simplifying, we obtain where R 2 and R 3 are defined by In order to make use of the results in Luhar et al. [9], we next consider the hypothetical case, A (1) = S (1) = 0, in which case equations (2.15) become These are in a form comparable to the SIR model previously described [9], written here as ð2:18Þ The two sets of equations (2.17, 2.18) become self-consistent if we take and 1 a þ s Thus, consistency implies the following dependence for the variables C 1 and C 2: where ρ j , Λ j , j = 1,2, are constants. The subsequent analysis will assume that most of the infection process is via exposure to asymptotic individuals, hence we will take C 1 (I ) >> C 2 (I).
In order to proceed, we will make use of the following results. Let e A and Q be defined in terms of A and D via Then, the functions e A and D satisfy the two ODEs,

Large time asymptotic analysis
With e A and D satisfying equations (2.22) and (2.23), C 1 given by (2.20) and taking C 2 < <C 1 , one can show that D approaches asymptotically the constant value D ∞ , and the dependence of σ on t for large t is not exponential The leading behaviour is given by the following equations: Àn=ðnÀ1Þ t À1=ðnÀ1Þ , m ¼ r 1 r 1 , Àn=ðnÀ1Þ t Àn=ðnÀ1Þ and e A (n À 1) Àn=ðnÀ1Þ Àn=ðnÀ1Þ t Àn=ðnÀ1Þ , t ! 1: Using in this equation, the expression of C 1 from (2.20) we obtain and 1 þ vm þ e L 1 (e s (1) ) 1=n1 þ me s ¼ 0: ð2:32Þ Hence, the requirement that σ decays implies

QED
In the subsequent section, we will fit the full equations to the experimental data with goal to determine the various parameters.
We proceed with the assumption C 2 = 0. First, we will consider fitting the data to the algebraic birational formula in (1.3).

Bidirectional long short-term memory network
Bidirectional long short-term memory (BiLSTM) is a powerful generalization of recurrent neural networks that can capture long-term dependencies while at the same time avoiding the problem of vanishing/exploding gradients. The BiLSTM networks are well suited for time series prediction and can potentially completely capture the contextual information of the time series. The BiLSTM network employed in this work was introduced to predict the number of infected for different countries in [13].

An algebraic formula and deep learning
We used data for the cumulative number of deceased, D(t), and the total number of hospitalized patients, H c (t), from the first wave of COVID-19 in Portugal dating from March to September 2020. During this period, Portugal experienced relatively low numbers of infected, although a larger number of deaths compared with other countries with similar populations, such as Greece. The corresponding data were fitted using the two formulae ð3:1Þ as well as the analogous formulae for the birational model. It can be shown that the following useful formula is valid: and R H denotes the number of people that recovered after they were hospitalized. We used H c , as this was the more accurate data reported. The respective parameters, shown in table 1, were determined using data up to the middle of August and the simplex optimization algorithm. We note that the above formulae provide an optimal fit, given that the use of the efficient deep learning algorithm, BiLSTM network [14], provided similar results. Figure 3 shows the performance of BiLSTM against the logistic and rational models for the number of diseased and hospitalized. Both BiLSTM and the rational had a similar strong correlation to the measured data (R 2~0 .999), with the logistic expression resulting into a smaller value R 2~0 .98.

The solution of the inverse problem
We solved the inverse problem by employing our new mechanistic model together with both the rational and the birational formulae.
We used the fact that

ð3:3Þ
Hence, We used equation (3.1) to match (3.5), and thus, we obtained with the aid of the simplex algorithm, the values for parameters R 3 , d and g, as shown in table 2.
To determine the remaining parameters, we first introduce the notation Equations (3.7), with the aid of the simplex algorithm, yield the values of the remaining constant parameters shown in table 3.
As expected, the value for n satisfies the constraint n > 2.

Numerical verification
Finally, we verified that the numerical solution of our mechanistic model matches, for large t, the rational formula. Similar considerations are valid for the birational model. For this purpose, we used the values of the various constants (tables 2 and 3), obtained from the experimental data. In this connection, it is interesting to note that from the whole set of equations (2.1)-(2.7), we only need to solve the following smaller set of four differential equations for the variables N, Q, e A, M: ð3:9Þ We will take C 2 = 0 and      royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 10: 230858 The function σ can be determined via which yields The solution of this equation with σ(t 0 ) = D(t 0 )/D ∞ − 1 is depicted in blue in figure 5. By plotting in the same figure, the equation depicted in red, we find It becomes clear that for large t, the function σ(t) approaches the r.h.s. of (3.17). The first of equations (2.33) implies the constraint vμ + 1 = 0. Using in this equation the expressions for v and μ given in (2.26) and (2.28) we find the equation Employing the values obtained via the solution of the inverse problem, it can be verified that equation (3.18) is indeed approximately satisfied: using for the constants appearing in the r.h.s. of (3.18) and the numerical values of the constants obtained in tables 2 and 3, we find that the r.h.s. of (3.18) is −0.014 (instead of zero).

Discussion
The forecasting models [2,3] suggest that the dynamics of the first wave of COVID-19 approach an equilibrium state algebraically, as opposed to the exponential behaviour normally predicted by the standard SEIR mechanistic models. Motivated by these results, we introduced a novel mechanistic model which for large times does exhibit algebraic behaviour. This model takes into consideration an additional nonlinear mechanism first understood in the context of SIR models in [9,10]. Specifically, in comparison with the standard SEIR models, such as those previously analysed in [11,12], we allowed for the values of the transmission constants to increase as the infection wanes, which is reflecting the social behavioural tendency to be less cautious as the infection diminishes. The specific relationship used fits well the actual data. As a result of this modification, the asymptotic analysis of the equations defining the new model yields the algebraic behaviour predicted by equations (3.8).
It should be emphasized that several authors have previously published papers addressing the need to incorporate sub-exponential/power law dynamics. Some of these papers are mentioned below. The main novelty of our SEIR-type model, in comparison with the earlier very interesting papers, is the introduction of an epidemiologically motivated mechanistic model which allows to determine the royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 10: 230858 associated large t-dynamics. Indeed, as we hope will become clear in what follows, the previous models concentrated on the early times behaviour. In particular, Chowell et al. [15], addressed early reactive behaviour changes associated with the 2014-2015 Ebola epidemic in West Africa, which exhibited a slower than exponential growth pattern. This was modelled by the equation dc dt ¼ rc p 0 p 1, ð4:1Þ which can also be used to model the spread of a range of pathogens, including the epidemics of influenza, Ebola, foot-and-mouth disease, HIV/AIDS, plague, measles and smallpox. Equation (4.1) was also investigated by Viboud et al. [16]. An extended version of (4.1), namely, the equation was investigated in [17]. Several modifications of the standard SIR epidemic model to support early sub-exponential growth dynamics are discussed by Chowell et al. [15], including the following: (i) metapopulation models, namely characterizing populations in terms of subpopulations (e.g. based on age, vulnerability, etc.), (ii) incorporating reactive behaviour changes by modelling a time-dependent transmission rate, and (iii) incorporating inhomogeneous mixing through a power law scaling parameter. These considerations motivated the introduction of the SIR type model, This model can be simply interpreted as an SIR model with a time-varying reproduction number. The models mentioned above were used in the paper by Chowell et al. [18] in order to fit the early epidemiological data in a number of infections.
Another important contribution of our work is that it addresses a well-known difficulty associated with mechanistic models, namely the determination of the values of the various parameters that enter in the models. In this paper, we present a computationally efficient and robust approach for solving this inverse problem. Moreover, we show that knowledge of the numbers of deceased and hospitalized patients is sufficient to determine all the parameters of the model uniquely and accurately.
The numerical solution presented is useful for two reasons. First, the numerical results confirmed the validity of the asymptotic analysis. Indeed, the numerical solution matched well the asymptotic formula, σ ∼ t −m . Second, they confirmed the validity of the solution of the inverse problem. Indeed, the asymptotic analysis gives rise to the constraint νμ + 1 = 0. Since ν and μ can be expressed in terms of the model parameters determined via the solution of the inverse problem, this equation provides a check of the accuracy of the solution of the inverse problem.
Concluding, it is worth noting that in addition to SEIR models several other types of models have been introduced in connection with COVID-19 [19][20][21][22]. It would be interesting to investigate whether any of these models exhibit large t algebraic behaviour.