Next Article in Journal
Refined Green–Lindsay Model for the Response of Skin Tissue under a Ramp-Type Heating
Next Article in Special Issue
Reducing the Reality Gap Using Hybrid Data for Real-Time Autonomous Operations
Previous Article in Journal
The Global Property of Generic Conformally Flat Hypersurfaces in R4
Previous Article in Special Issue
Machine-Learning Methods on Noisy and Sparse Data
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Modeling COVID-19 Using a Modified SVIR Compartmental Model and LSTM-Estimated Parameters

Departamento de Ingeniería Geológica y Minera, Escuela Técnica Superior de Ingenieros de Minas y Energía, Universidad Politécnica de Madrid, Ríos Rosas, 21, 28003 Madrid, Spain
*
Author to whom correspondence should be addressed.
Mathematics 2023, 11(6), 1436; https://doi.org/10.3390/math11061436
Submission received: 28 February 2023 / Accepted: 8 March 2023 / Published: 16 March 2023
(This article belongs to the Special Issue Mathematical Modeling, Optimization and Machine Learning)

Abstract

:
This article presents a modified version of the SVIR compartmental model for predicting the evolution of the COVID-19 pandemic, which incorporates vaccination and a saturated incidence rate, as well as piece-wise time-dependent parameters that enable self-regulation based on the epidemic trend. We have established the positivity of the ODE version of the model and explored its local stability. Artificial neural networks are used to estimate time-dependent parameters. Numerical simulations are conducted using a fourth-order Runge–Kutta numerical scheme, and the results are compared and validated against actual data from the Autonomous Communities of Spain. The modified model also includes explicit parameters to examine potential future scenarios. In addition, the modified SVIR model is transformed into a system of one-dimensional PDEs with diffusive terms, and solved using a finite volume framework with fifth-order WENO reconstruction in space and an RK3-TVD scheme for time integration. Overall, this work demonstrates the effectiveness of the modified SVIR model and its potential for improving our understanding of the COVID-19 pandemic and supporting decision-making in public health.

1. Introduction

On 31 December 2019, the Country Office in the People’s Republic of China of the World Health Organization (WHO) acquired media reports by Wuhan Municipal Health Commission on cases of viral pneumonia. On 10 January 2020, the first draft genome of the novel coronavirus was made public, this virus was denominated SARS-CoV-2, and it is the cause of Coronavirus disease 19 (COVID-19) [1]. In March 2020, due to the spread and severity of the infectious disease, COVID-19 was declared a pandemic [1]. Different stages have been observed throughout the pandemic, including developing and launching vaccines with high effectiveness, emerging new strains of concern, massive vaccination in developed countries, and significant lag in others. According to [2], up to 14 October 2022, there have been over 13.5 million reported cases and over 114 thousand deaths in Spain.
Compartmental models, such as the well-known SIR model, are commonly used in epidemic simulations to model the dynamics of Susceptible, Infected, and Recovered compartments. In previous studies, the SIR model with vaccination and immunization was proposed. The vaccination models include those with partial and waning immunity, pulse vaccination, and vaccination in newborns, which are modeled using two rates related to vaccination [3,4,5,6]. General vaccination models for homogeneous groups in the population were proposed in [7].
For more detailed vaccination policies and vaccine-related parameters, ref. [8] studied vaccine coverage, the average duration of immunity acquired by a vaccine, and vaccine leakiness. Additionally, ref. [9] proposed two models with delay, one with and without vaccination, to establish comparisons of the vaccination effort effects.
Overall, these models provide insight into the effects of vaccination on the spread and control of infectious diseases and can be used to inform public health policy.
The incidence rate is a critical characteristic of compartmental models that drives infections. The most frequently used incidence rates are the bilinear incidence rate and the standard incidence rate. Various incidence rates have been applied in different models and for different applications, such as the Holling-type incidence function [10], the saturated incidence rate [11,12], and other types, including those originally implemented in host–parasite population models, such as the saturated incident rate types in [13]. Some examples of applications for the different incident rate types are: [14] where bilinear incidence rate is used in a model considering vaccination, a Holling-type incidence function for a S E I R model in [10], a standard incidence rate in [15] where there are different coefficients considered for each defined infected sub-population (infected, asymptomatic, and dead infectious corpses), a saturated incidence rate in a S I R model used in [11]. Recent studies have applied compartmental models to simulate the COVID-19 outbreak in Spain and Italy, using nonlinear data-fitting approaches [16] and Runge–Kutta methods [17,18]. Many recent references dealing with mathematical models and numerical simulation of COVID-19 can be mentioned, such as [19,20,21,22]. In [16], a study on a modified SEIR model applied to the COVID-19 outbreak in Spain and Italy, solving the ODE system using a nonlinear data-fitting approach minimizing a least-squares error function.
Artificial Neural Networks (ANNs) have been used as efficient tools for estimating suitable values of parameters in epidemiological models. The first reference found with an application relating Artificial Neural Networks to epidemiological models is [23] and dates back to 1995. Even though other more recent applications [24,25,26] have been published, their use in this task is not widespread, therein lies the opportunity to propose their application as a parameter estimation technique. Several studies have applied ANNs to estimate the rate of infection and spread of an epidemic, including cooperative and supportive neural networks [24]. Some recent studies have used Long short-term memory (LSTM) ANNs to forecast time-varying parameters for a SIRD model for COVID-19 [26,27].
In this work, we modified a model that is a variant of the classical SIR model, including vaccination, called the SVIR model. The model’s formulation, equilibrium points, and reproduction number were studied, and the contact rate ( β ) was estimated using LSTM ANNs. The numerical simulations were conducted using the fourth-order Runge–Kutta scheme. A PDE model with diffusion was also included and studied, and numerical simulations were performed.
The results of the numerical simulations were discussed and compared with actual data for the ODE models, and insights from the PDE model implementation were presented. The study demonstrated the effectiveness of LSTM ANNs in estimating the contact rate and the suitability of compartmental models for simulating the COVID-19 pandemic.
In conclusion, this study modified a compartmental model, including vaccination, and applied ANNs to estimate the contact rate. The study also demonstrated the effectiveness of compartmental models and ANNs in simulating the COVID-19 pandemic. The results provide insights into the dynamics of the pandemic and suggest the importance of vaccination efforts in controlling the spread of the virus. This study’s proposed modifications to the SVIR model, including time-varying parameters, may offer additional insights into the dynamics of the pandemic.
The results of the numerical simulations for the PDE model with diffusion suggest that the diffusion term may play a role in reducing the spread of the virus, particularly at the early stages of the outbreak. This finding may have implications for policy decisions, such as the timing of lockdowns and social distancing measures, and highlights the importance of including spatial effects in epidemiological models.
Overall, this study’s use of ANNs and numerical simulations provides a valuable tool for understanding the COVID-19 pandemic’s spread and control. The application of these tools can help policymakers make informed decisions about vaccination and other measures to slow the virus’s spread, ultimately reducing the number of cases and deaths worldwide.

2. ODE SVIR Compartmental Model

In this section, we present the ODE SVIR compartmental model and its modifications. We begin with the model formulation and then discuss its positivity, equilibrium, and reproduction number. Next, we describe how the contact rate is estimated using LSTM, followed by a discussion on numerical simulations.

2.1. ODE Model Formulation

SVIR Compartmental models incorporate vaccinated individuals into the classical SIR model; this is relevant since vaccination greatly influences the evolution of an epidemic or pandemic. In the SVIR compartmental model, S is for individuals susceptible to the infection, V is for vaccinated, I is for infected, and R is for retired individuals (which in this case includes recovered and deceased). A simplified version of the SVIR compartmental model is considered a baseline for this work’s proposed techniques; this version is one of the cases in [7] where four compartments describe the dynamics of the epidemic in time, vital dynamics are not included, and reinfection is not considered. Even though the model shown in Equation (1) is a simplification with multiple assumptions, it is helpful in broadly studying the changes in population dynamics of epidemics involving vaccination.
d S d t = β S I N α S d V d t = α S d I d t = β S I N γ I d R d t = γ I
In system (1), α is the vaccination rate, β is the transmission rate and γ is the recovery rate. The model divides the total population (N) into these four compartments (S, V, I, R) such that S + V + I + R = N , and the values for these compartments are always positive ( S 0 , V 0 , I 0 , R 0 ).
It is important to state that several assumptions are considered in this model:
  • The contact rate ( β ), vaccination rate ( α ), and recovery rate ( γ ) are positive.
  • The recruitment rate of new susceptible individuals (such as newborns or migration into the established population) is equal to the death rate or migration outside of the established population and, therefore, not included in the model.
  • Once an individual is vaccinated with two dosages (these dosages are given to an individual with a delay, but the delay is not included in the model), the individual cannot be infected.
As mentioned in the introduction, Section 1, different versions of the incidence rate have been implemented in compartmental models. Modifications of the standard incidence rate are particularly relevant for the specific case of COVID-19 due to the impact of changes in restrictions and human behavior based on their perceived risk by implementing more personal cautions, which affect the rate of infections. Without explicitly defining parameters representing governmental restrictions and changes in human behavior, an alternative is to modify the incidence rate via parameters that regulate the model based on prior infection trends.
This work proposes modifications to the base S V I R model shown in system (1); the modified model is as follows:
d S d t = ( 1 ρ 2 ) ρ 1 β 1 + α 1 S + ( 1 ρ 1 ) β + ρ 2 α 2 β S I N α S d V d t = α S d I d t = ( 1 ρ 2 ) ρ 1 β 1 + α 1 S + ( 1 ρ 1 ) β + ρ 2 α 2 β S I N γ I d R d t = γ I
This system can be simplified by expressing the incidence rate as a function of S, denoted κ ( S ) :
κ ( S ) = ( 1 ρ 2 ) ρ 1 β 1 + α 1 S + ( 1 ρ 1 ) β + ρ 2 α 2 β
This modified model combines a saturated incidence rate (such as one mentioned in [12]) with a modified standard incidence rate to generate a mechanism that can decrease or increase the infection rate based on two parameters acting as logical gates: ρ 1 and ρ 2 . In an epidemiological context, these parameters can be considered as reflecting varying restrictions or changes in the spread of the disease due to the appearance of new variants that vary in contagiousness. Both scenarios have been observed at different points throughout the COVID-19 pandemic. The modification introduces two additional parameters, α 1 and α 2 , for regulating the incidence rate: α 1 is used to decrease the incidence rate, while α 2 has the opposite effect.
The parameter ρ 1 is used to decrease the incidence rate and is defined as follows:
ρ 1 ( t ) = 1 , i f d β d t δ 1 ( t ) > 0 0 , i f d β d t δ 1 ( t ) 0
The threshold at which the decrease in the incidence rate takes effect is defined by Equation (4) and is time-dependent. This threshold is related to outlier detection, which involves considering the median and standard deviation of the estimated β parameter. These parameters display highly different behaviors at different points of the epidemic, particularly at the onset of the epidemic and when new waves of infection are triggered.
δ 1 ( t ) = m e d β a b s ˜ + a σ β a b s ˜ .
Equation (5) introduces the time-dependent threshold δ 1 ( t ) , which determines the point at which the decrease in the incidence rate takes effect. This threshold is related to the detection of outliers in the estimated β parameter. Specifically, we calculate the median and standard deviation of the absolute daily changes in β , denoted by β a b s ˜ = d β d t . The constant a can be customized for each Autonomous Community, and we provide the determined values in Table A1 in Appendix A.
Similarly, ρ 2 serves as a mechanism to increase the incidence rate, and it is regulated by the recent changes in the infection trend d I d t . Its definition is as follows:
ρ 2 ( t ) = 1 , i f 1 i j = 1 i d I d t | n j δ 2 > 0 0 , i f 1 i j = 1 i d I d t | n j δ 2 0
Equation (5) sets a time-dependent threshold for the decrease in incidence rate, with β a b s ˜ = d β d t representing the absolute value of the daily changes for the estimated β parameter. The median ( m e d β a b s ˜ ) and standard deviation ( σ β a b s ˜ ) of this series are calculated, with the constant parameter a specified for each Autonomous Community. Table A1 in Appendix A includes the determined values.
Similarly, ρ 2 is regulated by the most recent changes in the infection trend d I d t , with its definition provided by Equation (6). This mechanism is designed to increase the incidence rate and depends on the i time steps preceding the current one, where i is an integer. The rate of infected individuals at time t n j is denoted by d I d t | n j , and the average at each particular time step is calculated by division with i. The parameter δ 2 sets a threshold of action for the ρ 2 parameter. Appendix B in Table A2 includes the specific values per Autonomous Community determined in this study.
In Equation (2) for the specific case where ρ 2 = 0 for any t, the dynamic system is simplified as:
d S d t = ρ 1 β 1 + α 1 S + ( 1 ρ 1 ) β S I N α S d V d t = α S d I d t = ρ 1 β 1 + α 1 S + ( 1 ρ 1 ) β S I N γ I d R d t = γ I
Equation (7) provides a simplified version that only includes a mechanism for reducing the incidence rate, which is often the most relevant scenario. Numerical simulations in Section 2.6 support this assumption by showing that this version is typically sufficient.
To analyze the model’s behavior at the equilibrium, the model is extended to include a recruitment rate for the Susceptible compartment, ϕ , and a natural death rate ω for all compartments.
d S d t = ϕ ( 1 ρ 2 ) ρ 1 β 1 + α 1 S + ( 1 ρ 1 ) β + ρ 2 α 2 β S I N α S ω S d V d t = α S ω V d I d t = ( 1 ρ 2 ) ρ 1 β 1 + α 1 S + ( 1 ρ 1 ) β + ρ 2 α 2 β S I N γ I ω I d R d t = γ I ω R
where d N d t = ϕ ω N ( t ) , and ω and ϕ are positive constants.

2.2. Positivity of the ODE Model

Theorem 1. 
All solutions S ( t ) , V ( t ) , I ( t ) , R ( t ) for the model in Equation (8) with initial conditions S ( 0 ) 0 , V ( 0 ) 0 , I ( 0 ) 0 , R ( 0 ) 0 , are non-negative for all t 0 .
Proof. 
To demonstrate the positivity of the ODE model we follow the strategy put forward in [28]. The model is given by the expression
d S ( t ) d t = ϕ S ( t ) κ ( S ) I N + α + ω ,
where κ ( S ) is given by Equation (3).
We first multiply Equation (9) by E ( t ) = exp 0 t κ ( S ) I N + α + ω d τ to get
d S ( t ) d t E ( t ) + S ( t ) κ ( S ) I N + α + ω E ( t ) = ϕ E ( t ) ,
so we have
d d t S ( t ) E ( t ) 0 t κ ( S ) I N + α + ω d τ = ϕ E ( t ) .
Integrating Equation (11) in [ 0 , t ]
S ( t ) exp 0 t κ ( S ) I N + α + ω d τ S ( 0 ) = 0 t ϕ E ( t ) d t .
Therefore we obtain
S ( t ) = S ( 0 ) exp 0 t κ ( S ) I N + α + ω d τ + exp 0 t κ ( S ) I N + α + ω d τ 0 t ϕ exp 0 t κ ( S ) I N + α + ω d σ d τ > 0 .
Similarly it can be established that V ( t ) 0 , I ( t ) 0 , and R ( t ) 0 , for all t 0 .    □

2.3. ODE Model Disease-Free Equilibrium and Reproduction Number

For Equation (8), at the disease-free equilibrium (DFE), the rate of change of each compartment is zero, that is, d S d t = 0 , d V d t = 0 , d I d t = 0 , d R d t = 0 . Equations (4) and (6) imply that ρ 1 and ρ 2 are both 0 at the DFE. At the DFE, I 0 = 0 . For S, from d S d t = 0 , we obtain S 0 = ϕ α + ω . For V, from d V d t = 0 , we obtain V 0 = α ϕ ω ( α + ω ) , and R 0 = 0 . Therefore, the DFE is E 0 = ϕ α + ω , α ϕ ω ( α + ω ) , 0 , 0 .
To calculate the Basic Reproduction Number R 0 for the proposed SVIR model in Equation (8), we use the next-generation matrix method to derive an expression for R 0 based on the I compartment. We linearize the differential equation for I as d x d t = ( F 1 + V 1 ) , where x = I , and calculate the matrices F 1 and V 1 .
F 1 = 1 ρ 2 ρ 1 β α 1 S + 1 + 1 ρ 1 β + ρ 2 α 2 β S I
V 1 = [ γ I + ω I ] .
These expressions are differentiated with respect to I and evaluated with the values at the DFE ( E 0 ). Therefore:
F = 1 ρ 2 ρ 1 β α 1 ϕ α + ω + 1 + 1 ρ 1 β + ρ 2 α 2 β ϕ α + ω
V = γ + ω
The inverse of V is needed:
V 1 = 1 ω + γ
An expression for the Reproduction Number is obtained with F V 1 :
R 0 = β ϕ ( α + ω ) ( ω + γ )
This equation provides a useful measure of the potential for disease spread in the proposed SVIR model, taking into account key parameters such as the contact rate, recovery rate, and vaccination rate.

2.4. ODE Model Local Stability

In the system (8), it is observed that the equation for d R d t is independent of the rest. Taking this into consideration, for the purposes of the local stability analysis, the system is reduced to:
d S d t = ϕ ( 1 ρ 2 ) ρ 1 β 1 + α 1 S + ( 1 ρ 1 ) β + ρ 2 α 2 β S I N α S ω S d V d t = α S ω V d I d t = ( 1 ρ 2 ) ρ 1 β 1 + α 1 S + ( 1 ρ 1 ) β + ρ 2 α 2 β S I N γ I ω I
Theorem 2. 
The disease-free equilibrium of the system defined in Equation (20) is locally unstable when R 0 < 1 .
Proof. 
The Jacobian matrix is calculated for the system in Equation (20) and evaluated for at the DFE ( E 0 ), with the following outcome:
J E 0 = α ω 0 β ϕ α + ω N α ω 0 0 0 β ϕ α + ω N γ ω
The characteristic equation is given by det ( λ I J E 0 ) = 0 , where I is the identity matrix, and its solution is a cubic equation of the form λ 3 + a 2 λ 2 + a 1 λ + a 0 = 0 , where
a 2 = R 0 1 γ + R 0 3 ω α a 1 = 2 R 0 3 ω 2 + R 0 2 α + 2 R 0 1 γ ω + α γ R 0 1 a 0 = ω + γ R 0 1 ω α + ω
The terms a 0 , a 1 , a 2 , in Equation (22), will always be negative when R 0 < 1 , there are no conditions in which the Routh–Hurwitz criterion is met. Therefore when R 0 < 1 , the disease-free equilibrium of the system in Equation (20) is locally unstable.    □
Theorem 3. 
The system defined in Equation (20) is a saddle point at the endemic equilibrium point  E * .
Proof. 
The Jacobian evaluated at E * is as follows:
J E * = j 11 j 12 j 13 j 21 j 22 j 23 j 31 j 32 j 33
where,
j 11 = 1 ρ 2 ρ 1 β α 1 S * I * S * α 1 + 1 2 N 1 ρ 2 ρ 1 β S * α 1 + 1 + 1 ρ 1 β + ρ 2 α 2 β I * N α ω j 12 = 0 j 13 = 1 ρ 2 ρ 1 β S * α 1 + 1 + 1 ρ 1 β + ρ 2 α 2 β S * N j 21 = α j 22 = ω j 23 = 0 j 31 = 1 ρ 2 ρ 1 β α 1 S * I * S * α 1 + 1 2 N + 1 ρ 2 ρ 1 β S * α 1 + 1 + 1 ρ 1 β + ρ 2 α 2 β I * N j 32 = 0 j 33 = 1 ρ 2 ρ 1 β S * α 1 + 1 + 1 ρ 1 β + ρ 2 α 2 β S * N γ ω
The first eigenvalue from the Jacobian matrix at E * is ω , which is negative. In order to further analyze this case, we computed the trace of the Jacobian matrix at E * and assessed the possible scenarios for ρ 1 and ρ 2 . To obtain a negative value in the trace, a condition relating I and S was determined:
I * > S * ( α 1 + 1 ) ( ( α + γ + 3 ω ) ( S * α 1 + 1 ) N S * β ) β
Subsequently, the determinant of the Jacobian matrix evaluated at E * , and there are two cases depending on the possible scenarios of ρ 1 and ρ 2 , as can be seen in the following expressions:
det J E 1 * = N ω 2 + ( ( α + γ ) N β α 2 ( S * I * ) ω + N α γ β α 2 I * γ + S * α ) ω N
A = N S * α 1 + 1 2 ω 2 B = S * α 1 + 1 2 α + γ N β S * 2 α 1 I * + S * ω C = α γ S * α 1 + 1 2 N D = β S * 2 α α 1 I * γ + S * α det J E 2 * = ( A + B + C D ) ω N S * α 1 + 1 2
Given the negative sign associated with the expressions for the determinant, we cannot ensure that the determinant will be positive, therefore it is established that the system defined in Equation (20), at E * is a saddle point.    □

2.5. ODE Model Contact Rate ( β ) Estimation with LSTM Neural Networks

The data for the SVIR compartmental model were collected from reliable sources, including the National Center of Epidemiology of Spain (CNE, Centro Nacional de Epidemiología), the Ministry of Health of Spain (Ministerio de Sanidad), the Ministry of Transport, Mobility and Urban Agenda (Ministerio de Transportes, Movilidad y Agenda Urbana), and the National Institute of Statistics of Spain (INE, Instituto Nacional de Estadística), as specified in Table 1. To address shortcomings in data availability, assumptions were made to simplify the model, such as the required granularity or a complete lack of information. This work is focused on modeling the SARS-CoV-2 epidemic for the Autonomous Communities of Spain, which have undergone a specific set of restrictions.
To handle irregularities in the data and improve the training process for artificial neural networks, data processing techniques were implemented. These techniques are commonly used in data-driven applications. Automated data collection (where applicable), data processing, artificial neural networks, and numerical simulations were coded in Python.
To calculate the system dynamics, we assumed that the epidemics would follow the established system from Equation (1), resulting in a system with time-varying parameters instead of constant values. The model dynamics were then calculated with available data, and these dynamics were used to train an LSTM neural network for estimating the contact rate ( β ). LSTM networks were chosen for their ability to handle sequential data, including time series, and their capacity to preserve trends in such data, which is desirable when estimating parameters for dynamical systems that describe the evolution of an epidemic in time. Additionally, LSTM networks can handle multi-dimensional input, allowing for consideration of multiple variables, such as mobility and other underlying contributing factors. The input for the LSTM network in this work is a multi-dimensional array with terms limited to the dynamics established by the system and the contributing factors. We implemented the LSTM neural network in Python using the torch library and set the following hyperparameters: learning rate of 0.1, 1 layer, 1000 training epochs, a sequence length of 28 data points, 80% of the data used for training, and 20% for testing. We trained a separate network for each Autonomous Community to estimate both the contact ( β ) and recovery ( γ ) rates. The input matrix for the LSTM networks consists of 28 time observations for the variables I, d I / d t , mobility percentage, d V / d t , and V, with the output variable being d R / d t . This output is subsequently used to calculate the parameters γ and β . An example of the input and output of the LSTM network for Madrid is displayed in Figure 1 and Figure 2, respectively.
Figure 2 shows the values of d R / d t used to train and test the LSTM network with data for the Community of Madrid and the output value of the network after training.
To train the LSTM network for each Autonomous Community, we used the data specific to that region, which resulted in different errors due to variations in the input data. For each community, we calculated three different errors: one for the training dataset, one for the test dataset, and a consolidated error combining both datasets. The root mean square error (RMSE) for all of the Autonomous Communities is presented in the heatmap shown in Figure 3. We observed that the errors were highest in the training dataset and lowest in the test dataset, which indicates that the LSTM is capable of generalizing well to unseen data. Specific error values for each community are included in Table A1 in Appendix A.
The trained LSTM networks output the time series for the parameter β in the system described by Equation (1). To obtain this parameter, the value of d R / d t from the LSTM is substituted into the system. However, since the LSTM output has an associated estimation error, as illustrated in Figure 3 and detailed in Table A1 in Appendix A, the calculated parameters are also affected by these errors.

2.6. ODE Model Numerical Simulations

This study employed numerical simulations using a fourth-order Runge–Kutta method with LSTM-estimated parameters to solve initial value problems (IVPs) for ordinary differential equation (ODE) systems. IVPs involve finding the solution of a differential equation at a specific initial time, given the initial conditions. The fourth-order Runge–Kutta method is widely used in solving IVPs of the form given in Equation (28), which is a system with time-varying parameters.
( I V P ) d U ( t ) d t = f ( t , U ( t ) ) U ( 0 ) = U ( 0 ) ,
The fourth-order Runge–Kutta method is an explicit multi-stage one-step method that requires only information from the previous time-step and is efficient for ODE systems. Here, U is the vector of unknowns, f is the right-hand side of Equation (2), and U ( 0 ) is the initial condition. Predictor-corrector methods are also efficient solvers for initial value problems but were not used in this study.
The fourth-order Runge–Kutta method is expressed as a set of equations, which include four intermediate values ( K 1 , K 2 , K 3 , K 4 ) and the new state vector, U n + 1 , at the next time step. The method is shown below:
K 1 = Δ t f ( t n , U n ) K 2 = Δ t f ( t n + Δ t 2 , U n + K 1 2 ) K 3 = Δ t f ( t n + Δ t 2 , U n + K 2 2 ) K 4 = Δ t f ( t n + Δ t , U n + K 3 ) U n + 1 = U n + 1 6 K 1 + 2 ( K 2 + K 3 ) + K 4 ,
Here, t n + 1 is the time for which the solution is computed, and t n is the previous time. U n is the vector of unknowns for time t n , U n + 1 is the vector of unknowns for time t n + 1 , and Δ t is the size of the time step considered.
The study conducted different simulations using the ODE models described in Section 2.1. Simulation A followed the dynamic from Equation (2), with specific values of α 1 and α 2 that fit best to each Autonomous Community. Simulation B followed the same dynamic, but with a different set of parameters for each Autonomous Community, as in Simulation A. Simulation C followed the same dynamic from Equation (2), but with the same set of parameters for all Autonomous Communities. Simulation D maintained the dynamic from Equation (1) but involved outlier handling for the LSTM-estimated parameters, removing them and imputing values through linear interpolation. Simulation E established the baseline, following the dynamic from Equation (1).
To compare the simulation results, the Root Mean Square Error (RMSE) was used as the error measure, computed as the difference between the actual data and the simulation results. The RMSE values for the different simulations and Autonomous Communities are displayed in Figure 4, and detailed results can be found in Table A3 in Appendix C. The largest differences in RMSE values were observed between Simulation E and the other simulations, except for the Autonomous City of Ceuta and the Autonomous Community of Canarias. Simulation C was largely ineffective, with gaps in the results, mainly in three Autonomous Communities (Andalucía, Aragón, and Extremadura), and in some cases, also for Simulations A and B. These results are discussed further in Section 4.

3. PDE SVIR Compartmental Model with Diffusion

This section extends the previous compartmental models to incorporate Diffusion, which allows the study of epidemic or pandemic transmission across physical space. The proposed system of one-dimensional PDEs is based on the modified SVIR Compartmental Model, and is expressed in Equations (30) and (31).
To solve the proposed PDE system, a finite volume framework is used with fifth-order WENO reconstruction in space, allowing for accurate and stable representations of sharp gradients in the solution. The time integration is performed using an RK3-TVD scheme, which is a third-order, total variation diminishing scheme that is highly effective in handling numerical oscillations and preserving the monotonicity of the solution. The accuracy of the numerical scheme is validated through the Manufactured Solutions methodology, yielding accurate numerical results compared to the exact solution. However, as the model is not based on known diffusion coefficients, the simulations cannot be directly compared to real data. While the current study uses arbitrary values for the diffusion parameters, the proposed model can be further studied and modified to explore more realistic scenarios.

3.1. PDE Model Formulation

To incorporate spatial variation in one dimension, we introduce Diffusion in the SVIR model described by Equation (2), which considers the case where ρ 2 = 0 for any t. The resulting PDE system is as follows:
S ( x , t ) t = d S x S S x ρ 1 β 1 + α 1 S + ( 1 ρ 1 ) β S I N α S V ( x , t ) t = d V x V V x + α S I ( x , t ) t = d I x I I x + ρ 1 β 1 + α 1 S + ( 1 ρ 1 ) β S I N γ I R ( x , t ) t = d R x V R x + γ I ,
In Equation (30), the diffusion constants d S , d V , d I , and d R are introduced to consider diffusion–dispersion mechanisms in the model. This inclusion of diffusion is important because it provides the model with an element that may increase its capacity to represent a more realistic situation, especially in determining the spatial propagation of the different compartments due to population mobility in one dimension. Related works, such as [29,30], have proposed incorporating diffusion in epidemic models to improve their realism and applicability. For instance, [30] proposes a reaction–diffusion SVIR epidemic model that considers spatial heterogeneity, while [29] introduced diffusion into a SIR model.
Similarly, when introducing Diffusion into the ODE version with both mechanisms ρ 1 and ρ 2 , i.e., Equation (2), the resulting model is:
S ( x , t ) t = d S x S S x κ ( S ) S I N α S V ( x , t ) t = d V x V V x + α S I ( x , t ) t = d I x I I x + κ ( S ) S I N γ I R ( x , t ) t = d R x R R x + γ I ,
where κ ( S ) = ( 1 ρ 2 ) ρ 1 β 1 + α 1 S + ( 1 ρ 1 ) β + ρ 2 α 2 β .
To further enhance the proposed diffusion model, the recovery and vaccination mechanisms have been modified. It is assumed that recovered individuals become susceptible again after an established constant rate, and vaccinated individuals also become susceptible again due to waning immunity at a fixed rate based on averages from the literature. By incorporating these factors into the model, a more realistic representation of the transmission dynamics can be achieved.
S ( x , t ) t = d S x S S x κ ( S ) S I N α S + η R R + η V V V ( x , t ) t = d V x V V x + α S η V V I ( x , t ) t = d I x I I x + κ ( S ) S I N γ I R ( x , t ) t = d R x R R x + γ I η R R ,
For Equations (30)–(32), the control mechanisms are similar to those of the ODE version of the modified model. The parameter β is only time-dependent; therefore, ρ 1 also follows Equation (4). Regarding ρ 2 , an average number of infected individuals I ¯ t is calculated for the number of cells n x at each time step. Moreover, an average of the last m time steps is considered, where a value A is calculated (as shown in Algorithm 1) to compare against the threshold δ 2 . This comparison establishes whether the control parameter ρ 2 will be active with a value of ρ 2 = 1 or inactive as ρ 2 = 0 . The algorithm is:
Algorithm 1: Steps to calculate ρ 2
At each time step t
I t ¯ 1 n x i = 1 n x I i
Considering j as the current time step
A 1 m t = j m j I t ¯ I ¯ t 1 t t t t 1
if   A δ 2 > 0   then
     ρ 2 1
else
     ρ 2 0
end if

3.2. PDE Model Numerical Simulations

To study the diffusion of the epidemic, we use a partial differential equation (PDE) approach, which can be represented more compactly as follows:
U t = F x + Q ,
where the vector of unknowns is U = ( S , V , I , R ) T . The vector of fluxes, F , is given by
F = d S S S x d V V V x d I I I x d R R R x .
and the source term, Q , is given by
Q = κ ( S ) S I N α S α S κ ( S ) S I N γ I γ I .
where d S , d V , d I , and d R are the diffusion coefficients for each compartment. The diffusion terms in Equation (33) account for the spatial variation of the epidemic dynamics, which can be important in modeling the spread of the disease across regions.
To discretize the system in (33) in space, we utilize the finite volume method, where we divide the spatial domain into discrete control volumes. Denoting the i t h control volume as T i = [ x i 1 / 2 , x i + 1 / 2 ] , we integrate the system over T i using cell averages, resulting in the following equation:
U i t = 1 Δ x i ( F i + 1 / 2 F i 1 / 2 ) + Q i ,
where the cell average of the solution is
U i = 1 Δ x i x i 1 / 2 x i + 1 / 2 U ( x , t ) d x ,
and the cell average of the source term reads
Q i = 1 Δ x i x i 1 / 2 x i + 1 / 2 Q ( x , t ) d x .
To obtain values and derivatives where needed, we use the Weighted Essentially Non Oscillatory (WENO) reconstruction method based on cell averages ([31,32,33]). The WENO method reconstructs the values based on r candidate stencils for reconstruction, each of which consists of r control volumes and a polynomial of degree r 1 . The reconstructed values are obtained as a convex combination of the values of the r possible polynomials, p r ( x ) , weighted with nonlinear weights.
The nonlinear weights are defined as
ω k = α k j = 0 r 1 α j , ( k = 0 , , r 1 )
where α k = d k ϵ + β k , ϵ is a small value introduced to avoid division by zero, and β k are smoothness indicators. In this work, we take ϵ = 10 6 . The smoothness indicators are computed according to the formulas
β k = m = 0 r 1 x i 1 / 2 x i + 1 / 2 d m d x m p k ( x ) 2 Δ x i 2 m 1 d x , ( k = 0 , , r 1 ) ,
For r = 3 , which corresponds to second degree polynomials, we have β 0 , β 1 , and β 2 given by the expressions shown in the paragraph above.
β 0 = 13 12 U i 2 U i + 1 + U i + 2 2 + 1 4 3 U i 4 U i + 1 + U i + 2 2 β 1 = 13 12 U i 1 2 U i + U i + 1 2 + 1 4 3 U i 1 + U i + 1 2 β 2 = 13 12 U i 2 2 U i 1 + U i 2 + 1 4 3 U i 2 4 U i 1 + 3 U i 2 ,
After discretizing the system in space using the finite volume method, we obtain a system of ODEs given by (42).
d d t U i = L i U , d d x U ,
To solve this system, we use a third-order Runge–Kutta TVD numerical scheme [34], which is designed to maintain the positivity of the solutions. The system takes the form of d d t U i = L i ( U , d d x U ) , where L i ( U , d d x U ) = 1 Δ x i ( F i + 1 / 2 F i 1 / 2 ) + Q i .
In the third-order Runge–Kutta TVD numerical scheme, U i n represents the numerical solution at time t n , and U i n + 1 represents the numerical solution at time t n + 1 = t n + Δ t . In this scheme, the intermediate states U i k , 1 and U i k , 2 are computed, followed by the final update step to obtain U i n + 1 . The third-order Runge–Kutta TVD numerical scheme is well suited for solving hyperbolic PDEs, as it maintains stability and accuracy even in the presence of discontinuities in the solution.
The manufactured solutions method is a validation technique used to test the accuracy of numerical schemes. In this method, an exact solution is artificially created by specifying a set of functions that satisfy the governing equations. For this study, we use the manufactured solutions method to verify the accuracy of the numerical scheme for the diffusion case. We consider Equation (43), where the small value of the derivatives can be considered zero for the numerical scheme, and the solution is considered to have Neumann boundary conditions. This system of equations describes a diffusion process for the S, V, I, and R variables.
S = 4 e x 2 t V = e 2 x 2 t I = 1.5 e 4 x 2 t R = 2.5 e 0.8 x 2 t
Figure 5 visually compares the exact solution obtained using Equation (43) and the numerical result, which serves as a validation of the accuracy of the numerical method.
The resulting order of accuracy for the numerical scheme is 3, as shown in Table 2. This was achieved with 640 cells. The lower order seen in Cells 160 can be disregarded as it trends towards three after the increase to 320 cells, which is expected since the scheme includes time integration with the third-order Runge–Kutta TVD method. Moreover, the order trending upwards of 3 is also expected because the space integration is performed via a 1-D WENO reconstruction with a higher order. However, the lower order of the time integration is causing the order reduction.
Incorporating the PDE model proposed in this study, Equation (31), with the model from [30] yields the following combined model:
S ( x , t ) t = d S x S S x κ ( S ) S I N α S + ϕ ω S V ( x , t ) t = d V x V V x + α S β 2 V I γ 2 V ω V I ( x , t ) t = d I x I I x + κ ( S ) S I N + β 2 V I γ I ω I R ( x , t ) t = d R x R R x + γ I + γ 2 V ω R ,
We simulate the PDE model using the parameter values ϕ = 0.392465 , β 2 = 0.0016 , γ 2 = 0.001832 , and ω = 3.9 × 10 5 , along with a time-dependent β estimated from data for the Autonomous Community of Madrid using the LSTM methodology described in Section 2.5. The initial values are set to ( S ( 0 ) , V ( 0 ) , I ( 0 ) , R ( 0 ) ) = ( 60 , 50 , 30 , 0 ) . Since no real data are available for comparison, we compare the behavior of this simulation to that observed in the reference paper [30]. The waves observed in our simulation are not uniformly distributed (Figure 6), unlike the comparison example, and exhibit more realistic behavior that could be attributed to changes in behavior or policies throughout the epidemic evolution. These variations are related to the estimated β and the control parameters underlying the mechanism defined by κ .

4. Results and Discussion

In this section, we discuss the performance of the proposed model for simulating the spread of COVID-19 in Spain. The model was trained using LSTM neural networks to estimate the parameters of the SIR epidemiological model, and it was validated against real data. We start by presenting the results of the baseline simulation and then discuss the effectiveness of the proposed outlier handling techniques.
To improve the performance of the baseline simulation, which is referred to as simulation E and is based on the base model from Equation (1), we utilized linear interpolation as a technique for outlier handling. This approach aimed to correct the under and overestimations introduced by the LSTM-estimated parameters. While this process improved the RSME values, significant miss-estimations remained, as observed in Figure 7 for the case of Madrid. Specifically, this figure shows that the infected ratio was underestimated in the first waves and vastly overestimated in the last wave included in the study.
Figure 4 displays the RMSE values for all Autonomous Communities for simulation E, highlighting the significant deviation of the simulation from the actual values. Thus, a pattern of mis-estimations emerges, and the need for further improvements is evident.
The neural network estimated parameters can yield simulations where the evolution of infected individuals is over or underestimated, depending on the input data and other outside factors not included in the input data used to train the neural network. To correct these potential behaviors observed with estimated parameters, a modified version of the SVIR compartmental model is introduced. This model is defined by Equation (2), which includes explicit terms controlled by established thresholds for the different Autonomous Communities of Spain.
Using the same parameters for the model described by Equation (2), the simulation is ineffective for nine communities, and poor results are obtained for four of the communities where the simulation can be carried out. These results indicate significant differences in behavior between the different communities, pointing to the need to determine specific values more adequate to the particular series/community. When utilizing specific parameter values for each Autonomous Community, the results obtained by the simulation are much closer to actual infected values, with lower RSME values, for the majority of Autonomous Communities. However, there were a few exceptions involving three communities for which Simulations A, B, and C were ineffective. The improvement is visually observable in the case of Madrid, shown in Figure 7, particularly in the simulated period’s last wave.
Assessing the specific values of RSME and visual analysis of the simulations, the following can be stated: without modifying the SVIR system from Equation (1) or the LSTM estimated parameters, there are significant errors in the numerical simulation compared to actual data. As a general solution, post-processing LSTM estimated parameters with well-established techniques can improve the simulation’s result, as observed by comparing Simulation D in Figure 4 vs. Simulation E. However, outlier handling is still insufficient, leading to the proposed modified SVIR model, for which not only the RSME improved but the deviation of the behavior in the waves from the simulation to actual data is reduced, highlighting the relevance of the proposed modified SVIR model.
We compared the relative L 2 error for 16 Autonomous communities (or autonomous cities in the case of Ceuta and Melilla) with a previous study by Bousquet et al. [27], who used Deep Learning to forecast values of the β parameter. Their study reported relative L 2 errors ranging from 3.13 × 10 3 to 6.29 × 10 2 . In contrast, our study achieved an average relative L 2 error of 1.10 × 10 3 , with a range from a minimum of 8.91 × 10 5 to a maximum of 3.92 × 10 3 . Table A4 lists the computed errors. Our findings suggest that our predictions for the β parameter in COVID-19 are more accurate and precise. Although the two studies used different compartmental models, the β parameter was the contact parameter in both. We did not compare the results for other parameters reported by Bousquet et al. as our study focused solely on the β parameter.
The PDE system from Equations (30) and (31), based on the modified SVIR Compartmental Model, can be used to study scenarios including Diffusion with the bi-parameter ( ρ 1 , ρ 2 ) control system. To validate the simulation mechanism, the Manufactured Solutions methodology was implemented, obtaining accurate numerical results compared to the exact solution. Although this is not a completely realistic system and does not have specified-known diffusion coefficients, reasonable constants were selected and applied to carry out the Manufactured Solutions methodology. These coefficients should be further studied and modified if the intention is to explore more realistic situations. In this study, the system was proposed as an extension of the ODE system to increase the potential to resemble a more realistic situation. So far, the study has validated the model and its corresponding order. Further studies could focus on comparing the model’s simulation when changing the currently used Diffusion coefficients.

5. Conclusions

This study presents a detailed simulation of COVID-19 in the Autonomous Communities of Spain using LSTM-estimated parameters and the modified SVIR model. The model includes regulating mechanisms based on the behavior of previous time steps. We have established the positivity of the ODE version of the model and explored its local stability. However, further analysis is needed to confirm the consistency and global stability of the modified model.
The proposed mechanisms show significant improvements in simulation accuracy, as evidenced by lower RSME and closer waves to actual data. Although LSTM proved to be efficient for estimating parameters for the compartmental model. The proposed modified SVIR model and LSTM neural networks approach is a promising tool for forecasting the spread of COVID-19, improving the accuracy of predictions for key parameters such as β , with a relatively low relative L 2 error, as shown by the results of the comparison with a previous study.
The proposed modified SVIR model has demonstrated its value as a tool for studying epidemic behavior in a short period, and its potential can be extended to other countries or more granular levels in Spain. However, the SVIR and proposed modified SVIR models have limitations, particularly with regard to their assumptions, such as granting permanent immunity to vaccinated and recovered individuals, which is not the case. To improve the accuracy and applicability of these models, future studies should explore strategies to incorporate waning immunity and recurrent vaccination, which could be crucial in accurately simulating the long-term impact of the COVID-19 pandemic.
This work presents an extension of the SVIR compartmental model in the form of a one-dimensional PDE system with diffusion. Although the stability analysis of the PDE system has not been carried out, the numerical scheme for the system is validated using the manufactured solutions methodology. The PDE system is solved using a fifth-order WENO reconstruction in space and an RK3-TVD scheme for time integration. While it is not possible to compare the results of the PDE system with actual data, this extension of the SVIR model has the potential to provide additional insights into the behavior of epidemics. Further studies are needed to fully explore the potential of this extension of the SVIR model.
In summary, this work highlights the relevance of developing and improving models to provide more accurate and comprehensive insights into epidemic behavior, and demonstrates the value of incorporating advanced numerical techniques such as LSTM and PDEs with diffusion into these models. Further studies can build upon the proposed methodology to enhance our understanding of epidemics and support decision-making in the public health sector.

Author Contributions

Conceptualization, A.W. and A.H.; methodology, A.W. and A.H.; software, A.W. and A.H.; validation, A.W. and A.H.; formal analysis, A.W.; investigation, A.W.; data curation, A.W.; writing—original draft preparation, A.W.; writing—review and editing, A.H.; visualization, A.W.; supervision, A.H. All authors have read and agreed to the published version of the manuscript.

Funding

The research of the second author is partially supported by the project PID2020-112517GB-I00 of Ministerio de Ciencia e Innovación (Spain).

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:
ANNArtificial Neural Networks
LSTMLong Short-Term Memory
RMSERoot Mean Square Error

Appendix A. LSTM Errors

Table A1 details the errors, RSME, obtained for each Autonomous Community, for the Training and Test datasets, as well as the Overall dataset, which includes both training and test.
Table A1. LSTM errors (RSME) per Autonomous Community.
Table A1. LSTM errors (RSME) per Autonomous Community.
Autonomous CommunityOverallTrainingTest
Andalucía 7.7988 × 10 6 6.3157 × 10 6 1.4880 × 10 5
Aragón 1.1609 × 10 5 8.0570 × 10 6 2.5895 × 10 5
Asturias 6.9153 × 10 6 7.2418 × 10 6 3.5171 × 10 6
Baleares 3.0618 × 10 5 1.0651 × 10 5 8.5197 × 10 5
C. Valenciana 2.3124 × 10 5 9.9428 × 10 6 6.2305 × 10 5
Canarias 2.6032 × 10 6 2.2519 × 10 6 4.4566 × 10 6
Cantabria 6.6877 × 10 6 6.3992 × 10 6 8.5846 × 10 6
Castilla La Mancha 4.6016 × 10 5 2.4329 × 10 5 1.1756 × 10 4
Castilla y Leon 9.8641 × 10 6 9.1941 × 10 6 1.3974 × 10 5
Cataluña 1.9340 × 10 5 7.4910 × 10 6 5.3035 × 10 5
Ceuta 1.9893 × 10 5 1.4680 × 10 5 4.2168 × 10 5
Extremadura 4.0817 × 10 5 9.5022 × 10 6 1.1728 × 10 4
Galicia 5.0179 × 10 6 4.9919 × 10 6 5.2295 × 10 6
La Rioja 1.4697 × 10 5 1.0096 × 10 5 3.3032 × 10 5
Madrid 1.8866 × 10 5 8.8077 × 10 6 4.9913 × 10 5
Melilla 2.0307 × 10 5 2.0553 × 10 5 1.8314 × 10 5
Murcia 1.0300 × 10 4 6.6730 × 10 6 3.0274 × 10 4
Navarra 1.2086 × 10 5 1.2723 × 10 5 4.9761 × 10 6
País Vasco 8.5856 × 10 6 8.3311 × 10 6 1.0332 × 10 5

Appendix B. Parameters for the Modified ODE SVIR Compartmental Model

Parameters were determined for the modified SVIR model. As can be seen in Table A2, most parameters are equal for most of the Autonomous Communities. In the case of a, found in Equation (5) part of the mechanism of the model from Equation (2), the parameter ranges from 2 to 9; these are related to the variations of the estimated β parameters which have different behaviours for the communities. One of the parameters remains constant for all simulations and all Autonomous Communities, α 1 = 10 .
Table A2. Parameters for the proposed modified SVIR model per Autonomous Community.
Table A2. Parameters for the proposed modified SVIR model per Autonomous Community.
Autonomous Communitya δ 2 α 2 i
AndalucíaN/AN/AN/AN/A
AragónN/AN/AN/AN/A
Asturias315001.53
Baleares615001.56
C. Valenciana315001.53
Canarias615001.59
Cantabria315001.53
Castilla La Mancha315001.53
Castilla y Leon615001.53
Cataluña215001.53
Ceuta315001.013
ExtremaduraN/AN/AN/AN/A
Galicia315001.53
La Rioja915001.59
Madrid63001.13
Melilla315001.53
Murcia315001.56
Navarra315001.53
País Vasco915001.53
Default 1315001.013
1 General parameter value used if the Autonomous Community, Simulation C.

Appendix C. ODE Numerical Simulation Errors

Table A3 details the errors, RSME, obtained for each Autonomous Community when performing numerical simulations. As mentioned in Section 2.6, some results are not available due to the nature of the data, and more adequate parameters were not determined during this study.
Table A3. Numerical simulation errors (RSME) per Autonomous Community.
Table A3. Numerical simulation errors (RSME) per Autonomous Community.
Autonomous CommunityABCDE
AndalucíaN/AN/AN/A 2.4203 × 10 4 1.3161 × 10 2
AragónN/AN/AN/A 7.4293 × 10 4 2.2352 × 10 2
Asturias 8.7013 × 10 4 8.7013 × 10 4 6.7207 × 10 3 6.9233 × 10 4 2.1149 × 10 2
Baleares 1.1484 × 10 3 1.1484 × 10 3 1.0802 × 10 2 2.5237 × 10 3 6.1029 × 10 3
C. Valenciana 1.7533 × 10 4 1.7533 × 10 4 N/A 1.2169 × 10 4 7.9763 × 10 3
Canarias 1.7698 × 10 4 1.7698 × 10 4 6.2609 × 10 4 2.9668 × 10 3 2.5351 × 10 3
Cantabria 2.5861 × 10 4 2.5861 × 10 4 3.1585 × 10 3 2.8453 × 10 4 1.6875 × 10 2
Castilla La Mancha 1.2475 × 10 3 1.2475 × 10 3 N/A 1.2316 × 10 3 4.5537 × 10 2
Castilla y Leon 5.7199 × 10 4 5.7199 × 10 4 1.7050 × 10 2 1.1459 × 10 3 1.8104 × 10 2
Cataluña 1.3750 × 10 3 1.3750 × 10 3 N/A 2.0126 × 10 4 9.1348 × 10 3
Ceuta 2.2848 × 10 5 2.2848 × 10 5 6.1686 × 10 4 7.0457 × 10 5 1.1597 × 10 3
ExtremaduraN/AN/AN/A 8.8996 × 10 4 9.2140 × 10 3
Galicia 7.9528 × 10 5 7.9528 × 10 5 8.2338 × 10 5 2.2975 × 10 3 8.1395 × 10 3
La Rioja 1.0924 × 10 4 1.0924 × 10 4 1.7713 × 10 3 1.1966 × 10 3 3.4083 × 10 2
Madrid 5.4030 × 10 4 5.4030 × 10 4 3.2903 × 10 4 8.9100 × 10 4 1.2344 × 10 2
Melilla 1.9580 × 10 4 1.9580 × 10 4 3.2850 × 10 4 3.2719 × 10 3 2.1193 × 10 2
Murcia 4.6199 × 10 4 4.6199 × 10 4 N/A 1.7584 × 10 4 1.7155 × 10 2
Navarra 2.0010 × 10 4 2.0010 × 10 4 N/A 5.5822 × 10 4 2.1931 × 10 2
País Vasco 1.2436 × 10 3 1.2436 × 10 3 N/A 4.1042 × 10 4 4.5271 × 10 2
Table A4. Relative L 2 error for β .
Table A4. Relative L 2 error for β .
Autonomous CommunityRelative L 2 Error
AndalucíaN/A
AragónN/A
Asturias 3.9174 × 10 3
Baleares 1.0725 × 10 3
C. Valenciana 6.3781 × 10 4
Canarias 3.1516 × 10 4
Cantabria 5.4328 × 10 4
Castilla La Mancha 1.6311 × 10 3
Castilla y Leon 7.2528 × 10 4
Cataluña 1.0131 × 10 3
Ceuta 5.8867 × 10 4
ExtremaduraN/A
Galicia 1.8901 × 10 4
La Rioja 8.9091 × 10 5
Madrid 5.0681 × 10 4
Melilla 3.3019 × 10 3
Murcia 5.1324 × 10 4
Navarra 2.5385 × 10 4
País Vasco 1.7748 × 10 3

References

  1. World Health Organization. Listings of WHO’s Response to COVID-19; World Health Organization: Geneva, Switzerland, 2020. [Google Scholar]
  2. Ministerio de Sanidad. Enfermedad por Nuevo Coronavirus, COVID-19 Situación Actual; Paseo del Prado: Madrid, Spain, 2022; Available online: https://www.sanidad.gob.es/profesionales/saludPublica/ccayes/alertasActual/nCov/situacionActual.htm (accessed on 24 October 2022).
  3. Hethcote, H.W.; Tudor, D.W. Integral equation models for endemic infectious diseases. J. Math. Biol. 1980, 9, 37–47. [Google Scholar] [CrossRef]
  4. Alexander, M.E.; Bowman, C.; Moghadas, S.M.; Summers, R.; Gumel, A.B.; Sahai, B.M. A Vaccination Model for Transmission Dynamics of Influenza. SIAM J. Appl. Dyn. Syst. 2004, 3, 503–524. [Google Scholar] [CrossRef]
  5. Meng, X.; Chen, L. The dynamics of a new SIR epidemic model concerning pulse vaccination strategy. Appl. Math. Comput. 2008, 197, 582–597. [Google Scholar] [CrossRef]
  6. Arino, J.; McCluskey, C.C.; van den Driessche, P. Global Results for an Epidemic Model with Vaccination That Exhibits Backward Bifurcation. SIAM J. Appl. Math. 2003, 64, 260–276. [Google Scholar] [CrossRef] [Green Version]
  7. Hethcote, H.W. An immunization model for a heterogeneous population. Theor. Popul. Biol. 1978, 14, 338–349. [Google Scholar] [CrossRef] [PubMed]
  8. Kribs-Zaleta, C.M.; Velasco-Hernández, J.X. A simple vaccination model with multiple endemic states. Math. Biosci. 2000, 164, 183–201. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  9. Kaymakamzade, B.; Hincal, E. Delay epidemic model with and without vaccine. AIP Conf. Proc. 2018, 1997, 020025. [Google Scholar] [CrossRef]
  10. Safi, M.A.; Garba, S.M. Global Stability Analysis of SEIR Model with Holling Type II Incidence Function. Comput. Math. Methods Med. 2012, 2012, 826052–826058. [Google Scholar] [CrossRef] [Green Version]
  11. Kaddar, A. Stability analysis in a delayed SIR epidemic model with a saturated incidence rate. Nonlinear Anal. 2010, 15, 299–306. [Google Scholar] [CrossRef]
  12. Hou, J.; Teng, Z. Continuous and impulsive vaccination of SEIR epidemic models with saturation incidence rates. Math. Comput. Simul. 2009, 79, 3038–3054. [Google Scholar] [CrossRef]
  13. Anderson, R.M.; May, R.M. Regulation and Stability of Host-Parasite Population Interactions: I. Regulatory Processes. J. Anim. Ecol. 1978, 47, 219–247. [Google Scholar] [CrossRef]
  14. Brauer, F. Backward bifurcations in simple vaccination models. J. Math. Anal. Appl. 2004, 298, 418–431. [Google Scholar] [CrossRef] [Green Version]
  15. Nistal, R.; De la Sen, M.; Alonso-Quesada, S.; Ibeas, A. On a New Discrete SEIADR Model with Mixed Controls: Study of Its Properties. Mathematics 2019, 7, 18. [Google Scholar] [CrossRef] [Green Version]
  16. López, L.; Rodó, X. A modified SEIR model to predict the COVID-19 outbreak in Spain and Italy: Simulating control scenarios and multi-scale epidemics. Results Phys. 2021, 21, 103746. [Google Scholar] [CrossRef] [PubMed]
  17. Marinov, T.T.; Marinova, R.S. Inverse problem for adaptive SIR model: Application to COVID-19 in Latin America. Infect. Dis. Model. 2022, 7, 134–148. [Google Scholar] [CrossRef] [PubMed]
  18. Kopfová, J.; Nábělková, P.; Rachinskii, D.; Rouf, S.C. Dynamics of SIR model with vaccination and heterogeneous behavioral response of individuals modeled by the Preisach operator. J. Math. Biol. 2021, 83. [Google Scholar] [CrossRef] [PubMed]
  19. Aràndiga, F.; Baeza, A.; Cordero-Carrión, I.; Donat, R.; Martí, M.C.; Mulet, P.; Yáñez, D.F. A Spatial-Temporal Model for the Evolution of the COVID-19 Pandemic in Spain Including Mobility. Mathematics 2020, 8, 1677. [Google Scholar] [CrossRef]
  20. Bertaglia, G.; Pareschi, L. Hyperbolic compartmental models for epidemic spread on networks with uncertain data: Application to the emergence of Covid-19 in Italy. Math. Model. Methods Appl. Sci. 2021, 31, 2495–2531. [Google Scholar] [CrossRef]
  21. Ramos, A.; Vela-Pérez, M.; Ferrández, M.; Kubik, A.; Ivorra, B. Modeling the impact of SARS-CoV-2 variants and vaccines on the spread of COVID-19. Commun. Nonlinear Sci. Numer. Simul. 2021, 102, 105937. [Google Scholar] [CrossRef] [PubMed]
  22. Bertaglia, G.; Boscheri, W.; Dimarco, G.; Pareschi, L. Spatial spread of COVID-19 outbreak in Italy using multiscale kinetic transport equations with uncertainty. Math. Biosci. Eng. 2021, 18, 7028–7059. [Google Scholar] [CrossRef]
  23. Antoniou George, E.; Mentzelopoulou, S. Neural Networks: An Application to the Epidemics. In Proceedings of the Neural, Parallel and Scientific Computations, Atlanta, GA, USA, 28–31 May 1995; Volume 1, pp. 18–21. [Google Scholar]
  24. Rao, V.S.H.; Kumar, M.N. Estimation of the parameters of an infectious disease model using neural networks. Nonlinear Anal. Real World Appl. 2010, 11, 1810–1818. [Google Scholar] [CrossRef] [Green Version]
  25. Atencia, M.; García-Garaluz, E.; Arazoza, H.D.; Joya, G. Estimation of parameters based on artificial neural networks and threshold of HIV/AIDS epidemic system in Cuba. Math. Comput. Model. 2013, 57, 2971–2983. [Google Scholar] [CrossRef]
  26. Tessmer, H.L.; Ito, K.; Omori, R. Can Machines Learn Respiratory Virus Epidemiology?: A Comparative Study of Likelihood-Free Methods for the Estimation of Epidemiological Dynamics. Front. Microbiol. 2018, 9, 343. [Google Scholar] [CrossRef]
  27. Bousquet, A.; Conrad, W.H.; Sadat, S.O.; Vardanyan, N.; Hong, Y. Deep learning forecasting using time-varying parameters of the SIRD model for COVID-19. Sci. Rep. 2022, 12, 3030. [Google Scholar] [CrossRef] [PubMed]
  28. Huo, H.F.; Zou, M.X. Modelling effects of treatment at home on tuberculosis transmission dynamics. Appl. Math. Model. 2016, 40, 9474–9484. [Google Scholar] [CrossRef]
  29. Gai, C.; Iron, D.; Kolokolnikov, T. Localized outbreaks in an S-I-R model with diffusion. J. Math. Biol. 2020, 80, 1389–1411. [Google Scholar] [CrossRef]
  30. Zhang, C.; Gao, J.; Sun, H.; Wang, J. Dynamics of a reaction–diffusion SVIR model in a spatial heterogeneous environment. Phys. A 2019, 533, 122049. [Google Scholar] [CrossRef]
  31. Jiang, G.S.; Shu, C.W. Efficient implementation of weighted ENO schemes. J. Comput. Phys. 1996, 126, 202–228. [Google Scholar] [CrossRef] [Green Version]
  32. Liu, X.D.; Osher, S.; Chan, T. Weighted Essentially Non-oscillatory Schemes. J. Comput. Phys. 1994, 115, 200–212. [Google Scholar] [CrossRef] [Green Version]
  33. Titarev, V.; Toro, E. Finite-volume WENO schemes for three-dimensional conservation laws. J. Comput. Phys. 2004, 201, 238–260. [Google Scholar] [CrossRef]
  34. Gottlieb, S.; Shu, C.W. Total Variation Diminishing Runge-Kutta schemes. Math. Comput. 1998, 67, 73–85. [Google Scholar] [CrossRef] [Green Version]
Figure 1. Input variables used for the LSTM network for Madrid.
Figure 1. Input variables used for the LSTM network for Madrid.
Mathematics 11 01436 g001
Figure 2. Output from the LSTM network for Madrid.
Figure 2. Output from the LSTM network for Madrid.
Mathematics 11 01436 g002
Figure 3. RSME for LSTM networks per Autonomous Community.
Figure 3. RSME for LSTM networks per Autonomous Community.
Mathematics 11 01436 g003
Figure 4. RSME values for the numerical simulations per Autonomous Community.
Figure 4. RSME values for the numerical simulations per Autonomous Community.
Mathematics 11 01436 g004
Figure 5. SVIR with Diffusion Manufactured Solutions for Fortran code implementation Exact vs Numeric Solution.
Figure 5. SVIR with Diffusion Manufactured Solutions for Fortran code implementation Exact vs Numeric Solution.
Mathematics 11 01436 g005
Figure 6. PDE Numerical simulation of infected individuals (I), with Equation (44).
Figure 6. PDE Numerical simulation of infected individuals (I), with Equation (44).
Mathematics 11 01436 g006
Figure 7. Numerical simulations of infected individuals (I) for Madrid.
Figure 7. Numerical simulations of infected individuals (I) for Madrid.
Mathematics 11 01436 g007
Table 1. Data sources, assumptions, and pre-processing considerations.
Table 1. Data sources, assumptions, and pre-processing considerations.
VariableData SourceAssumptionsPre-Processing
InfectedCNE 1, SiViES 2Only reported casesMA-7 3
RecoveredN/A 4Based on infectedMA-7
VaccinatedMinistry of Health of Spain (Ministerio de Sanidad)MA-7
MobilityMinistry of Transport, Mobility and Urban Agenda (Ministerio de Transportes, Movilidad y Agenda UrbanaSet as 1 after data were no longer available 5MA-7
PopulationINECalculated the variation from available points in time and divided the variation in daily changes for smoothnessMA-7
1 CNE: National Center of Epidemiology of Spain, Centro Nacional de Epidemiología. 2 SiViES: a platform of the CNE, Sistema de Vigilancia de España. 3 MA-7: defined as moving averages of 7 days. 4 N/A: reliable data source is unavailable. 5 Mobility data are only available from 29 February 2020, to 9 May 2021; they are the percentage of mobility compared to a reference period (established from 14–20 February 2020). Manually tabulated from the published dashboards since processed data were not made available for download.
Table 2. SVIR with Diffusion Manufactured Solutions for Fortran code Error Norms and corresponding Order, for t = 0.25.
Table 2. SVIR with Diffusion Manufactured Solutions for Fortran code Error Norms and corresponding Order, for t = 0.25.
Susceptible
Cells L 2 Order L 1 Order L Order
40 1.2000 × 10 1 1.9140 × 10 1 1.0900 × 10 1
80 8.8400 × 10 3 3.76 1.3700 × 10 2 3.80 8.1730 × 10 3 3.74
160 5.4800 × 10 3 0.69 8.6100 × 10 3 0.67 4.9680 × 10 3 0.72
320 8.1200 × 10 4 2.75 1.2690 × 10 3 2.76 7.4320 × 10 4 2.74
640 7.1200 × 10 5 3.51 1.1655 × 10 4 3.44 5.8970 × 10 5 3.66
Vaccinated
Cells L 2 Order L 1 Order L Order
40 2.5030 × 10 2 3.3300 × 10 2 2.6520 × 10 2
80 1.7500 × 10 3 3.84 2.3000 × 10 3 3.86 1.8660 × 10 3 3.83
160 1.1260 × 10 3 0.64 1.4920 × 10 3 0.62 1.1980 × 10 3 0.64
320 1.6415 × 10 4 2.78 2.1670 × 10 4 2.78 1.7494 × 10 4 2.78
640 1.6572 × 10 5 3.31 2.2510 × 10 5 3.27 1.7450 × 10 5 3.33
Infected
Cells L 2 Order L 1 Order L Order
40 3.1750 × 10 2 3.5810 × 10 2 3.8650 × 10 2
80 2.2700 × 10 3 3.81 2.5500 × 10 3 3.81 2.5990 × 10 3 3.89
160 1.4350 × 10 3 0.66 1.6100 × 10 3 0.66 1.7498 × 10 3 0.57
320 2.1092 × 10 4 2.77 2.3680 × 10 4 2.77 2.5086 × 10 4 2.80
640 2.0140 × 10 5 3.39 2.2130 × 10 5 3.42 2.9039 × 10 5 3.11
Recovered
Cells L 2 Order L 1 Order L Order
40 7.9300 × 10 2 1.3300 × 10 1 6.7040 × 10 2
80 5.6110 × 10 3 3.82 9.3380 × 10 3 3.83 4.7760 × 10 3 3.81
160 3.5660 × 10 3 0.65 5.9560 × 10 3 0.65 3.0220 × 10 3 0.66
320 5.2250 × 10 4 2.77 8.7110 × 10 4 2.77 4.4377 × 10 4 2.77
640 5.0478 × 10 5 3.37 8.5370 × 10 5 3.35 4.2030 × 10 5 3.40
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Wyss, A.; Hidalgo, A. Modeling COVID-19 Using a Modified SVIR Compartmental Model and LSTM-Estimated Parameters. Mathematics 2023, 11, 1436. https://doi.org/10.3390/math11061436

AMA Style

Wyss A, Hidalgo A. Modeling COVID-19 Using a Modified SVIR Compartmental Model and LSTM-Estimated Parameters. Mathematics. 2023; 11(6):1436. https://doi.org/10.3390/math11061436

Chicago/Turabian Style

Wyss, Alejandra, and Arturo Hidalgo. 2023. "Modeling COVID-19 Using a Modified SVIR Compartmental Model and LSTM-Estimated Parameters" Mathematics 11, no. 6: 1436. https://doi.org/10.3390/math11061436

APA Style

Wyss, A., & Hidalgo, A. (2023). Modeling COVID-19 Using a Modified SVIR Compartmental Model and LSTM-Estimated Parameters. Mathematics, 11(6), 1436. https://doi.org/10.3390/math11061436

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop