key: cord-0313444-e05oyl8r
authors: Paul, S.; Lorin, E.
title: A hybrid approach to predict COVID-19 cases using neural networks and inverse problem
date: 2022-05-17
journal: nan
DOI: 10.1101/2022.05.17.22275205
sha: 72e4e5193e2265d9ed8f50e6c233c4117b3374b7
doc_id: 313444
cord_uid: e05oyl8r

We derive a novel hybrid approach, a combination of neural networks and inverse problem, in order to forecast COVID-19 cases, and more generally any infectious disease. For this purpose, we extract a second order nonlinear differential equation for the total confirmed cases from a SIR-like model. That differential equation is the key factor of the present study. The neural network and inverse problems are used to compute the trial functions for total cases and the model parameters, respectively. The number of suspected and infected individuals can be found using the trial function of total confirmed cases. We divide the time domain into two parts, training interval (first 365/395 days) and test interval (first 366 to 395/ 396 to 450 days), and train the neural networks on the preassigned training zones. To examine the efficiency and effectiveness, we apply the proposed method to Canada, and use the Canadian publicly available database to estimate the parameters of the trial function involved with total cases. The trial functions of model parameters show that the basic reproduction number was closed to unity over a wide range, the first from 100 to 365 days of the current pandemic in Canada. The proposed prediction models, based on influence of previous time and social economic policy, show excellent agreement with the data. The test results revel that the single path prediction can forecast a period of 30 days, and forecasting using previous social and economical situation can forecast a range of 55 days.

correspond to a particular social economic condition. Now, if that particular social economic condition occurs in the future, we 48 can consider the values of the model parameters to be β (t 1 ) and δ (t 1 ). This assists us to develop a prediction model. 49 Finally, we develop the prediction models using Taylor series; there are three terms in the prediction model, previous total 50 cases, influence of the preceding time and effect of social economic policy. According to the choice of prior data/situation, we 51 obtain two types of prediction models. The proposed forecasting model can capture the social economic policy and can use it to 52 forecast the total cases. The estimated values P est T of the parameters, involved in the NN of the trial function of the total cases T a , vary in a 57 wide range of −100 to 300 (Figure 1(a) ), and the trial function T a (t, P est T ) has excellent agreement with the database ( Figure   58 1(b))during the training period, first 365 days (January 27 2020 to January 25, 2021) of the pandemic in Canada. We analytically 59 compute the first (Figure 1 (c)) and second (Figure 1(d) ) order derivatives where d dt (T a (t, P est T )) is non negative, as T a (t, P est T ) is 60 a monotonic increasing function. However d 2 dt 2 (T a (t, P est T )) is an oscillating function and can be negative as well as positive.

Using the trial function of total cases and its first and second order derivatives and the second order differential Eq. (5) for total cases, we obtain trial functions β a (t, P est β ) and δ a (t, P est δ ) for the model parameters β (t) and δ (t) (Figures 2(a) and 2(b)) during the training period of the first 365 days of the current pandemic in Canada. P est β and P est δ are the estimated values of the parameters P β and P δ , involved in the trial functions of β a and δ a , respectively. The range of β a (t, P est β ) for that period is 0.65 to 0.85. However, over a long interval, first 65th day to 365th day, the value of β a (t, P est β ) ranges from 0.65 to 0.7. A similar phenomenon occurs for the trial function δ a (t, P est δ ); the minimum value of δ a (t, P est δ ) during the training period is 0.4, but over a long interval the range of δ a (t, P est δ ) is 0.6 to 0.7. The trial functions E a 1 and E a 2 of the error functions E I and E S range from 1 to 15 (Figure 2(c) ). Using the trial functions β a (t, P est β ) and δ a (t, P est δ ), we can calculate the time dependent basic reproduction number R 0 (t) for the first 365 days of the current pandemic where R 0 (t) = β a (t, P est β ) δ a (t, P est δ )

.

During the interval, first 100th day to 365th day, the value of R 0 (t) was closed to unity (Figure 2(d) ). 62 From Eq. 4 in the Methods section, it follows that the estimation of the infected population can be expressed by I a (t) = 1 δ a (t, P est δ ) d dt T a (t, P est T ) + E a 1 (t, P est E I ) , and the number of susceptible individuals can be obtained using the expression S a (t) = N − I a (t) − T a (t, P est T )

where N is the total population size as mentioned in the Methods section. The number of susceptible individuals decreased 63 rapidly after the first 200 days (Figure 3(a) ), and the curve of infected individuals shows two peaks (Figure 3 (b)), the first peak 64 2/17 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted May 17, 2022. ; Figure 1 . Total number of confirmed coronavirus cases: (a) Estimated values P est T of the NN parameters P T , involved in the trial function of total cases T a (t). The number of hidden unit is 20 so i = 1, 2, · · · , 20. (b) Trial function of total cases T a (t, P est T ) compared to the available data 40 . Blue dots indicate the results obtained from the NN, and red circles are the publicly available data 40 . (c) First order derivatives of the trial function T a (t, P est T ). (d) Second order derivatives of the trial function T a (t, P est T ).

. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted May 17, 2022. ; Figure 2 . Outcomes of NNs and Inverse problem: (a) Trial function β a (t, P est β of the model parameter β (t). P est β is the estimated value of the NN parameters P β . (b) Trial function δ a (t, P est δ of the model parameter δ (t). P est δ is the estimated value of the NN parameters P δ . (c) Trial functions of the error corrections E I and E S . P est E I and P est E S are the estimated values of the NN parameters P E I and P E S involved with the trial functions of E I and E S , respectively. (d) Time dependent basic reproduction number R 0 (t) for the first 365 days of the current pandemic in Canada.

. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted May 17, 2022. ; Figure 3 . Outcomes of NNs and Inverse problem: (a) Number of susceptible individuals obtained from the publicly available data 40 of the total confirmed cases. (b) Number of infected individuals obtained from the publicly available data 40 of the total confirmed cases. (c) The trajectory (daily cases, β a (t, P est β )). (d) The trajectory (daily cases, δ a (t, P est δ )).

. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) is the prediction starting at the 360th day; red circles are the publicly available data 40 . Time domain has been divided into two parts, testing (green + pink region) and test (light blue region). Narrow pink region is to check the validation of the prediction model. indicates average over 10 path predictions. (e) T est paths denotes a bundle of 30 path predictions starting points 365th, 364th, · · · , 336th days of the current pandemic in Canada. T est av30 indicates average over 30 path predictions. (f) T est paths denotes a bundle of 50 path predictions starting points 365th, 364th, · · · , 316th days of the current pandemic in Canada. T est av50 indicates average over 50 path predictions.

. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted May 17, 2022. ; T a presents trial function of total cases; T est SE is the prediction using previous social economic condition; red circles are the publicly available data 40 . Time domain has been divided into two parts, training (green + pink region) and prediction (light blue region). Narrow pink region is for check the validation of the prediction model. (c) Time scale is the first 350 to 450 days of the current pandemic in Canada, focus on the validation (pink) as well as test (light blue) regions.

. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted May 17, 2022. ; https://doi.org/10.1101/2022.05.17.22275205 doi: medRxiv preprint due to α-variant and the second peak due to β -variant. The trajectories (β a (t, P est β ), daily cases) (Figure 3 (c)) and (δ a (t, P est δ ), daily cases) (Figure 3 (d)) show that after the first 250 days the number of daily cases increase even though the trial functions 66 β a (t, P est β ) and δ a (t, P est δ ) vary in a small range. This phenomenon occurred due to larger number of infected individuals after the first 250 days.

approach 'parameters at the fixed point'; a detail is provided in the Methods section. The single path prediction T est 360 is started at 72 the 360th day of the current pandemic and is extended up to 450th day (Figure 4 (a)) along with the trial function of total cases 73 T a . We obtain excellent agreement between T est 360 and the database in the validation region (Figure 4(b) ). In the initial part of 74 the test domain, the first 365 to 370 days, both curves, T est 360 and T a , have good agreement with the database (Figure 4 (c)), and 75 after the initial part, beyond the first 370 days, T a does not agree with the data. However, we observe an excellent agreement 76 between T est 360 and the data for a wide range of 30 days, the first 366 to 396 days, and at 396th day the difference between T a and 77 the data is 20,000 which is significantly large. Now, instead of single path prediction, we present the outcome of the multiple The COVID-19 timeline is a combination of several brands of various social economic policies, such as reopening, lockdown, curfew, released curfew, etc. A small part of that timeline ( Figure 5 (a)), August 2020 to February 2021, shows that due to the increase in public gatherings and the reopening of schools, new COVID-19 cases increased, and to control the situation, the government reimposed various restrictions. After several months in which the situation was under control, the government released those restrictions again. One can observe that pattern multiple times during the current pandemic. It seems that the situation on 382th day of the pandemic was similar to that on 228th day (154 days ago), reopening and releasing restrictions, and we can expect that the social economic policy just after 382th day would be the same as the social economic policy just after 228th days. Therefore, we can assume that the model parameters β (t) and δ (t) satisfy the conditions:

where t ≥ 382; a detailed is provided in the Methods section. In the calculation we use trial functions β a (t, P est β ) and δ a (t, P est δ ) 86 instead of the model parameters β (t) and δ (t), respectively. We divide the time domain into three parts; training region, the 87 first 395 days, where we train the trial functions T a (t), β a (t) and δ a (t); validation region, a small interval of the first 382 88 to 395 days, where we can justify the prediction model T est SE ; prediction region, an interval of the first 396 days to 450 days.

Here, we describe the outcomes T est SE employing the approach 'parameters map previous social economic situation'; a detail 90 is provided in the methods section. We obtain an excellent agreement between T est 360 and the database in the validation region 91 ( Figure 5 (b)). We observe outstanding agreement between T est SE and the database in the validation region ( Figure 5 (c)). We 92 observe an outstanding agreement ( Figure 5 (c)) between T est SE and the data for a wide range of 54 days, the first 396 to 450 days, 93 and at 450th day there is a large discrepancy between T a and the data.

The function F β , δ ,

, the third term of the prediction model, is equal to d 2 dt 2 T a (t, P est T ), which can be 96 either positive or negative (Figure 1(d) similar to those of the common flu that confused the population. As a consequence, the value of the function δ a (t, P est δ ) was corresponds to a certain socioeconomic condition; in other word, the trial functions of the model parameters β (t) and δ (t) can capture the social economic policy.

The 'prediction using previous data' T est s (t) is based on the socioeconomic conditions of the starting time s and forecasting 109 for the time domain t > s. The prediction shows the expected outcomes if the socioeconomic conditions remain the same as 110 at the start time; in single path prediction, we consider β a (t) = β a (s) and δ a (t) = δ a (s) for t > s. Since β a (t) and δ a (t) are The lockdown/socioeconomic restriction and reopening/release of the socioeconomic restriction happened and are happening 117 one after another throughout the current pandemic. The repetition of the events is the principle concept to use the previous 118 situation to predict the continuing scenario. Naturally, the proceeding situation will not be exactly the same as the previous 119 situation. However, the earlier situation can help to predict the growing situation. In this article, we proposed a hybrid approach to forecast the COVID-19 cases using previous database of total cases. For this 122 purpose, we derived a second order nonlinear differential equation, involved the time dependent model parameters, for the total 123 cases; this differential equation is the foundation of the present study. The trial function of total cases has been generated using 124 a single layer NN and publicly available database. The first and second order derivatives of the total cases have been evaluated 125 analytically using that trial function. Using inverse problem technique on the second order nonlinear differential equation is free to solve any system of differential equations because of that the proposed hybrid method is computationally simple. In 131 future work, we will consider COVID-19 dynamics for two neighbouring states with moving population. One can study the 132 recovery and death toll rates using the second order differential equations mentioned in the Methods section. The approach 133 used in this work could be applied in the study of other infectious disease transmission dynamics as well.

In this section, we introduce a compartment based infectious disease model, extending the standard SIR model, including a total Thereafter, we calculate the derivatives of that trial function analytically.

•

Step 3 : Solving the inverse problem for the differential equation, from Step 1, using the trial function and its derivatives,

Step 2, we obtain trial functions for the time dependent model parameters, from Step 1.

•

Step 4 : Finally, we develop several prediction models for total cases using Taylor's expansion. The variety of the 147 prediction models is based on the choice of the time dependent model parameters. . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted May 17, 2022. ; https://doi.org/10.1101/2022.05.17.22275205 doi: medRxiv preprint • The total population is constant (neglecting the migrations, births and unrelated deaths) and initially every individual is 156 assumed susceptible to contract the disease.

• The disease is spread through the direct (face-to-face meeting) or indirect (through air current, common used or delivery 158 items like door handles, grocery products) contact of susceptible individuals with the infected individuals. 

In the model, we assume that there is no overlap between these two compartments, infected(I) and total cases (T ). In other 170 words, tested corona-positive individuals are assumed to be unable to substantially spread the disease due to isolation and are 171 immune to re-infection after recovery 41 . The time-dependent model is the following set of coupled delay differential equations:

where β (t), δ (t), λ (t) and µ(t) are real positive parameters, respectively, modeling the rate of infection, the rate of tested corona-positive, the rate of recovery and the rate of death toll. It follows from (1), that for any time t

where N (constant) is the total population size. Using the fact, total cases T = C + R + D the system of equations (1) can be written as

It follows from third equation of (3)

where E I is the model correction associated with I. The cause of model correction is the inherent simplifications and approximations that are made in postulating a mathematical model of a physical phenomenon, as even high-fidelity models 10/17 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted May 17, 2022. ; 

. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted May 17, 2022. ; https://doi.org/10.1101/2022.05.17.22275205 doi: medRxiv preprint simply approximate reality. These simplifications may arise from lack of knowledge about the system and/or computational constraints. Eliminating S from (3) and using (4), we obtain a second order nonlinear ordinary differential equation for T involving the time dependent parameters β (t) and δ (t) as

and E S is the model correction associated with S.

It follows from the last three equations of (1) and third equation of (3):

and

where E R , E D and E C are the model corrections. The above two Eqs. (7) and (8) can be used to calculate the recovery rate λ (t) 174 and the rate of death toll µ(t). However, the present work is devoted for forecasting the COVID-19 cases so we merely focus 175 on (5), the second order differential equation for total cases.

The trial function 177 In this section, we construct a trial function for T (t), and the parameters of this function are estimated using publicly available database. The trial function for T (t), satisfied the Eq. (5), may be written in the following form:

above function has two parts, the first part satisfies the initial condition (at t = t 0 ) and contains no adjustable parameters, the second part involves a NN N T (t, P T ) (Figure 6(b) ) containing adjustable parameters P T . Here, T d (t 0 ) indicates the observed data, total cases at time t 0 , obtained from the publicly available database. The NN N T (t, P T ) that takes a single input value of time t is defined as follows:

where the sigmoid activation function

and

The NN N T (t, P T ) has one hidden layer with H T activation functions and a linear activation function in output unit. Here, w j is 178 a weight parameter from input layer to the jth hidden layer, b j is a bias for the jth hidden layer, p j is a weight parameter from 179 the jth hidden layer to the output layer. Here, we choose the sigmoid activation function, presented in Eq. (11), because T (t) is 180 a monotonic increasing function.

The NN should be trained using the publicly available database. Therefore, the time domain is divided into two parts which are training, t 0 ≤ t ≤ t 1 , and test, t 1 < t ≤ t 2 (Figure 6(c) ). Using the publicly available data {T d i } M i=1 corresponding the M discrete points in time domain [t 0 ,t 1 ] with step length h, the NN is trained by finding the optimal values of the weight and bias parameters (p j , w j and b j ) which provide the minimum value of the root mean square error function. The error function E T (P T ) has the following form:

where 

In this section, we derive the trial functions for β (t) and δ (t), and the parameters of these trial functions are estimated using inverse problem techniques. After computing the trial function T a (t, P est T ), we calculate the first and second order derivative of T a (t, P est T ) analytically for the time domain [t 0 ,t 1 ]. The trial function for β (t) defined in the time domain [t 0 ,t 1 ], involved in the Eq. (5), may be written in the following form:

above NN N β (t, P β ) contains adjustable parameters P β . The NN N β (t, P β ) that takes a single input value of time t is defined as follows:

where the hyperbolic activation function

and

The NN N β (t, P β ) has one hidden layer with H β activation functions and a linear activation function in output unit. Here, u j is a weight parameter from input layer to the jth hidden layer, d j is a bias for the jth hidden layer, q j is a weight parameter from the jth hidden layer to the output layer. Similar type of NNs can be used as the trial functions for δ (t), E I and E S , involved in the Eq. (5). The error function associated with the inverse problem that must minimize, has the following form:

where P β ,P δ ,P E I and P E S are the parameter sets associate with the trial functions β a , δ a , E a I and E a S , respectively. The 186 discrete points y i for i = 1, 2, · · · , J can be expressed as y i = t 0 + (i − 1)∆t and y 1 = t 0 , y J = t 1 . Here, error function 187 E IP (P β , P δ , P E I , P E S ) for J equally spaced points with step length ∆t inside the interval [t 0 ,t 1 ] is trained. The sufficient 188 small values of E IP (P β , P δ , P E I , P E S ) indicates the estimated values of P β ,P δ ,P E I and P E S denoted by P est β ,P est δ ,P est E I and P est E S , 189 respectively.

In this section, we develop prediction models for the total number of cases using Taylor series. Suppose that total confirmed 192 cases at time t 1 i.e., T (t 1 ) is given, and we want to predict T (t 1 + ∆t) i.e., total confirmed cases at time t 1 + ∆t. In this context, 193 we can use Taylor series method for T (t 1 + ∆t).

194

where T (t 1 ) and T (n) (t 1 ) are the first and n-th order derivatives. Here, we use the fact T (n+2) = F (n) , obtained from the Eq. (5).

The right hand side of the final expression for T (t 1 + ∆t), presented in Eq. (19), has three terms ( Figure 6(d) ). For simplicity we neglect all the derivatives of the function F(t 1 ) and replacing T (t 1 ) by dT a dt (t 1 ), and obtain an estimated value for T (t 1 + ∆t), denoted by T est (t 1 + ∆t), as follows

where at time t 1 T est (t 1 ) = T (t 1 ) and dT est dt = dT a dt . Again using Taylor series expansion, we can express

Using the recurrence relations, Eqs. (19) and (20), we can predict the total corona positive cases for the time domain t ≥ t 1 + ∆t.

According to the choice of the trial functions β a and δ a , we can derive two types of prediction model.

• To consider the fixed values of β a and δ a , i.e., for the prediction domain t ≥ t 1 + ∆t,

and δ a (t) = δ a (s) where s belongs to the training domain.

• To consider the previous time evolution for β a and δ a , i.e., for the prediction domain t ≥ t 1 + ∆t,

and δ a (t) = δ a (t − τ) where t − τ belongs to the training domain. Here, β a and δ a correspond to the previous social economic situation. Here, previous data means that we consider the same social economic condition as the starting time s through out the calculation. As a consequence in Eqs. (19) and (20), the functions β a (t, P est β ) and δ a (t, P est δ ) are replaced by β a (s, P est β ) and δ a (s, P est δ ), respectively, for t 0 < s ≤ t 1 . We can obtain an estimated value for T est s , as follows 

where s ≤ t ≤ t 2 . Here, we neglect the time evolution of the functions E a I and E a S .

The starting points for the estimation can be t 1 − k∆t, k ≥ 0, and the prediction starting through the point t 1 − k∆t can be defined as k-th path and denoted by T est t 1 −k∆t (t) for t 1 < t ≤ t 2 . The first K paths are as follows:

The bundle of the first K paths can be denoted by T est paths . Taking average over first K paths, we obtain estimation of total cases T over K paths as follows:

where t 1 < t ≤ t 2 .

. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) In this section, we demonstrate the prediction model when the trial functions β a and δ a correspond to the previous Social Economic (SE) situation. Suppose the time domains for training and test are [t 0 ,t 1 ] and [t 1 + ∆t,t 2 ] (Figure 6(c) ), respectively, and we use the Eqs. (19) and (20) and a previous SE situation to estimate the values of T , defined as T est SE , for the second domain. Here, we assume that the social economic condition of the period [t 1 + ∆t,t 2 ] is similar with the social economic condition of τ days ago. As a consequence in Eqs. (19) and (20), the functions β a (t, P est β ) and δ a (t, P est δ ) are replaced by β a (t − τ) and δ a (t − τ), respectively, for t 1 + ∆t ≤ t ≤ t 2 and t 0 ≤ t − τ ≤ t 1 . We can obtain an estimated value for T est SE , as follows

where at time t 1 , T est SE (t 1 ) = T (t 1 ) and

. Again using Taylor series expansion, we can express

where t 1 < t ≤ t 2 . Here, we neglect the time evolution of the functions E a I and E a S . length ∆t = 0.1. We obtain an estimated value P est β of P β , P est δ of P δ , P est E I of P E I , P est E S of P E S for corresponding the value of the 219 error function E IP (P est β , P est δ , P est E I , P est E S ) = 61. For the prediction of total confirmed cases using Taylor series expansion method, 220 we consider the step length ∆t = 0.1.

one 17, e0262708 (2022). . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The data used to estimate the model parameters are publicly available and are available with the code. 

Rolling updates on coronavirus disease (COVID-19)

Modeling and forecasting the early evolution of the COVID-19 pandemic in Brazil

Simulation of coronavirus disease 2019 (covid-19) scenarios with possibility of reinfection

Evaluating ideologies of coronacrisis-related self-isolation and frontiers 247 closing by sir compartmental epidemiological model

Distribution of incubation periods of COVID-19 in the

Estimation of covid-19 recovery and decease periods in canada using delay model

Analysis of a vector-borne diseases model with a two-lag delay differential equation. The 253 North Carol

New approximations, and policy implications

Stability analysis of an age-structured seirs model with time delay

Elementary time-delay dynamics of COVID-19 disease. medRxiv

A time delay dynamic system with external source for the local outbreak of 259 2019-ncov

Solvable delay model for epidemic spreading: the case of COVID-19 in Italy

Effects of age and sex on recovery from COVID-19: Analysis of 5769 israeli 262 patients

COVID-19 pandemic and its recovery time of patients in India: A 264 pilot study

COVID-19 in a designated infectious diseases hospital outside Hubei Province

Comparisons of viral shedding time of SARS-CoV-2 of different samples in ICU and non-ICU patients

Prolonged presence of SARS-CoV-2 viral RNA in faecal samples

Epidemiology and transmission of COVID-19 in 391 cases and 1286 of their close contacts in Shenzhen, 272 China: a retrospective cohort study

Predictors of the prolonged recovery period in COVID-19 patients: a cross-sectional study

Inferred duration of infectious period of SARS-CoV-2: rapid scoping review and analysis of available 276 evidence for asymptomatic and symptomatic COVID-19 cases

Increase in covid-19 cases and case-fatality and case-recovery rates in europe: A cross-temporal meta-analysis

Covid-19 forecasting based on an improved interior search algorithm and multilayer 280 feed-forward neural network

Deep learning via lstm models for covid-19 infection forecasting in india. PloS 36. Pacheco, C. & de Lacerda, C. Function estimation and regularization in the

This research is supported by the Natural Sciences and Engineering Research Council of Canada (NSERC).