Analysis of a COVID-19 compartmental model: a mathematical and computational approach

1.   Introduction

Mathematical models for the transmission dynamics of infectious diseases have been a powerful tool to understand and control epidemics. An important amount of these mathematical models are given by systems of ordinary differential equations (ODE's), considering a continuous-time framework. In the last decades, a great number of compartmental models has been proposed and many of them share some of the assumptions of the well known SIR model, first proposed by Kermack-McKendrick in 1927 ^[5] and also given by a system of ODE's.

The importance of compartmental models given by systems of ODE's has been even more highlighted since the beginning of the COVID-19 pandemic. In fact, SIR, SEIR, SEIRD-type models, among many others, have been used to analyze, predict and control the spread of SARS-CoV-2 virus worldwide, in what follows we refer to some of this models applied to COVID-19 given by systems of ODE's. The effect of the lockdown on the spread of COVID-19 was analyzed in ^[4] by considering a mathematical model that assesses the imposition of the lockdown in Nigeria. Besides lockdown measures, in ^[6] quarantine, and hospitalization of COVID-19 infected individuals is analyzed. A SIDARTHE compartmental model, is proposed in ^[10] which discriminates between infected individuals depending on whether they have been diagnosed and on the severity of their symptoms and is fitted to the COVID-19 epidemic in Italy. The authors claim that restrictive social-distancing measures will need to be combined with widespread testing and contact tracing to end the ongoing COVID-19 pandemic. The spread of COVID-19 in Portugal is modeled in ^[14], fitting the model to real data of simultaneously the number of active infected and hospitalized individuals, SEIR types models applied to COVID-19 are proposed in, for example, ^[12,16,21], were the models are fitted to real data. Quarantine and lockdown measures are considered in ^[24] with SEIQR type models. Fractional SEIR type models are also important in modeling and predicting the spread of infectious diseases. Their advantages in comparison to other type of models are highlighted in, for example, ^[1] were COVID-19 epidemic in Pakistan is analyzed through a fractional SEIR type model using the operator of Atangana-Baleanu. The Atangana-Baleanu derivative is also considered in ^[15]. A nabla discrete ABC-fractional order COVID-19 model is analyzed in ^[13]. In ^[2] a fractional model is proposed to study the first COVID-19 outbreak in Wuhan, China. The outbreak in Wuhan is also considered in ^[7,19]. A Caputo fractional model is proposed in ^[18] and fitted to COVID-19 data of Galicia, Spain and Portugal. A Bats-Hosts-Reservoir-People transmission fractional-order COVID-19 model for simulating the potential transmission with the thought of individual response and control measures by the government is analyzed in ^[23]. Stochastic epidemic models for COVID-19 spread and control are proposed in, e.g., ^[8,11]. A new class of distributions are applied to the generalized log-exponential transformation of Gumbel Type-II and implemented using real data of COVID-19, in ^[27]. The exponential transformation of Gumbel Type-II distribution for modeling COVID-19 data is used to analyze the number of deaths due to COVID-19 for Europe and China, in ^[28].

Among other important issues, in the mathematical analysis of compartmental models given by systems of ODE's, we may emphasize the stability analysis of the equilibrium points, the basic reproduction number and the model fitting to real data. Although the difficulty of this analysis depends on the complexity of the model under study, part of it is common to the majority of the models and may become simplified if we use adequate mathematical software. In this paper, we show how to use free and open-source software for the mathematical analysis of the models. We focus on the mathematical software SageMath (version 9.2) ^[22] and $\texttt{Python}$ programming language ^[20], version 3.8, and the Python Libraries: Numpy (version 1.18.5), Pandas (version 0.25.3), Scipy (version 1.4.1) and Matplotlib (version 2.0.0).

In this note, we consider a SAIRP model, given by a system of five ODE's, for the transmission of SARS-CoV-2, first proposed in ^[26] and after generalized to piecewise constant parameters and complex networks model in ^[25]. We provide all the codes that allow us to compute the equilibrium points, basic reproduction number, visualize the global stability of the equilibrium points and by considering piecewise constant parameters. Moreover, we provide the Python code that allows to estimate some of the piecewise constant parameters and fit the model to the real data of COVID-19 transmission in Portugal, from March 2, 2020 until April 15, 2021. It is important to note, that all the codes are elementary and thus can be easily adapted to other compartmental models, see e.g., ^{[3,6,9,17,24,30]}.

The paper is organized as follows. In Section 2, the SAIRP model, given by a system of ODE's, is described. The equilibrium points of the SAIRP model are computed in Section 3 and their global stability is illustrated using SageMath (version 9.2) ^[22] and Python (version 3.8) ^[20]. The generalized SAIRP model with piecewise constant parameters is presented and the Python code for the model fitting to COVID-19 Portuguese real data is given, in Section 4.

2.   Mathematical model

In this section, we consider the compartmental SAIRP mathematical model proposed in ^[25,26] for the transmission dynamics of COVID-19. We start by recalling the assumptions of the model. The total population $N(t)$, with $t \in [0, T]$ (in days) and $T > 0$, is subdivided into five classes: susceptible individuals ($S$); asymptomatic infected individuals ($A$); active infected individuals ($I$); removed (including recovered individuals and COVID-19 induced deaths) ($R$); and protected individuals ($P$). Therefore, $N(t) = S(t) + A(t) + I(t) + R(t) + P(t)$, considering a continuous time framework, with $t \in [0, T]$. The total population is homogeneous and has a variable size, with constant recruitment rate, $\Lambda$, and natural death rate, $\mu > 0$. The susceptible individuals $S$ become infected by contact with active infected $I$ and asymptomatic infected $A$ individuals, at a rate of infection $\beta \frac{\left(\theta A + I \right)}{N}$, where $\theta$ represents a modification parameter for the infectiousness of the asymptomatic infected individuals $A$ and $\beta$ represents the transmission rate. Only a fraction $q$ of asymptomatic infected individuals $A$ develop symptoms and are detected, at a rate $v$. Active infected individuals $I$ are transferred to the recovered/removed individuals $R$, at a rate $\delta$, by recovery from the disease or by COVID-19 induced death. A fraction $p$, with $0 < p < 1$, is protected (without permanent immunity) from infection, and is transferred to the class of protected individuals $P$, at a rate $\phi$. A fraction $m$ of protected individuals $P$ returns to the susceptible class $S$, at a rate $w$. Let $\nu = v q$ and $\omega = w m$. The previous assumptions are described by the following system of ordinary differential equations:

The compartments and parameters descriptions and notations of system (2.1) are resumed in Table 1.

Consider the compact invariant region

The model (2.1) is biologically and mathematically well-posed, that is for any initial condition $x_0 = (S_0, \, A_0, \, I_0, \, R_0, \, P_0)^T \in \Omega$, the system (2.1) admits a unique solution defined on $[0, \, \infty)$, whose components are non-negative. Furthermore, the region $\Omega$ defined by (2.2) is positively invariant ^[25].

3.   Equilibrium points and stability analysis

In this section, we use the free and open-source mathematical software SageMath (version 9.2) ^[22] to help us

compute the equilibrium points and the basic reproduction number of model (2.1). The global stability analysis of the equilibrium points is illustrated through numerical simulations developed using Python (version 3.8) ^[20].

3.1.   Computing equilibrium points and basic reproduction number in SageMath

The model (2.1) has two equilibrium points:

● disease-free equilibrium, denoted by $\Sigma_0$, given by

● endemic equilibrium, $\Sigma_+$, whenever $R_0 > 1$, given by

with

where the basic reproduction number, $R_0$, is given by

The SageMath code reads as follows:

var('S', 'A', 'I', 'R', 'P', 'phi', 'mu', 'nu', 'delta', \

'omega', 'Lambda', 'theta', 'beta', 'p', 'lambdat');

a0 = p*phi+mu; a1 = nu+mu; a2 = delta+mu; a3 = omega+mu;

lambdat = beta*(A*theta+I)*(1-p);

N = S+A+I+R+P;

eqS = Lambda + omega*P-lambdat*S/N-a0*S;

eqA = lambdat*S/N-a1*A;

eqI = A*nu-I*a2;

eqR = I*delta-R*mu;

eqP = S*p*phi-P*a3;

pretty_print((eqS+eqA+eqI+eqP+eqR).full_simplify())

sistema = [eqR, eqA = = 0, eqI = = 0, eqP = = 0, eqS = = 0];

sol = solve(sistema, A, I, P, R, S); pretty_print(sol)

The output $\texttt{sol}$ gives the two biological meaningful equilibrium points $\Sigma_0$ and $\Sigma_+$ given by models (3.1) and (3.2), respectively.

To compute the basic reproduction number $R_0$, (3.3), we follow the approach presented in ^[29].

3.2.   Stability analysis

The following theorems hold and are numerically illustrated in Figure 1.

Theorem 1 (Local stability of the DFE, ^[25]). The disease-free equilibrium, $\Sigma_0$, is locally asymptotically stable whenever $R_0 < 1$.

Theorem 2 (Global stability of the DFE, ^[25]). If $R_0 < 1$, then the disease-free equilibrium, $\Sigma_0$, is globally asymptotically stable in $\Omega$.

Theorem 3 (Global stability of the EE, ^[25]). The compact region $\Gamma$ defined by

is positively invariant under the flow induced by system $(2.1)$. It contains the disease-free equilibrium, $\Sigma_0$, and the endemic equilibrium, $\Sigma_+$, if $R_0 > 1$. Furthermore, if $R_0 > 1$, then the endemic equilibrium $\Sigma_+$ is globally asymptotically stable in $\Gamma$.

The Python code to generate Figure 1 is given below. We first import several Python scientific libraries: $\texttt{numpy}$ contains standard routines for numerical computations; $\texttt{scipy}$ contains the function $\texttt{odeint}$, which implements the Runge-Kutta method for integrating ODE's systems; $\texttt{matplotlib}$ contains useful functions for producing figures.

Next, we define the parameters for the SAIRP model (2.1).

Using expression (3.3), we easily compute the basic reproduction number $R_0$.

Afterwards, we define the SAIRP model given by system (2.1).

Finally, we integrate the SAIRP model (2.1) with several randomly chosen initial conditions and we produce the 3D Figure 1.

4.   Mathematical model with piecewise constant parameters

In this section, we consider the SAIRP model with piecewise constant parameters, proposed in ^[25], which allows to model the impact of public health policies and the human behavior in the dynamics of the COVID-19 epidemic.

For the sake of simplicity, the equations of the SAIRP model (2.1) can be rewritten as

with $x = (S, \, A, \, I, \, R, \, P)^T \in \mathbb{R}^5$ and $\alpha = (\Lambda, \, \mu, \, \beta, \, p, \, \theta, \, \phi, \, \omega, \, \nu, \, \delta)^T \in \mathbb{R}^9$, where the non-linear operator $f$ is defined in $\mathbb{R}^5 \times \mathbb{R}^9$ by

The time line $[0, \, T_{\rm{end}}]$ is subdivided into a finite number of $n+1$ intervals

with disjoint unions, and we introduce a piecewise constant function $\alpha$ defined on each time interval as

with $T_0 = 0$, $T_{n+1} = T_{\rm{end}}$ and $\alpha_i \in \mathbb{R}^9$. Next, consider the sequence of Cauchy problems defined for each initial condition $x_0 \in \Omega$ by

The following theorem establishes the well-posedness of system (4.3).

Theorem 4. For any initial condition $x_0 \in \Omega$, the sequence of Cauchy problems given by model (4.3) admits a unique global solution, denoted again by $x(t, \, x_0)$, whose components are non-negative. Furthermore, the region $\Omega$ is positively invariant ^[25].

We recall that the solutions of problem (4.3) are continuous on the time interval $[T_0, \, T_{\rm{end}}]$, but may not be of class $\mathscr{C}^1$ at $t = T_i$, $0 \leq i \leq n-1$. From the modeling point of view, each change of parameters occurring at time $t = T_i$ ($1 \leq i \leq n-1$) corresponds, for example, to a public announcement of confinement/lift of confinement or prohibition of displacement ^[25].

In what follows, we sub-divide the time interval $[0, 410]$ days into 9 sub-intervals and consider a set of piecewise parameters which are estimated in order to fit the real data of COVID-19 spread in Portugal, since the first confirmed case on March 2, 2020, until April 15, 2021.

Model fitting to real data with Python In order to fit the real data of active infected individuals by SARS-CoV-2 (detected by test) daily provided by the health authorities in Portugal ^[31], we use the model (4.3) with piecewise constant parameters and use the programming language Python (version 3.8) ^[20], and the Python Libraries: Numpy (version 1.18.5), Pandas (version 0.25.3), Scipy (version 1.4.1) and Matplotlib (version 2.0.0). The real data are available in ^[31] or, for example, in the data repository site, see e.g. ^[32]. In what follows we explain the goal of each code block.

First we need as previously to import several $\texttt{Python}$ scientific libraries.

We consider the real data from ^[31,32], here denominated data-pt-covid19 with extension .xlsx. To read the real data one can use the following code.

Next, we subdivide the time window of the total 410 days into 9 subintervals. This specific subdivision is related to the increase/decrease of the transmission of the virus in the community, the public health measures implemented by the Portuguese authorities and the human behavior, that changes over time.

Next, we draw the curve of active infected individuals with COVID-19 in Portugal (real data), from the first confirmed case day, March 2, 2020, until April 15, 2021.

In order to solve the model (4.1) we need to define the initial conditions, for $t = 0$, that correspond to the number of individuals in each class on March 2, 2020.

Then, we define the values of the parameters that take constant values for all time $t \in [0, 410]$.

The parameters $\beta$, $p$ and $m$ are assumed to be piecewise constant. Moreover, we estimate the parameter $\beta$ for all $t \in [0, 410]$ days and the parameters $m$ and $p$ in some of the sub-intervals, using $\texttt{scipy.optimize.curve_fit}$ from Python, see ^[33] for more details, that uses non-linear least squares method to fit a function to data.

Finally, we make the plot with the real data and the model solution $I(t)$, for $t \in [0, 410]$.

The graphic that results from the previous $\texttt{Python}$ code is given in Figure 2.

We recall, that the mathematical method and numerical code can be applied to other models given by systems of ODE's that describe the transmission dynamics of different virus or bacteria in human, animals or cells, for example.

5.   Conclusions

In this note, we provided a complete toolkit for performing both a symbolic and numerical analysis of a recent epidemic model, presented in ^[25], introduced in order to study the spreading of the COVID-19 pandemic.

This innovative compartmental model takes into account the possible transmission of the virus by asymptomatic individuals, as well as the possibility to protect a fraction of the affected population by public health strategies, such as confinement or quarantine. Since the public health policies have a strong influence on the spreading of the epidemic, we have also considered a piecewise constant parameters extension of the initial autonomous system, which have proved its ability to fit with real data. Once the mathematical analysis of such a complex compartmental model can be tedious, we have presented a computational approach, which is intended to support the theoretical analysis:

● using the free and open-source software Sagemath, we have first computed the symbolic expressions of the basic reproduction number $R_0$ (Eq (3.3)), of the disease-free equilibrium $\Sigma_0$ (Eq (3.1)) and of the endemic equilibrium $\Sigma_+$ (Eq (3.2)), which highlight the role of each parameter of the model;

● using the scientific libraries numpy, $\texttt{scipy}$ and $\texttt{matplotlib}$ of the free and open-source language Python, we have integrated the epidemic model, in order to illustrate the theoretical stability statements (see Figure 1);

● finally, we fitted the piecewise constant parameters extension of the initial model with recent real-world data (see Figure 2).

Overall, the programs presented in this note have been written in a sufficiently general manner, so that they can easily be adapted to a great number of other epidemic models.

In a near future, we aim to apply our computational approach to an improved version of our complex compartmental model, so as to consider the effects of the mutations of the virus and the benefits of vaccination.

Acknowledgments

This research is partially supported by the Portuguese Foundation for Science and Technology (FCT) by the project UIDB/04106/2020 (CIDMA). Cristiana J. Silva is also supported by FCT via the FCT Researcher Program CEEC Individual 2018 with reference CEECIND/00564/2018.

Conflict of interest

The authors declare that they have no conflict of interest.