Log In Sign Up

Generic probabilistic modelling and non-homogeneity issues for the UK epidemic of COVID-19

by   Anatoly Zhigljavsky, et al.

Coronavirus COVID-19 spreads through the population mostly based on social contact. To gauge the potential for widespread contagion, to cope with associated uncertainty and to inform its mitigation, more accurate and robust modelling is centrally important for policy making. We provide a flexible modelling approach that increases the accuracy with which insights can be made. We use this to analyse different scenarios relevant to the COVID-19 situation in the UK. We present a stochastic model that captures the inherently probabilistic nature of contagion between population members. The computational nature of our model means that spatial constraints (e.g., communities and regions), the susceptibility of different age groups and other factors such as medical pre-histories can be incorporated with ease. We analyse different possible scenarios of the COVID-19 situation in the UK. Our model is robust to small changes in the parameters and is flexible in being able to deal with different scenarios. This approach goes beyond the convention of representing the spread of an epidemic through a fixed cycle of susceptibility, infection and recovery (SIR). It is important to emphasise that standard SIR-type models, unlike our model, are not flexible enough and are also not stochastic and hence should be used with extreme caution. Our model allows both heterogeneity and inherent uncertainty to be incorporated. Due to the scarcity of verified data, we draw insights by calibrating our model using parameters from other relevant sources, including agreement on average (mean field) with parameters in SIR-based models.


page 1

page 2

page 3

page 4


Mean field game for modeling of COVID-19 spread

The paper presents the one of possible approach to model the epidemic pr...

Epidemics on Hypergraphs: Spectral Thresholds for Extinction

Epidemic spreading is well understood when a disease propagates around a...

Clustering of countries based on the associated social contact patterns in epidemiological modelling

Mathematical models have been used to understand the spread patterns of ...

A stochastic SIR model for the analysis of the COVID-19 Italian epidemic

We propose a stochastic SIR model, specified as a system of stochastic d...

Dynamic Survival Analysis for non-Markovian Epidemic Models

We present a new method for analyzing stochastic epidemic models under m...

Can a latent Hawkes process be used for epidemiological modelling?

Understanding the spread of COVID-19 has been the subject of numerous st...

Age-stratified epidemic model using a latent marked Hawkes process

We extend the unstructured homogeneously mixing epidemic model introduce...

1 Introduction

We model the development of the COVID-19 epidemic in the UK under different scenarios of handling the epidemic. There are many standard epidemiological models for modelling epidemics, see e.g. [1, 2]. In this research, we use a more generic simulation model which has the following features:

  • it can take into account spatial heterogeneity of the population and heterogeneity of development of epidemic in different areas;

  • it allows the use of time-dependent strategies for analysing the epidemics;

  • it allows taking into account special characteristics of particular groups of people, especially people with specific medical pre-histories and elderly.

Standard epidemiological models such as SIR and many of its modifications do not possess these properties and hence are not applicable for any study that requires the use of the features above. In particular, In particular, for influenza the mortality changes less significantly with age in comparison to coronavirus; hence common influenza models do not give much insight in modelling the COVID-19 epidemic. The report [3], which is widely considered as the main document specifying the current epidemic in the UK and in the world, is almost entirely based on the use of standard epidemiological models and hence the conclusions of [3] seem to be lacking specifics related to important issues of the study such as heterogeneity of development of epidemic at different locations and even flexible use of different death rates across different ages. Moreover, as noted in [4], the basic assumptions for the model of [3] were written more than 13 years ago and based on the specific dynamics of a flu pandemic and hence the model was not calibrated for COVID-19. Unreliability of COVID-19 data, including the numbers of COVID cases and COVID deaths, is a serious problem. It is discussed, in particular, by G. Antes, in [5].

Despite the SIR models including the model of [3] have been heavily criticised [4, 6], when choosing the main parameters of our model we calibrate it so that its mean field version approximately reproduces the same output as SIR-based models of [3, 7] and we use the parameter values consistent with  [3, 7]. Also, the main notation (which sometimes does not look natural from a statistician’s point of view) is taken from [7].

The primary objective of our work is construction of a reliable, robust and interpretable model describing the epidemic under different control regimes. In this paper, we make the first step towards this objective. Due to lack of time and unreliability of available data on COVID-19, most of our results serve as an illustration of the role of different control options and we hope that even as illustration they are useful. However, we try to stay close to the COVID-19 epidemic scenario and hence we use appropriate recommendations about the choice of parameters and models of the virus behaviour we found from the studies based on the use of standard models.

See conclusions in Sections 3,4,5 for the main conclusions of this work.

Acknowledgement. We are grateful to several colleagues in Univ. Cambridge and Univ. Warwick for useful comments on drafts of this paper and to M. Hairer (ICL) and J. Ball (Univ. Nottingham) for valuable comments and discussions.

2 The model

Variables and parameters

  • - time (in days)

  • - intervention time (e.g., the time when self-isolation starts)

  • - population size

  • - number of infected at time

  • - number of immune (recovered or dead) at time (this makes full sense from a modelling perspective, but we appreciate that this may look kind of strange when read by non-scientists)

  • - number of susceptible at time

  • - number of groups selected for the study

  • - size of -th group ()

  • - number of infected at time in group

  • - number of immune at time in group

  • - number of susceptible at time in group

  • - average number of transmissions of the virus per unit of time with no intervention

  • - average infectious period

  • - shape parameter of the Erlang distribution

  • - reproductive number (average number of people who will capture the disease from one contagious person)

  • - reproductive number after intervention

  • - probability to recover for an infected person from the

    -th group

Values of parameters and generic model

The reproductive number is the main parameter defining the speed of development of an epidemic. There is no true value for as it varies in different parts of the UK (and the world). In particular, in rural areas one would expect a considerably lower value of than in London. Authors of [3] suggest and as typical; the authors of [7] use values for in the range [2.25, 2.75]. We shall use the value as typical which may be a slightly pessimistic choice overall but could be an adequate choice for the mega-cities where the epidemics develop faster and may lead to more causalities. In rural areas, in small towns, and everywhere else where social contacts are less intense, the epidemic is milder.

We assume that the person becomes infected days after catching the virus, where

has Poisson distribution with mean of 1 week. To model the time to recover (or die) we use

Erlang distribution with shape parameter and rate parameter so that the mean of the distribution is (in simulations, we discretise the numbers to their nearest integers). This implies that we assume that the average longevity of the period of time while the infected person is contagious is 21 days, in line with the current knowledge, see e.g. [8, 9, 10]

. Standard deviation of the chosen Erlang distribution is approximately 12, which is rather large and reflects the uncertainty we currently have about the period of time a person needs to recover (or die) from COVID-19. An increase in

would prolong the epidemic and smaller values of would make it to cause people to be contagious for less time. The use of Erlang distribution is standard for modelling similar events in reliability and queuing theories, which have much in common with epidemiology. We have considered the sensitivity of the model in this study with respect to the choice of parameters and but more has to be done in cooperation with epidemiologists. As there are currently many outbreaks epidemics, new knowledge about the distribution of the period of infection by COVID-19 can emerge soon.

The model varies depending on purposes of the study. The main ingredients are: is the birth-and-death processes and is the associated pure birth processes. The process of transmission is the Poisson process with intensity (time to next transmission has the exponential density ). After the intervention (for ), the Poisson process of transmission for -the group has intensity .

We treat this model as purely stochastic despite parts of it can be written it terms of systems of stochastic differential equations. Despite running pure simulation models taking longer than running combined models, they are simpler and less prone to certain errors.

The number of infected at time as the main quantity of interest

To start with, in Figure 1 we consider an uninterrupted run of an epidemic with in a homogeneous (one-group) population. The starting time of an epidemic is unknown and is even hard to define as the first transmissions of the virus take random and perhaps long times. We start plotting the curves after 0.5% of the population are infected.

In red colour, in Figure 1 and all plots below we plot values of , the proportion of people non-infected by time . In blue, we plot , the proportion of people recovered from the disease (or dead) and in green we plot our main quantity of interest which is , the proportion of infected people at time . The values of do not play any part in modelling and are plotted for information only.

Figure 1: An uninterrupted run of a COVID-19 epidemic in homogeneous conditions

The duration of time while a person with a virus is infected is modelled by Erlang distribution with shape parameter and mean 21. The simulations are flexible and we can easily change these values and plot similar curves with updated data. Despite the full recovery taking slightly longer,

days is a good estimation for the period of time when people with severe infection may require intensive care, ventilator etc. Also, we believe that this is a suitable distribution for the period of time when such people may die.

From the value , we can estimate the distribution of the number of deaths, at time , as follows: first, since on average an infected person is considered as infected for 21 days, we compute . At each particular day, any person can die with the probability which is typical for the chosen population or sub-population. If we consider the whole population, then we apply the mortality rate for the population. However, if refers to a particular group only, then we should apply the corresponding coefficient for the group. In any case, the distribution for the number of deaths at time can be roughly considered as binomial with parameters and the probability of ‘success’ , where is the mortality rate. Different authors disagree on the values of mortality for the COVID-19; see, for example, [11]. UK’s experts believe  [3], WHO sets the world-wide mortality rate at 0.034, the authors of [12] believe is very small and could be close to , an Israeli expert D. Yamin sets , see [13].

Example 1. Assume that the epidemic in Birmingham (with population size m) was running uninterrupted according to the scenario depicted in Figure 1. The maximum value of is giving the maximum expected daily death toll of assuming we use UK experts value .

Although simple to implement, the estimation in Example 1 is likely to be wrong, as it would be for any heterogeneous city (e.g. London). Critically, as will be discussed in the next section, the maximum death toll is likely to be significantly lower than the above calculations suggest.

Summarizing, the expected daily mortality curve is simply a suitably scaled version of of Figure 1. The same is applied to the curves representing expected number of hospital beds and ventilators. The variability for the number deaths is naturally considerably higher than the variability for the numbers of required hospital beds and ventilators as the latter ones are correlated in view of the fact that each person occupies a bed (requires a ventilator) for a few days in a row but dies only once.

Unlike , the number of deaths at each particular day is (approximately) known and potentially can be used for estimating the first and second derivatives of and hence the stage of the epidemic. It is, however, difficult to do for the following reasons: randomness in the values , see above, and, more importantly, heterogeneity of large sub-populations, see next section.

One of the main targets of decision-makers for dealing with epidemics like COVID-19 can be referred to as ‘flattening the curve’, where ‘the curve’ is or any of its equivalents and ’flattening’ roughly means ‘suppressing the maximum’. We shall consider this in the next sections.

3 Spatial heterogeneity of the population and heterogeneity of epidemic development in different areas

The purpose of this section is to demonstrate that heterogeneity of epidemic developments for different sub-populations has significant effect on ‘flattening the curve’.

Already when this paper was completed, the authors learned about a preprint [14] which uses SIR modelling to produce somewhat similar conclusions with more specifics for COVID-19 in the UK. However, SIR modelling using numerous age-classes and many equations for each age-class different across different parts of the country resulting in an astronomical number of parameters; this makes any sensitivity analysis to various parameters uncheckable. We believe that a combination of stochastic models like the present one with standard SIR-based models can help in significantly reducing the number of parameters while also adding more flexibility to a hybrid model.

Consider the following situation. Assume that we have a population consisting of sub-populations (groups) with similar demographic and social characteristics and these sub-populations are subject to the same epidemic which has started at slightly different times. Let the sizes of sub-population be with . In all cases, the curves are shifted in time versions of the curve of Figure 1.

In Figures 3 and 3, and the second epidemic started 50 days after the first one (incubation period is set to be 0). In Figures 5 and 5, and each next epidemic cycle has started 7 days after the previous one. In all cases, we evidence the significant ‘flattening the curve’ phenomenon. In Figures 35, the green line is used for the resulting curves . The maximal values of in the examples depicted in Figures 35 are 0.124, 0.136, 0.134 and 0.132, respectively. This is significantly lower than the value 0.163 for the original curve. In the assumptions of Example 1, these would respectively lead to the maximum expected number of deaths 57.7, 63.3, 62.9 and 62.7.

Extra heterogeneity of different epidemics caused by the social and demographic heterogeneity of the sub-populations would further ‘flatten the curve’. The following conclusions are in line with what the Israeli expert D. Yamin states [13] and agree with the main conclusions of the paper [15] which is rather critical towards standard epidemiological models.

Figure 2: , .
Figure 3: ,
Figure 4: , ()
Figure 5: , , ()

Conclusions. (a) Isolation of sub-populations at initial stages of an epidemic is very important for preserving heterogeneity of the epidemic and ‘flattenning the curve’; this flattening can be very significant. (b) Epidemiological models based on the assumption of a homogeneous population but applied for populations consisting of heterogeneous sub-populations may give completely misleading results.

4 An epidemic with intervention

In this section, we assume that we make an intervention to the epidemic by introducing an isolation at a certain stage. Moreover, we assume that isolation may be different for 2 different groups. We define a special group of more vulnerable people consisting, for example, of all people aged . We define , where is the total population size and is the size of this special group .


be the moment of time when the isolation occurs. It is natural to define

from the condition , where, for example, . For , the virus has been transmitted to people uniformly so that, conditionally a virus is transmitted, the probability that it reaches a person from group is . For , the virus has been transmitted to people in such a way that, conditionally a virus is transmitted, the probability that it reaches a person from is with . Moreover, at time the value of may change to in view of self-isolation. The parameters in this model are:

  • : relative size of the group ;

  • defines the start of the isolation strategy;

  • defines the strength of isolation of the group ;

  • (initial);

  • : the reproductive number after intervention for ; it defines the strength of the overall isolation.

We have run a large number of scenarios coded using a combination of R and Julia [16]; the code is provided in Appendix. In Figures 611, we illustrate a few of these scenarios. In all these scenarios, we have chosen and . The results are robust towards values of these parameters (subject to faster or slower rate of the epidemic in dependence on ). In Section 5 we use to illustrate some specific results; can be considered either as a generic value or, in the contents of Section 5, as the relative size of the group of vulnerable people which is larger than the group of people aged .

We would like to emphasise that there is still a lot of uncertainty on who is vulnerable. It is very possible that the long term effect of the virus might cause many extra morbidities in the coming years coming from severe cases who recover. It is also possible that the virus will mutate and the new strain might affect younger populations more than the current strain. These, and many more possibilities, have a non-negligible probability of occurring, and they have devastating effects. In future models we will try to address these possibilities.

To estimate the number of COVID cases at a given time is difficult [5]. This implies that at the time of making an intervention only very rough guesses about the value of , which is crucial for the future development of the epidemic, can be made.

We distinguish different strengths of isolation by values of parameters and :

The values and mean that under the condition that a virus is infecting a new person, the probabilities that this new person belongs to are and respectively.

Measuring the level of compliance in the population and converting this to simple epidemiological measures and is hugely complex problem which is beyond the scope of this paper.

The lines and respective colours in Figures 611 are as follows.

  • Solid blue line: , the proportion of people recovered from the disease (or dead) in case of no intervention, as in Figure 1.

  • Dashed blue line: in case of intervention.

  • Solid black line: the proportion of people from group recovered from the disease (or dead) in case of intervention.

  • Solid red line: , the proportion of people non-infected by time in case of no intervention, as in Figure 1.

  • Solid orange line: , the proportion of people susceptible to the virus from group non-infected by time in case of intervention.

  • Solid green line: , the proportion of infected people at time with no intervention, as in Figure 1.

  • Dashed green line: , the proportion of infected people at time with intervention.

  • Solid purple line: , the proportion of infected people from group at time with intervention.

Figure 6: , ,

Considerably important, in view of the discussions of the next section, is Figure 6. In this figure, for the initial data of Figure 1, we strongly separate group at the time when 10% of the population is infected. We make no call to the general public for social distancing. In this scenario, the intervention does not considerably change the proportion of infected in the total population but it significantly ‘flattens the curve’ for the group : compare the purple and green lines. The maximum of is 0.086 which is 2.5 times lower than 0.216, the maximum of . Note also the fact that the curve is rather flat for long time, during the main stage of development of the epidemic. Just after the peak of the epidemic in the whole population, peaks; it is caused by the presence of very large number of infected people from the general population.

Figure 7: , , , no incubation period

The scenario which led to Figure 6 is important and hence it has been re-run for a few variations of the model. Figure 7 shows the results of the simulations where we have removed the incubation period. The maximum of is now 0.075 which is more than twice lower than 1.63, the maximum of . Note that in the scenario with no incubation period the whole epidemic is milder as the virus lives longer.

In Figure 8 we use similar scenario as for Figure 6 but we separate group only mildly at the time when 10% of the population is infected. We make no call to general public for social distancing. The curve is ‘flattened’ for the group but in a considerably smaller degree than in the case of strong separation of .

Figure 8: , ,

Figures 9 and 10 illustrate the situation with mild and strong separation of people from complemented with introduction of the social distancing for general public. Interestingly enough, the effect of social distancing for general public gives less benefit than even mild isolation of people from group .

Figure 11 illustrates the scenario with no call to the general public for social distancing but with strong separation of the people from at at a later stage of epidemic, when 20% of the population is infected. The effect is similar to the one observed in Figures 6 and 7 except for the fact that the call for isolating the group came slightly late.

Figure 9: , ,
Figure 10: , ,
Figure 11: , ,

Conclusion. By considering a number of scenarios we have observed an extreme sensitivity of the negative consequences of the epidemic to the degree of separation of vulnerable people. Sensitivity to the parameter measuring the degree of self-isolation for the whole population is less apparent although it is much costlier.

This conclusion is very much in line with recommendations of leading German epidemiologists [17].

5 Consequences for the mortality and impact on NHS

Results of the type presented in Figures 611 can be translated into the language of the expected number of death and expected number of beds required. To do this we extend the observations we made at the end of Section 2 concerning translation of the curve into the curves for the expected number of death and expected number of hospital beds in the UK. We use the common split of the UK population into following age groups:

and corresponding numbers (; in millions) taken from [18]

The death probabilities are given from Table 1 in [3] and replicated many times by the BBC and other news agencies are:

Unfortunately, these numbers do not match the other key number given in [3]: the UK average mortality rate which is estimated to be about As we feel the value of the UK average mortality rate is more important, we have multiplied all probabilities above by 0.732 to get the average mortality rate to be

Defining the group as a union of groups and , we have m and . We then compute the mortality rate in the group and for the rest of population by

For estimating the average number of death in the group and for the rest of population we can the use the formulas

respectively. The average numbers of hospital beds required for two different groups are proportional to these numbers.

We have run a series of scenarios for an UK epidemic without taking into account spatial and social heterogeneity of the society assuming that we would have separated (mildly and strongly) the group of people at the time when 10% of the population were infected with no call to general public for social distancing. We marked as March 23, 2020.

The two scenarios (with and ) we have used for illustrating this technique are different from the scenarios used for plotting Figures 6 and 7 only by the value of . For Figures 6 and 7, we have chosen but for the group of old people in the UK we have .

In Figures 13 and 13 we plot (using the same colours as in Figures 6 and 7) the following curves: (solid green, no intervention on March 23), (dashed green, intervention on March 23) and (purple, intervention on March 23).

In Figures 15 and 15 we plot the curves for the estimated average number of deaths: (solid green, no intervention on March 23), (purple, intervention on March 23), (dashed green, intervention on March 23) and combined (black, intervention on March 23).

We can deduce from Figure 15 (taking into account the extra factor of hospital beds availability) that strong separation of old people alone would have reduced the expected number of death by at least the factor of 2. Another feature of the scenario with is a roughly 50/50 split between the number of deaths in the and groups.

The y-axis in Figures 15 and 15, after multiplication by 140, can be roughly interpreted as hundreds in the London epidemic assuming , on March 23 and homogeneity of the epidemic (as there is about 9m people in London out of 64.1m). As mentioned in Section 3, deaths numbers in other regions of the UK should be expected to be lower.

Figure 12:
Figure 13:
Figure 14:
Figure 15:

We then have run through the scenarios with partial lockdown on March 23 reducing to and with a return, 30 days later, to with and . Results are plotted in Figures 17-19. The style of Figure is exactly the same as for Figures 13-15. If the value (proportion of non-infected people on March 23) happens to be larger than then the second wave of epidemic should be expected to be (perhaps, considerably) larger. If then the second wave will be less pronounced.

All these figures are given for illustration only without any claim on accurate predictions as there is an uncertainty of the outputs towards the choice of several parameters describing the virus characteristics but the sensitivity of the model towards the choice of these parameters is not yet adequately assessed. Figures 21-21 illustrate sensitivity of the scenario of Figures 17-17 with respect to the choice of , the shape parameter of the Erlang distribution defined in Section 2. In these figures, we plot 100 trajectories similar to the trajectories of Figures 17-17 with the values of taken at random in the interval .

Figure 16:
Figure 17:
Figure 18:
Figure 19:
Figure 20: ,
Figure 21: ,

Conclusion. The outputs of the model developed above can be translated into the language of the expected number of death due to the epidemic. This language amplifies the findings of the previous section stating that isolation of a relatively small percentage of population may hugely reduce the death toll of the epidemic.

Some of the next steps and open problems

  • Run many different scenarios to get better understanding of the current situation with the epidemic and what can be done to effectively control it.

  • Continue calibrating the model against existing epidemiological models and new data as it emerges.

  • By using methods of stochastic global optimization [19, 20], teach the model to learn from the data emerging daily, for adapting values of parameters describing the virus and hence the course of the epidemic.

  • Incorporate into the model algorithms for quickest on-line detection of changes [21, 22] to learn about temporal and spatial heterogeneity of the development of the epidemic.

  • Understanding the reproductive number in dependence on the population density in local areas and hence the mixing distribution for needed to combine sub-populations.

  • Understanding

    and other virus parameters as random variables due to mutations.

  • Understanding the risk as a function depending on different factors such as age and social groups.

  • Use of OR models for studying NHS capacity issues under different scenarios and correlated to the percentage of health providers getting infected by the virus.

  • Properly quantify uncertainty for model predictions.

  • Full-scale sensitivity analysis to different parameters.

  • With better understanding of the role of parameters, formulate inverse problems like finding the stage of an epidemic by its early development.

  • Use game theory

    [23] to produce tools for optimal decision making with respect to slowing down the process of epidemic by enforcing spatial isolation and isolation of different sub-populations.

  • Develop a more sophisticated model combining stochastic differential calculus [24], fractional diffusion [25] and elements of direct simulation (the models of fractional diffusion will be used to understand the effects of slowing down of the epidemic).

  • Understand and quantify the differences in epidemic developments across different countries such as UK, USA, Italy and Spain.

  • Studying long-term effects of coronovirus and future mortality from Covid-19.

  • Studying mortality from other diseases due to the capacity/fragility of the health system during the epidemic.

  • Modelling epidemic control strategies based on testing, tracing and isolation.

Appendix: Julia/R code used for computing the scenarios of Section 4

1#### Main script with R and Julia ####
3## Running code will require installation of Julia
6## In the below, specify where Julia is located
7julia <- julia_setup(JULIA_HOME= "/Applications/")
8## Install required packages in Julia and import into R
12## Main simulation
14## Set total time to run simulation
17## Below is the main procedure. Refer to paper for description of what parameters are.
19julia_command("function Main_procedure(Time,R,R_prime,N,c,quantile,alpha,sigma)
20              Time=Int(Time)
21              N=Int(N)
23              ## Set parameters for study
24              beta = R*sigma              ## Transmission coefficient
25              k=3                         ## parameter for Erlang Distribution
26              alpha_prime = c*alpha       ## New alpha
27              beta1 =R_prime*sigma        ## New transmission coefficient if R0 changes
28              Infectious_parameter=1/7    ## Average of incubation period
30              ## indexes of columns in the matrix Population
31              person=1
32              age=2
33              Group=3
34              Susceptible=4
35              Infected=5
36              Recovered=6
37              Time_til_infected=7
38              time_of_next_inf=8
39              time_of_recovery=9
40              previous_infection=10
42              ## Create initial population
44              G1_split = Int(N*(1-alpha))   ## Number of people in first group
45              G2_split = N-G1_split         ## Number of people in second group
47              Population =zeros(Int32,N,10)
48              Population[:,person]=[1:N;]
49              Population[:,Susceptible].=1
50              Population[1:G1_split,Group].=1
51              Population[G1_split+1:N,Group].=2
53              Num_of_ran = Int(round(N*0.005))
54              ran=randperm(N)[1:Num_of_ran]
55              qq=(1:Num_of_ran)/(Num_of_ran+1)
56              temp1=floor.(cquantile.(Exponential(1/Infectious_parameter), qq))[randperm(Num_of_ran)]
57              temp1[temp1.==0].=1
58              Population[ran,Time_til_infected]=temp1  ## Time they will become infected and can pass on virus
59              #Population[ran,]$Infected=0
60              Population[ran,Susceptible].=1
61              temp = floor.(cquantile.(Exponential(1/beta), qq))[randperm(Num_of_ran)]
62              temp[temp.==0].=1
63              Population[ran,time_of_next_inf] = Population[ran,Time_til_infected] + temp  ## Time they pass on virus to someone
64              temp2 = floor.(cquantile.(Gamma(k,1/(k*sigma)), qq))[randperm(Num_of_ran)]
65              Population[ran,time_of_recovery] = Population[ran,Time_til_infected] + temp2 # Time this person recovers when infected
69              ## Create variables for tracking
71              Susceptible_track       = zeros(Time)     ## vector created to track those susceptible in whole population
72              Infected_track          = zeros(Time)     ## vector created to track those infected in whole population
73              Recovered_track         = zeros(Time)     ## vector created to track those recovered in whole population
74              Second_group_track      = zeros(Time)     ## vector created to track those recovered in vulnerable population
75              Second_group_track_susc = zeros(Time)     ## vector created to track those susceptible in vulnerable population
76              Second_group_track_Inf  = zeros(Time)     ## vector created to track those infected in vulnerable population
78              ## Begin simulation
80              time_to_start_isolation=Time
81              for t in 1:Time
82              ##Isolation code. Only initialise isolation when proportion of suceptible people = quantile
84              if (sum(Population[:,Susceptible])/N < quantile) && (t< time_to_start_isolation)
85              time_to_start_isolation=t
86              end
88              ## Index the current people who have the chance to infect others.
90              Current = Population[(Population[:,Time_til_infected].<t) .& (Population[:,Time_til_infected].>Int(0)) .& (Population[:,Recovered].==Int(0)),1]
92              for i in Current
93              Population[i,Infected] = 1               ## Person now becomes infected
94              Population[i,Susceptible] = 0            ## Remove susceptible flag
95              Population[i,previous_infection] = 1     ## Initialise previous infection flag
97              ## If it is time to infect and not time to recover proceed
99              if (Population[i,time_of_next_inf]==t & t<Population[i,time_of_recovery] )
101              if ( Population[i,Group]==1  )
103              ## Before isolation. We take alpha and not alpha prime
104              if (t<time_to_start_isolation)
105              X = ifelse(rand()<alpha,2,1)   ## Sample from which group to infect
106              ## Randomly select one person from that group
107              Y = Int(ifelse(X==1,ceil(G1_split*rand()),G1_split+ceil(G2_split*rand())))
108              ## Check if person has virus previously. Otherwise give virus
109              if ( Population[Y,previous_infection]==0 & Population[Y,Time_til_infected]==0 )
110              temp = floor(rand(Exponential(1/beta)))
111              Population[Y,Time_til_infected]= t+ floor(rand(Exponential(1/Infectious_parameter)))
112              Population[Y,time_of_next_inf] = Population[Y,Time_til_infected]+ ifelse(temp==0,1,temp)
113              Population[Y,time_of_recovery] = Population[Y,Time_til_infected]+ floor(rand(Gamma(k,1/(k*sigma))))
114              end
115              ## isolation initialised. We take alpha prime here
116              else
117              X = ifelse(rand()<alpha_prime,2,1)   ## Sample from which group to infect
118              ## Randomly select one person from that group
119              Y = Int(ifelse(X==1,ceil(G1_split*rand()),G1_split+ceil(G2_split*rand())))
120              ## Check if person has virus previously. Otherwise give virus
121              if ( Population[Y,previous_infection]==0& Population[Y,Time_til_infected]==0 )
122              temp = floor(rand(Exponential(1/beta1)))
123              Population[Y,Time_til_infected]= t+ floor(rand(Exponential(1/Infectious_parameter)))
124              Population[Y,time_of_next_inf] = Population[Y,Time_til_infected]+ ifelse(temp==0,1,temp)
125              Population[Y,time_of_recovery] = Population[Y,Time_til_infected]+ floor(rand(Gamma(k,1/(k*sigma))))
126              end
127              end
130              ## Repeat the same above for Group 1 but for Group 2 (Vulnerable group)
132              elseif ( Population[i,Group]==2  )
134              if (t<time_to_start_isolation)
135              X = ifelse(rand()<alpha,2,1)
136              Y = Int(ifelse(X==1,ceil(G1_split*rand()),G1_split+ceil(G2_split*rand())))
137              if ( Population[Y,previous_infection]==0& Population[Y,Time_til_infected]==0  )
138              temp = floor(rand(Exponential(1/beta)))
139              Population[Y,Time_til_infected]= t+ floor(rand(Exponential(1/Infectious_parameter)))
140              Population[Y,time_of_next_inf] = Population[Y,Time_til_infected]+ ifelse(temp==0,1,temp)
141              Population[Y,time_of_recovery] = Population[Y,Time_til_infected]+ floor(rand(Gamma(k,1/(k*sigma))))
142              end
143              else
144              X = ifelse(rand()<alpha_prime,2,1)
145              Y = Int(ifelse(X==1,ceil(G1_split*rand()),G1_split+ceil(G2_split*rand())))
146              if ( Population[Y,previous_infection]==0 & Population[Y,Time_til_infected]==0  )
147              temp = floor(rand(Exponential(1/beta1)))
148              Population[Y,Time_til_infected]= t+ floor(rand(Exponential(1/Infectious_parameter)))
149              Population[Y,time_of_next_inf] = Population[Y,Time_til_infected]+ ifelse(temp==0,1,temp)
150              Population[Y,time_of_recovery] = Population[Y,Time_til_infected]+ floor(rand(Gamma(k,1/(k*sigma))))
151              end
152              end
153              end
156              ## If it is not time for the person infected to recover, find next time the person will transmit
157              if (t<time_to_start_isolation)
158              temp = floor(rand(Exponential(1/beta)))
159              Population[i,time_of_next_inf]=Population[i,time_of_next_inf]+ ifelse(temp==0,1,temp)
160              else
161              temp = floor(rand(Exponential(1/beta1)))
162              Population[i,time_of_next_inf]=Population[i,time_of_next_inf]+ ifelse(temp==0,1,temp)
164              end
166              ## If it is time to recover, set recovery flag =1 and remove infected flag.
168              elseif ( Population[i,time_of_recovery]<t )
169              Population[i,Recovered]=1
170              Population[i,Infected]=0
171              end
172              end # for i in Current
174              ## Update tracking variables
176              index = t
178              Susceptible_track[index]   = sum(Population[:,Susceptible])/N       ## Keep track of susceptible rate
179              Recovered_track[index]     = sum(Population[:,Recovered])/N         ## Keep track of recovered rate
180              Infected_track[index]      = sum(Population[:,Infected])/N          ## Keep track of infected rate
181              Second_group_track[index]  = sum(Population[G1_split+1:N,Recovered])/G2_split   ## Keep track of recovered rate for vulnerable
182              Second_group_track_susc[index] = sum(Population[G1_split+1:N,Susceptible])/G2_split ## Keep track of susceptible rate for vulnerable
183              Second_group_track_Inf[index]  = sum(Population[G1_split+1:N,Infected])/G2_split  ## Keep track of Infected rate for vulnerable
185              #global t=t+1
186              end # for t
188              Output     =  zeros(Time,8)
189              Output[:,1] = Susceptible_track
190              Output[:,2] = Recovered_track
191              Output[:,3] = Infected_track
192              Output[:,4] = Second_group_track
193              Output[:,5] = Second_group_track_susc
194              Output[:,6] = Second_group_track_Inf
195              Output[:,7] .= time_to_start_isolation
196              Output
197              end")
199# Set the parameters for the study
209Scenario_1  =  julia_call("Main_procedure",Total_time,R,R_prime,N,c,quantile,alpha,sigma)
210## Change c for different scenarios
212Scenario_2  =  julia_call("Main_procedure",Total_time,R,R_prime,N,c,quantile,alpha,sigma)
217##  Plots for scenario 1
219lines(seq(0,Total_time-1,1),Scenario_1[2,],type=’l’,col=’blue’,lwd=3 ,lty=1  )
220lines(seq(0,Total_time-1,1),Scenario_1[3,],type=’l’,col=’green’,lwd=3 ,lty=1  )
221abline(v= Scenario_1[7,],lwd=2)
223##  Lines for scenario 2
227lines(seq(0,Total_time-1,1),Scenario_2[4,],type=’l’,col=’black’,lwd=3,lty=1  )
228lines(seq(0,Total_time-1,1),Scenario_2[5,],type=’l’,col=’orange’,lwd=3,lty=1  )
229lines(seq(0,Total_time-1,1),Scenario_2[6,],type=’l’,col=’purple’,lwd=3,lty=1  )