A mathematical model for simulating the spread of a disease through a country divided into geographical regions with different population densities.

The SIR (susceptible-infectious-recovered) model is a well known method for predicting the number of people (or animals) in a population who become infected by and then recover from a disease. Modifications can include categories such people who have been exposed to the disease but are not yet infectious or those who die from the disease. However, the model has nearly always been applied to the entire population of a country or state but there is considerable observational evidence that diseases can spread at different rates in densely populated urban regions and sparsely populated rural areas. This work presents a new approach that applies a SIR type model to a country or state that has been divided into a number of geographical regions, and uses different infection rates in each region which depend on the population density in that region. Further, the model contains a simple matrix based method for simulating the movement of people between different regions. The model is applied to the spread of disease in the United Kingdom and the state of Rio Grande do Sul in Brazil.

Introduction

The global spread of infections with the recently emerging coronavirus SARS-CoV-2, henceforth denoted COVID-19, was declared a pandemic by the World Health Organisation on March 11, 2020. The risk of infection is strongly dependent on individual behaviour and can be reduced by following simple rules such as distance keeping recommendations, hygiene with alcohol gel and wearing of a face mask. In addition to the rules there are some factors that play a further role in the spread of the infection such as the regional population density distribution and circumstances such as living conditions. Contact with other individuals, both in the private sphere with family members and friends, and in the professional environment can lead to an increased risk of transmission of the infection both in- and outdoors, even over a distance larger than . Depending on the country, national and/or regional political decisions for limiting the spread of COVID-19 were based on risk assessments obtained by considering the number of and trends in reported cases and in accordance with the appropriate national Infection Protection Acts. A further source of information for decision making came from monitoring the proportions of the populations with mild, severe and fatal outcome after becoming infected with the disease and evaluating long-term consequences of the pandemic.

Statistics and findings from model calculations provided some information relevant for the spread of the disease, such as the incubation time being between 0 and 14 days and that the transmission time was found to be in an interval between 5 and 7 days. However, a method for determining what are the conditions for a mild or severe outcome of the infection has remained elusive (Famulare 2020; Ikejezie 2019). One of the criteria for measuring the spread of the disease established in the literature is the number of people infected by a previously infected person, known as the base reproduction number R. This estimates how fast the disease spreads and whether political decisions and measures are sufficient to limit its growth of the disease (Coburn et al. 2009; Delamater et al. 2019; Egger et al. 2017; Ferguson et al. 2006; Guerra Fiona et al. 2017; Kucharski et al. 2015; Qun et al. 2020; Lipsitch et al. 2003; Milligan and Barrett 2015; Nishiura and Chowell 2014; Julien et al. 2020; Steven et al. 2020; Thompson et al. 2019; Truelove et al. 2020; Wallinga and Teunis 2004; Wu et al. 2020). For the situation turns supercritical and the total number of infected individuals rises exponentially and without control. In order to be able to predict the time evolution of the pandemic we can make use of mathematical models that were developed for other epidemics, and the following describes some of the more common models. The SI model considers the spread of a contagion without recovery, the SIS model takes into account the spread of an infectious disease without a build-up of immunity, the SIR model considers the spread of the infectious disease together with immunity response and the SEIR model simulates the spread of the contagious disease an immunity response and an initial period where infected people are not immediately infectious. It is noteworthy that most models are simply based on the number of individual who are infected, are immune or have died without considering further details such as differences in transmission probabilities within the regions due to different population densities in different regions, or the movement of individuals between regions.

The model presented in this paper addresses these issues by developing a mathematical model of how a disease spreads through a population which is distributed over multiple geographical regions and with a different population density in each region. There is evidence in the literature that the population density influences how quickly a disease can spread through a population. Alirol et al. (2011) discuss how urbanisation affects the spread of an infectious disease as the global population becomes more concentrated in large cities. Carl-Johan (2015) discussed the challenges presented by the more rapid spread of infectious diseases as the global population becomes more urbanised and are further discussed by Reyes et al. (2013). These, and more recent studies, indicate that any model of an infectious disease through a population must take the differences in population density between urban and rural regions into account.

Mathematical models for simulating the spread of a disease through a population in a single geographical region have been developed since the early 20th century (Kermack et al. 1927). Most of these have been based on the SIR (susceptible-infectious-recovered) models which utilise a system of differential equations to describe the number of individuals in each SIR category. Variations of the SIR model can simulate phenomena such as deaths from the disease, or the number of people who are exposed to the disease but not yet infectious (See Busenberg and Vandendriessche 1990; Allen 1994; Cha et al. 1998; Earn et al. 2000; Greenhalgh et al. 2016; Zhang and Zhou 2019 for example). A summary of the SIR model and its variations is given in Hethcote Herbert (2000).

Models of the spread of a disease through a country or state divided into a number of different regions have been developed in recent years. Liang and Bian (2010) present a statistical model of the spread of the influenza virus through an urban environment which accounted for the way in which individuals move through the urban area considered. Rakowski et al. (2010) presents an individual based statistical model of the spread of influenza in Poland. Lau et al. (2017) present a spacial-temporal model for simulating the spread of the Ebola virus in West Africa, and critically compare their model to an SEIR (susceptible-exposed-infectious-recovered) model. However, none of these models are based on an SIR differential equation model of the spread of a disease.

Yin et al. (2020) present a SIR model which includes the movement of people between a number of cities and which includes the spread of disease from one city to another. However, this paper only considers movement of people between different cities which have high population densities but does not include the surrounding rural areas which have low population densities, and there is observational evidence that a disease will spread faster in dense urban population when compared to the spread though a sparse rural population. A stochastic SIRS (susceptible-infectious-recovered-susceptible) model which includes the mobility of the population has been developed by Wanduku (2017).

This paper presents a mathematical model for simulating the spread of a disease through a country or state that is divided into a number of geographical regions. We derive a system of differential equations that can be solved for the number of people in each category in each region where the infection rate is inversely proportional to the area of each region. This means a region with a small area will have a larger infection rate than one with a large area, but if both regions initially have approximately the same number of susceptible and infectious people then the disease will spread quicker in the smaller region due to its larger population density. The model will also simulate the movement of people between the different regions by including a matrix term in the differential equations which will account for the proportion of the population in one region which move to each of the other regions. This matrix can be set up such that the total population of each region remains the same, and that the proportion of people who move to an adjacent region is much greater than the proportion who move to a region which is further away.

Mathematical modelling

Population kinetics model

Consider a single region when an individual within the population can be classified as susceptible (never had the disease), infected with the disease, recovered from the disease (and who is assumed to be immune) or having died from the disease. Let S, I, R, and D to denote the total number individuals in the population who are susceptible, infected, recovered or who have died respectively, then Heathcote (2000)where is the infection rate, is the rate at which infectious people recover and is the rate at which infectious people die. We note that the sum of the right-hand sides of the differential equations in (1) is zero, indicating that there will be no change in the size of the population.

However, this model is not applicable to diseases where people may be infected with the disease but have no symptoms and so are unaware that they are infected. Here, we will refer to these people as carriers. Hence, the spread of a disease may be more rapid and widespread in a population than is indicated by the number of infected people. This can be incorporated into the basic SIRD model with the inclusion of an additional class, C, of carriers who are either infected with the disease but do not have any symptoms or who have very mild symptoms and so do not get tested for the disease. It is important to note and emphasize here that both infected individuals and carriers can infect other people as both are infected with the disease and the only difference is that carriers are unaware that they are infected. Including the carrier category leads to the modified system of differential equationswhere is the rate at which the carriers recover from the disease and cease to be infectious, is the death rate for carriers and is the proportion of the population who become infected with the disease, have symptoms and are diagnosed as having the disease. It follows is the proportion of the population who become infected with the disease but who are not diagnosed as having the disease, either because they do not have symptoms, they have mild symptoms and mistake the disease under consideration for another disease, or they simply refuse to get a diagnosis.

Systems of differential equations, such as (1) and (2) can be applied to model the evolution of a disease through an entire population. However, these models do not make distinction between rural and urban areas which can have widely differing population densities and which may in turn affect how quickly the disease is transmitted. A disease which is spread by person-to-person contact will spread much more quickly in an urban area as each individual in an urban area will come into close contact with many more other people than someone who is living in a sparsely populated rural area (see Alirol et al. 2011; Reyes et al. 2013; Carl-Johan 2015 for example). In order to consider this aspect we propose a new model which describes the spread of a disease through different geographical regions of a country. An example which shows the importance of considering population density rather than just the population is given in Sect. 3 which discusses the numerical results.

Regional kinetics model

A large population may be divided into a number of smaller populations according to the geographical region in which the people live. In this case it is likely that the rates at which people recover or die from a disease is the same in every region but the infection rate may be different in each region and may depend on quantities such as population density.

The model will also need to include the effect of people moving between the different regions. Here we will assume that the movement of the people does not produce any change in the population of any one region. That is, the number of people who leave region A will be the same as the number of people who enter region A.

Assume that the country or state under consideration is divided into N geographical regions, where it is also assumed that if the country is divided into small enough regions any inhomogeneities in the population distribution within each region can be neglected. Let , , , and denote the number of people who are susceptible, infected, carriers, recovered or who have died in the ith region respectively; and let , , , and denote the vectors which list the values of the corresponding category in each region. For simplicity of notation, introduce the vector valued function defined byThe differential equations (2) for a single region can be extended to the multi-region case and can be expressed aswhere T is a matrix which models how the population moves between the different regions. Note that since the people who have died cannot move between regions there is no transport term in the last equation.

In the model presented here the matrix T is constructed so that the total population of living people in each region remains the same, and the total number of people in each category is not changed by the transport terms. Let M be a symmetric matrix with zeros on the diagonal and where the off diagonal element gives the proportion of people who move from region i to region j. In the work presented herewhere is the geographical distance between regions i and j; is the maximum distance that people move from their original location and is a scaling parameter. However, there are some exceptions to this as there can be large numbers of people moving between cities which are a long way apart. In such cases we can simulate the greater movement of people by setting the distance between the appropriate regions to be smaller than the geographical distance. The converse is also true. There can be parts of a country where fewer people move between the different regions. For example, there may be islands which are not served by a ferry every day and which have no airport or landing strip. In these cases we can reduce the simulated movement of people by setting the distance between the appropriate regions to be bigger than the geographical distance. We note that generally only moderate values of should be used in (4). In most countries and states only a small proportion of the population of a region moves into an adjoining region and even fewer move to regions that are further away.

Let P be the diagonal matrixThen the matrix T appearing in the differential equations (3) is defined aswhere the notation denotes the (i, j) element of the matrix product PM and the normalisation is the sum over the N regions. It is important to note at this point that the matrix T is not constant in time since it is formed from the matrix P which will change as the number of people in each category change.

However, a problem with the system of differential equations (3) is when two regions have similar populations but significantly different areas. For example, if a state has two regions with the same population where one region is twice the area of the other we would expect the disease to spread more rapidly in the region with the smaller area as the population density is bigger. A better approach to modelling the spread of a disease in a multi-region state is to consider the population densities rather than the populations.

Let be the area of the ith region and letbe the densities of the people who are susceptible, infected, carriers, recovered and who have died in each of the N regions. Assume that the population densities rather than the actual populations will satisfy a system of differential equations similar to (3). That iswhere A is the diagonal matrix with areas of the regions on the diagonal and . Replacing by , by and so on in (5) givesMultiplying both sides of the equations in (6) by A givesIt is worth noting that the system of equations (7) is the same as the system (3) except that the infection rate in each region is now inversely proportional to the area of the region. Hence the spread of the infection will be slower in a region with a smaller population density.

Time integration

The system of ordinary differential equations (7) have to be integrated through time using a numerical method. Let , , , , and denote the quantities at the nth time-step. We will assume that the initial values of the quantities (denoted by subscript 0) are known, although the exact values will depend on the country under consideration and the scenario for how the spread of the disease starts.

Here we apply the (Crank and Nicholson, 1947) to solve the system of differential equations (7) which uses the known quantities at the current time step (which have subscript n) to calculate the unknown quantities at the new time-step (which have subscript ). However, this leads to an implicit set of equations and so we adopt the iterative predictor-corrector schemewhere h is the time-step and the superscript [i] on the unknown quantities denotes the corrector iteration.

The iterative process is stopped whenwhere(and the other terms on the numerator of (8) are similarly defined) and is some predetermined accuracy level. Additionally, the iterations will also be stopped when the maximum number of iterations is reached. The choice of norm to use in (8) is not significant and in the the results presented here we have used the uniform or -norm with . For further details of predictor-corrector schemes see any text on numerical methods, such as (Atkinson 1989).

The accuracy of the calculated solution is further controlled by using an error checking time-stepping algorithm. The solution at each new time-step is calculated using a sequence of time-steps until the solutions calculated using two consecutive values of h agree to a predetermined accuracy. In the examples presented in this work we use (giving an initial time-step that is equivalent to one hour) and . The difference in the solution using two different time-steps is measured using (8) where the index i in (8) is now the same index as the one used to denote the different values of h and when considering the differences in the solution for different values of h we used as the stopping criteria.

The initial conditions used here are that on day 0 the total population and the number of infected people in each region are known, and that initially no one has recovered or died from the disease. The initial number of carriers is computed from the initial number of infected people usingwhere denotes the nearest integer. The initial number of susceptible can then be calculation in each region usingwhere is the vector of the total population in each region.

Numerical results

In this section we present the results of using our model to simulate the spread of a disease through country or region. Most of the results are for the United Kingdom and the state of Rio Grande do Sul in Brazil, although the first example is an artificial case which illustrates the importance of using the population density in the equations rather than just the population. In all of the numerical results presented here we have assumed that and that .

Initial example

This first example considered here illustrates the importance of using the population density in the modelling process. Consider the spread of a disease through the four unconnected regions detailed in Table 1. As the regions are unconnected the matrix T in (7) is the zero matrix meaning that there is no movement of people between the regions. If we only consider the size of the population in the model then results for Regions A and D (and Regions B and C) would be identical as they have the same population. Figure 1 shows that the percentage of people infected with the disease over time in each region if the population density is considered instead of the population. In this case, the results for Regions A and C are the same as these have the same population densities. Region B has an earlier and higher peak number of infections as the population density of this region is higher although the population is the same as for Region C. The peak number of infections for Region D is smaller and the peak occurs later as this region has a smaller population density.

Table 1

Example regions for illustrating the importance of considering population density

		Initial		Population
Region	Population	Infected	Area (km^2)	Density (km^-2)
A	1,000,000	10	1000	1000
B	2,000,000	20	1000	2000
C	2,000,000	20	2000	1000
D	1,000,000	10	2000	500

Fig. 1

The percentage of the population infected with a disease in each of the regions. Note that the results for Regions A and C are superimposed

United Kingdom

The first location considered in this work is the United Kingdom (UK) excluding Northern Ireland. The UK is divided into 143 regions which represent the local government authorities except in the large cites (such as London, Manchester and the West Midlands for example) where the small local authorities have been merged to form single metropolitan areas. For the UK we assumed that the initial conditions are that there are two people infected with the disease and they are located in London.

The first set of results explore how varying the parameters , and affects the number of infections and number of deaths and the results are summarised in Table 2. When considering the UK we used , which corresponds to it taking 14 days for an individual to recover from the disease, and corresponding to a death rate of 1% of those who are infected. Finally, for the UK we used in equation (4) for determining the proportion of people moving between the different regions.

Table 2

The simulated maximum number of infected people and total deaths for the United Kingdom using different values of the parameters

β	α	λ	Peak	Day of peak	Died
			Infected	Infected	Day 400
0.05	0.00	4×10^-5	132,676	84	52,352
		6×10^-5	194,554	51	54,556
		8×10^-5	234,478	37	54,951
	0.05	4×10^-5	190,778	112	213,851
		6×10^-5	294,870	69	271,762
		8×10^-5	375,682	53	301,174
	0.10	4×10^-5	218,424	136	221,022
		6×10^-5	367,931	76	280,718
		8×10^-5	470,074	55	309,754
	0.25	4×10^-5	174,314	235	197,634
		6×10^-5	435,442	99	288,841
		8×10^-5	602,362	65	324,306
0.10	0.00	4×10^-5	264,309	89	104,583
		6×10^-5	388,158	53	109,081
		8×10^-5	467,941	39	109,894
	0.05	4×10^-5	379,269	117	422,326
		6×10^-5	587,626	72	539,801
		8×10^-5	749,180	55	599,518
	0.10	4×10^-5	433,199	143	437,479
		6×10^-5	732,748	79	558,295
		8×10^-5	937,288	58	617,102
	0.25	4×10^-5	343,042	249	388,073
		6×10^-5	865,607	103	574,721
		8×10^-5	1,200,474	67	646,436
0.25	0.00	4×10^-5	653,825	95	260,510
		6×10^-5	966,603	57	272,458
		8×10^-5	1,169,077	41	274,672
	0.05	4×10^-5	930,723	125	1,020,882
		6×10^-5	1,453,214	75	1,323,300
		8×10^-5	1,856,550	58	1,478,300
	0.10	4×10^-5	1,055,557	154	1,061,052
		6×10^-5	1,807,695	83	1,371,879
		8×10^-5	2,321,189	60	1,524,244
	0.25	4×10^-5	815,705	271	919,957
		6×10^-5	2,125,866	110	1,414,053
		8×10^-5	2,968,683	71	1,599,267

The results in Table 2 shows that as increases the peak in the number of people infected and the number of people who have died from the disease after 400 days also increases. This is as expected as is the infection rate and so increasing should produce an increase in the number of infections. We can also see that increasing without changing any of the the other parameters causes the day on which the peak in the number of people who are infected to become earlier in the epidemic.

The results in Table 2 also show that as increase so does the peak in the number of infected people and the long-term number of deaths. This is to be expected as increasing means that a larger proportion of the population is being diagnosed as being infected with the disease rather than being just an undiagnosed carrier. Increasing delays the day on which the number of people who are infected reaches its peak.

These results show that when , which controls the proportion of the population of one region which moves to another region, is zero then the number of deaths after 400 days from the disease is greatly reduced although the reduction in the maximum number of infected people is not so large. However, the results also show that the biggest change in the number of deaths is when is increased from 0 to 0.05 and that further increases in does not produce such large changes in the number of deaths, and in some cases the number of deaths after 400 days decreases as increases.

The percentage of the population in each region who are infected or who are carriers on the day on which the number of people who are infected or are carriers reaches its peak is shown in Fig. 2 for different values of when and . The day refers to the day on which the number of people who are infected or who are carriers reaches its maximum value. As expected, when (Fig. 2a) London is the only region with infected people (or carriers) since in this case none of the population moves between the different regions. When (Fig. 2b) the infection spreads to the other metropolitan area since these have good transport links with each other and it was assumed that a proportionally large number of people travel between these regions compared to the more rural areas. As we further increase (Fig. 2c, d) the infection spreads more into the rural areas. It is worth noting that in these cases the proportion of the population who are infected or who are carriers in the urban areas decreases as increases.

Fig. 2

The percentage of people in each region who are infected or who are carriers on the day when the number of people who are infected or who are carriers reaches its maximum for each value of . Here and

Rio Grande do Sul, Brazil

The other geographical location considered in this paper is the State of Rio Grande do Sul in Brazil. For the simulations the state is divided into its 35 microregions. The initial conditions were that there were 3 cases in Porto Alegre, and depending on the value of the corresponding number of carriers in Porto Alegre has been calculated using (9). It is assumed that both the number of cases and the number of carriers are zero for all the other microregions in Rio Grande do Sul. When considering Rio Grande do Sul we used in equation (4) for determining the proportion of people moving between the different regions.

The first set of results that we present investigates how the parameter , which gives the proportion of infected people who have been diagnosed as having Covid-19, affects the values of the other parameters when approximately fitting the model to the observed data. Table 3 gives the values of the parameters in the different cases considered where we used in (4) to simulate a large proportion of the population moving between the different microregions. Figure 3 compares the predicted number of cases for each value of with the observed number of cases in Rio Grande do Sul, and Fig. 4 gives the corresponding comparison for the number of deaths. These figures show that the results of the simulation are broadly the same as the observed data.

Table 3

Values of the fitted parameters for different values of

Parameter	β=0.001	β=0.01	β=0.1
λ	3.50×10^-4	6.30×10^-4	1.52×10^-3
μ_IR=μ_CR	0.020	0.190	0.742
μ_ID	0.0017	0.0020	0.0020

Fig. 3

The observed and simulated number of cases for Rio Grande do Sul using different values of the parameter

Fig. 4

The observed and simulated number of deaths for Rio Grande do Sul using different values of the parameter

It can be seen that as the parameter increases, both the infection rate given by and the recovery rate given by and also increase. There is observational evidence that says that the recovery time for Covid-19 is around 14 days, which would give a recovery rate of 0.0714. Clearly the values of which do this are between 0.01 and 0.001. This is implying that less than one percent of people who have Covid-19 are actually being diagnosed with the disease.

The second set of results investigates how the spread of the disease is affected by the mobility of people between the different microregions. The proportion of people who move between different microregions is proportional to the parameter in (4). Here has been used since the results discussed above demonstrate that the simulated numbers of infections and deaths are closest to the observed values for this . The other parameters that were use are , , and .

Figure 5 shows the percentage of the population in each microregion that are infected on days 45, 90, 135 and 180 for different values of . These results clearly show that as the proportion of people who move between microregions increases then the disease is spread over a larger geographical area, as expected. The results also show that increasing the mobility of the people can delay when the maximum occurs. The first column of maps in Fig. 5 show how the geographical spread of the disease is slower when the population moves around less. The results in the first column of Fig. 5 show that the microregion Campanha Ocidental (the most westerly microregion), which is a rural region with the 6th lowest population density in Rio Grande du Sol, has only a small percentage of infections and carriers showing that the disease is slow to spread in regions with a low population density, as expected.

Fig. 5

The simulated percentage of people who are infected or who are carriers in each microregion on days 45, 90, 135 and 180

Conclusions

Motivated by the actual pandemic, in the present work we have developed a new mathematical model to simulate the spread of a disease taking into account realistic geographical domains with subdomains and their associated population densities. The work presented here is different from most existing models which describe members of a population as being divided into four groups (susceptible, infected, recovered or dead) and which cover the whole country or state. The present approach is based on the population being divided into local regional populations, and each of these are further divided into five local groups, namely the numbers of susceptible, infected, recovered, carriers and dead in each local region. Further, unlike previous models which have used a single infection rate for the whole country or state, this present model uses a different infection rate in each local region that depends on the geographical area of that local region. In each local region the infection statistics are computed using kinetics, whereas the spread of the disease between regions was implemented by a detailed balance diffusive transport matrix, with a mobility parameter which controls the number of individuals traveling between regions and a distance rule defining the probability of an individual traveling between two regions. One academic example was analysed in order to justify our reasoning with respect to the model structure, two realistic scenarios were simulated, the first one for the UK and the second one for the southernmost Brazilian state Rio Grande do Sul.

The results for the small example with four isolated regions show that it is important to use the population densities in the calculations rather than just the populations. If there are two regions with the same size of population then the disease will spread more slowly through the region with the larger area, whereas using a model based purely on population size would predict that the same rate would happen in both regions. Further, if one of the regions becomes big enough then the disease will simply decline in that region since the interactions between the susceptible and infected populations becomes small. Further, the results for both the UK and Rio Grande du Sol show that the population density does play a significant role in how rapidly a disease spreads through a population. The results for the UK shown in Fig. 2 show that on the day when the number of people who are infected or carriers peaks the largest percentages of infected people and carriers are in the large cities, such as London, Birmingham and Manchester, whilst the percentages are much lower in the more sparsely populated rural areas. The results for Rio Grand du Sol given in Fig. 5 show that the disease spreads quickly though the densely populated region around Porto Alegre with a large percentage of the population becoming infected or carriers whilst the percentage of infected people is much lower in the more rural areas.

The results for both the UK and Rio-Grande do Sul show that the key parameter that affects the magnitude of both the number of people who are infected and the number of deaths is the proportion of infected people that are diagnosed with the disease, although we have made the assumption that no carriers die and that they all recover. Increasing the infection rate causes the maximum number of people who are infected to increase and also the number of long term deaths to increase. Also, as the infection rate increases the day on which the number infected people reaches its peak becomes earlier. As the parameter (which determines the proportion of the population in each region who travel to another region) increases then generally the peak in the number of people who are infected increases and the day on which the maximum number of infected people occurs is later.