The Multi-Compartment SI(RD) Model with Regime Switching: An Application to COVID-19 Pandemic

Abstract

We study—with existence and unicity results—a variant of the SIR model for an infectious disease incorporating both the possibility of a death outcome—in a short period of time—and a regime switch that can account for the mitigation measures used to control the spreading of the infections, such as a total lockdown. This model is parametrised by three parameters: the basic reproduction number, the mortality rate of the infected, and the duration of the disease. We discuss a particular example of application to Portuguese COVID-19 data in two short periods just after the start of the epidemic in 4 March 2020, with the first two cases dated that day. We propose a simple and effective method for the estimation of the main parameters of the disease, namely, the basic reproduction number and the mortality rate of the infected. We correct these estimated values to take into account the asymptomatic non-diagnosed members of the population. We compare the outcome of the model in the cases of the existence, or not, of a regime switch, and under three different scenarios, with a remarkable agreement between model and data deaths in the case of our basis scenario. In a final short remark, we deal with the existence of symmetries for the proposed model.

1. Introduction

In this paper, we study a variant of the SIR model for an infectious disease, incorporating the possibility of a death outcome, and focus on a short period of time in order to describe the difference in the evolution of subpopulations with, or without, mitigating measures such as a lockdown. The SIR model—SIR for Susceptible, Infected, Recovered—was introduced in [1] (see also [2] for a present day account) and still is at the root of the most used models for epidemics. It is similar to the system of ordinary differential equations given in Formula (1), the difference being that in the SIR model there is no equation for the proportion of those in the population that have died once infected; that is, the equation for the variable d. The introduction of a death outcome in a SIR model has not been a widely studied subject, maybe because of the difficulties in the biological interpretation of parameters or parameter estimation (as is clearly stated in [3] (p. 34)). In the specific case of the data analysis presented here, it is possible to satisfactorily estimate a death rate parameter, and the biological interpretation is simple, within the bounds of the following caveat that pervades the data set we use: we do not distinguish deaths classified under the two categories—death caused by the infection and death caused by some other disease than COVID-19 on a person with a positive diagnostic of the infection. In the real data application we perform, our goal is to analyse a disease in a period of less than 70 days, and so the demographic fluctuations are not significative, since the lethality rate is known to be small, at this date; moreover, as it follows from [1] (p. 110): “Since the course of an epidemic is short compared with the life of an individual, the population may be considered as remaining constant, except in as far as it is modified by deaths due to the epidemic disease itself”. Nevertheless, if the disease has a rapid progression with an outcome of possible death, the death toll of the disease has a socially significant impact. Moreover, the death toll of such an epidemic can be determined daily; as such, if the model has, as one of its computed results, the deaths of the infected in the population, an assessment of the quality of the model can be performed; we attempt such an assessment in Remark 11. In this work, we need only the most elementary concepts on mathematical epidemic modelling. A very thorough introduction to epidemic modelling can be read in [2]. Some other introductory expositions of the mathematics of epidemics are, for instance, [3] (pp. 15–53) or [4] (pp. 345–410).

The modelling of the current epidemic has been the subject of many studies, and it is impossible to mention even a significant part; we chose to recall some references that we find to be closely related to our approach. For instance, in [5], the authors consider eight compartments, and they estimate the parameters using data between day 1 and day 46; the details of the estimation procedures are not given. In [6], the authors use—among other models—a compartment model of SEIR type, with an additional group of exposed patients for which there is a delay between being exposed and being infected; they refer to the use of maximum likelihood estimation with a Poisson likelihood as an estimation procedure. The work [7] also considers a SIR type model but with natural deaths, and considers explicitly a discretisation procedure to solve the differential equations. Let us stress the following author’s very important remark that justifies the care needed to consider prospective results of these kinds of models. “…Predicting the behaviour of the epidemic with high accuracy is nearly impossible due to there being many unknown factors, e.g., …and parameters of the epidemic”. We must point out that different models have different advantages and weaknesses. For instance, in [8], a review of 72 articles on the mathematical models of the COVID-19 pandemic leads the authors to the conclusion that, although epidemiological models take into account the dynamics of the spread of the infectious disease in a population, they are also highly dependent on the estimation of parameters, and do not take into account the variability of parameters during the course of the epidemic. This observation of May 2021 shows the opportunity of studying epidemiological models with regime switching, which is a major contribution of the present work.

The introduction of birth and death in compartmental SIR type models (see, for instance [9] (p. 30)) aims at describing the evolution of a disease in a long term scenario, of several years, in which demographic fluctuations have to be taken into account. We must presently stress the following caveat for these models. Relying on SIR type epidemiological models to produce predictions on the course of an epidemic can be deceiving, as the work [10] clearly shows; in it, the authors consider a SIR model, with scenarios defined by parameters adjusted for qualitative information—such as the real number of infected—and then compare the evolution given by the model with the observed data. The main conclusion is that: ...The constructed SIR model failed to forecast the actual spread and pattern of the epidemic in the long term. The study of SIR type compartmental models with a subpopulation of the deaths caused by infection, which is to be applied in shorter periods for which the mortality and birth rates of the population can be considered as irrelevant, has already been considered by several authors.

In [11], not related to the COVID-19 pandemic, there is a thorough study—with equilibrium points and stability—of a SEIR model with two additional compartments—asymptomatic infectious and dead once infected. The model has at least eleven primary parameters and a simulation study with parameters obtained from real data from a study of Ebola disease. The qualitative analysis of a system of ODE of the SIR type model can shed light on the characteristic features of the evolution of the variables, such as susceptible, infected, and recovered. Such a study is developed in [12] for a SIR model, in which the proportion of deaths caused by disease is not a variable and, having as additional parameters, the protection rate and the isolation rate; these parameters being non-estimable; due to the lack of data, the model is exploited with plausible given values for the parameters. The models proposed in [13] consider, besides the usual subpopulations in a SEIR type model “...the number of individuals that pass away due to the disease, divided by the total population”. The author considers that all variables are measurable with the exception of the exposed, and considers eight parameters for the models. The estimation procedure considered is of constrained weighted least squares type and simulation results are presented for different sets of parameters corresponding to typified disease cases such as recurrent epidemic and short lockdown period. A remarkable modelling work on the current COVID-19 pandemic is presented in [14]. The authors consider an ODE model with ten variables—susceptible, exposed, asymptomatic, quarantined asymptomatic, severely infected, hospitalised, recovered, recovered asymptomatic, deceased, and protective susceptible. They calibrate the parameters employing the least square method to the proposed model to obtain the best-fit parameters for Ghana and Egypt; they perform a sensitivity analysis and study the optimal control using two control variables that are incorporated into the ODE initial system to reduce the disease transmission and the quarantine for asymptomatic infected individuals. The work in [15] aims at describing the COVID-19 epidemics in Wuhan, Toronto, and Italy. It is a network model, as the “...individuals (are classified) by their average number of contacts in a typical day (time unit for the modelling)”, and to each of these numbers, there corresponds a set of classes of a SEIR type model with additional subpopulations of asymptomatically infected, symptomatically infected, hospitalised and dead. The main results are the theoretical derivations of the “...epidemic threshold defined as the average number of secondary cases produced by an infected individual in a completely susceptible population and of the final epidemic size”. These results are exploited in a simulation study for which the authors consider phases such as before and after lockdown. In [16], the author considers an alternative for the regime switching—of rough type—considered here, but for a SEIR type model in which the coefficients are supposed to take two (or three values) on two (or three) complementary intervals of the time domain. A study of a regime switching models for the evolution of an infected population that influenced this work is [17]. A recent study of a SI(RD) model with regime switching with a logistic function is presented in [18]; the authors use a neural network to detect the regime change points, and a grid search for the estimation of the usual parameters. A more recently published study of a SIR model with deaths and with changing parameters applied to Spanish data is given in [19]; the SI(RD) model studied in this article was, independently, being discussed by us at least since November 2020 (see [20,21]). In this work, we consider data from the first wave of the COVID-19 pandemic in Portugal. Other very recent works that also consider data from Portugal are [22,23].

We now state the main results of this work.

• We study the model proposed—a susceptible-infected-(recovered-dead) model with regime switching—proving that there are two ways to achieve the regime switching (the rough way and the smooth way), and that these are equivalent for all practical purposes;
• Given a bivariate set of data—(daily cases-lethality)—we propose a new simple method to estimate the parameters of the model specially suited to the jagged character of the observations using a binomial model; this amounts to hypothesise that there is a mechanism that corrupts, by noise, the observed trajectories of the ordinary differential equations of the model.
• We show that, taking into account important qualitative information, there are secondary parameters—that modify the primary parameters—that can be set according to reasonable scenarios, such that in one of these scenarios, the lethality rate can be recovered with a very small error.

2. The SI(RD) Model

In this section, we detail the main concepts and specify the notations used next. Some of the concepts presented in the following are classical concepts in the infectious disease modelling literature. Variants of the SIR model were considered by one of the authors in [24] and, for the COVID-19 epidemic context, in [25,26]. In the following, we will focus on the main differences between previous models and the one we study in this work. In the SI(RD) model—which is a variant of the SIR model with an additional compartment represented by D in the acronym SI(RD) for the deceased elements of the population—the members of a population are categorised into one of four groups: S, those who are susceptible to being infected; I, those who have been infected and are able to spread the disease to susceptible individuals; R, those who have recovered from the disease and are immune to subsequent re-infection and D, those who died while infected. The movement of individuals is one-way only according to the scheme: $S\to I\to (R,D)$. The parameters of the model are:

• $\alpha$ the daily infection rate of the susceptible population;
• $\beta$ the recovery rate is a rate which controls how fast members progress into the I and R groups, respectively;
• $\mu$ is the death rate among those infected;
• $\gamma =\alpha /\beta$ is a composite parameter, which—in case $\mu =0$—is often used, and is referred to as the contact number.

The SI(RD) model is described by a system of ordinary differential equations (ODE) for the functions $s,i,r,d$—given in formulas (1)—depicting the time evolution of the proportions of the population in each group, that is, s in S, i in I, r in R and d in D; as so, $s,i,r,d$ refer, respectively, to the fraction of the population in the susceptible, infected, recovered and dead subpopulation groups. The ODE system in the Cauchy problem of Formula (1) can also be found in [19] (p. 2).

$$\{\begin{array}{l}{s}^{\prime }\left(t\right)=-\alpha i\left(t\right)s\left(t\right),s\left(t_1\right)=s_0\\ i^{\prime }\left(t\right)=i\left(t\right)(-\beta -\mu +\alpha s\left(t\right)),i\left(t_1\right)=i_0\\ r^{\prime }\left(t\right)=\beta i\left(t\right),r\left(t_1\right)=r_0\\ d^{\prime }\left(t\right)=\mu i\left(t\right),d\left(t_1\right)=d_0.\end{array}$$

(1)

When initial conditions for these groups—susceptible, infected, recovered and dead—are specified, the change in size of these groups may be plotted over time by solving numerically the system of four ordinary differential equations which, for constant or regular time varying parameters, has a unique regular solution.

The inspection of these equations will reveal details of the dynamics of the disease in a population (see, for a related approach of this discussion, [27] (pp. 44, 45)).

• If the fraction of the population in the infected group is initially increasing (i.e., $i^{\prime }\left(0\right)=i_0>0$), an epidemic has begun. For any $t>0$ in a neighbourhood of zero—such that $i\left(t\right)>0$—and if almost all population is susceptible—that is, $s_0\approx 1$—we have that:

  $$i^{\prime }\left(t\right)>0\iff i\left(t\right)(-\beta -\mu +\alpha s\left(t\right))>0\iff \frac{\alpha }{\beta +\mu }>\frac{1}{s\left(t\right)}\approx \frac{1}{s_0}\approx 1.$$

  (2)
  As a consequence, for the SI(RD) model, we define the composite parameter $\tilde{\gamma }$, corresponding to $\gamma :=\alpha /\beta$ in the SIR model—since in SIR model, there is no $\mu$ parameter—by:

  $$\tilde{\gamma }=\frac{\alpha }{\beta +\mu }.$$
• At the peak point of an epidemic—with $i\left(t\right)\ne 0$—we should have that $i^{\prime }\left(t\right)=0$, by Formula (2), and so, the peak point of an epidemic occurs at a time $t_{\pi }$ such that:

  $$s\left(t_{\pi }\right)=\frac{\beta +\mu }{\alpha }=\frac{1}{\tilde{\gamma }}.$$

Remark 1

(Basic reproduction number).The basic reproduction number ${\mathcal{R}}_0$ is the expected number of secondary cases produced by a single infected individual in a completely susceptible population. Usually, it is defined as,

$${\mathcal{R}}_0=T·C_a·D_u,$$

with:

• T—the transmissibility which is the probability an individual infecting another given there was contact between them;
• $C_a$—the average rate of contact between susceptible and infected individuals;
• $D_u$—the duration of the infection in individuals; in the SI(RD) model, the infection ends with two possible outcomes: either recovered or dead.

We may identify $\alpha =T·C_a$ and $D_u={(\beta +\mu )}^{-1}$, and doing so, we have the identification:

$${\mathcal{R}}_0=\frac{\alpha }{\beta +\mu }=\tilde{\gamma }.$$

(3)

The contact number $\tilde{\gamma }$—alias, the basic reproduction number—may be understood as the average number of susceptible members of the population; an infected individual spreads the disease to, while that individual is in the infected group.

Remark 2

Let us notice that, at some time $t_1$, we have $i^{\prime }\left(t_1\right)>0$ if, and only if, with $i\left(t_1\right)>0$ we have,

$$s\left(t_1\right)·\frac{\alpha }{\beta +\mu }=s\left(t_1\right)·\tilde{\gamma }>1.$$

(4)

Remark 3

(On the positive and bounded solutions of ODE system in Formula (1)). We observe that the sum of the righthand terms of the four equations in Formula (1) is zero. As a consequence of this observation, we have that, for all t and by integration, the sum $s\left(t\right)+i\left(t\right)+r\left(t\right)+d\left(t\right)$ is constant; since at time $t_1$, we want to deal with proportions we have that $s_0+i_0+r_0+d_0=1$ and so, we may conclude that for all t, we have that $s\left(t\right)+i\left(t\right)+r\left(t\right)+d\left(t\right)=1$. Furthermore, if the solutions are all non-negative, the relation $s\left(t\right)+i\left(t\right)+r\left(t\right)+d\left(t\right)=1$, valid for all t, this implies that all the solutions are bounded. The proof to show that the solutions are bounded for the SIR type models may be done in several ways; the most efective proceeds by considering one equation at a time, by expressing the equation in integral form in such a way that an exponential function appears, guaranteeing the positivity of solutions (see for instance [28]). By dealing with one equation after another, it is possible to take into account relations such as the one observed in Formula (4) above.

3. SI(RD) Models with Regime Switching (Two Different Values for the Parameter ${\mathcal{R}}_{\mathbf{0}}$)

For the modelling purposes of our work, we are interested in considering that the parameters of the SI(RD) model, namely, the basic reproduction number ${\mathcal{R}}_0$ alias $\tilde{\gamma }$, the duration of the infection $D_u$ alias $\delta$ and the mortality rate $\mu$—that have to be estimated or calibrated from the observations—may vary in disjoint and complementary intervals of time. So, firstly, we are going to select the main parameters of the model. For that, we first observe that the parameters $\alpha$ and $\beta$ can be expressed as functions of the primary parameters ${\mathcal{R}}_0$ alias $\tilde{\gamma }$, $D_u$ alias $\delta$ and $\mu$. By using Formula (3) and the fact that $\alpha =T·C_a$ and $D_u={(\beta +\mu )}^{-1}$, we have that:

$$\alpha (\tilde{\gamma },\delta ):=\frac{\tilde{\gamma }}{\delta }\mathrm{and} \beta (\delta ,\mu ):=\frac{1}{\delta }-\mu ,$$

and so, the system of ODE in (1) can be rewritten as:

$$\{\begin{array}{l}{s}^{\prime }\left(t\right)=-\alpha (\tilde{\gamma },\delta )i\left(t\right)s\left(t\right),s\left(t_1\right)=s_0\\ i^{\prime }\left(t\right)=i\left(t\right)(-\beta (\delta ,\mu )-\mu +\alpha (\tilde{\gamma },\delta )s\left(t\right)),i\left(t_1\right)=i_0\\ r^{\prime }\left(t\right)=\beta (\delta ,\mu )i\left(t\right),r\left(t_1\right)=r_0\\ d^{\prime }\left(t\right)=\mu i\left(t\right),d\left(t_1\right)=d_0,\end{array}$$

(5)

a formulation that only uses the parameters $\tilde{\gamma }$, the contact number, $\delta$, the duration of the disease and $\mu$, the mortality rate among the infected.

We now present the main ideas for the study of the SI(RD) model that we have been considering, but with regime switching. Consider a phenomena modelled by a system of ODE—such as the one in Formulas (5)—denoted for simplicity by $y^{\prime }\left(t\right)=F(t,y\left(t\right),{\theta }_1)$ with ${\theta }_1\in \Theta \in {\mathbb{R}}^d$ a chosen fixed parameter and $t\in [0,T]$ and $F:[0,T]\times {\mathbb{R}}^4\times \Theta ⟶{\mathbb{R}}^4$ a regular function of the first two variables say, twice differentiable. The change in the dynamics of phenomena may be rendered, in the ODE of the model, by a change in the right-hand side of the ODE, that is, by considering $y^{\prime }\left(t\right)=H(t,y\left(t\right),{\theta }_1)$ with a function H different from the function F for, let us say as an example, $t\in [T_a,T_b]\subset [0,T]$; this may be called a pure functional regime switching. Another basic rendition amounts to just change the parameter, that is considering $y^{\prime }\left(t\right)=F(t,y\left(t\right),{\theta }_2)$ for $t\in [T_a,T_b]\subset [0,T]$ and ${\theta }_2\in \Theta \in {\mathbb{R}}^d$ another chosen fixed parameter different from ${\theta }_1$; this type of regime switching may be called a pure parametric regime switching. A more general mixed regime switching can be considered by changing both the righthand side of the ODE and the parameter in some subintervals of the original time domain $[0,T]$. In this work we consider, and compare, two different forms of regime switching with the same purpose, to wit, the rough (R) regime switching and the smooth (S) regime switching. Let us consider a first example of pure parametric regime switching, the rough regime switching in the context of a standard existence and unicity theorem for ODE for instance, the global Cauchy-Lipschitz theorem (see [29] (pp. 152–154) and [30] (pp. 119–121)) or with the denomination Picard-Lindelöf theorem [31] (p. 8)).

Definition 1

(A rough regime switching for ODE). Suppose that the following conditions are verified for $\theta \in \{{\theta }_1,{\theta }_2\}\subset \Theta \subset {\mathbb{R}}^d$.

(j)

The function $F(t,y,\theta ):[0,T]\times {\mathbb{R}}^4⟶{\mathbb{R}}^4$ is a continuous function.

(jj)

For some constant $M=M\left(\theta \right)>0$, the uniform Lipschitz condition in the variable t given by,

$$∥F(t,x_1,\theta )-F(t,x_2,\theta )∥\le M∥x_1-x_2∥,$$

is verified for all values of the variables $(t,x)\in [0,T]\times {\mathbb{R}}^4$.

Let the first regime solution be $y_1\left(t\right)$ for $t\in [0,T_a]$, the unique solution of the ODE $y^{\prime }\left(t\right)=F(t,y\left(t\right),{\theta }_1)$ for $t\in [0,T_a]\subset [0,T]$ with a given set of initial conditions, which exists and is unique by the Cauchy–Lipschitz theorem. Let the second regime solution be $y_2\left(t\right)$ for $t\in [T_a,T]$ the unique solution of $y^{\prime }\left(t\right)=F(t,y\left(t\right),{\theta }_2)$ for $t\in [T_a,T]\subset [0,T]$ taking as initial conditions $y_1\left(T_a\right)$, that is, the terminal value of the solution in the first regime. Then, the continuous function defined by,

$$y_{1,2}\left(t\right)=\{\begin{array}{l}{y}_1\left(t\right),t\in [0,T_a]\\ y_2\left(t\right),t\in [T_a,T],\end{array}$$

(6)

is the a rough regime switching function for the ODE given by $y^{\prime }\left(t\right)=F(t,y\left(t\right),\theta )$ for $\theta \in \{{\theta }_1,{\theta }_2\}$.

We observe that $y_{1,2}\left(t\right):[0,T]⟶{\mathbb{R}}^4$ is obtained in $t\in [0,T]$ by glueing the solution $y_1\left(t\right)$ for $t\in [0,T_a]$ and the solution $y_2\left(t\right)$ for $t\in [T_a,T]$ and so $y_{1,2}\left(t\right)$ is the continuous function in $[0,T]$ that coincides with $y_1\left(t\right)$ for $t\in [0,T_a]$ and with $y_2\left(t\right)$ for $t\in [T_a,T]$.

Remark 4

(On the rough regime switching as a solution of an ODE). A question that arises naturally may be formulated in the following way: is the rough regime switching function for the ODE given by $y^{\prime }\left(t\right)=F(t,y\left(t\right),\theta )$ for $\theta \in \{{\theta }_1,{\theta }_2\}$ a bona fide solution of an ODE? The answer is no in the usual sense, because of the following. The rough regime switching here described by $y_{1,2}\left(t\right)$ is not differentiable at the point $T_a$ in general. In fact, supposing that $F(T_a,y_1\left(T_a\right),{\theta }_1)\ne F(T_a,y_1\left(T_a\right),{\theta }_2)$ and observing that by construction, under the hypothesis of the continuity of the function $F(t,y,\theta )$, we have that:

$$\begin{array}{cc}\underset{t\to T_a,t<T_a}{lim}{y}_{1,2}^{\prime }\left(t\right)& =\underset{t\to T_a,t<T_a}{lim}{y}_1^{\prime }\left(t\right)=F(T_a,y_1\left(T_a\right),{\theta }_1)\ne F(T_a,y_1\left(T_a\right),{\theta }_2)=\hfill \\ & =\underset{t\to T_a,t>T_a}{lim}{y}_2^{\prime }\left(t\right)=\underset{t\to T_a,t>T_a}{lim}{y}_{1,2}^{\prime }\left(t\right)\hfill \end{array}$$

despite the fact that $y_{1,2}\left(t\right)$ is differentiable in $[0,T_a[\cup ]T_a,T]$ as $y_1\left(t\right)$ is differentiable in $[0,T_a]$ and $y_2\left(t\right)$ is differentiable in $[T_a,T]$. Being so, the function $y_{1,2}\left(t\right)$ may not be considered as a usual solution of a Cauchy type problem. In Theorem 2 of Section 6, we consider extended solutions of an ODE in the sense of Carathéodory, and we show that the rough regime switching function described above can, nevertheless, be considered as a solution of an ODE, although only in an extended sense.

We now define a second form of regime switching that we will also use in this work, namely the smooth (S) regime switching.

Definition 2

(A smooth regime switching with transition function G). Let $t_0,t_1\in [0,T]$, with $t_0<t_1$, be some fixed values. Let the function $G(t,\lambda )$ defined for all $t\in [0,+\infty [$ and for every value of the parameter $\lambda \in \mathbb{R}$, be a twice differentiable function of its two variables. Suppose that there exists a value of the parameter ${\lambda }_{01}=\lambda (t_0,t_1)$ such that $G(0,{\lambda }_{01})=0$, $G(t,{\lambda }_{01})\approx 0$ for $t\in [0,t_0]\subset [0,T]$ and $G(t,{\lambda }_{01})\approx 1$ for $t\in [t_1,T]\subset [0,T]$. The smooth regime switching solution associated to the two Cauchy problems with ODE given by:

$$y^{\prime }\left(t\right)=F(t,y\left(t\right),\theta ),$$

(7)

for the parameters $\theta \in \{{\theta }_1,{\theta }_2\}\subset \Theta \in {\mathbb{R}}^d$ and with transition function G, is the solution of the correspondent ODE given by:

$$y^{\prime }\left(t\right)=F(t,y\left(t\right),{\theta }_1)(1-G(t,{\lambda }_{01}))+F(t,y\left(t\right),{\theta }_2)G(t,{\lambda }_{01}).$$

(8)

An example of a function G of Definition 2 is the function given by Formula (9) and depicted in Figure 1. In this example, the parameter $\lambda$ controls the extension of the regime changing interval $[t_0,t_1]\subset [0,T]$.

Remark 5

(On the existence of the smooth regime switching with transition function G as in Definition 2). It is clear that, due to the regularity of the function G, if the two Cauchy problems—for the values ${\theta }_1$ and ${\theta }_2$ of the parameters—given by the ODE in Formula (7) and arbitrary initial conditions have unique solutions, then any Cauchy problem given by the ODE in Formula (8) also has a unique solution.

This second form of the regime switching here considered, namely the smooth (S) regime switching, provides a solution for the problem of the change in the dynamics of the phenomena in our ODE model representation, given a choice of function G a smooth solution in $[0,T]$, and given a choice of fixed times $t_0,t_1\in [0,T]$.

In the following, we will consider a smooth regime switching model, as just described with a function such as the one in Figure 1, that is:

$$G(t,\nu ,\eta ,\theta ,\tau ):={\left(1-e^{(\eta -1)\theta (\nu +t)}\right)}^{\frac{1}{1-\tau }},$$

(9)

with a multidimensional parameter $\mathbf{\lambda }:=(\nu ,\eta ,\theta ,\tau )$ controlling the main characteristics of the function. Our use of the specific function G in Formula (9), in our modelling, stems from the possibility of performing a stability analysis similar to the one presented in [17]. With the following definitions: ${\mathbf{\alpha }}_1=\alpha ({\tilde{\gamma }}_1,{\delta }_1)$, ${\mathbf{\alpha }}_2=\alpha ({\tilde{\gamma }}_2,{\delta }_2)$, ${\mathbf{\beta }}_1:=-\beta ({\delta }_1,{\mu }_1)$, and ${\mathbf{\beta }}_2:=-\beta ({\delta }_2,{\mu }_2)$, the ODE that we consider is obtained from ODE in Formulas (5) by means of the procedure that gave the ODE in Formula (8), that is,

$$\{\begin{array}{l}{s}^{\prime }\left(t\right)=\left[-{\mathbf{\alpha }}_1(1-G(t,\mathbf{\lambda }))-{\mathbf{\alpha }}_{2}G(t,\mathbf{\lambda })\right]i\left(t\right)s\left(t\right),s\left(t_1\right)=s_0\\ i^{\prime }\left(t\right)=i\left(t\right)\left[(-{\mathbf{\beta }}_1-{\mu }_1+{\mathbf{\alpha }}_{1}s\left(t\right))(1-G(t,\mathbf{\lambda }))+(-{\mathbf{\beta }}_2-{\mu }_2+{\mathbf{\alpha }}_{2}s\left(t\right))G(t,\mathbf{\lambda })\right],i\left(t_1\right)=i_0\\ r^{\prime }\left(t\right)=\left[{\mathbf{\beta }}_1(1-G(t,\mathbf{\lambda }))+{\mathbf{\beta }}_{2}G(t,\mathbf{\lambda })\right]i\left(t\right),r\left(t_1\right)=r_0\\ d^{\prime }\left(t\right)=\left[{\mu }_1(1-G(t,\mathbf{\lambda }))+{\mu }_{2}G(t,\mathbf{\lambda })\right]i\left(t\right),d\left(t_1\right)=d_0.\end{array}$$

(10)

In order to obtain the results of the model, for the smooth regime switching case—results presented in Section 5—the system of ODE in Formula (10) is to be integrated numerically. The software used was Mathematica ™ (see [32]). The existence and uniqueness of solutions for such a model pose no problem, as pointed out in Remark 5.

4. Data and Parameter Estimation

The data considered in this study—represented in

Table 1. Data: 208 daily infection numbers; 4 March to 27 September 2020, Portugal.

$2,3,4,8,9,9,2,18,19,34,57,76,86,117,194,143,235,260,320,460,302,633,549,724,$
$902,792,446,1035,808,783,852,638,754,452,712,699,815,1516,515,598,349,514,$
$643,750,181,663,521,657,516,603,371,444,595,472,163,295,183,551,295,161,92,$
$242,178,480,533,553,138,175,98,234,219,187,264,227,226,173,223,228,252,288,$
$271,152,165,219,285,304,350,257,297,200,195,366,331,377,382,342,192,421,294,$
$310,270,283,227,346,300,336,417,375,377,292,259,345,367,311,451,323,457,266,$
$229,313,328,374,413,328,232,287,443,418,602,342,291,306,233,375,339,312,313,$
$246,135,127,252,229,313,263,209,135,111,203,255,204,238,153,106,112,167,213,$
$290,186,131,157,120,278,325,235,198,121,132,214,253,291,219,241,145,123,192,$
$362,399,401,374,320,244,231,390,418,406,486,315,249,388,646,585,687,497,673,$
$613,425,605,770,780,849,552,623,463,793,700,919,864,665.$

and

Table 2. Data: 208 death daily numbers; 4 March to 27 September 2020, Portugal.

$0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,2,3,6,3,8,10,10,17,16,24,19,21,20,27,22,37,20,29,16,34,35,$
$29,26,35,34,31,32,32,30,28,30,27,21,27,23,35,34,26,23,25,20,25,16,18,16,20,20,11,15,16,$
$9,12,9,9,19,12,9,6,13,15,13,16,16,14,12,13,14,14,12,14,13,14,13,14,14,12,11,8,10,9,5,6,$
$7,5,7,1,7,5,3,2,1,1,3,1,2,4,6,3,6,6,6,3,4,8,3,8,11,7,9,6,9,2,13,2,8,6,2,6,8,3,3,2,5,2,6,$
$5,3,7,4,1,2,3,3,2,8,2,1,0,1,1,3,3,4,6,3,2,3,6,2,3,3,1,5,2,2,4,2,2,5,4,2,2,6,3,1,3,2,3,2,$
$4,5,2,3,3,3,3,3,5,7,4,4,3,10,6,5,13,8,5,3,3,5,8,9$

and Figure 2—are composed of two daily reported time series; the series of the number of infected and the series of the number of deaths attributed to the infection.

The data used were extracted from a public available source transcribing information from the Portuguese General Direction of Health (DGS) and can be retrieved at https://pt.wikipedia.org/wiki/Predefinição:Tabela_Covid-19_Portugal_região/dia_(DGS) (the website was accessed on on 31 October 2020). The precise numbers do not really matter for our work, only the general trends matter; the methods introduced can be applied to other sets of data with the same characteristics. The number of infected is the number of individuals with positive outcome on an infection test; most commonly a PCR test.

As far as we know, with obvious exceptions, the reported deaths are of those persons who died having had a positive outcome on an infection test. There are natural questions arising from this data description; these questions do not have known answers because, for the moment, there is no practical and affordable way of getting the answers. What is the real number of infections in the whole population? What is the real lethality rate—alias, the case fatality rate—of the disease? This uncertainty is reflected on the outcomes of any model relying on parameters connected to these unknown numbers. In Section 4.5, we will consider the use of some secondary parameters to overcome this issue.

We now address the problem of estimating the parameters of the model from the data. We first observe, as it is usually observed in the usual graphic curves of epidemics evolution, e.g., Figure 2, that we may distinguish several periods the most remarkable being the initial growth period, starting 4 March 2020, with the first two cases being acknowledged and ending in 18 March (day 15)—with the emergency state declaration, effective 22 March—and the following period, ending in 9 May, day 67, as a consequence of the lockdown measures. We will concentrate our analysis in these first two periods as they correspond to a regime switch—derived from the emergency status applied from 22 March on (day 19)—for a determined set of viral clades; it is already known (see the genomic analysis of Instituto Ricardo Jorge in Portugal—https://insaflu.insa.pt/covid19/(accessed on 30 June 2020)—and, in France [33]) that the two following pairs of increasing-decreasing infection numbers corresponding to sets of different viral clads and strains—as shown in Figure 3—and so, they may require a completely different set of parameters; if not, a different model.

The identification of these periods can be done, for instance, with a moving average that smoothes the rough oscillations. In Figure 2, it can be observed that there is a clear maximum in the 12th day moving average of the infected, and there is a sort of plateau of maximum values for the death numbers. This allows us to define the starting and ending dates of the two periods under consideration for computational purposes.

From the data on the infected, we can estimate the parameter ${\mathcal{R}}_0\equiv \tilde{\gamma }$ and get some information on the duration of the infection, the $D_u\equiv \delta$ parameter. From the mortality data, we can estimate the mortality rate $\mu$ and get some additional information on the $\delta$ parameter.

4.1. On the Assumed Value for $\delta$ the Duration of the Disease

It is an acknowledged fact (see [34]) that the viral load for infections of this type of disease becomes undetectable after a few days, between 6 and 14. In the particular case of this disease, it was observed that, in the vast majority of cases, after 5 days (see [35]), the viral load becomes undetectable. As a consequence, we will assume that our basis scenario corresponds to $\delta =7.5$ days, and we will consider two extreme scenarios with parameter values $\delta =5$ days and $\delta =10$ days. Further on—see Section 4.4—we will see that the lethality rate estimates, with a three and eleven days delay, define an interval $[3,11]$ consistent with the proposed duration day interval for the disease $[5,10]$, as the main contagious period of the infection.

4.2. Determination of the Initial Date of the Decreasing Trend Period

Given a whole set of time labeled data points $i_1,i_2,\cdots ,i_n,\cdots$ of the number of infected the determination of the date where the regime switch may be taken to happen can be performed with a moving average sequence ${\overline{M}}_n\left(k\right)$, the average of the k day values prior to the current date, defined by:

$${\overline{M}}_n\left(k\right):=\frac{1}{k}\sum _{j=1}^k{i}_{n-k+j},\mathrm{for} n\ge k.$$

Intuitively, the sequence ${\left({\overline{M}}_n\left(k\right)\right)}_{n\ge k}$, for $k=1,2,\cdots$ is closer to ${\left(i_n\right)}_{n\ge k}$ the smaller the value of k. Moreover, for a sequence ${\left(i_n\right)}_{n\ge 1}$ with a rising trend, we should always have ${\overline{M}}_n\left(k\right)<i_{n+1}$, and so, the first time label n such that ${\overline{M}}_n\left(k\right)\ge i_{n+1}$ may be considered the time stamp of rising trend break. The same applies mutatis mutandis to a decreasing trend. For our data we know that the regime switch must have occurred between days 15 and 30 on account of the lockdown determination on 18 March 2020 (day 15). For $k=4,5,\cdots ,12$ we had $n=26$ for the first index n such that ${\overline{M}}_n\left(k\right)\ge i_{n+1}$ and so we take the date $n=26$ as the starting date of the decreasing trend in the data that will be used for estimation of the contact number.

4.3. Estimation of $\tilde{\gamma }$ the Contact Number

We propose the following method that is suited for the type of data considered with well defined growth characteristics: data defining a jagged line such that a regularisation obtained by a moving average that smoothes the rough oscillations—see Figure 2—is increasing in the first considered sub-period and decreasing in the second sub-period. To account for the jagged characteristic of the line corresponding to the data—a characteristic that is possibly due noise contamination originated from delays in collecting and reporting the observed numbers—we will consider a binomial type model, that is, the following hypothesis for each of the sub-periods of dates $\{1,\cdots ,67\}$.

(1)

The length of the period $\tau$ is determined and we have ${\left(i_n\right)}_{0\le n\le n_{\tau }}$ a sequence of daily observations of the number of infected. There is a sequence of random variables ${\left(I_n\right)}_{1\le n\le n_{\tau }}$ such that the initial $i_0$ is arbitrary but the remaining observed data ${\left(i_n\right)}_{1\le n\le n_{\tau }}$ is a realisation of this sample.

(2)

There exists a random variable $Z^{\tau }$, taking two values $u,d$—with u representing the magnitude of the upward jump and d representing the magnitude of a downward jump—such that $\mathbb{P}[Z^{\tau }=u]=p_u$ and $\mathbb{P}[Z^{\tau }=d]=p_d$ with $p_u+p_d=1$ and $p_u,p_d>0$. For ${\left(Z_n^{\tau }\right)}_{n\ge 1}$ a sample of $Z^{\tau }$ we have, for $n\ge 1$, that:

$$I_{n+1}=Z_{n+1}^{\tau }{I}_n,I_0\equiv i_0.$$

(11)

The values $u,p_u,d,p_d$ are the parameters to be estimated and they depend on the period considered.

As a consequence of these hypothesis we have the following immediate conclusions.

(1)

The maximum likelihood estimator of $p_u$, as a consequence of the distribution of $Z^{\tau }$, is given by:

$$\widehat{p_u}\left(N\right):=\frac{\#\left\{1\le n\le N:I_{n+1}>I_n\right\}}{N}.$$

(12)

(2)

By the law of large numbers we have that for $N\ge 1$ large enough,

$${\widehat{s}}_N:=\frac{1}{N}\sum _{n=1}^N\frac{i_{n+1}}{i_n}\approx {\mathbb{E}}^{\mathbb{P}}\left[Z^{\tau }\right],\widehat{s_N^2}:=\frac{1}{N-1}\sum _{n=1}^N{\left(\frac{i_{n+1}}{i_n}-{\widehat{s}}_N\right)}^2\approx {\mathbb{V}}^{\mathbb{P}}\left[Z^{\tau }\right],$$

(13)

that is, ${\widehat{s}}_N$ may be taken as a good approximation of ${\mathbb{E}}^{\mathbb{P}}\left[Z^{\tau }\right]$—the expectation of $Z^{\tau }$—and $\widehat{s_N^2}$ is a good approximation of ${\mathbb{V}}^{\mathbb{P}}\left[Z^{\tau }\right]$, the variance of $Z^{\tau }$.

(3)

By the method of moments, we can determine the couple of estimators $(\widehat{u},\widehat{d})$ for the parameter couple $(u,d)$ as the solutions of the equations derived from the law of $Z^{\tau }$ given by,

$$\left\{\begin{array}{c}{p}_{u}u+(1-p_u)d=\mathbb{E}\left[Z^{\tau }\right]\hfill \\ p_u{u}^2+(1-p_u)d^2-{(p_{u}u+(1-p_u)d)}^2=\mathbb{V}\left[Z^{\tau }\right],\hfill \end{array}\right.$$

considering the estimator $\widehat{p_u}\left(N\right)$ of $p_u$ given by Formula (12) and the estimator of $\mathbb{E}\left[Z^{\tau }\right]$ (respectively $\mathbb{V}\left[Z^{\tau }\right]$) given by Formula (13). There are two sets of solutions—Formulas (14) and (15)—and we have to choose the set that corresponds to the period characteristic, that is, in the first period we should have the greatest u with the highest probability and in the second period the greatest d with the higher probability.

$$\left\{\begin{array}{c}\widehat{u_1}:={\widehat{s}}_N-\frac{\sqrt{(1-\widehat{p_u}\left(N\right))\widehat{p_u}\left(N\right)\widehat{s_N^2}}}{\widehat{p_u}\left(N\right)}\hfill \\ \\ \widehat{d_1}:=\frac{{\widehat{s}}_N(1-\widehat{p_u}\left(N\right))+\sqrt{(1-\widehat{p_u}\left(N\right))\widehat{p_u}\left(N\right)\widehat{s_N^2}}}{1-\widehat{p_u}\left(N\right)}\hfill \end{array}\right.$$

(14)

$$\left\{\begin{array}{c}\widehat{u_2}:={\widehat{s}}_N+\frac{\sqrt{(1-\widehat{p_u}\left(N\right))\widehat{p_u}\left(N\right)\widehat{s_N^2}}}{\widehat{p_u}\left(N\right)}\hfill \\ \\ \widehat{d_2}:=\frac{{\widehat{s}}_N(1-\widehat{p_u}\left(N\right))-\sqrt{(1-\widehat{p_u}\left(N\right))\widehat{p_u}\left(N\right)\widehat{s_N^2}}}{1-\widehat{p_u}\left(N\right)}\hfill \end{array}\right.$$

(15)

Finally, we will take $\tilde{\gamma }=\widehat{u}$ for the first period and $\tilde{\gamma }=\widehat{d}$ for the second period, and the computed results are presented in

Table 3. Estimates of $p_u,p_d$ and of ${\mathcal{R}}_0$ for the two successive sub-periods.

Periods	$\widehat{{\mathit{p}}_{\mathit{u}}}\left(\mathit{N}\right)$	u	$\widehat{{\mathit{p}}_{\mathit{d}}}\left(\mathit{N}\right)$	d	$\widehat{\tilde{\mathit{\gamma }}}=\widehat{{\mathcal{R}}_0}$
$[1,25]$	0.833333	2.33474	0.166667	−2.04242	2.33474
$[26,67]$	0.487805	1.94555	0.512195	0.411811	0.411811

.

Remark 6

(Comments on the estimates of Table 3). In the first period, days 1 to 25, that is, 4 March to 28 March, the estimated value of $p_u$ is perfectly compatible with the strongly increasing character of the observations, and the estimated values of ${\mathcal{R}}_0$— highlighted in boldface—are similar to other estimations of this parameter; for instance ${\mathcal{R}}_0=2.6$ a value used in the Imperial College study paper, see [36] or [37] and ${\mathcal{R}}_0=2.28$ for a very particular subpopulation of a cruise ship. In the second period, days 25 to 67, that is, 28 March to 9 May, we observe that there is almost equal probability for either an upward or a downward probability and the downward parameter estimate gives us a lesser but still high estimate for ${\mathcal{R}}_0$.

4.4. Estimation of $\mu$ the Lethality Rate

The estimation of the lethality rate must be accomplished by establishing some functional relation with the number of infected. Among several different models tried, including ARMA models, we finally observed that the simplest model, that is the linear fit, was the adequate model for the two periods. We then considered, for the two periods independently, a linear fit with vanishing constant, of the daily deaths $d{D}_t$ given as a function of the daily infected $d{I}_t$ but with varying time delays. Based on both the lowest AIC—Akaike Information Criterion (see [38] for both AIC e BIC thorough developments)—and the highest R-squared the best model was a linear fit with vanishing constant with a delay of three days in the first period and a delay of eleven days in the second period. That is,

$$d{D}_{t+\Delta }=\mu d{I}_t+{ϵ}_t,t\ge 1,\Delta =3,11.$$

(16)

We recall that the p-value in

Table 4. Estimates $\widehat{\mu }$ of the lethality $\mu$ for the two successive sub-periods.

Periods	$\mathbf{\Delta }$	$\widehat{\mathit{\mu }}$	St. Error	t-Stat	p-Value	Adj. ${\mathit{R}}^2$	BIC	AIC
$[1,25]$	3	0.0306118	0.00155902	19.6353	$2.71833\times {10}^{-16}$	0.938957	123.008	120.57
$[26,67]$	11	0.0350414	0.00210826	16.621	$8.02743\times {10}^{-20}$	0.867616	302.68	299.21

is the probability of observing a t-statistic at least as far from zero as the one obtained from the data used in the estimation; given that the p-value in the table the model is adequate.

Remark 7

(Comments on the estimates of Table 4). Despite the remarkable statistical properties of the estimates obtained in the two periods under analysis, the values obtained for the lethality rate $\widehat{\mu }$—highlighted in boldface—are excessive if we consider the whole population of infected including also the undiagnosed and/or asymptomatic cases. In the following Section 4.5 we put forward an explanation and a correction to this excess of deaths.

Remark 8

(On the nature of the most adequate ODE model). The lethality model we obtained as the best fit in Formula (16) suggests that the ODE model should comprise a delay in the differential equation corresponding to the lethality. The use of delay differential equations (DDE) in epidemiology is not new, see for instance [39]. Also, in [40] there is an application to the multi-group SIR model with mass action incidence and with discrete time delay based on the construction of a Lyapunov functional for the DDE model. As it is well known, the DDE models require the study of stability that can be quite demanding (see, for instance, [41] (p. 36)).

4.5. Further Hypothesis on the Results of the Estimation Procedure

As early as 8 May 2020, it was a common assumption that the total number of infections was much higher than the officially reported (see [42]). In fact, already in [43], it is stated that the total number of infected was between 3 and 20 times the number of reported infection cases. Moreover, recently in 5 October 2020, Dr. Michael Ryan, head of emergencies at the WHO said, that… “best estimates indicate that roughly 1 in 10 people worldwide may have been infected by the coronavirus—more than 20 times the number of confirmed cases” (see https://www.cnbc.com/2020/10/05/who-10percent-of-worlds-people-may-have-been-infected-with-virus-.html (accessed on 6 October 2020)). As a consequence, since our model applies to the whole population, we will further assume the following:

• We will consider that the real number of infected is 20 times higher than the number reported. As a consequence, the model parameters $\tilde{\gamma }$ and $\mu$ are to be modified to take this assumption in consideration.
• The estimated value for the parameter $\tilde{\gamma }$ has to be replaced by the value ${\tilde{\gamma }}_c$ given by,

  $${\tilde{\gamma }}_c=\tilde{\gamma }\frac{1+19\zeta }{20},$$

  with $\zeta \in ]0,1[$, a secondary parameter, being the correction needed to take into account that the non diagnosed infected cases must correspond to non symptomatic cases and, as such, with a contagion force inferior to those infected with symptoms and tested with a positive outcome. Having no prior information available, we will consider that $\zeta =0.5$ correspond to our basis scenario, and we will base our main conclusions on the observations of the outcomes of the model under this particular scenario.
• In the same way, the estimated value for the parameter $\mu$ should be replaced by the value ${\mu }_c$ given by,

  $${\mu }_c=\mu \frac{1+19\theta }{20},$$

  (17)

  with $\theta \in ]0,1[$, another secondary parameter, being the correction needed to take into account the mortality rate due to the infection in the asymptomatic and not tested with a positive outcome. Having no prior information available, we will suppose that $\theta =0$ correspond to our basis scenario and $\theta =0.25$ and $\theta =0.5$ to alternative scenarios.
• For the initially immune $I_i$ we consider that in the first period these are supposed to be 2.5% of the population. In the second period the initially immune are those who either recovered or died at the end of the preceding period.
• For the smooth regime switching model drive function G, shown in Figure 1, the parameter choice ${\mathbf{\lambda }}_0:=(0,\eta ,\theta ,\tau )=(0,0.85,1.39,0.994)$ ensures that for $t_0=16$ we have $G(t_0,{\mathbf{\lambda }}_0)=0.00238659\approx 0$, for $t_1=22$ we have $G(t_1,{\mathbf{\lambda }}_0)=0.181604$, for $t_2=26$ we gave $G(t_2,{\mathbf{\lambda }}_0)=0.477708\approx 0.5$ and that for $t_3=40$ we gave $G(t_3,{\mathbf{\lambda }}_0)=0.960981\approx 1$ and so, the smooth regime switching unfolds, essentially, between day 16, the declaration of starting time of the lockdown, and day 40.

5. Results of the Model

We now apply the model given by the ODE systems in Formula (5) and in Formula (10) to the two different sub-periods considering the different set of parameters indicated above. For the second sub-period and the ODE systems in Formula (5), the initial conditions are the terminal conditions of the preceding sub-period. In the numerical integration of the rough regime switching ODE we take day 20 that is, 23 March the following day to the emergency state declaration, in Portugal, to be the beginning of the second sub-period with a different regime.

Let us give a summary of the assumptions for the basis scenario by interpreting the parameters in

Table 5. Parameters for basis scenario and extreme scenarios I and II.

Parameters	Extreme Scenario I	Basis Scenario	Extreme Scenario II
$\zeta$	0.3	0.5	0.7
$\theta$	0.25	0	0.5
$\delta$	5	7.5	10

. In this scenario $\zeta =0.5$, that is, those infected but not diagnosed as such have a contagion force ${\gamma }_c=0.525\gamma$ with $\gamma$ being the contagion force of those with an positive diagnostic infection; $\theta =0$ and so ${\mu }_c=\mu /20$ meaning that mortality rate for the asymptomatic infected and not tested with a positive outcome is one twentieth of the mortality for those infected and with a positive outcome on a test; and, finally, $\delta =7.5$ states that the duration of the infection in the basis scenario is $7.5$ days. We will first consider the results for the basis scenario and after, the results for the extreme scenarios I and II defined for specific choices of the secondary parameters $\zeta$, $\theta$ and $\delta$ as in Table 5 (see Section 4.5 for the need of establishing scenarios for the secondary parameters).

5.1. Results of the Model in the Basis Scenario

In Figure 4 we have the numerical solutions of the ODE systems for the two periods considered, for the constant parameter case—that is, with no regime switching—and for both the regime switching cases the smooth and the rough case in the case of the base scenario (see Table 5 for the values of the secondary parameters $\zeta$, $\theta$ and $\delta$ in the basis scenario). The yellow constant line at height one—explained in Remark 3—is just the sum of proportions and is drawn for control purposes. The terminal values, at date 67 corresponding to May 9th, for the proportions in each subpopulations—susceptible, infected, recovered and dead—are given in

Table 6. Numerical results of the ODE model of Figure 4 and Figure 5 at date 67.

Regime	$\mathit{\zeta }$	$\mathit{\theta }$	$\mathit{\delta }$	Susceptible	Infected	Recovered	Deaths
Constant	0.5	0	7.5	0.588081	0.0147572	0.368176	0.00398521
Smooth RS	0.5	0	7.5	0.71313	0.000537805	0.258525	0.00280666
Rough RS	0.5	0	7.5	0.746942	0.000271394	0.225357	0.00243014

.

Remark 9

(Comments on the ODE solutions displayed in Figure 4). With, susceptible in blue, infected in red, recovered in green, dead in magenta, and with the sum of proportions in yellow, the ODE model without regimes, on the left-hand side, and the ODE models with regimes are clearly different. In the case of the rough regime switching the effect of the regime change is clear both in the susceptible and the infected. At first sight, the final values of the proportions between the rough regime switching and the smooth regime switching; on the right-hand side, if they differ, they differ slightly (see Table 6).

Remark 10

(Comments on the differences displayed in Figure 5 and on Table 6). As expected, the regime switch has the following effects: the increase of the susceptible, the sharp fall of the infected, the decline of the recovered and the most remarkable diminution of the deaths of the infected. In each instance, the change between the constant parameter evolution and the regime switching evolution was more pronounced for the rough regime switching, as expected.

Remark 11

(On the agreement between the death results of the model and the death numbers of the data). The introduction of a death compartment in the model can be used to ascertain the quality of the model as well as of the parameter estimation procedures. We will now compare the results of the model with the real data. The lowest death rate—for the whole infected population corresponding to the rough regime switching—given by the model in Table 6, at day 67, is 2430.14 deaths per million; by reversing the procedure we applied to the estimated lethality rate in order to have a rate applied to the whole infected population in Formula (17) we have that the number of deaths of infected in regard to the diagnosed population is 2430.14/20 = 121.507 deaths per million. The data number of deaths caused by infection corresponding to the pandemic status at day 67, the end of the second period under analysis, should be taken at date 78 in Table 2, corresponding to 20 May 2020, due to the 11-day time lag between the deaths and the infected, and this value is 1263; in [44], we get, for 20 May, 1247 deaths because they did not count the number of deaths which occurred on 20 May, which was 37. For coherence with our computations, we will consider the death toll as 1263. Taking into account the estimate of the Portuguese population of 10,188,013 (see https://www.worldometers.info/world-population/portugal-population/ (accessed on 15 October 2020), we have the number of deaths of infected per million equal to 123.969. So, the ratio of the data number/model number is 123.969/121.507 = 1.0203, that is, an excess of 2% of the data over the model. This can be considered as a good agreement between model and data given that we used an estimate for the total population number in Portugal. We nevertheless observe that the true death toll of the epidemic may be larger than the accounted (see [45], for an early suggestion of this possibility) and that there is also the possibility that the model studied does not cover some important aspects of the epidemic evolution. This problem may be a feature common to the compartmental models; for instance in [5] (p. 857a) the authors predict a 0.06% percentage of deaths in a time horizon of 350 days—a very high value for a death of infected estimate of the whole population—which is a result comparable to the rough regime switching evolution for the extreme scenario I in

Table 7. Numerical results of the ODE model for extreme scenario I at date 67.

Regime	$\mathit{\zeta }$	$\mathit{\theta }$	$\mathit{\delta }$	Susceptible	Infected	Recovered	Deaths
Constant	0.3	0.25	5	0.801685	0.000413846	0.107069	0.0658329
e Smooth RS	0.3	0.25	5	0.815603	$6.82517\times {10}^{-6}$	0.0981462	0.0612436
Rough RS	0.3	0.25	5	0.822341	$3.1775\times {10}^{-6}$	0.0937749	0.058881

ahead.

5.2. Results of the Model in the Extreme Scenarios I and II

The consideration of alternative scenarios—with a different set of secondary parameters—allows for a better overview of the descriptive power of the model. As proposed in Section 4.5, we can consider several different values for the secondary parameters $\zeta$, $\theta$ and $\delta$, allowing the definition of two extreme scenarios in Table 5. The extreme scenario I has the lowest values for the $\zeta$ secondary parameter—and, as a consequence, also for the contact number $\tilde{\gamma }$— and also for the $\delta$ parameter, which is the duration of the disease. The extreme scenario II, with the numerical results in

Table 8. Numerical results of the ODE model for extreme scenario II at date 67.

Regime	$\mathit{\zeta }$	$\mathit{\theta }$	$\mathit{\delta }$	Susceptible	Infected	Recovered	Deaths
Constant	0.7	0.5	10	0.361501	0.0420978	0.0849807	0.48642
Smooth RS	0.7	0.5	10	0.604511	0.00413745	0.0279379	0.338413
Rough RS	0.7	0.5	10	0.67301	0.00232243	0.0214171	0.27825

, is the opposite.

Remark 12

(On the results of the model for the extreme scenarios I and II and consequent conclusions about the basis scenario assumptions). With respect to the basis scenario, the extreme scenario I has a larger death toll with an increasing proportion in the susceptible population while the proportions of infected and recovered are much smaller in the extreme scenario I when compared with the same proportions in the basis scenario. In the extreme scenario II (see Table 8), the proportion of susceptible is smaller but the proportion of recovered is comparable while the proportion of deaths is larger almost 114-fold for the regime switching model. As so and given the agreement between data and model numbers for the death toll in Remark 11, we must conclude that the assumptions made for the basis scenario are the most adequate. We stress that this fact gives us important information about the disease stated in the following conclusions. For the asymptomatic and/or undiagnosed the most plausible value of the contact number is 50% of the contact number of the diagnosed. The most plausible mortality rate, caused by the disease, for the asymptomatic and/or undiagnosed is zero and the most plausible value for the duration of the disease is 7.5 days.

6. On the Regime Switching SI(RD) ODE Models: Existence and Unicity Results

An allegedly alternative to the rough regime switching considered in this work—that it seems being used by some authors (see [16], for instance)—is to consider that the coefficients of the SIR type model—for instance, SEIR or SI(RD) models—are not constant but take two (or three) different values in two (or three) complementary intervals of the time domain. In fact, we will observe in Remark 14 that the two approaches to the rough regime switching are essentially the same. For the alternative rough regime switching, the usual Lipschitz type existence theorem or even the Peano existence theorem for ODE do not apply as the function in the righthand term of the ODE—with such discontinuous time varying parameters—is not continuous. Existence and unicity of solutions for such a model is a crucial question. A different approach is needed; following [46] (pp. 41–44), we consider the definition of an extended solution of a differential equation,

Definition 3

(Extended solution of an ODE). Let a differential equation be given by,

$${\mathit{{\rm Y}}}^{\prime }\left(t\right)=f(t,\mathit{{\rm Y}}\left(t\right)),\mathit{{\rm Y}}\left(0\right)={\mathit{y}}_0\in {\mathbb{R}}^4,$$

(18)

or in the equivalent integral form with the appropriate Lebesgue measure $du$ in ${\mathbb{R}}^4$,

$$\mathit{{\rm Y}}\left(t\right)={\mathit{y}}_0+{\int }_0^{t}f(u,\mathit{{\rm Y}}\left(u\right))du,$$

(19)

for $f(t,\mathit{y}):I\times \mathcal{D}\to {\mathbb{R}}^4$ anonnecessarily continuous function, with $I\subset [0,+\infty [$ and $\mathcal{D}\subset {\mathbb{R}}^4$. An extended solution $\mathit{{\rm Y}}\left(t\right)$ of the ODE in Formula (18) is an absolutely continuous function $\mathit{{\rm Y}}\left(t\right)$ (see for a precise definition [47] (pp. 144–150)) such that $f(t,\mathit{{\rm Y}}(t\left)\right)\in \mathcal{D}$ for $t\in I$ and Formula (18), or equivalently, Formula (19), is verified for all $t\in I$, possibly with the exception of a set of null Lebesgue measure.

We now recall Caratheodory’s existence theorem (see [46] (p. 43)), in the context of the model we are studying, a theorem that ensures the existence of an extended solution under general conditions.

Theorem 1

(Caratheodory’s existence theorem). Suppose that $f(t,\mathit{y}):I\times \mathcal{D}\to {\mathbb{R}}^4$ verifies:

(i)

$f(t,\mathit{y})$ is measurable in the variable t, for fixed $\mathit{y}$, and continuous in the variable $\mathit{y}$, for fixed t.

(ii)

there exists a Lebesgue integrable function $m\left(t\right)$, defined on a neighbourhood of the initial time, let us say I, such that $\left|f(t,\mathit{y})\right|\le m\left(t\right)$ for $(t,\mathit{y})\in I\times \mathcal{D}$.

Then, for a given initial condition, given in a neighbourhood of the initial time, there exists anextended solutionaccording to Definition 3.

We now state and prove an important result for our work, showing that the rough regime switching function of Definition 1 is, essentially, an extended solution of the ODE of the SI(RD) model in Formula (1), in agreement with what was announced in Remark 4.

Theorem 2

(Existence and unicity in the SI(RD) model with rough regime switching). Consider the SI(RD) model given by the Cauchy problem n Formulas (1), that is:

$$\{\begin{array}{l}{s}^{\prime }\left(t\right)=-\alpha (\tilde{\gamma },\delta )i\left(t\right)s\left(t\right),s\left(t_1\right)=s_0\\ i^{\prime }\left(t\right)=i\left(t\right)(-\beta (\delta ,\mu )-\mu +\alpha (\tilde{\gamma },\delta )s\left(t\right)),i\left(t_1\right)=i_0\\ r^{\prime }\left(t\right)=\beta (\delta ,\mu )i\left(t\right),r\left(t_1\right)=r_0\\ d^{\prime }\left(t\right)=\mu i\left(t\right),d\left(t_1\right)=d_0,\end{array}$$

with a time dependent vectorial parameter $\theta =\theta \left(t\right)=(\tilde{\gamma },\delta ,\mu )$, such that for some $t_0\in ]0,T[$:

$$\theta \left(t\right)=\{\begin{array}{l}{\theta }_1,t\in [0,t_0[\\ {\theta }_2,t\in [t_0,T].\end{array}$$

Then, by Theorem 1, this Cauchy problem has an extended solution in the sense of Definition 3. Moreover, any two different solutions differ possibly only on a set of zero Lebesgue measure, and so the extended solution is unique.

Proof. 

In order to address the problem of the existence of extended solutions for the rough regime switching models, we observe that the ODE in Formulas (5)—with time changing, discontinuous, parameters $\tilde{\gamma }\left(t\right)$, $\delta \left(t\right)$ and $\mu \left(t\right)$—may be represented as:

$${\mathit{{\rm Y}}}^{\prime }\left(t\right)=\mathit{M}(\Theta \left(t\right),\mathit{{\rm Y}}\left(t\right))·\mathit{{\rm Y}}\left(t\right),$$

(20)

with $\mathit{{\rm Y}}\left(t\right):={(s\left(t\right),i\left(t\right),r\left(t\right),d\left(t\right))}^t$, $\Theta \left(t\right):=(\alpha (\tilde{\gamma }\left(t\right),\delta \left(t\right)),\beta (\delta \left(t\right),\mu \left(t\right)),\mu \left(t\right))$ and:

$$\begin{array}{cc}& \mathit{M}(\Theta (t),\mathit{{\rm Y}}(t\left)\right):=\hfill \\ & :=\left(\begin{array}{cccc}-\alpha (\tilde{\gamma }\left(t\right),\delta \left(t\right))i\left(t\right)& 0& 0& 0\\ 0& -\beta (\delta \left(t\right),\mu \left(t\right))-\mu \left(t\right)+\alpha (\tilde{\gamma }\left(t\right),\delta \left(t\right))s\left(t\right)& 0& 0\\ 0& \beta \left(\delta \right(t),\mu (t\left)\right)& 0& 0\\ 0& \mu \left(t\right)& 0& 0\end{array}\right).\hfill \end{array}$$

(21)

As a consequence, the ODE in Formula (20) can be represented as in Formula (18) with $f(t,\mathit{y}):[0,+\infty [\times {\mathbb{R}}^4\to {\mathbb{R}}^4$ such that $(t,\mathit{y})\mapsto \mathit{M}(\Theta (t),\mathit{y})·\mathit{y}$. We have, for the Euclidean norm $∥·∥$ in ${\mathbb{R}}^4$, with $m_{ij}(t,\mathit{y})$ the matrix components of $\mathit{M}(\Theta (t),\mathit{y})$ and using the fact that the Frobenius norm dominates the matrix norm associated with the Euclidean norm,

$$\begin{array}{cc}∥f(t,\mathit{y})∥& =∥\mathit{M}(\Theta (t),\mathit{y})·\mathit{y}∥\le ∥\mathit{M}(\Theta (t),\mathit{y})∥∥\mathit{y}∥\le {\left(\sum _{i,j=1}^4{\left|m_{ij}(t,\mathit{y})\right|}^2\right)}^{\frac{1}{2}}∥\mathit{y}∥\le \hfill \\ & \le N\left(t\right)L\left(∥\mathit{y}∥\right),\hfill \end{array}$$

(22)

with $L\left(s\right):=s$ and since $s\left(t\right),i\left(t\right)\le 1$, with N(t) given by the expression:

$${\left[{\alpha }^2(\tilde{\gamma }\left(t\right),\delta \left(t\right))+{\left(\beta (\delta \left(t\right),\mu \left(t\right))+\mu \left(t\right)+\alpha (\tilde{\gamma }\left(t\right),\delta \left(t\right)\right)}^2+{\beta }^2(\delta \left(t\right),\mu \left(t\right))+{\mu }^2\left(t\right)\right]}^{\frac{1}{2}}.$$

(23)

By Caratheodory’s existence theorem (see Theorem 1 above), the extended solution for equation given in Formula (18)—for a given initial condition and in the case of the rough regime switching models here considered—exists as a consequence of the estimate in Formula (22) by observing that the components of $\mathit{y}$ are proportions, and so $L\left(∥\mathit{y}∥\right)$ is bounded, and that $N\left(t\right)$ in Formula (23) is bounded and measurable, and so it is Lebesgue-integrable in any compact interval of the time variable.

We now address the unicity of solutions for the regime switching ODE problems already mentioned in this section. For that purpose, we will prove that the function $f(t,\mathit{y})$ satisfies a Lipschitz condition in the variable $\mathit{y}$ for each fixed t. By resorting to Formulas (5) and with $\beta :=\beta \left(\delta \right(t),\mu (t\left)\right)$, $\alpha :=\alpha (\tilde{\gamma }\left(t\right),\delta \left(t\right))$ and $\mu :=\mu \left(t\right)$, for notation simplification, we have that,

$$f(t,{\mathit{y}}_1)-f(t,{\mathit{y}}_2)=\{\begin{array}{l}-\alpha (i_1\left(t\right)s_1\left(t\right)-i_2\left(t\right)s_2\left(t\right))\\ (i_1\left(t\right)-i_2\left(t\right))(-\beta -\mu )+\alpha (i_1\left(t\right)s_1\left(t\right)-i_2\left(t\right)s_2\left(t\right))\\ \beta (i_1\left(t\right)-i_2\left(t\right))\\ \mu (i_1\left(t\right)-i_2\left(t\right)).\end{array}$$

(24)

We now use the fact that $i_1\left(t\right)s_1\left(t\right)-i_2\left(t\right)s_2\left(t\right)=i_1\left(t\right)(s_1\left(t\right)-s_2\left(t\right))+s_2\left(t\right)(i_1\left(t\right)-i_2\left(t\right))$ to obtain that,

$$f(t,{\mathit{y}}_1)-f(t,{\mathit{y}}_2)={\mathit{M}}_{\Delta }(\Theta \left(t\right),{\mathit{y}}_1,{\mathit{y}}_2)·({\mathit{y}}_1-{\mathit{y}}_2),$$

where, with the same simplified notations for $\alpha$, $\beta$ and $\mu$, we have:

$${\mathit{M}}_{\Delta }(\Theta \left(t\right),{\mathit{y}}_1,{\mathit{y}}_2):=\left(\begin{array}{cccc}-\alpha i_1\left(t\right)& -\alpha s_2\left(t\right)& 0& 0\\ \alpha i_1\left(t\right)& -\beta -\mu +s_2\left(t\right)& 0& 0\\ 0& \beta & 0& 0\\ 0& \mu & 0& 0\end{array}\right).$$

(25)

Now, using the Euclidean norm $∥·∥$ in ${\mathbb{R}}^4$, the fact that the Frobenius norm dominates the matrix norm associated with the Euclidean norm, the fact that all components of $\mathit{y}$ are non-negative and bounded by one, as proportions always are, we have that, if the parameters $\alpha$, $\beta$ and $\mu$ are bounded measurable functions, then $∥{\mathit{M}}_{\Delta }(\Theta \left(t\right),{\mathit{y}}_1,{\mathit{y}}_2)∥$ is bounded, by a bounded measurable function $k\left(t\right)$ given by:

$$k\left(t\right):={\left[3{\alpha }^2(\tilde{\gamma }\left(t\right),\delta \left(t\right))+{(\beta (\delta \left(t\right),\mu \left(t\right))+\mu \left(t\right)+1)}^2+{\beta }^2(\delta \left(t\right),\mu \left(t\right))+{\mu }^2\left(t\right)\right]}^{\frac{1}{2}}$$

and then we have:

$$∥f(t,{\mathit{y}}_1)-f(t,{\mathit{y}}_2)∥=∥{\mathit{M}}_{\Delta }(\Theta \left(t\right),{\mathit{y}}_1,{\mathit{y}}_2)∥∥{\mathit{y}}_1-{\mathit{y}}_2∥\le k\left(t\right)∥{\mathit{y}}_1-{\mathit{y}}_2∥,$$

(26)

thus ensuring a boundedness condition that implies the Lipschitz character of $f(t,\mathit{y})$. Now, either directly using theorem 18.4.13 in [48] (p. 337) or using Osgood’s uniqueness theorem—as presented for instance, in [49] (p. 58) or in [50] (pp. 149–151)—we may conclude that the extended solution, that we know to exist, is unique, in the sense that two solutions may only differ on a set of Lebesgue measures equal to zero. □

Remark 13

(Existence alternative proof). We could also quote Wintner’s theorem as referred by [51], stating that if $∥f(t,\mathit{y})∥\le N\left(t\right)L\left(∥\mathit{y}∥\right)$ with N and L piecewise continuous, positive and L nondecreasing, such that for some $c>0$

$${\int }_c^{+\infty }\frac{1}{L\left(s\right)}ds=+\infty ,$$

then the ODE in Formula (18) has a solution for a given initial condition, by observing that the quoted theorem is valid, under the assumption that $f(t,\mathit{y})$ is continuous with the possible exception of points of a null Lebesgue set of the time variable, by considering extended solutions instead of usual solutions, which have a continuous derivative, and as $L\left(t\right)$ and $N\left(t\right)$ both satisfy the hypotheses of Wintner’s theorem for parameters $\tilde{\gamma }\left(t\right)$, $\delta \left(t\right)$ and $\mu \left(t\right)$, taking each two (or three) distinct values in two (or three) complementary intervals of the time domain.

We observe that these existence and uniqueness results are essential for a numerical integration of the ODE, but that no result of convergence numerical is implied—in the existence and uniqueness results above—for the regime switching ODE with discontinuous coefficients. Nevertheless, the Lipschitz condition with respect to the $\mathit{y}$ variable—such as the one in Formula (26)—is sufficient for the convergence of the Euler method (see [52] (p. 74)).

Remark 14

(On the existence of essentially only one rough regime switching model). We also observe that our definition of a rough regime switching as the continuous function obtained by glueing the usual solutions of ODE in complementary time intervals—as described in Section 3—is equal almost everywhere in the time variable, with respect to the Lebesgue measure, to any extended solution of the alternative rough regime switching obtained by solving an ODE with discontinuous parameters, which is the best we can expect for extended solutions.

Remark 15

(Existence and unicity of solutions for the smooth regime switching). It is clear that the existence and unicity of a solution for the smooth regime switching model in Section 3 may be dealt either by the usual Cauchy-Lipschitz existence and unicity theorem or by the results of this section, with the necessary adaptations of the definitions of the matrices $\mathit{M}(\Theta (t),\mathit{{\rm Y}}(t\left)\right)$ and ${\mathit{M}}_{\Delta }(\Theta \left(t\right),{\mathit{y}}_1,{\mathit{y}}_2)$. An example of smooth regime switching was studied in [17] with two coupled ODE, having in each equation two evolution regimes in such a way that the transition in time, between these regimes, is achieved by a coupling function similar to the function in Formula (9).

Remark 16

(On the Lie symmetries for ODE models with regime switching). There are at least two ways of considering symmetries on the model developed in this work. The first and more simple one is the usual reversibility in time for ODE. In fact it is clear from the expression of an ODE in integral form in Formula (19)—using $t=-s$ as the variable change—that the generalised solution of a ODE for a regime switching model proved using Carathéodory’s theorem is reversible in time. The second and profound one is related to the Lie symmetries for the differential model in Formula (5) with added regime switching as studied in this work. Lie type symmetries for ODE in general are the object of several reference works, for instance in [53]. For the SIR and the SI(RD) models, Lie symmetries were presented in [54,55], following the general proposals of using Lie symmetries in epidemiology developed in [56,57]. The study of Lie symmetries for ODE with regime switching remains to be done. A natural approach is to use the ideas of bifurcation theory, as masterfully developed in [58,59].

7. Discussion, Conclusions and Further Work

We introduced, studied, and applied a SIR type model with a compartment for deaths—occurring by force of the infection—and with regime switching to model pandemic mitigation measures such as a lockdown. The model depends on three main parameters, to wit, the basic reproduction number, the lethality rate and the duration of the disease. The statistical methodologies introduced to estimate the main parameters of the model—using real data of the initial period of the pandemia in Portugal—are natural, and provide estimates in line with values reported in the literature. Nevertheless, these values have to be resized in order to account for the true dimension of the infected population; this resizing creates secondary parameters which are impossible to estimate, due to the lack of proper data. When fed with the resized parameters—for some scenario chosen values of the secondary parameters—the model gives results in line with what is to be expected, for instance, the lowering of the death toll of mortality when a mitigating measure of lockdown type is present and, for the lethality also, a remarkable agreement within an excess of 2% of the observed data over the model. We must conclude that the assumptions made for the basis scenario are the most adequate of the three scenarios proposed; we observe that this fact gives us important information, about the disease, stated in the following. For the asymptomatic and/or undiagnosed the most plausible value of the contact number is 50% of the contact number of the diagnosed. The most plausible mortality rate—caused by the disease—for the asymptomatic and/or undiagnosed is zero and the most plausible value for the duration of the disease is 7.5 days. The methodology introduced is applicable whenever there are clear regime switchings in the evolution and correspondent estimates of the relevant parameters occurring by force of control measures such as lockdown. The main advantage of the SI(RD) model studied in this work is that, in its simplicity, it allowed us to recover the past behaviour of the lethality which is the main social nefarious consequence of an epidemic. One line of investigation that our work opens up to is to replace parameters by parametric functions depending on parametric variables that may capture the effects of our secondary parameters; this line of investigation seems justified by the importance of the secondary parameters in the conclusions we attained. An alternative to this line of investigation is to randomise the parameters; with this approach we would have to prescribe probability distributions for the parameters and then use Monte Carlo simulation to get probability distributions for the proportions of susceptible, recovered and dead, at the end of the period studied.