Discrete spread model for COVID-19: the case of Lebanon

Abstract

Background: Mathematical models are essential to predict the likely outcome of an epidemic. Various models have been proposed in the literature for disease spreads. Some are individual based models and others are compartmental models. In this study, discrete mathematical models are developed for the spread of the coronavirus disease 2019 (COVID-19).

Methods: The proposed models take into account the known special characteristics of this disease such as the latency and incubation periods, and the different social and infectiousness conditions of infected people. In particular, they include a novel approach that considers the social structure, the fraction of detected cases over the real total infected cases, the influx of undetected infected people from outside the borders, as well as contact-tracing and quarantine period for travelers. The first model is a simplified model and the second is a complete model.

Results: From a numerical point of view, the particular case of Lebanon has been studied and its reported data have been used to estimate the complete discrete model parameters using optimization techniques. Moreover, a parameter analysis and several prediction scenarios are presented in order to better understand the role of the parameters.

Conclusions: Understanding the role of the parameters involved in the models help policy makers in deciding the appropriate mitigation measures. Also, the proposed approach paves the way for models that take into account societal factors and complex human behavior without an extensive process of data collection.

Graphical abstract

Keywords

discrete stochastic modeling / COVID-19 / numerical simulation

1 INTRODUCTION

The world is in a battle with the coronavirus disease 2019 (COVID-19), which is caused by a highly virulent virus SARS-Cov-2 that targets the human respiratory system. The World Health Organization (WHO) declared the situation as a pandemic in March 2020 [1,2].

With the uncertainties about the COVID-19, predictive mathematical models are essential in order to discover the likely outcome of the epidemic, to inform healthcare needs and to minimize its socioeconomic consequences. Various models have been proposed in the literature for disease spreads. They can be categorized into agent-based models (ABM) [3,4] and compartmental models [5-7]. Compartmental models are built on differential equations and assume that the population is perfectly mixed with people moving between compartments such as susceptible (S), infected (I) and recovered (R) [8-11]. These models revealed the threshold nature of epidemics and triumphed in explaining “herd immunity”. However, they fail to capture complex social networks and the behavior of individuals who may adapt depending on disease prevalence. On the other hand, agent-based models can capture irrational behavior and are used to simulate the interactions of autonomous agents that can be either individuals or collective entities [12-16]. The ABM approach needs realistic data, typically obtained from a census, and important assumptions and data collection to set their structural parameters [17], but this is not always possible at the early stages of the outbreak. Several models have been developed for the COVID-19 pandemic. Anastassopoulou et al. included dead individuals in a discrete-time SIRD model and provided estimates of the main epidemiological parameters for Hubei (China) [18]. Casella derived a simplified control-oriented model comparing the outcomes of different policies [19] and Wu et al. inferred clinical severity estimates using transmission dynamics [20]. Giordano et al. proposed a mean-field epidemiological model extending the classical SIR model for the COVID-19 epidemic in Italy [21]. Sameni also proposed an SEIRP including compartments for asymptomatic exposed individuals and passed-away population [22] while Goel and Sharma suggested a mobility-based SIR model for COVID-19 pandemic [23]. Stochastic transmission models have been proposed in [24,25] and agent-based models have been used for a computational simulation of the pandemic in Australia by Chang et al. in [26] and for recommending universal masking by De Kai et al. in [16]. Moreover, in [27], Varotsos and Krapivin developed the COVID-19 decision making system to study disease transmission. Since the list of models developed so far for COVID-19 is much longer, and the aforementioned list is non-exhaustive, we also refer the reader to [28-31] for some insight on parameter calibration and inverse modeling.

The above compartmental models consider human populations as homogeneous. However, a realistic model should take into account that there are many heterogeneities in societies that affect disease transmission. In the present work, we propose a model that shows how such population heterogeneity can lead to substantial heterogeneity among the infected population (that appear in clusters). The consequent impact of such modeling strategy is on exit policies intending to minimize the risk of future infection spikes. We are trying to understand herein the effect of the social activity level of the individual on the spread of the disease. A common way to do this is by means of network models [32]. Since we are interested in a layer that lies between the micro-scale modeling represented by agent-based approach and the macro-scale modeling represented by compartmental models (such as SIR, SEIR or SEIRD among others), we propose two stochastic discrete models (a simplified and a complete) for the spread of COVID-19. We opt for discrete models in view of the daily reporting by countries of the infection indicators and we choose a non-deterministic approach in view of the randomness in the transmission of COVID-19. The models take into account the characteristics of the virus by attributing to them probability distributions. We mention, for instance, the incubation period, infectiousness, and testing sensitivity. They also incorporate probability distributions for social and individual conditions such as family size within the same household, number of people contacted per day and their subdivision into known versus unknown contacts. Up to the authors’ knowledge, the only works that were based on discrete modeling of COVID-19 epidemic are [33] where a discrete deterministic model is proposed for forecasting the temporal evolution of the epidemic from day to day and [34] where a discrete stochastic compartmental model is proposed. An important discrepancy between [34] and our work is that we use two discrete time variables one for the time of infection and the second one is for the daily update. This is crucial in modeling the COVID-19 epidemic due to the prevalence of pre-symptomatic transmission. Indeed, the distinction of two time variables permits the tracking, in the future, of people who get infected on a particular day and how long they stay actively infecting others before they get detected either by appearance of symptoms or by contact tracing combined with polymerase chain reaction (PCR) testing. We also distinguish the transmission between family members, known contacts and unknown contacts, hence including the social network in the model.

The remainder of the paper is laid out in the following way. Firstly, we describe the main parameters that must be accounted for in modeling epidemics such as COVID-19. Secondly, we present the discrete stochastic model in two steps. First, we describe a simplified version that does not include contact-tracing and quarantine period and second, we present the complete model. Numerical simulations for the validation, the prediction and the parameter analysis of the complete discrete stochastic model are presented then. Finally, we conclude with some remarks.

2 DESCRIPTION OF KEY PARAMETERS OF THE MODEL

To comprehend the dynamics of COVID-19 epidemic, several realistic parameters must be integrated in the model. Specifically, the chronological characteristics of the virus, the context of contact between individuals and the influx of infected undetected individuals from outside the borders are fundamental to be able to forecast the evolution of the epidemic from day to day. We shed light herein on the rationale for selecting keys parameters as follows:

- The times of infection and detection: In COVID-19 epidemic, many infected people may be infecting others in the community without being detected. So, it is important to reduce the delay between the infection time and the detection time [35]. In the proposed model, we use, for the flow of individuals from one category to another, two time variables to distinguish infection and detection times and to keep track in the future of the people infected on a particular day. This will be clarified in the next section.

- The probability of infection (or the transmission coefficient of the virus): This parameter depends on several factors such as mask wearing, adherence to hygiene measures and social distancing in public. This probability differs also depending on the nature of contact between the infected person and his contacts where we differentiate family members in the same household from contacts outside the house such as at work or in public transportation. The probability of infecting family members sharing the same house is assumed to be higher during the period extending from the instant of infection to the time of detection.

- The incubation period (time from exposure to illness onset): The general incubation period for COVID-19 ranges from 2 to 14 days depending on each individual. According to several recent studies, the average incubation time is about 5 days with the 95th percentile of the distribution at 12.5 days and the log-normal distribution provides the best fit to the data for incubation period estimates [36, 37]. We use herein a discrete random variable corresponding to the incubation period obtained from the suggested log-normal distribution in [36, 37].

- The serial interval (duration from onset of symptoms in an infector (a primary-case patient) to the onset of symptoms in an infectee (a secondary-case patient)): The mean serial interval for COVID-19 was estimated as 7.5 days in [36] and 6.9 days in [38] but more recent studies [39, 40] suggest that it is around 4 days. Being shorter than the mean incubation period, pre-symptomatic transmission is likely and may be more frequent than symptomatic transmission [21, 39]. In [39, 40] a normal distribution was suggested for the serial interval estimates. The knowledge of the distribution of the serial interval will help us estimate the proportion of infected people who will be detected according to the appearance of symptoms opposed to those who will be isolated according to contact tracing.

- The connections: These are the contacts of a person that we classify into family members, known contacts such as coworkers, neighbors and unknown contacts such as people encountered in public transportation or in social events. To account for family members we consider a discrete random variable obtained from a log-normal probability distribution function for the average family size or household size (see Fig.1 for the case of Lebanon). For known contacts, we consider an exponential distribution for the daily new encounters. So if a person was in contact with the individual on the first day, he is not counted in the new encounters of the following days. In other words, the first encounter is considered for the possible infection of the person (see Fig.1 for the case of Lebanon). Finally, the unknown contacts are chosen to be uniformly distributed.

- Contact-tracing: Identifying the source of infection, tracing and isolating its contacts are crucial for breaking the chain of transmission and for the control of the epidemic [24]. Very high levels of contact tracing are required in the case of presymptomatic infectiousness [24]. In the present model, we account for contact tracing efficiency as detailed in Section of Discrete Mathematical Model.

- The effect of border opening: In order to account for the effect of airport and border opening, we split infected individuals into travelers and locals. The probability distribution function for infected travelers is based on data collected from local health authorities. Furthermore, we distinguish between travelers having a positive screening result and those with a false negative result. Those with a positive test result are assumed to be isolated and are no longer infectious (see Fig.1 where based on data collected over 188 days from the beginning of the epidemic, we show the number of days where a certain number of travelers are screened positive for COVID-19, e.g. on 51 out of 188 days, there were no positively screened travelers and on 16 days there was only 1 infected traveler.). The number of travelers with false negative test result is deduced from the number of travelers with positive test result based on the ratio of false negative screening results in the country. Moreover, a quarantine period is imposed on travelers and it changes according to policies adopted by local authorities. However, the rate of compliance to this measure varies from one individual to another. On the other hand, the local infected residents are subdivided into different categories: the symptomatic individuals with positive test result versus the individuals screened due to their contact with an infected person before onset of symptoms.

Some of the aforementioned parameters are calibrated to fit the real data by numerical optimization (e.g. probability of infection) or obtained from different studies and publications (e.g. incubation period distribution).

3 DISCRETE MATHEMATICAL MODEL

Daily infection being reported by countries worldwide, the unit of time adopted in the model is a single day, denoted herein by n.

3.1 Simplified model

We start from a fully susceptible population. After the introduction of the virus, a susceptible individual may be infected and then removed. By removed individuals, we designate those that have been detected and therefore isolated. The infected members of the population are split into four categories: \{P, N,F,C \}, so that the cumulative number I(n) of infected individuals up to a day n is given by

(1)\begin{array}{l}I(n)=P(n)+N(n)+F(n)+C(n),\end{array}

where

P:\mathbb{N}\to {\mathbb{R}}^{+}, the cumulative number of travelers with positive PCR result,

N:\mathbb{N}\to {\mathbb{R}}^{+}, the cumulative number of infected travelers with false negative PCR,

F:\mathbb{N}\to {\mathbb{R}}^{+}, the cumulative number of infected family members, and

C:\mathbb{N}\to {\mathbb{R}}^{+}, the cumulative number of infected contacts.

Let \mathrm{\Delta }I: \mathbb{N}\to {\mathbb{R}}^{+} denote the daily number of new infected individuals and \mathrm{\Delta } P:\mathbb{N}\to {\mathbb{R}}^{+} the daily number of new travelers with positive PCR result which is a source term to the model and it is equal to zero when borders are closed. Moreover, we introduce the following functions: \mathrm{\Delta } N(\cdot , \cdot ), \mathrm{\Delta } F(\cdot , \cdot ), and \mathrm{\Delta } C(\cdot , \cdot ) defined from \mathbb{N}\times \mathbb{N} into {\mathbb{R}}^{+}, where \mathrm{\Delta }X(n,n ) denotes the daily new number of infected individuals from category X\in \{N, F,C\} that have been infected on day n and \mathrm{\Delta }X(k,n ), for k>n, denotes the remaining number out of \mathrm{\Delta } X(n, n) that are still active on day k. By convention, we write \mathrm{\Delta } X(0, 0)=0.

The cumulative number of cases from category X, X\in \{N, F,C\}, is given by

(2)\begin{array}{l}X(n)=X(n-1)+\mathrm{\Delta }X(n,n).\end{array}

Furthermore, the daily new number of infected cases, \mathrm{\Delta } I(n), is modeled by the following formula

(3)\begin{array}{l}\mathrm{\Delta }I(n)=\mathrm{\Delta }P(n)+\mathrm{\Delta }N(n,n)+\mathrm{\Delta }F(n,n)+\mathrm{\Delta }C(n,n)\end{array}

where \mathrm{\Delta }N(n,n ), the daily number of new infected travelers with false negative PCR test, is a source term that is related to the number \mathrm{\Delta } P(n). For instance, in Lebanon, the PCR test gives false negative result for about 30\mathrm{\%} of the infected travelers.

Remark 1. Studies of false-negative PCR results for COVID-19 are variable demonstrating false-negative rates ranging from 1% to 30% [41,42]. False-negative test results are attributed to several reasons including sub-optimal specimen collection, testing too early in the disease process, low analytic sensitivity, low viral load, or variability in viral shedding [43-48].

In order to take into consideration the decrease of the number of susceptible local individuals, we distinguish between the local infected cases and the infected travelers. So if we denote by N_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}} the total number of the local population, by S(n) the number of susceptible local cases on day n and I_l(n) the cumulative number of infected local individuals, then we have

(4)\begin{array}{rl}{I}_l(n)=& I_l(n-1)+\mathrm{\Delta }F(n,n)+\mathrm{\Delta }C(n,n)\\ =& \sum _{k=1}^n\left[\mathrm{\Delta }F(k,k)+\mathrm{\Delta }C(k,k)\right]\end{array}

and the following conservation equation as follows

(5)S(n)+V(n) +I_l(n) =N_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}},

where we denote by V(n) the total number of vaccinated people on day n, and it changes daily according to the rate of vaccination in the population. Thus, with this conservation equation, the immunity of both infected and vaccinated cases is assumed to be lifelong. However, waning immunity can be considered by adding to the number S(n) of susceptible cases on day n, the number of new infected and vaccinated cases on day n-k where some period of k previous days (for instance after 6 months, k is about 180 days).

In this model, we assume that the number of non-infected travelers is the same as the number of susceptible locals who leave the country.

The infection process of the susceptible population will evolve as follows. New infections in the category F on day n are introduced either from an infected contact or from a traveler (with a false negative test result) who transmit the virus to their household members. We assume that household members are infected only on the next day of the infection time of their infector.

Similarly, new infections in the category C on day n result from family members infecting their new contacts from the time of their infection until day n-1 or similarly from travelers or other contacts (see Fig.2 for a schematic of the infection process between the three categories N, F and C). These assumptions are represented by the following relationships

(6)\begin{array}{rl}\mathrm{\Delta }F(n,n)=& \alpha p_S(n)[dF(1)\mathrm{\Delta }C(n-1,n-e)\\ & +dF(1)\mathrm{\Delta }N(n-1,n-e)],\end{array}

(7)\begin{array}{rl}\mathrm{\Delta }C(n,n)=& \beta p_S(n)[\sum _{i=e}^{n-1}dC(i)\mathrm{\Delta }F(n-1,n-i)\\ & +\sum _{i=e}^{n-1}dC(i)\mathrm{\Delta }C(n-1,n-i)\\ & +\sum _{i=e}^{n-1}dC(i)\mathrm{\Delta }N(n-1,n-i)],\end{array}

where \alpha and \beta are the infection probability coefficients within and outside the same household respectively and that will be estimated according to the real data, p_S(n )={\displaystyle \frac{S(n-1)}{N_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}}} is the proportion of susceptible individuals within the total population, dF is the random variable of the family size whose probability distribution for the case of Lebanon is given in Fig.1, and dC is a time dependent random variable for the number of new contacts met per day starting from the infection time. We assume that dC, as a function of time, has an exponential form that depends on both social habits and measures imposed by authorities. On the first day after the infection time, dC is maximal and it decreases to zero as time increases. We denoted by e the latency period (the period of time in which newly infected individuals are asymptomatic and noninfectious). We consider the latency period e as a constant time delay in our model and we assume that it is fixed and equal to 3 days.

Remark 2. In [49], the latent period for COVID-19 was estimated to have a mean of 4.19 days while in [50], the authors used two differential equations’ models to study the impact of this period. In the first model, they included a compartment for individuals exposed but not yet infectious and in the second model they included a time delay to model the latency period. They also concluded that the second model gives a better prediction than the first. However, the latency period was estimated to range between 6 and 12 hours.

On the other hand, as for the removal process, we denote by RX(n) the number of removed cases from category X on day n. We have the following relationship

(8)RX(n )=\sum _{k=1}^{n-1}\left[\mathrm{\Delta }X(n-1,k)-\mathrm{\Delta }X(n,k)\right].

In a first step, we assume that the removal of infected individuals takes place according to the appearance of symptoms which leads to isolation or to the cure of asymptomatic cases so that they are no longer infectious. Thus, there is no contact tracing in this simplified model. Also, we do not take into account the quarantine period imposed on the travelers by the government.

We denote by dS the discrete random variable of the incubation period. Then we have the following formula

(9)\mathrm{\Delta }X(n+i,n)=\mu (i+1 )\mathrm{\Delta } X(n, n),

where , with g(k)= \mathbb{P}[dS= k] is the probability that the incubation period is equal to k. Notice that \mu (1)=1 and \underset{i\to +\mathrm{\infty }}{lim}\mu ( i)=0.

Accordingly, we can prove that

(10)RX(n )=\sum _{i=1}^{n-1}g( n-i)\mathrm{\Delta }X(i,i ).

So, now the daily new number of removed cases is given by

(11)\mathrm{\Delta }R(n)= \mathrm{\Delta }P(n) +RN( n)+R F(n)+RC (n),

and the total number of removed cases on day n is

(12)R(n)=R(n-1 )+\mathrm{\Delta }R(n).

Therefore, the current number of active infected cases (i.e. that can infect others) is given by

(13)I_a(n)= I(n)-R(n) .

As for the vaccinated people, we assume that vaccination takes place with a constant rate, denoted \tau, out of the number of the current susceptible individuals. So, we propose the following equation

(14)V(n)=V(n-1 )+\tau S(n),\mathrm{f}\mathrm{o}\mathrm{r}n⩾n_0

where n_0 denotes the starting time for vaccination.

For a better understanding of the discrete model, we aimed to solely express \mathrm{\Delta }F(n,n ) and \mathrm{\Delta } C(n, n) (given in Eqs. (6) and (7)) in terms of \mathrm{\Delta }F(k,k ), \mathrm{\Delta } C(k, k) and \mathrm{\Delta } N(k, k) that can be viewed as functions of one time variable k. The resulting system can be considered well fitting into the discrete time series models. Moreover, this way of expressions turns out to be helpful in the derivation of a continuous model.

Now, for \mathrm{\Delta }F(n,n ) and \mathrm{\Delta } C(n, n) given as in Eqs. (6) and (7) respectively, we have the following system

(15)\begin{array}{rl}{p}_S(n)=& 1-{\displaystyle \frac{I_l(n-1)+V(n-1)}{N_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}}}\\ =& 1-{\displaystyle \frac{1}{N_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}}}\sum _{k=1}^{n-1}(\mathrm{\Delta }F(k,k)+\mathrm{\Delta }C(k,k))\end{array}

(16)\begin{array}{rl}\mathrm{\Delta }F(n,n)=& \alpha p_S(n)[dF(1)\mathrm{\Delta }C(n-e,n-e)\\ & +dF(1)\mathrm{\Delta }N(n-e,n-e)]\end{array}

(17)\begin{array}{rl}\mathrm{\Delta }C(n,n)=& \beta p_S(n)\sum _{i=1}^{n-e}[d\lambda (n-i)\mathrm{\Delta }F(i,i)\\ & +d\lambda (n-i)\mathrm{\Delta }C(i,i)+d\lambda (n-i)\mathrm{\Delta }N(i,i)]\end{array}

where d\lambda (i)=d C(i)\mu ( i).

Remark 3. System (15)-(16)-(17) shows that the daily new number of infected family members \mathrm{\Delta } F(n, n) and contacts \mathrm{\Delta }C(n,n) have a similar formulation to a nonlinear autoregressive model with a delayed input source term \mathrm{\Delta } N(n, n) (nonlinear ARX model) but whose coefficients are random variables. These coefficients vanish after some interval of time related to the maximum incubation period.

3.2 Complete model

In this model, contact tracing and quarantine period are explicitly involved. In addition to the categories P,N,F and C, we introduce a new category denoted U. Now the category C denotes the contacts that are infected by individuals they know (such as cousins, neighbors, colleagues, friends) which facilitates the contact tracing, and U is concerned with individuals that are infected from an unknown source of infection (such as public transportation, social events) so that contact tracing is impossible for these cases. Thus the daily new number of infected cases, \mathrm{\Delta } I(n), will be given by

(18)\begin{array}{rl}\mathrm{\Delta }I(n)=& \mathrm{\Delta }P(n)+\mathrm{\Delta }N(n,n)+\mathrm{\Delta }F(n,n)\\ & +\mathrm{\Delta }C(n,n)+\mathrm{\Delta }U(n,n),\end{array}

the cumulative number of local infected cases verifies

(19)\begin{array}{l}{I}_l(n)=I_l(n-1)+\mathrm{\Delta }F(n,n)+\mathrm{\Delta }C(n,n)+\mathrm{\Delta }U(n,n),\end{array}

and the daily new number of removed cases is given by

(20)\begin{array}{l}\mathrm{\Delta }R(n)=\mathrm{\Delta }P(n)+RN(n)+RF(n)+RC(n)+RU(n)\end{array}

where each RX(n ) is given by Eq. (8) for X\in \{ N,F,C ,U\}.

Moreover, we split the categories F and C into two categories each: F is split into F_S and F_T, and C is split into C_S and C_T. For instance, F_S is the subset of F containing those who are infected from category F and removed later according to symptoms or cured after a while, whereas F_T corresponds to the subset of F containing those who are removed before symptoms occur and isolated according to contact tracing when their infectors have been detected. This idea of splitting is valid only for F and C and it is not valid for U and N because these last two categories contain the individuals whose source of infection is unknown (see Fig.3 for a schematic of the transmission process of each category). We assume that the subsets F_S and C_S form a certain proportion p of the categories F and C respectively, with . Therefore, we have the following relations

(21)\begin{array}{l}\begin{array}{r}\mathrm{\Delta }{F}_S(n,n)=p\mathrm{\Delta }F(n,n),\mathrm{\Delta }{F}_T(n,n)=(1-p)\mathrm{\Delta }F(n,n),\\ \mathrm{\Delta }{C}_S(n,n)=p\mathrm{\Delta }C(n,n),\mathrm{\Delta }{C}_T(n,n)=(1-p)\mathrm{\Delta }C(n,n).\end{array}\end{array}

For the cases that are detected after the appearance of symptoms or after healing, we use the same formulation that we developed in the simplified model. Indeed, the numbers \mathrm{\Delta }{F}_S(n,n), \mathrm{\Delta } C_S(n,n), \mathrm{\Delta }N(n,n ) and \mathrm{\Delta } U(n, n) of new infected people on day n, will change during the next days according to Eq. (9).

In order to take into account the quarantine period imposed by the government on travelers, we modify the term that accounts for the number of people infected by travelers with false negative test and have been quarantined for some period. Indeed, we propose the following system

(22)\begin{array}{rl}\mathrm{\Delta }F(n,n)=& \mathit{\alpha }{p}_S(n)[dF(1)\mathrm{\Delta }C(n-1,n-e)\\ & +dF(1)\mathrm{\Delta }U(n-1,n-e)\\ & +\sum _{i=e}^{n-1}{\delta }_F(n-1,n-i)\mathrm{\Delta }N(n-1,n-i)]\end{array}

(23)\begin{array}{rl}\mathrm{\Delta }C(n,n)=& \mathit{\beta }{p}_S(n)[\sum _{i=e}^{n-1}dC(i)\mathrm{\Delta }F(n-1,n-i)\\ & +\sum _{i=e}^{n-1}dC(i)\mathrm{\Delta }C(n-1,n-i)\\ & +\sum _{i=e}^{n-1}dC(i)\mathrm{\Delta }U(n-1,n-i)\\ & +\sum _{i=e}^{n-1}{\delta }_C(n-1,n-i)\mathrm{\Delta }N(n-1,n-i)]\end{array}

(24)\begin{array}{rl}\mathrm{\Delta }U(n,n)=& \mathit{\gamma }{p}_S(n)[\sum _{i=e}^{n-1}dU(i)\mathrm{\Delta }F(n-1,n-i)\\ & +\sum _{i=e}^{n-1}dU(i)\mathrm{\Delta }C(n-1,n-i)\\ & +\sum _{i=e}^{n-1}dU(i)\mathrm{\Delta }U(n-1,n-i)]\end{array}

where q denotes the quarantine period and, {\delta }_F(n+q,n)= dF(1) and {\delta }_F(k,n)=0 for k\ne n+q. Similarly, {\delta }_C( n+q+k ,n)= dC(k +1) for k⩾0 and {\delta }_C(k,n)= 0 for .

Since not all the travelers respect this period, we assume that the number q is uniformly distributed over an interval [1,q_{max}] where q_{max} is the quarantine period imposed by the government. The real number \gamma denotes the infection probability coefficient related to the unknown contacts and dU is a random variable for the number of unknown contacts met per day which is supposed to have a uniform probability distribution that also depends on social habits and measures imposed by authorities.

Now, regarding the subsets F_T and C_T their corresponding individuals will be isolated as a result of contact tracing. By contact tracing we mean that when an infector has been detected then those who were infected from the categories F_T and C_T will be isolated. We can express this statement by the following formulation

(25)\begin{array}{rl}R{F}_T(n)=& (1-p)\alpha p_S(n)[dF(1)RC(n-1)\\ & +dF(1)RU(n-1)+dF(1)RN(n-1)].\end{array}

Using Eq. (8) for all the terms in Eq. (25), we obtain, on one hand,

(26)R{F}_T(n) =\sum _{k=1}^{n-1}\left[\mathrm{\Delta }{F}_T(n-1,k)-\mathrm{\Delta }{F}_T(n,k)\right]

and on the other hand,

(27)\begin{array}{rl}& R{F}_T(n)\\ & =(1-p)\alpha p_S(n)dF(1)\sum _{k=1}^{n-2}\left[\mathrm{\Delta }C(n-2,k)-\mathrm{\Delta }C(n-1,k)\right]\\ & +(1-p)\alpha p_S(n)dF(1)\sum _{k=1}^{n-2}\left[\mathrm{\Delta }U(n-2,k)-\mathrm{\Delta }U(n-1,k)\right]\\ & +(1-p)\alpha p_S(n)dF(1)\sum _{k=1}^{n-2}\left[\mathrm{\Delta }N(n-2,k)-\mathrm{\Delta }N(n-1,k)\right].\end{array}

By identification between Eqs. (26) and (27), one gets the updates \mathrm{\Delta }{F}_T(k, n) for k⩾n+1 as follows:

(28)\{\begin{array}{rl}& \mathrm{\Delta }{F}_T(n+i,n)=\mathrm{\Delta }{F}_T(n,n),1⩽i⩽e\\ & \mathrm{\Delta }{F}_T(k,n)=(1-p)\alpha p_S(n)[dF(1)\mathrm{\Delta }C(k-1,n)\\ & +dF(1)\mathrm{\Delta }U(k-1,n)\\ & +dF(1)\mathrm{\Delta }N(k-1,n)],\mathrm{f}\mathrm{o}\mathrm{r}k⩾n+e+1.\end{array}

We proceed in the same way to obtain the updates \mathrm{\Delta } C_T(k,n) for k⩾n+1. The only difference with \mathrm{\Delta } F_T(k,n) is in the terms that model the number of people met from family members or known contacts. We obtain the following formula

(29)\{\begin{array}{rl}& \mathrm{\Delta }{C}_T(n+i,n)=\mathrm{\Delta }{C}_T(n,n),1⩽i⩽e\\ & \mathrm{\Delta }{C}_T(k,n)=(1-p)\beta p_S(n)[dC(k-n)\mathrm{\Delta }F(k-1,n)\\ & +dC(k-n)\mathrm{\Delta }C(k-1,n)\\ & +dC(k-n)\mathrm{\Delta }U(k-1,n)\\ & +dC(k-n)\mathrm{\Delta }N(k-1,n)],\mathrm{f}\mathrm{o}\mathrm{r}k⩾n+e+1.\end{array}

Moreover, as in the simplified model, we were able to express the updates \mathrm{\Delta }F(n+k ,n) and \mathrm{\Delta }C(n+k ,n) in terms of \mathrm{\Delta }{F}_S(n,n), \mathrm{\Delta }{F}_T(n,n), \mathrm{\Delta }{C}_S(n,n), \mathrm{\Delta }{C}_T(n,n), \mathrm{\Delta }U(n,n ), and \mathrm{\Delta } N(n, n). Again, we obtain a sort of nonlinear ARX model for the numbers of daily new cases in all of the categories F_S,F_T,C_S,C_T and U where N and P are source terms or delayed inputs.

4 NUMERICAL SIMULATIONS: CASE OF LEBANON

The first confirmed case of COVID-19 in Lebanon was detected on February 21, 2020 and by February 29, 2020, there was a total of 7 confirmed cases going back from travel trips. Consequently, educational institutions were closed, and gradual measures were introduced until the declaration of general mobilization and state of emergency by mid-March where public gatherings were banned, cultural venues closed and social distancing measures imposed in public. Furthermore, the airport, the land borders and the seaports were closed as of March 19. On April 5, new arrivals were intermittently allowed through the airport. By the end of May, authorities in Lebanon have been gradually easing restrictions. Public transportation has resumed, with social-distancing measures. Government institutions and certain private companies, including various shops and stores, were permitted to return to normal operations as of June 1st and the lockdown was lifted on June 7. Starting July 1st, flights were resumed and the airport started operating up to 10 \mathrm{\%} of its capacity [51]. Travelers with negative PCR results upon arrival were isolated for a maximum of 3 days. However, people were not abiding by the preventive measures, and travelers were not respecting the isolation period [52]. The virus rebounded to reach a total number of 4205 by July 29 and the authorities reinstated lock-down from July 30 to August 3 and from August 6 to August 10. However, due to the explosion in Beirut, the sanitary situation deteriorated in the country and the total number of cases reached 16,870 by August 31. In this context, it is obvious that there is an urgent need for evidence-based decisions. Furthermore, the number of cases in December shockingly increased indicating that something new was happening during the first two weeks of December. Most probably, this may be related to the new UK variant of the virus that has been discovered for the first time in Lebanon by mid-December. This variant is known for its fast spread and currently most of the infected cases in Lebanon are due to this variant.

4.1 Calibration and parameters’ estimation

In order to obtain a realistic model for the case of Lebanon, we used realistic data about the various inputs of our model. For instance, we used for the incubation period a discrete distribution obtained from [36, 37]. On the other hand, due to the lack of census in Lebanon, we used the reports [53, 54] to deduce a distribution for family size in the same household (see Fig.1). In the absence of official reports, the distribution of new people met per day was also deduced from our familiarity with the Lebanese context where public transport is not prevalent and the number of extended family members met per day is high.

However, as for the coefficients that are involved in the model, we used the real data published on the website of the ministry of public health of Lebanon (https://coronanews-lb.com/). Indeed, these coefficients have been estimated using numerical optimization techniques in order to fit the real observed data. The estimated parameters are displayed in Tab.1.

These parameters change depending on the period of time that corresponds to different measures in the country. For the whole period, the probability of infection for the same household was estimated as \alpha =0.15 and the probability p related to contact tracing is p=0.1.

The model being stochastic, we use 60 runs for each simulation and we display their average along with the corresponding standard deviation. The result of the simulation along with the real data are shown in Fig.5.

The results of the simulations greatly overlap with the real data. The discrepancy in the latest stages has been intentionally added and it is attributed to the fact that the real number of daily infected people is not fully reported. Indeed, many asymptomatic people are never detected because contact tracing is not 100\mathrm{\%} efficient while others may develop minor to mild symptoms but do not get tested for several reasons (some of them are financial). It is crucial to have mass testing to detect such cases and better contain the epidemic.

Remark 1. The error bars in the simulations are about 30\mathrm{\%} from the mean. This could be explained by the choice of the probability distributions of the family size and the known contacts where the standard deviation is about 30\mathrm{\%} from the average. Indeed, reducing the standard deviation in those distributions also reduces the error bars for the simulations. On the other hand, after running the simulation 60 times, we noticed that the results have a kind of normal distribution and we considered their average as the representative result.

For a deeper insight on the impacts of the main parameters involved in our model, we studied by numerical simulations the effect of the change of the parameters for the case study of Lebanon. Mainly we focused here on the role of the two main parameters \beta and \gamma, that model the infection coefficients for known and unknown contacts respectively, since these are the parameters that are mostly tied to measures. We considered four different scenarios. In scenario (I), we reduced beta by 2% compared to the scenario with constant coefficients (where we assume that following February 8, \beta = 0.040 and \gamma = 0.035). In (II), we reduced beta by 2% and gamma by 10%, in (III), we reduced beta by 5% and in (IV) we reduced gamma by 10%, see Fig.6.

We observe from these simulations that, for the case of Lebanon, the parameter \beta is the most effective in the control of the spread of the epidemic. Since \beta is the infection coefficient for known contacts, this result suggests that the epidemic can be efficiently controlled when strict measures are imposed in the context of workplaces, universities, schools, family gatherings or shops in villages where people know each other. On the other hand, if such measures are not taken in those places, authorities must impose much stricter measures on public transportation, malls, restaurants in cities (to name a few places where people get in contact with strangers).

We can also confirm this observation when we look at the number of infected people in each of the categories F, C and U, see Fig.7. Indeed, the curves of the number of infected cases in the category C are the highest however they are the lowest for the category U. These observations hold for the case of Lebanon and these are not necessarily valid for the case of other countries where the category U may have an important role in transmitting the virus mainly when for instance public transportation is essential in the society.

4.2 Prediction

Because the parameters are strongly related to people habits and to government actions regarding lockdown or ease of measures, it will be hard to predict a realistic evolution of the spread of the virus. However, we propose two scenarios that help us better understand how the virus may evolve.

First, we assume that the parameters hold constant after the last used observed data (February 8, 2021), see Fig.7. The simulations are done till the end of July, 2021.

Second, we considered variable coefficients \beta and \gamma. Their variation takes place according to a significant increase or decrease in the reported number of cases. In other words, when the number of the daily new reported cases highly increases, then the coefficients decrease and conversely, when the number of the daily new reported cases highly decreases, then the coefficients increase, see Fig.8. Such a scenario represents the population’s reaction to an increase in daily reported cases (then they tend to respect the sanitary measures) or a decrease in those cases (then they tend to disregard the sanitary measures). The simulations are done till the end of February, 2025.

Analyzing Fig.7 and Fig.8, we clearly see that with constant coefficients the spread of the virus will slow down significantly, however, with variable coefficients, which is a more realistic scenario since it models the reaction of humans to measures and to the spread speed, several waves of the epidemic are expected to hit the population. The waves’ amplitudes become smaller with time. This means that there will be a high risk that the virus will sustain for several years or probably forever if new variants with faster spreading will appear in the future.

Vaccination may help avoid this problem. But this must be done rapidly before the appearance of new virus variants that are not covered by the existing vaccines. For instance, the vaccination campaign started in Lebanon on February 21, 2021. Up to April 27, there has been 155,812 fully vaccinated people. The vaccination rate is about 2,500 people per day. In Fig.9, we show the simulation results for a scenario where we assume that the immunity against the virus lasts only 6 months and we consider variable coefficients. It is clear that under this scenario the epidemic will decline as long as no new variants appear. Moreover, comparing the trend of the daily new cases in this simulation to the real data till June 2, 2021, we clearly see that the simulation results greatly overlap with the reported cases by the ministry of public health in Lebanon as shown in Fig.4.

5 CONCLUSION

In this work, we proposed two stochastic discrete models for the spread of COVID-19. These models may be generalized to any other epidemic by considering appropriate characteristics of the pathogen. In these models, we include several social indicators related to family structure and habits by splitting the infected individuals into several categories and we take into account contact-tracing and quarantine period. We used two discrete time variables to permit tracing infected people from the time of infection until the time of detection and account for delayed effects. We were able to express the daily new cases at a given day in terms of the daily new cases of previous days. The resulting model is similar to a nonlinear auto-regressive discrete time series model with delayed input corresponding to the travelers’ inflow to the country. Furthermore, we validated the efficiency of the complete discrete model using the available data for the case of Lebanon. Indeed, we were able to reproduce the pattern of this data by estimating the parameters involved in the model using optimization techniques. We deduced that in the context of Lebanon, to control the epidemic without imposing a strict lock-down, local authorities must insure that social distancing and hygiene measures are respected within, for instance, workplaces and educational institutions. In addition, familial gatherings (in the context of extended families) must be banned. As for the prediction, we considered three scenarios: the first one with constant parameters over time, the second with variable parameters and the third one taking into account the effect of vaccination. The scenario with variable parameters seems to be more realistic since it captures the reaction of individuals in response to a significant increase or decrease in the reported number of cases. On the other hand, it is clear that vaccination plays a central role in reducing the number of new infections and avoiding the new infection spikes observed in the variable parameters scenario.

Although the proposed model was applied to COVID-19 epidemic in Lebanon, the same approach may be applied to other countries by changing the distributions of household size, known contacts and unknown contacts. Moreover, a continuous integro-differential model can be obtained from the proposed discrete model and in this case, one may add the spatial distribution of the population to the model and take into account the effect of mobility between different regions on the spread of the disease. In summary, we consider that this approach paves the way for models that take into account societal factors and complex human behavior without an extensive process of data collection.

References

[1] WorldHealth Organization. Statement on the second meeting of the international health regulations. Emergency committee regarding the outbreak of novel coronavirus (2019-ncov). Accessed: December 1, 2020

[2] WorldHealth Organization. Director-general’s opening remarks at the media briefing on covid-19—11 march 2020. Accessed: December 1, 2020

[3] BonabeauE.. ( 2002) Agent-based modeling: Methods and techniques for simulating human systems. In: Proceedings of the National Academy of Sciences, 99(suppl 3), 7280‒ 7287

[4] Ajelli,M., alves,B., Balcan,D., Colizza,V., Hu,H., Ramasco,J. J., Merler,S. ( 2010). Comparing large-scale computational approaches to epidemic modeling: agent-based versus structured metapopulation models. BMC Infect. Dis., 10 : 190

[5] Hethcote,H. ( 2000). The mathematics of infectious diseases. SIAM Review, 42 : 599– 653

[6] DiekmannO., HeesterbeekH.. ( 2012) Mathematical Tools for Understanding Infectious Disease Dynamics, volume 7. Princeton University Press

[7] BrauerF.. ( 2012) Mathematical Models in Population Biology and Epidemiology, volume 2. Springer

[8] KermackW. O. McKendrickA.. ( 1927) A contribution to the mathematical theory of epidemics. In: Proceedings of the royal society of london. Series A, Containing papers of a mathematical and physical character, 115, 700‒ 721

[9] sellI., ( 1996) The quasi-stationary distribution of the closed endemic sis model. Adv. Appl. Probab., 895‒ 932

[10] Hurley,M., Jacobs,G. ( 2006). The basic SI model. In: New Directions for Teaching and Learning, 2006 : 11– 22

[11] JinY. WangW.. ( 2007) An sirs model with a nonlinear incidence rate. Chaos, Solitons, Fract., 34, 1482‒ 1497

[12] Epstein,J. ( 2009). Modelling to contain pandemics. Nature, 460 : 687– 687

[13] Hunter,E., Namee,B. M. ( 2017). A taxonomy for agent-based models in human infectious disease epidemiology. J. Artif. Soc. Social Simul., 20 : 2

[14] Hunter,E., Mac Namee,B. ( 2018). An open-data-driven agent-based model to simulate infectious disease outbreaks. PLoS One, 13 : e0208775

[15] Tracy,M., Cerda,M. Keyes,K. ( 2018). Agent-based modeling in public health: current applications and future directions. Annu. Rev. Publ. Health, 39 : 77– 94

[16] De KaiG. MorgunovA., NangaliaV.. ( 2020) Universal masking is urgent in the COVID-19 pandemic: Seir and agent based models, empirical validation, policy recommendations. arXiv, 2004.13553

[17] Epstein,J. M., Goedecke,D. M., Yu,F., Morris,R. J., Wagener,D. K. Bobashev,G. ( 2007). Controlling pandemic flu: the value of international air travel restrictions. PLoS One, 2 : e401

[18] Anastassopoulou,C., Russo,L., Tsakris,A. ( 2020). Data-based analysis, modelling and forecasting of the COVID-19 outbreak. PLoS One, 15 : e0230405

[19] CasellaF.. ( 2020) Can the COVID-19 epidemic be managed on the basis of daily data? arXiv, 2003.06967

[20] Wu,J. T., Leung,K., Bushman,M., Kishore,N., Niehus,R., de Salazar,P. M., Cowling,B. J., Lipsitch,M. Leung,G. ( 2020). Estimating clinical severity of COVID-19 from the transmission dynamics in Wuhan, China. Nat. Med, 26 : 506– 510

[21] Giordano,G., Blanchini,F., Bruno,R., Colaneri,P., Di Filippo,A., Di Matteo,A. ( 2020). Modelling the COVID-19 epidemic and implementation of population-wide interventions in Italy. Nat. Med, 26 : 855– 860

[22] SameniR.. ( 2020) Mathematical modeling of epidemic diseases; a case study of the COVID-19 coronavirus. arXiv, 2003.11371

[23] Goel,R. ( 2020). Mobility based SIR model for pandemics−with case study of COVID-19. In: 2020 IEEE/ACM International Conference on Advances in Social Networks Analysis and Mining (ASONAM), pp. 110– 117

[24] Hellewell,J., Abbott,S., Gimma,A., Bosse,N. I., Jarvis,C. I., Russell,T. W., Munday,J. D., Kucharski,A. J., Edmunds,W. J., Funk,S. . ( 2020). Feasibility of controlling COVID-19 outbreaks by isolation of cases and contacts. Lancet Glob. Health, 8 : e488– e496

[25] Kucharski,A. J., Russell,T. W., Diamond,C., Liu,Y., Edmunds,J., Funk,S., Eggo,R. M., Sun,F., Jit,M., Munday,J. D. . ( 2020). Early dynamics of transmission and control of COVID-19: a mathematical modelling study. Lancet Infect. Dis., 20 : 553– 558

[26] Chang,S. Harding,N., Zachreson,C., Cliff,O. ( 2020). Modelling transmission and control of the COVID-19 pandemic in Australia. Nat. Commun., 11 : 5710

[27] Varotsos,C. Krapivin,V. ( 2020). A new model for the spread of COVID-19 and the improvement of safety. Saf. Sci., 132 : 104962

[28] Comunian,A., Gaburro,R. ( 2020). Inversion of a SIR-based model: A critical analysis about the application to COVID-19 epidemic. Physica D, 413 : 132674

[29] Calafiore,G. C., Novara,C. ( 2020). A modified SIR model for the COVID-19 contagion in Italy. In: 2020 59th IEEE Conference on Decision and Control (CDC), pp. 3889– 3894

[30] Calvetti,D., Hoover,A., Rose,J. ( 2020). Bayesian dynamical estimation of the parameters of an SE(A)IR COVID-19 spread model. arXiv, 2005.04365

[31] Rouabah,M. T., Tounsi,A. Belaloui,N. ( 2020). Early dynamics of COVID-19 in Algeria: a model-based study. ArXiv, 2005.13516

[32] Pastor-Satorras,R. ( 2001). Epidemic spreading in scale-free networks. Phys. Rev. Lett., 86 : 3200– 3203

[33] Boulmezaoud,T. ( 2020). A discrete epidemic model and a zigzag strategy for curbing the COVID-19 outbreak and for lifting the lockdown. Math. Model. Nat. Phenom., 15 : 75

[34] He,S., Tang,S. ( 2020). A discrete stochastic model of the COVID-19 outbreak: Forecast and control. Math. Biosci. Eng, 17 : 2792– 2804

[35] Klinkenberg,D., Fraser,C. ( 2006). The effectiveness of contact tracing in emerging epidemics. PLoS One, 1 : e12

[36] Li,Q., Guan,X., Wu,P., Wang,X., Zhou,L., Tong,Y., Ren,R., Leung,K. S. M., Lau,E. H. Y., Wong,J. Y. . ( 2020). Early transmission dynamics in Wuhan, China, of novel coronavirus-infected pneumonia. N. Engl. J. Med., 382 : 1199– 1207

[37] Linton,N. M., Kobayashi,T., Yang,Y., Hayashi,K., Akhmetzhanov,A. R., Jung,S. M., Yuan,B., Kinoshita,R. ( 2020). Incubation period and other epidemiological characteristics of 2019 novel coronavirus infections with right truncation: A statistical analysis of publicly available case data. J. Clin. Med., 9 : 538

[38] Lavezzo,E., Franchin,E., Ciavarella,C., Cuomo-Dannenburg,G., Barzon,L., Del Vecchio,C., Rossi,L., Manganelli,R., Loregian,A., Navarin,N. . ( 2020). Suppression of a SARS-CoV-2 outbreak in the Italian municipality of Vo’. Nature, 584 : 425– 429

[39] Nishiura,H., Natalie,M. Akhmetzhanov,A. ( 2020). Serial interval of novel coronavirus (COVID-19) infections. Int. J. Infect. Dis. 93, 284– 286

[40] Du,Z., Xu,X., Wu,Y., Wang,L., Cowling,B. J. Meyers,L. ( 2020). Serial interval of COVID-19 among publicly reported confirmed cases. Emerg. Infect. Dis., 26 : 1341– 1343

[41] LongD. GombarS., HoganC. GreningerA. OReillyS. V., Bryson-CahnC., StevensB., RustagiA., JeromeK. KongC.. ( 2020) Occurrence and timing of subsequent SARS-CoV-2 RT-PCR positivity among initially negative patients. medRxiv, 20089151

[42] Arevalo-Rodriguez,I., Buitrago-Garcia,D., Simancas-Racines,D., Zambrano-Achig,P., Del Campo,R., Ciapponi,A., Sued,O., Rutjes,A. W., Low,N. . ( 2020). False-negative results of initial RT-PCR assays for COVID-19: A systematic review. PLoS One, 15 : e0242958

[43] Kucirka,L. Lauer,S. Laeyendecker,O., Boon,D., ( 2020). Variation in false-negative rate of reverse transcriptase polymerase chain reaction-based SARS-CoV-2 tests by time since exposure. Ann. Intern. Med, 173 : 262– 267

[44] Yang,Y., Yang,M., Shen,C., Wang,F., Yuan,J., Li,J., Zhang,M., Wang,Z., Xing,L., Wei,J. . ( 2020). Laboratory diagnosis and monitoring the viral shedding of 2019-ncov infections. Innovation (N Y), 1 : 100061

[45] Pan,Y., Long,L., Zhang,D., Yuan,T., Cui,S., Yang,P., Wang,Q. ( 2020). Potential false-negative nucleic acid testing results for severe acute respiratory syndrome coronavirus 2 from thermal inactivation of samples with low viral loads. Clin. Chem., 66 : 794– 801

[46] He,X., Lau,E. Wu,P., Deng,X., Wang,J., Hao,X., Lau,Y. Wong,J. Guan,Y., Tan,X. . ( 2020). Temporal dynamics in viral shedding and transmissibility of COVID-19. Nat. Med., 26 : 672– 675

[47] Dugdale,C. Anahtar,M. Chiosi,J. Lazarus,J. McCluskey,S. Ciaranello,A. Gogakos,T., Little,B. Branda,J. Shenoy,E. . ( 2020). Clinical, laboratory, and radiologic characteristics of patients with initial false-negative severe acute respiratory syndrome coronavirus 2 nucleic acid amplification test results. Open Forum Infect Dis, 8 : ofaa559

[48] Kinloch,N. Ritchie,G., Brumme,C. Dong,W., Dong,W., Lawson,T., Jones,R. Montaner,J. Leung,V., Romney,M. . ( 2020). Suboptimal biological sampling as a probable cause of false-negative COVID-19 diagnostic test results. J. Infect. Dis., 222 : 899– 902

[49] Zhu,H., Li,Y., Jin,X., Huang,J., Liu,X., Qian,Y. ( 2021). Transmission dynamics and control methodology of COVID-19: a modeling study. Appl. Math. Model, 89 : 1983– 1998

[50] Liu,Z., Magal,P., Seydi,O. ( 2020). A COVID-19 epidemic model with latency period. Infect. Dis. Model, 5 : 323– 337

[51] Ministryof Public Health. Accessed: December 1, 2020

[52] NajiaHOUSSARI. Lebanon reinstates lockdown measures after virus rebound. Accessed: December 1, 2020

[53] NaylaMadi Masri. The development and state of the art of adult learning and education (ale): National report of Lebanon. National Committee for Illiteracy and Adult Education, Ministry of Social Affairs, 2008. Accessed: December 1, 2020

[54] NajwaYaacoub LaraBadre. Population & housing in Lebanon, 2012. Accessed: December 1, 2020

RIGHTS & PERMISSIONS

The Author(s). Published by Higher Education Press.