Deterministic and Stochastic Dynamics of COVID-19: The Case Study of Italy and Spain

In December 2019, a severe respiratory syndrome (COVID-19) caused by a new coronavirus (SARS-CoV-2) was identi ﬁ ed in China and spread rapidly around the globe. COVID-19 was declared a pandemic by the World Health Organization (WHO) in March 2020. With eventually substantial global underestimation, more than 225 million cases were con ﬁ rmed by the end of August 2021, counting more than 4.5 million deaths. COVID-19 symptoms range from mild (or no symptoms) to severe illness, with disease severity and death occurring according to a hierarchy of risks, with age and preexisting health conditions enhancing the risks of disease severity manifestation. In this paper, a mathematical model for COVID-19 transmission is proposed and analyzed. The model strati ﬁ es the studied population into two groups, older and younger. Applied to the COVID-19 outbreaks in Spain and in Italy, we ﬁ nd the disease-free equilibrium and the basic reproduction number for each case study. A sensitivity analysis to identify the key parameters which in ﬂ uence the basic reproduction number, and hence regulate the transmission dynamics of COVID-19, is also performed. Finally, the model is extended to its stochastic counterpart to encapsulate the variation or uncertainty found in the transmissibility of the disease. We observe the variability of the infectious population ﬁ nding its distribution at a given time, demonstrating that for small populations, stochasticity will play an important role.


Introduction
The coronavirus pandemic and its emerging variants are currently a major global public health threat. Globalisation has speeded up the spread of infections over a short period of time. This has an impact on the public healthcare system and is also detrimental to the economic development of many countries. Since the outbreaks began, more than 215 million cases were confirmed by mid of August 2021, with more than 4 million deaths [1].
COVID-19 symptoms can range from mild (or no symptoms) to severe illness, with disease severity and death occurring according to a hierarchy of risks, with age and preexisting health conditions enhancing risks of disease severity [2]. Vaccines against COVID-19 have been developed in record time and are now globally distributed [3][4][5][6][7][8]. The analysis of the impact of different vaccine administration is ongoing, mostly using previous research experiences applied to other infectious diseases [9][10][11][12][13][14][15][16][17] and will be discussed in detail in our forthcoming publications [18,19].
As an example of the impact of the pandemic in Europe, Spain has reported, up to date, more than 4.7 million cases with around 83 thousand deaths [20], while in Italy, although the total number of cases are similar, 4.5 million, a higher mortality rate was reported, counting more than 128 thousand deaths [21]. It is important to notice that mortality is higher in older than younger individuals in all populations [22,23].
As the COVID-19 pandemic progressed, research on mathematical modeling became imperative and very influential to understand the epidemiological dynamics of disease spreading. Task forces were created to assist public health managers and governments, with many research papers being recently published. Applied to the outbreaks in the Basque Country, Spain, a flexible framework was developed within the so-called COVID-19 Basque Modeling Task Force (BMTF). As an extension of the basic SHAR (Susceptible-Hospitalized-Asymptomatic-Recovered) model, the SHARUCD models were parameterized and validated with epidemiological data continuously collected and provided by the Basque Health Department and the Basque Health Service (Osakidetza) and have been used (up until now) to monitor COVID-19 spreading and control over the course of the pandemic [18,19,[24][25][26][27][28]. Modeling refinements and results on the evolution of the epidemic in the Basque Country are regularly updated and publicly available on the "SHARUCD Dashboard" [29].
Over the course of the pandemic, a broad spectrum of research has been produced, e.g., statistical work using two common approaches, the SIR model and a log-linear model, analyzing the available empirical data and estimating the reproduction values for Spain and Italy countries [30]. Deterministic and stochastic models were developed [31,32] to study the effects of facial masks and hospitalization of symptomatic people and asymptomatic quarantine of people on the prevalence of the coronavirus outbreak in India [31] and the transmission dynamics of the COVID-19 in Wuhan, China [32]. And still, many questions remain to be investigated.
In this paper, we present a mathematical model to describe COVID-19 dynamics in Spain and Italy. The population is divided into two groups, young and old. Considering simple mass action type incidence, the model is formulated by assuming that the course COVID-19 infection leads to a different outcome for the elderly population as compared to the younger population. This paper is organized as follows. Section 2 describes the model formulation, followed by the model analysis in Section 3, showing the existence of equilibrium and basic reproduction number. Section 4 presents the results for the sensitivity analysis for the parameters involved in reproduction number. Section 5 describes impact of different parameter on COVID-19 prevalence. Section 6 presents the stochastic modeling approach, and in Section 7, the simulation results are shown. Finally, in Section 8, we conclude this work, with a discussion on both modeling approaches.

The Model
This model is a refinement of the model proposed by Srivastav et al. [33]. A new parameter ε is introduced to differentiate the infectivity of young infected individuals ðI 1 Þ with respect to the baseline infectivity of elderly individuals I 2 ðt Þ in a population of N individuals.
The value of ε can be tuned to reflect different situations: a value of ε > 1 reflects the fact that young individuals, which are likely to develop mild disease and higher mobility, have larger infectivity, than elderly individuals, which are at higher risk of developing severe disease and are more likely to be detected and isolated [24]. On the other hand, ε < 1 indicates that young individuals have smaller infectivity than elderly individuals, and that could be justified due to lower or higher viral load during the infection, which is correlated with disease symptoms. Here, the assumption relies on the epidemiological observation of young population developing mild or no symptoms with lower viral load, affecting disease transmissibility, versus severe disease and higher viral load, mostly observed in older ages.
The total human population NðtÞ is divided into eight compartments, stratified into two age classes, namely, young and old: susceptible S 1 ðtÞ and S 2 ðtÞ, exposed E 1 ðtÞ and E 2 ðtÞ, and infected I 1 ðtÞ and I 2 ðtÞ for the young and for the old, respectively. Two extra classes to accommodate individuals from both age groups are also considered: quarantined/hospitalized HðtÞ for those identified as COVID-19-positive case with symptoms and needing medical assistance and finally the recovered individual class RðtÞ.
We assume that the transitions/movements, from one disease related class to another, are different between the elderly and the young individuals. However, the HðtÞ class includes both groups. For the mathematical modeling framework development, we make the following assumptions: (1) Total population N is constant (2) The susceptible young individuals S 1 become exposed to the infection and join the young exposed class E 1 on effective contacts with infectious human population I 1 , I 2 and H at rates β 1 and β 3 , respectively, with β 3 < β 1 (3) Exposed people E 1 will move to I 1 with rate of η 1 . If young exposed people E 1 will get contact with I 1 , I 2 unknowingly, then exposed people E 1 will move fast into I 1 class with rate γ 1 (4) The susceptible old individuals S 2 become exposed to the infection and join the old exposed class E 2 on effective contacts with infectious human population I 1 , I 2 and ðHÞ at the rates β 2 and β 4 , respectively, with β 4 < β 2 (5) Exposed people E 2 will move to I 2 with rate of η 2 ; if old exposed people E 2 will get contact with I 1 , I 2 unknowingly, then exposed people E 2 will move fast in I 2 with rate γ 2 (6) Infected, young I 1 and old I 2 , will move to hospitalized class ðHÞ with rates ν 1 and ν 2 , respectively (7) Hospitalized people ðHÞ will get recovered from COVID-19 and join recovered class ðRÞ with rate α The schematic diagram of our proposed model is shown in Figure 1, and the mathematical model is given as follows: 2 Computational and Mathematical Methods Here, all the parameters are positive and the description of these parameters is given in Table 1. Parameter values given in Tables 2 and 3 are already defined in [33].

Analysis of the Model
We consider the system (1) and find the disease-free equilibrium. For our model, we have disease-free equilibrium as We find the basic reproduction number R 0 by following the next-generation matrix method described in [38]. Following the same notations as in [38], we find the vector F and V as follows:  Disease-related death rate in I 1 compartment δ 2 Disease-related death rate in I 2 compartment δ 3 Disease-related death rate in H compartment ν 1 Rate of detection/quarantine in I 1 compartment ν 2 Rate of detection/quarantine in I 2 compartment η 1 Rate of progression of individuals from E 1 to I 1 Recovery rate of quarantine/hospitalized people 3 Computational and Mathematical Methods and V = Jacobian of V at And it follows that where Three eigenvalues of the above matrix are zero and the remaining two are the roots of the following quadratic equation: So the basic reproduction number (R 0 ) is the positive root of the above quadratic and is given by where

Sensitivity Analysis
We also perform a sensitivity analysis for the parameters involved in reproduction number ðR 0 Þ, which reflects that the increase or decrease in these parameter values will lead to the increase or decrease of ðR 0 Þ. The sensitivity of R 0 to different parameters is shown in Figure 2. This analysis is performed to evaluate which parameters have the highest impact on R 0 and hence being targeted as the most effective intervention measures for each case study. The sensitivity indices allow to measure the relative change in a variable when parameter changes. For that, we use the forward sensitivity index of a variable with respect to a given parameter, which is defined as the ratio of the relative change in the variable to the relative change in the parameter. If the variable is varying with respect to a parameter, then the sensitivity index is defined using partial derivatives [39]. The normalized forward sensitivity index of R 0 , which is differentiable with respect to a given parameter p, is defined by The above formula can be used to compute the analytical expression for the sensitivity of R 0 to each parameter included in the system. From Figure 2, we can conclude that β i , ν i for i = 1, 2 and p, for both case study, Italy and Spain, are very sensitive parameters, with small changes in these parameters leading to a significant change in the value of R 0 .

Impact of Different Parameters on Prevalence of COVID-19
From Figures 3 and 4, we observe that the increase of transmission rate β 1 , from infected young and old individuals to  Computational and Mathematical Methods and 6 show that the infections in old populations I 2 will increase as β 2 , the transmission rate from infected young and old individuals to susceptible old individuals, increases. Nevertheless, the increase of ν 2 , the detection rate of infected old individuals, will decrease the number of infections in the old population, but also increasing the overall hospitalizations. From these figures, we can conclude that controlling the transmission between human to human and increasing the detection rates, providing better treatment for the infected hospitalized people, will be of a major importance to control the epidemics in Italy and in Spain, with a significant decrease of number of infections in the population.
Finally, Figure 7 shows the effect of the variation of the parameter ν 1 on the number of infections. As ν 1 increases, the overall infected population decreases. Note that variations of the parameter ν 2 show similar effects on the number of overall infections.

Stochastic Model
As all natural systems are prone to stochastic perturbations, we extended our deterministic model, equation system (1), to the corresponding stochastic model. The derivation of the stochastic model and its analysis are important when populations are small and hence with the dynamics being severely affected by small changes in the parameter values. Thus, for the initial phase of the disease outbreak, such as COVID-19, the stochastic model setup is the most appropriate modeling approach to be used for a local epidemiological evaluation.
Using the approach of [26,41,44], we construct a matrix M such that Ω = MM T , where M is an 8 × 10 matrix given as where Then, the Itô stochastic differential equation system has the form with initial condition

Stochastic Simulation Results
To emphasize the impact of stochasticity in the system, we simulate the stochastic model shown in system (28) by using Euler-Maruyama method [45]. For this purpose, we use the following set of parameter values given in Tables 2 and 3. We perform simulations of system (28) for 120 days and 100 simulation runs. First, we compare the mean of 100 runs of stochastic model simulation with the results of corresponding deterministic model and plot the time series of all the variables; see Figure 8 for Italy and Figure 9 for Spain. From these figures, we observe that the mean of 100 runs of stochastic simulation is very close to the simulation results of the deterministic model for the susceptible population ðS 1 Þ and ðS 2 Þ, whereas for the other populations a small deviation between stochastic and deterministic simulations results is observed.
In order to understand these results better, we also plot the distributions of the infected young ðI 1 Þ and old ðI 2 Þ populations at the 50th, 80th, and 120th days. These can be seen in Figures 10 and 11 for Italy and Spain, respectively. One can easily see the change in distribution of the population as time progresses.

Results and Discussion
With an unprecedented global health burden, the COVID-19 pandemic has changed the behavior of societies around the globe which were significantly affected by the extreme measures implemented to control disease transmission.
In this paper, a mathematical model for COVID-19 transmission is formulated and analyzed. We computed the disease-free equilibrium and the basic reproduction number R 0 . Sensitivity analysis was performed to find the key parameters that are most sensitive to basic reproduction number R 0 . The deterministic model was extended to its In this work, we have introduced a new parameter ε that plays a very important role to reduce the infectivity of the young population. While ε > 1 represents higher infectivity of young population than the infectivity of older population, the assumption of ε < 1 represents a lower infectivity of young population than the infectivity of older population. A detailed analysis of this specific parameter is ongoing.
Results presented here support the reduction of the transmission rate between human to human, by the proper use of nonpharmaceutical intervention and vaccination, to affect the most the dynamics of the pandemic. Moreover, a fast detection of the infected individuals would lead to a better treatment, decreasing the disease mortality rate in the population. We would like to emphasize that this is an 13 Computational and Mathematical Methods ongoing work and models will be refined and extended, using the present model as a baseline for future research. As an example, applied to the Basque Country, Spain, scenario, the incorporation of asymptomatic classes for both young and old groups is proposed to evaluate its impact on the control measures implemented over the pandemic time.

Data Availability
Previously reported data on COVID-19 cases for Italy and Spain were used to support this study and are available at doi:10.1002/mma.7344. These prior studies (and datasets) are cited at relevant places within the text as references [20,21].

Conflicts of Interest
The authors declare no potential conflict of interests. 14 Computational and Mathematical Methods