A Modified SIRD Model to Study the Evolution of the COVID-19 Pandemic in Spain

symmetry S S Article A Modified SIRD Model to Study the Evolution of the COVID-19 Pandemic in Spain Vicente Martínez   Citation: Martínez, V. A Modified SIRD Model to Study the Evolution of the COVID-19 Pandemic in Spain. Symmetry 2021,13, 723. https:// doi.org/10.3390/sym13040723 Academic Editor: Stephen Anco Received: 14 March 2021 Accepted: 14 April 2021 Published: 19 April 2021 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2021 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). Departament de Matemàtiques, Institut de Matemàtiques i Aplicacions de Castelló, Universitat Jaume I, E-12071 Castelló, Spain; [email protected] Abstract: In this paper, we use an SIRD model to analyze the evolution of the COVID-19 pandemic in Spain, caused by a new virus called SARS-CoV-2 from the coronavirus family. This model is governed by a nonlinear system of differential equations that allows us to detect trends in the pandemic and make reliable predictions of the evolution of the infection in the short term. This work shows this evolution of the infection in various changing stages throughout the period of maximum alert in Spain. It also shows a quick adaptation of the parameters that define the disease in several stages. In addition, the model confirms the effectiveness of quarantine to avoid the exponential expansion of the pandemic and reduce the number of deaths. The analysis shows good short-term predictions using the SIRD model, which are useful to influence the evolution of the epidemic and thus carry out actions that help reduce its harmful effects. Keywords: pandemic; numerical simulation; SIRD model; mathematical modeling 1. Introduction The detection of a new coronavirus (SARS-CoV-2) in December 2019 by the Municipal Health Commission in Wuhan (China) and its subsequent expansion throughout the world caused the World Health Organization (WHO) to declare a global pandemic of this infection on 11 March 2020. Its disease pattern has been named COVID-19. As of 2 May 2020, this disease has presented 3,272,202 confirmed cases and 230,104 deaths globally—with 25,264 in Spain [1]. In this work, we perform an analysis of the evolution of the epidemic in Spain using an SIRD model, which will be detailed below. The use of SIRD-type models to study epidemics has been very popular for decades [ 2 – 7 ]. They were proposed in 1927 [ 8 ] by two Scottish scientists: William O. Kermack, a biochemist, and Anderson G. McKendrick, a physicist and epidemiologist. For this reason, they are also known as Kermack–McKendrick models. These models show how an infection evolves when a certain number of sick people are detected within a healthy population. They are governed by systems of differential equations, which means that the solutions they provide are highly accurate from a mathematical point of view, but small variations in their parameters can cause significant changes in the solution. Another drawback is that accurate knowledge of the evolution of the epidemic requires data that are only known at the end of it. Therefore, at the moment when the epidemic is taking place, it is difficult to predict its evolution in the long term. However, the simplicity of the model allows us to make quick parameter adjustments and carry out instant simulations and measure trends in the evolution of the epidemic in the short term. In this paper, our aim is to analyze what has happened during the development of the epidemic and demonstrate how the SIRD model has great adaptability to predict the behavior of the COVID-19 pandemic in Spain during the course of the epidemic. Our analysis seeks to detect trends that allow short-term decisions to be made. Fortunately, in today’s society, the COVID-19 pandemic is not expected to be as devastating as the Spanish Flu Pandemic at the beginning of the last century, which took Symmetry 2021,13, 723. https://doi.org/10.3390/sym13040723 https://www.mdpi.com/journal/symmetry Symmetry 2021,13, 723 2 of 14 place between 1918 and 1920 [ 9 – 11 ]. That terrible pandemic resulted in more than 40 million victims worldwide and Spain was one of the most badly-affected countries, with 8 million people infected and 300,000 deaths. The healthcare and technological resources, as well as the research potential that exists today, are far superior to those of the early twentieth century; it is therefore not likely to cause as many victims as in 1918 (in proportion to the population). At that time, the entire population was exposed to the risk of infection and there were not as many means of treatment as nowadays. Therefore, it is obvious that in this case, the isolation measures (quarantine) taken can alter the number of people exposed to the virus. This is one of the most effective actions to prevent contagion [ 12 – 15 ] and, consequently, deaths [ 16 ]. In Spain, several royal decrees have been published to ensure quarantine [17–19]. In addition to the works already considered about SIRD models, interested readers can find a wide variety of models for the study of the COVID-19 pandemic in the scientific literature published in the last months: using Bayesian and stochastic techniques [ 20 – 22 ], including mobility [ 23 ], confinement and quarantine [ 15 , 24 ], fractional models [ 25 ], and logistic models [26], among others. 2. Model Description The simplest SIRD model considers three types of people: ( S ) Susceptible—the people who could become infected. ( I ) Infected—the people who are infected at that moment. ( R ) Recovered—the people who have had the disease and are now healthy. In the case of the COVID-19 pandemic, many deaths have taken place [ 1 , 12 , 27 ], so we introduce a new variable. In this regard, we consider the variable that represents deaths: ( D ) Deceased—the people who have died of the disease. If we assume that the population affected by the epidemic is N , satisfying N= S+I+R+D , the model, which we will call SIRD, it is governed by the following system of differential equations: S0(t) = −βS(t)I(t), (1) I0(t) = βS(t)I(t)−αI(t)−γI(t), (2) R0(t) = αI(t), (3) D0(t) = γI(t). (4) This model depends on three parameters: α (recovery rate per unit of time), β (infected rate per unit of time), and γ(death rate per unit of time). Figure 1shows the usual way in which the variables S , I , R , and D evolve. No ranges appear, because they depend on each particular epidemic. 0 Susceptible Infected Recovered Deceased Time Number of individuals Figure 1. Typical behavior of an epidemic’s evolution over time. Symmetry 2021,13, 723 3 of 14 On the other hand, the recovery rate α may change during the course of the epidemic, influenced by new medical treatments, which are becoming increasingly effective as research progresses. In addition, the value of α can also change when the healthcare system is on the verge of collapse and the mildly ill are discharged to complete their recovery at home. Actually, such patients are not fully recovered, but the official data do not count them as infected. Considering the presented reasons, during the evolution of the epidemic, the predictions that can be made will only be reliable in the short term. However, this does not mean that they are not useful, since knowledge of these predictions can help to make decisions to minimize harmful effects during the course of the epidemic. 2.1. Equations Governing the Model The SIRD epidemiological model that we are going to use to describe the evolution of the COVID-19 pandemic in Spain is a nonlinear ordinary differential equation system, presented in the previous section. It is governed by Equations (1)–(4), the formulation of which is described below. Usually, the parameter β , infected rate per unit of time, is considered the product of two quantities: β=p·m , where p is the probability that contact with an infected person results in an infection and m is the number of contacts made by an infected person. So, pmS is the number of people who have been infected by an infected person. Therefore, pmSI is the number of susceptible people who have been infected by all the infected people. Now, substituting the value of pm for β, we obtain Equation (1): S0=Variation of susceptible people per unit of time caused by all infected people =−βS I. On the other hand, Equations (3) and (4) are deduced from the definitions of α and γ: R0=Variation of recovered people per unit of time =αI, D0=Variation of deaths per unit of time =γI. Finally, from the two previous expressions, Equation (2) is deduced: I0=Variation of infected people per unit of time =βS I −αI−γI. Obviously, when we face an unknown epidemic, the constants to be determined to ensure adequate action are the parameters α , β , and γ . For this reason, estimating them is a fundamental task that we will address in the following sections. As discussed in the previous section, social isolation measures have been the most effective measures to control the spread of the epidemic and reduce the number of deaths [ 16 ]. In this way, the value of another important parameter has also been reduced. We refer to the number ρ0 of infected people who have been infected by an infected person (we will call this the reproduction number). In mathematical terms: ρ0(t) = β α+γS(t). If this parameter is greater than 1, then the expansion of the epidemic is exponential and therefore uncontrollable. However, when it is less than 1, the number of new infections decreases over time and the epidemic tends to subside. The reader can find more details of this parameter in [ 28 , 29 ]. For this reason, social isolation is effective, since, for a given epidemic with a given parameter β , the only way to decrease ρ0 is to decrease the number of people at risk of being infected, that is, S(0). 2.2. Estimation of Parameters αand γ Suppose that the infection has become stable. Therefore, the class of infected people will no longer increase, so the differential equation that governs this class is Symmetry 2021,13, 723 4 of 14 I0(t) = −αI(t)−γI(t),t≥0 I(0) = I0. This equation has the following solution: I(t) = I0e−(α+γ)t. In this way, the proportion of individuals who remain infected at time t≥ 0, or in other words, the probability that an individual remains infected at time t≥ 0, is given by the quotient I(t) I0 =e−(α+γ)t. Thus, the fraction of individuals who have left the infected class at time t≥ 0 will be F(t) = 1−e−(α+γ)t,t≥0, where F(t) takes the form of the distribution function of a random variable X , which is defined as the time it takes an individual to leave the infected class. In our case, this happens when an infected person has overcome the disease or has died. Calculating the density function of the variable X, we obtain f(t) = F0(t) =    (α+γ)e−(α+γ)t,t≥0, 0, t<0. Under these conditions, we can calculate the average time it takes for a sick person to leave the infected class. E[X] = Z+∞ −∞ t f (t)dt =Z+∞ 0t(α+γ)e−(α+γ)t=1 (α+γ). (5) According to the available data [ 30 , 31 ], it takes an average of 14 days for an ill person to recover, while deaths occur on average 42 days after onset of the illness. If the fraction of deaths in relation to discharges of the sick in Spain as of 14 April 2020 [ 1 ] is 18579 70853 ≈ 0.26, then E[X]≈ 0.74 × 14 + 0.26 × 42 = 21.28. From Equation (5), we have α≈ 0.0373 and γ≈0.0095 in the case of the COVID-19 pandemic in Spain. 2.3. Estimation of Parameter β In order to estimate the value of β , we consider that α and γ are known. Thus, system (1)–(4) is reduced to S0(t) = −βS(t)I(t), I0(t) = βS(t)I(t)−αI(t)−γI(t). Operating, we obtain I0 S0=βS I −(α+γ)I −βS I =−1+(α+γ) βS⇒I0=−1+(α+γ) βSS0. Therefore, I=−S+(α+γ) βln S+C, where C is the integration constant. Thus, the solutions of the system are given implicitly by I+S−(α+γ) βln S=C. (6) Symmetry 2021,13, 723 5 of 14 To continue with our estimate of β , we need to make some assumptions: S0=S( 0 ) , I0=I( 0 ) , S∞=lim t→+∞ S(t) , and I∞=lim t→+∞ I(t) = 0. Now, substituting in (6) to remove C, we obtain I0+S0−(α+γ) βln S0=S∞−(α+γ) βln S∞⇒β (α+γ)=ln(S0/S∞) I0+S0−S∞ . Finally, β= (α+γ)·ln(S0/S∞) I0+S0−S∞ . (7) Obviously, when the epidemic is evolving and has not yet ended, S∞ is unknown and it will be necessary to estimate it with the data available at each moment. 3. Analysis of the COVID-19 Pandemic in Spain During the Course of the Disease In order to have enough data for the estimates to have a significant value, we consider official data [ 1 , 27 ] provided by the Spanish Government from 11 March 2020, when the disease had been in the country for several weeks (or months, this is still unknown). We differentiate several stages in the evolution of the epidemic as the data have changed, either due to changes in case reporting criteria or due to changes in the trend of the disease. 3.1. First Stage: 11–21 March 2020 We consider data supplied by the Ministry of Health of the Spanish Government up to 21 March 2020. Table 1shows the status of the epidemic up to that date. During the first days, official data were completed with data from the GitHub Enterprise group [ 32 ]. To fit an SIRD model that approximates the known values, the following data have been taken: S0= 1, 550, 000, I0= 2039, and S∞= 5000. The rest of the parameters have been estimated as follows: • The discussion of the above section allows us to consider the Recovered/Deceased fraction as 55 183 ≈ 0.3, in accordance with the figures from 11 March 2020 (see Table 1 ). Taking 0.3 as the Recovered/Deceased fraction, we can estimate that E[X]≈ 0.7 × 14 +0.3 ×42 =22.4 days. • From Equation (5) and considering γ≈0.3 α, we obtain α=0.0343 and γ=0.0103. • From Equation (7), we estimate βonce αand γare known. So, β=1.6554 ×10−7. In this work, we have chosen to use the MATLAB Ode Solver ode45 [ 33 ] to solve the system (1)–(4). The results obtained are shown in Figure 2. In this first stage, when the epidemic is beginning to spread, the estimates show great uncertainty, although the forecasts are bleak: the number of infected people is around 800,000. It is surely this fact that led the Spanish authorities to take severe isolation measures [ 17 ]. It is true that there are more optimistic estimates, but also more pessimistic ones. The results obtained show a possible scenario. Figure 2shows how the evolution of the epidemic may have progressed if no action had been taken and if the efficacy of the treatments had not improved. The figures about those recovered and deceased do not yet show a clear trend and we will not comment on them. In Figure 2and below, the evolution of the number of susceptible people is not shown because it is an unknown figure introduced by us to adjust the other variables that are of most interest in terms of allowing us to influence the subsequent development of the epidemic. Symmetry 2021,13, 723 6 of 14 Table 1. Number of infected, recovered, and deceased people in Spain from 11–21 March 2020. Date Infected Recovered Deceased Total Cases 11 March 2039 183 55 2277 12 March 4906 193 133 5232 13 March 5679 517 195 6391 14 March 6992 517 289 7798 15 March 9070 530 342 9942 16 March 10,187 1028 533 11,748 17 March 12,206 1081 623 13,910 18 March 16,026 1107 830 17,963 19 March 17,779 1588 1043 20,410 20 March 21,874 2125 1375 25,374 21 March 24,421 2575 1772 28,768 0 50 100 150 ×105 0 2 4 6 8 10 12 I estimated R estimated D estimated I real R real D real Time (a) Number of individuals I,Rand D 0 10 20 30 40 50 60 70 ×105 0 1 2 3 4 5 6 I estimated R estimated D estimated I real R real D real Time (b) Figure 2. ( a ) Comparison of real data with the predictions of the SIRD model on 21 March 2020. (b) Zoomed-in graph of real data. 3.2. Second Stage: 22–31 March 2020 At this stage, we consider the data provided by the Ministry of Health of the Spanish Government up to 31 March 2020 (see Table 2). The values we have taken at this stage to adjust the aforementioned SIRD model have been as follows: S0= 500, 000, I0= 27, 552, and S∞= 20, 000. The estimation of the other parameters has been carried out analogously to the previous section: • We observed in the data that the Recovered/Deceased fraction follows the sequence tending towards 0.3 and E[X]≈20 days. • From Equation (5) and γ≈ 0.3 α , we obtain α= 0.0385 and γ= 0.0155. From Equation (7), we estimate the value of β once we know α and γ . So, β= 3.171 × 10 −7 . The results obtained using the MATLAB Ode Solver considered in the first stage are shown in Figure 3. It is observed that the influence of the quarantine imposed by the Spanish Government begins to be noticed and the figures, although still worrying, are no longer so alarming. Obviously, the approximations show trends that the model variables follow according to the aforementioned conditions. It is important to note that we are estimating the values of the variables that will later be officially counted according to the conditions described. In the event that these conditions change, the values of the variables would undergo changes. However, the information provided will continue to be useful for making medical and political decisions, since the observed trends are representative of the situation. In addition, the cases that are officially detected are those of people who have a serious evolution or an uncertain prognosis, for whom hospitalization is therefore required or at least require that a follow-up be carried out by Family Medicine Units. In summary, these are the cases that must be attended to urgently and followed up on. In the evolution of an epidemic caused by the action of a little-known virus, it is important to anticipate the ICU beds that may be Symmetry 2021,13, 723 7 of 14 needed, and more so in this case, where the evolution of those infected with severe cases of COVID-19 require ICU care. Table 2. Official data provided by the Spanish Government from 22–31 March 2020. Date Infected Recovered Deceased Total Cases 22 March 27,552 3355 2182 33,089 23 March 33,183 3794 2696 39,673 24 March 38,809 5367 3434 47,610 25 March 45,084 7015 4089 56,188 26 March 49,844 9357 4858 64,059 27 March 54,273 12,285 5690 72,248 28 March 57,560 14,709 6528 78,797 29 March 61,075 16,780 7340 85,195 30 March 66,969 19,259 8189 94,417 31 March 70,436 22,647 9053 102,136 0 50 100 150 ×105 0 0.5 1 1.5 2 2.5 I estimated R estimated D estimated I real R real D real Time Number of individuals I,Rand D Figure 3. Comparison of the real data with the estimations of the SIRD model during the second stage. 3.3. Third Stage: 1–13 April 2020 At this stage, we will consider the data in Table 3. Table 3. Official data provided by the Spanish Government from 1–13 April 2020. Date Infected Recovered Deceased Total Cases 1 April 73,492 26,743 10,003 110,238 2 April 76,262 30,513 10,935 117,710 3 April 78,773 34,219 11,744 124,736 4 April 80,261 38,080 12,418 130,759 5 April 81,540 40,437 13,055 135,032 6 April 83,504 43,208 13,798 140,510 7 April 84,114 48,021 14,555 146,690 8 April 85,043 52,165 15,238 152,446 9 April 85,511 55,668 15,843 157,022 10 April 86,390 59,109 16,353 161,852 11 April 86,656 62,391 16,972 166,019 12 April 87,280 64,727 17,489 169,496 13 April 86,981 67,504 18,056 172,541 The values that we have taken in this stage to adjust the SIRD model, described in previous paragraphs, have been S0= 170, 000, I0= 73, 492, and S∞= 20, 000. The other parameters have been estimated as follows: Symmetry 2021,13, 723 8 of 14 • The Recovered/Deceased fraction has been observed to follow the sequence: 0.65, . . . , 0.58, . . . , 0.31, . . . , 0.27, . . . . Consequently, we have estimated the value 0.22, which produces good approximations for the data that are already known. Taking 0.22 as the Recovered/Deceased fraction, we can estimate that E[X]≈20 days. • From Equation (5) and γ≈ 0.22 α , we obtain α= 0.041 and γ= 0.009. From Equation (7) and the values of αand γ, we estimate β. So, β=4.7878 ×10−7. The results obtained using the MATLAB Ode Solver considered in previous stages are shown in Figure 4. These results are conditioned by the number of tests carried out, which provide us with the number of infected people, the number of patients discharged (in many cases the patients continue their convalescence at home, due to the saturation of hospitals), and the effectiveness of treatments. 20 40 60 80 100 120 140 ×105 0 0.5 1 1.5 2 2.5 I estimated R estimated D estimated I real R real D real Time Number of individuals I,Rand D Figure 4. Comparison of the real data with the estimates of the SIRD model during the third stage. In the current part of the process in which we find ourselves, that is, fully immersed in the course of the epidemic, there are many doubts about it; if we take into account the number of hidden carriers of the disease, the data obtained would change shape very significantly [ 16 ]. However, the estimates that have been obtained are very useful for managing actions in the short term. There may be 7 million people infected in Spain, as indicated in the previous reference, but at the moment what is of interest is managing cases with severe symptoms. These patients who require a follow-up are those detected among the people who are cared for in health centers and hospitals, and their tests indicate that they are infected and are included in the official figures of those infected. SIRD models can also estimate the magnitude of the death toll. Figure 4indicates that the number of deceased people may be between 40,000 and 50,000 people with the data available on 13 April 2020. The number of deaths will be known more accurately once the epidemic has been overcome and comparative analyses of the data on deaths are carried out for the months during which the epidemic has remained active and the same months in previous years. However, the data provided by the models can help make decisions that modify the final figure. 3.4. Fourth Stage: 14–23 April 2020 We consider the data in Table 3, adding the values in Table 4. We carry out estimations with the parameters of the third stage. At this stage, the Spanish Government has reported an increase in the number of tests carried out [ 1 , 27 ], which largely disarrange the parameters of the model considered in previous stages. Symmetry 2021,13, 723 9 of 14 Figure 5shows the evolution of the pandemic during this stage. A more detailed analysis indicates that on 14 April, there was a slight mismatch; which was accentuated on 15 April, with an upward trend in real data with respect to those expected by the SIRD measure. On 16 April, the change in trend was already evident in the data; during the days after 16 April, the change in the trend of the data was confirmed. Table 4. Official data provided by the Spanish Government from 14–23 April 2020. Date Infected Recovered Deceased Total Cases 14 April 70,853 88,201 18,579 177,633 15 April 74,797 88,889 19,130 182,816 16 April 72,963 95,627 19,478 188,068 17 April 74,662 97,021 20,043 191,726 18 April 77,357 98,134 20,453 195,944 19 April 80,587 98,771 20,852 200,210 20 April 82,514 100,382 21,282 204,178 21 April 85,915 100,757 21,717 208,389 22 April 89,250 101,617 22,157 213,024 23 April 92,355 88,111 22,524 202,990 20 40 60 80 100 120 140 ×105 0 0.5 1 1.5 2 2.5 I estimated R estimated D estimated I real R real D real Time (a) Number of individuals I,Rand D 10 20 30 40 50 60 70 ×104 0 2 4 6 8 10 12 I estimated R estimated D estimated I real R real D real Time (b) Figure 5. ( a ) Comparison of the state of the real data with the predictions of the SIRD model on 23 April 2020. (b) Zoomed-in graph of real data. Given the circumstances, the model needs to be revised. The Spanish authorities have clarified [ 1 , 27 ] the way of accounting for the data and the types of tests being carried out. Figure 6shows the official data regarding the increase in the number of people infected and the types of tests performed. Symmetry 2021,1, 0 9 of 14 Figure 5shows the evolution of the pandemic during this stage. A more detailed analysis indicates that on 14 April, there was a slight mismatch; which was accentuated on 15 April, with an upward trend in real data with respect to those expected by the SIRD measure. On 16 April, the change in trend was already evident in the data; during the days after 16 April, the change in the trend of the data was confirmed. Table 4. Official data provided by the Spanish Government from 14–23 April 2020. Date Infected Recovered Deceased Total Cases 14 April 70,853 88,201 18,579 177,633 15 April 74,797 88,889 19,130 182,816 16 April 72,963 95,627 19,478 188,068 17 April 74,662 97,021 20,043 191,726 18 April 77,357 98,134 20,453 195,944 19 April 80,587 98,771 20,852 200,210 20 April 82,514 100,382 21,282 204,178 21 April 85,915 100,757 21,717 208,389 22 April 89,250 101,617 22,157 213,024 23 April 92,355 88,111 22,524 202,990 20 40 60 80 100 120 140 ×105 0 0.5 1 1.5 2 2.5 I estimated R estimated D estimated I real R real D real Time (a) Number of individuals I,Rand D 10 20 30 40 50 60 70 ×104 0 2 4 6 8 10 12 I estimated R estimated D estimated I real R real D real Time (b) Figure 5. ( a ) Comparison of the state of the real data with the predictions of the SIRD model on 23 April 2020. (b) Zoomed-in graph of real data. Given the circumstances, the model needs to be revised. The Spanish authorities have clarified [ 1 , 27 ] the way of accounting for the data and the types of tests being carried out. Figure 6shows the official data regarding the increase in the number of people infected and the types of tests performed. Daily increase PCR test Antibodies test Time Number of new cases daily % daily increase Figure 6. Types of tests carried out and the % increase in the number of infected people on 23 April 2020. Source: Spanish Government [1] on 24 April 2020. Figure 6. Types of tests carried out and the % increase in the number of infected people on 23 April 2020. Source: Spanish Government [1] on 24 April 2020.