key: cord-329357-ujh2nmh5
authors: Ben Miled, S.; Kebir, A.
title: Simulations of the spread of COVID-19 and control policies in Tunisia
date: 2020-05-06
journal: nan
DOI: 10.1101/2020.05.02.20088492
sha: 
doc_id: 329357
cord_uid: ujh2nmh5

We develop and analyze in this work an epidemiological model for COVID-19 using Tunisian data. Our aims are first to evaluate Tunisian control policies for COVID-19 and secondly to understand the effect of different screening, quarantine and containment strategies and the rule of the asymptomatic patients on the spread of the virus in the Tunisian population. With this work, we show that Tunisian control policies are efficient in screening infected and asymptomatic individuals and that if containment and curfew are maintained the epidemic will be quickly contained.

On March 11, 2020, WHO announced that the COVID-19 epidemic had passed the pandemic stage, indicating its autonomous spread over several continents. Since March 22, Tunisia has experienced a turning point and general health containment has begun. Tunisia's strategy of containment and targeted screening corresponds to the first WHO guidelines, the aim being to detect clusters by diagnosing only suspicious persons and then to trace the people who came into contact with the positive cases. This method is now showing its limitations. The mass screening carried out in some countries shows that asymptomatic patients or those who develop only a mild form of the disease may exist in significant numbers. So what is the rule of the asymptomatic patients on the spread of the virus in the Tunisian population and does containment and mass screening strategies are sufficient to control the spread of the virus in the Tunisian population?

In this work, a mathematical epidemiological model for COVID-19 is developed to study and predict the effect of different screening, quarantine, and containment strategies on the spread of the virus in the Tunisian population. This model is more detailed than the classical model (SIR) but it remains very simple in its structure. Indeed, all individuals are assumed to react on average in the same way to the infection, there are In what follows, we present the model and its assumptions. Then we calibrate different parameters of the model based on the Tunisian data and calculate the expression of the basic reproduction number R0 as a function of the model parameters. Finally, we carry out simulations of interventions and compare different strategies for suppressing and controlling the epidemic.

COVID-19 is a respiratory disease that spreads mainly through the respiratory droplets expelled by people who cough. So the transmission is usually direct from person to person. Infection is considered possible even when in contact with a person with mild symptoms. In fact, in the early stages of the disease, many people with the disease have only mild symptoms. It is, therefore, possible to contract COVID-19 through contact with a person who does not feel sick. Subsequently, in this work, we consider susceptible individuals, noted S, who are infected first go through a stage where they are infected but asymptomatic, noted As for unreported asymptomatic infectious. This stage appears to be particularly important in the spread of COVID-19. The individuals then develop symptoms and become symptomatic infectious, so either enter directly into a quarantine stage, noted Q, corresponds to reported symptomatic infectious individuals, or go through a moderate or severe infectious stage, noted I for unreported symptomatic infectious and then can go through the quarantine stage or not. Finally, the infection ends and the individuals are then immunized, denoted R or dead, denoted D. This life cycle can be represented using the following flow chart (1) followed by table 1 that lists the model parameters.

The quarantine is assumed to act on the As and I stages. Indeed, we assumed that the state can detect asymptomatic individuals by for example random screening, and then positive testing ones go into quarantine. Let τ1, respectively τ2, be the quarantine rate for As class, respectively I. We assumed that the asymptomatic individual, As, turn out to be I at a rate β. We further assumed that quarantined individuals, Q and infected individuals, I either die at a rate of µ per unit of time or become recovered/immune, R, at a rate of γ per unit of time.

Finally, we assume that each healthy individual is infected proportionally by As asymptomatic individuals, with a rate of α1 and by I infected individuals, with a rate of α2. As I state is constituted with the moderate or severe state, they are more contagious than the As state, therefore, α1 < α2. Therefore, our model consists of the following system of 6 ordinary differential equations:

With an initial condition at time t = t0 defined as following:

One of the advantages of the basic reproduction number R0 concept is that it can be calculated from the moment the life cycle of the infectious agent is known. We calculate the R0 for our model using the Next Generation Theorem [3] (see section A.1),

3 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted May 6, 2020. . https://doi.org/10.1101/2020.05.02.20088492 doi: medRxiv preprint

The first term of the expression of (2) corresponds to infections generated by asymptomatic types (healthy carrier to mild symptoms). The second term corresponds to secondary infections caused by moderate or severe symptomatic infection. With this expression, we can see that there are several ways to lower the R0 and thus control the epidemic. For example, we can reduce the number of susceptible people (decreasing α1 and α2) by confining the population, reducing contacts, and wearing masks. We can also reduce the rate of contact with an infected person by increasing quarantine rates (τ1 and τ2) by isolating asymptomatic or symptomatic infected persons through mass screening.

The estimation of the different parameters of the model is done in three steps (see section A.2). In the first step, we will estimate the start date of the epidemic, t0, the initial states As(t0) and I(t0) as well as the infection rates α1 and α2. In the second step, we estimate the mortality rate, µ, and the recovery rate, γ. In the third step, we evaluate the parameters τ1 and τ2 by an optimization method. The program is available for download 1 .

We used the Tunisian Health Commission 2 data-set of reported data to model the epidemic in Tunisia. It represents the daily new-cases, death, and recoveries in Tunisia. The first case was detected on March 2, 2020.

To estimate the initial conditions As(t0) and I(t0) and parameters α1 and α2, we fix S0 = 11694720, which corresponds to the total population of Tunisia and assume that the variation in S(t) is small during the period consider. We also fix the parameters β, γ, τ1 and τ2. For this estimation, we adapt the method developed by [4] in our case. Let CRt) the cumulative number of reported infectious cases at time t, defined by,

Let's assume that CR(t) = χ1 exp(χ2t) − χ3 with χ, χ2 and χ3 three positive parameters that we estimate using log-linear regression on cases data (see figure 2 and table 2). 

We obtain the model starting time of the epidemic t0 by assuming that CR(t0) = 0 and therefore equation (3) implies that:

1 https://github.com/MayaraLatrech/covid19_sasymodel.git 2 https://covid-19.tn/ 4 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted May 6, 2020. For now, we assume that α2 = f α1 with f a fixed parameter bigger than 1 and let's τ = τ2/τ1. Then, by following using the approach described in the step 1 of section A.2, we have:

In figure 3 , we plotted the number of individuals in quarantine without curfew predicted by the ODE model, Q(t) i.e. the number of diagnosed cases per day and compared to the Tunisians data. We observe that from the 20th day onwards, the simulated curve deviates from the observed data. This deviation is due to the epidemic control policies put in place between 20 and 25 March (closure of cafes shops and the introduction of a curfew).

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted May 6, 2020. figure 4 we study the effect of the curfew installed by the Tunisian state since March 20, 2020, on the number of infected people reported by the state. Figure 4 shows the effect of two curfew strategies on the dynamics of the epidemic: a 12-hour curfew (the chosen Tunisian policy) and an 18-hour curfew (a more restrictive policy). During the period of curfew, the rate of infection α1 is divided by 100. In Figure 4 (a), it can be seen that, for the chosen policy to maintains a 12-hour curfew for 100 days, the epidemiological peak in terms of the number of reported infected is reached after about 50 days with a value equal to 953. After the peak, we observe a slow decrease in the number of reported infected persons. On the other hand, in the more restrictive case of18 hours curfew, the peak would be reached more quickly after 27 days with a more rapid decrease. These values should be compared with the 774 cases given by [1] and the fact that the epidemiological peak was reached around April 29, 2020, after about 58 days with a reported infected number equal to 975.

We represent on figure 4(b) , the number of deaths by time, it appeared that the peak of the deaths is shifted to the peak of the infected for about 50 days, this shift corresponds to the hypothesis that we made on the duration between the beginning of the symptoms and the deaths ( i.e. 30 days). We note that the simulated death curve overestimates the observed curve. This may be due to the mortality of unreported infected and therefore the surplus corresponds to unreported COVID-19 mortality, which makes the optimization of the model's variables for deaths imprecise. Moreover, the model predicts 228 deaths at the end of the epidemic in the current case (12 hours curfew). In the case where the curfew was 18 hours, the number of deaths would be 84. This information should be taken with caution, because at the time the simulations were made the number of deaths was low.

Finally, we notice that the ratio between the undeclared cases (asymptomatic and symptomatic) represents between 45% at the beginning of the epidemic for less than 10% at the end of the epidemic (see figure 5 ).

We study in this section the effect of more intensive mass screening, i.e. by increasing τ1 and τ2, on the basic reproduction number, R0, and on the number of declared infected, Q (see figures 6, 7).

In figure 6 , we can see that no matter how intensively we screen, we have R0 > 1. Moreover, we note that to minimize the basic reproduction number, it is necessary to increase the massive screening of the class of 6 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted May 6, 2020. is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted May 6, 2020. . https://doi.org/10.1101/2020.05.02.20088492 doi: medRxiv preprint asymptomatic infections. Similarly in figure 7, we can see that the increase in the mass screening effort allows a more rapid decrease in the number of cases (see figure 7(a) ). This is probably since cases are detected earlier, they do not contribute to the contamination of the healthy ones.

Indeed, mass screening has an indirect effect on recruitment in the infect compartment. More specifically, we assume that β + τ1 does not vary when τ1 changes. Consequently, the effort of mass screening on asymptomatic patients cannot exceed this value, and then any additional effort beyond β + τ1 will be passed on to healthy patients and therefore will be useless.

It is observed that the calibration of the model on the Tunisian data using Metropolis-Hastings (MH) algorithm, gives a value of τ1 = 0.78, which represents a very important screening effort. This would prove that the screening strategy is very efficient. However, we didn't had access to the testing campaign methodology that would have allowed us to adjust our estimates.

Moreover, it is observed that the number of deaths at the end of the epidemic varies from 228 to 167, i.e. a 20% decrease in the number of cases (see figure 7 (b)). 

It can be noted that since April, 15, 2020, Tunisia has succeeded in slowing down the speed of propagation thanks to and containment and a curfew. With this work, we suggest that if containment and curfew are maintained, short-term projections could be more optimistic. The fact that the epidemic is quickly contained tends to show that the number of undeclared infected is low, which may suggest that our model is efficient for the evaluation of undeclared cases. In fact, we show that at the time of the epidemiological peak, the number of unreported infected persons constitutes at most 1/3 of the infected population. However, this will need to be confirmed by a field evaluation. Moreover, using Tunisian data, the optimization algorithm fixes the rate at which asymptomatic enter in 8 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted May 6, 2020. . https://doi.org/10.1101/2020.05.02.20088492 doi: medRxiv preprint quarantine, τ1, at a high value. This expresses the good performance of the control policy of the Tunisian government. Indeed, in Tunisia, the control policy consists of an intense isolation campaign targeting sick individuals and their relatives. An effort of testing was carried out in a targeted manner, similar to snowball sampling. This approach enabled to have a major screening effort on infected and asymptomatic individuals. Finally, the model was successfully able to predict the time of the peak at the end of April.

A.1 Computation of the basic reproduction number R 0

We use the next generation matrix to drive the basic reproduction number R0 [3] . In the system (1) we have two infected compartments represented by the second and third equations of the system. Therefore, at the infection-free steady state, i.e. for a small (As, I) and S = S0, the linear epidemic subsystem is :

If we set X = (As, I) T as the vector of infected, F is the matrix that represents the production of new infections and T the matrix of transfer into and out of the compartment by transmission, mortality, quarantine, and recovery, then the matrix form of the linear epidemic subsystem is:

Therefore, the next generation matrix is :

Knowing that the basic reproductive number R0 is the largest eigenvalue of the next-generation matrix, then:

Step 1: In this part, to estimate the initial conditions As(t0) and I(t0)

and parameters α1 and α2 we adapt the method developed by [4] . Let 9 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted May 6, 2020. . https://doi.org/10.1101/2020.05.02.20088492 doi: medRxiv preprint γ1 = γ + µ and CR(t) the cumulative number of reported infectious cases at time t, defined by,

Let's assume that CR(t) = χ1 exp(χ2t) − χ3 with χ, χ2 and χ3 three positive parameters.

By assuming that, CR(t0) = 0 equation (3) implies that:

and then t0 = 1 χ2 (ln(χ3) − ln(χ1)).

Using equation (3), we have also:

Let's note τ = τ 2 τ 1 and H(t) = As(t) + τ I(t). Then we have,

In order to simplify the calculus, we will use the normalized functions, As H and I H . We have:

Rewriting the third equation (1), with H variable, dI dt = βH − (τ2 + βτ + γ1)I

By assuming, that I(t) = I(t0) exp(χ2(t − t0)) and substituting in equation (15), we obtain: χ2I(t0) = βH(t0) − (τ2 + βτ + γ1)I(t0)

Equation (16), implies

By using equation (14) and (17), we obtain:

Let's assume that α2 = f α1 with f a fixed parameter bigger than 1. The parameter α1 is evaluated using As(t) = As(t0) exp(χ2(t − t0)) and the second equation of (1) at t0, we obtain:

χ2

As(t0) H(t0) = α1S0(f I(t0) H(t0) + As(t0) H(t0) ) − (β + τ1) As(t0) H(t0) ⇔ (19) α1 = χ2 + β + τ1

(f I(t 0 ) As(t 0 ) + 1)S0

(20)

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted May 6, 2020. . https://doi.org/10.1101/2020.05.02.20088492 doi: medRxiv preprint and therefore using equations (17) and (18), α1 = (χ2 + β + τ1) ( f β χ 2 +τ 2 +γ 1 + 1)S0 (21)

Step 2 We hereby propose to estimate γ and µ. We notice that, R(t) = γ µ D(t), for all t > 0. Let ρ = µ γ , ρ is estimate using dead and recoveries data.

Let p the fraction of infectious (quarantined or not) that become reported dead ( i.e. 1 − p become reported recovered). Thus ρ = pμ (1−p)γ , with 1/μ the average time to death and 1/γ the average time to recover. Therefore,

Step 3 Parameters τ1 et τ2 was estimated using Metropolis-Hastings (MH) algorithm developed in the pymcmcst python package [5] 

Estimation of Tunisia COVID-19 infected cases based on mortality rate

Infectious Diseases of Humans: Dynamics and Control

On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations

Understanding Unreported Cases in the COVID-19 Epidemic Outbreak in Wuhan, China, and the Importance of Major Public Health Interventions

A Python Package for Bayesian Inference Using Delayed Rejection Adaptive Metropolis