THE MONTE CARLO

DIFFUSION OF DISEASES

SIMULATION PROGRAM:

USER HANDBOOK

Daniel A. Griffith

Department of Geography

Syracuse University

Syracuse, NY 13244-1020

October, 1999

INTRODUCTION

The Monte Carlo diffusion of diseases simulation program permits concepts from geographic and spatial interaction theories to be explored with numerical examples. This type of experimentation focuses on the importance of the roles played by intervening distance and opportunities on spatially variant origin and destination size. The reader is referred to Hagerstrand (1965) and Forster (1973) for a formal discussion of this particular experimental approach to geographical analysis. More general discussions appear in Dalton (1972), Choguill (1973), Walford (1973), and Cole (1975).

Results from experiments employing this computer program will tend to be those that may be anticipated from geographic theory. The objective, then, simply is to permit a student to both demonstrate and observe numerical manifestations of the spread of selected diseases, within the context of either purely hierarchical flows or distance constrained gravity type flows. Consequently, in order to make the experimentation more understandable, first a demonstration execution should be inspected. Next, a sufficient number of time periods should be evaluated—preferably the disease should die out—for anexecuted experiment, and perhaps several simulation replications examined. Of note is that a space-time data set is generated for any given realization (i.e., replication) of a spreading disease.

Intervention strategies can be evaluated during an epidemic—e.g., vaccination, dispensing of antibiotics, and hospitalization. Costs of intervention, including deaths, are calculated. Because random numbers are generated from the seconds of the computer clock, a stochastic element is built into the analysis. Nevertheless, this version of the diffusion simulation game is based upon average tendencies for selected diseases.

This handbook is divided into five sections. The first section presents a brief history of the Hagerstrand Monte Carlo diffusion simulation program. The second section outlines principal assumptions of the model utilized here. The third section reviews some basic concepts of spatial diffusion. The fourth section summarizes how to use this computer program. Finally, the fifth section discusses the graphical and tabular output generated by this program.

HISTORY OF THE ALGORITHM

The first version of the Hagerstrand Monte Carlo simulation experiment program was proposed by Hagerstrand himself. An operationalized version was developed at The Pennsylvania State University in the late 1960s, and revised at Virginia Polytechnical Institute and State University in the early 1970s. Further modifications of this program occurred at Ryerson Polytechnical University in 1976, followed by refinements being made to it at SUNY/Buffalo in 1983 and Syracuse University in 1994. This current version is a variant that has been modified specifically for diffusion of diseases. All of these versions were written in FORTRAN.

PRINCIPAL ASSUMPTIONS

Five prominent assumptions are made in the simplest diffusion simulation work. First, the likelihood of a diffusion materializing is cast as being directly proportional to the population of the destinations (a hierarchical component), and inversely proportional to the extent of intervening distance (a contagion component). Hence, the type of interaction involved may be described by a gravity model.

The surface over which interaction occurs is isotropic, so that movement is equally easy in all directions over the planar surface in question. This assumption precludes the presence of channels funneling and/or accelerating diffusion, or barriers to diffusion.

The disease carriers (i.e., interactors) are assumed to make a single interpersonal contact during each time period. This contact is the result of a random process, in which every actor on the plain has the same probability of being interacted with. Therefore, each interactor is assumed to randomly select someone with which to interact. These single, equally likely interactions then are assigned to locations in proportion to the corresponding gravity model probabilities associated with the locations, where contact either is with an infected or an uninfected (e.g., susceptible) inhabitant. This procedure works the same way that a lottery does, where each ticket has the same chance of being selected, and individuals' probabilities are a function of the number of tickets they hold.

Finally, a saturation level less than the total population may be specified. Moreover, only a portion of the population may be susceptible to a given disease.

Basic Concepts of Spatial Diffusion

The model around which this computer program was built is based upon the notion of a mean information field (MIF). A MIF defines the probability of contact between a given origin and each destination in a geographic landscape. Let spatial interaction be characterized by the gravity model, such that

I_ij = G P_i^" P_j^ß exp(_ * d_ij) , (1)

where I_ij = the magnitude of spatial interaction between origin location i and destination location j,

P_i = the population of location i,

d_ij = the distance separating origin i and destination j, and

G, ", ß, * = parameters.

Then for any given origin i, the probability of interaction between this origin and any destination j is defined as

P_ij = G P_i^" P_j^ß exp(_ * d_ij)/[E G P_i^" P_j^ß exp(_ * d_ij)]

= G P_i^" P_j^ß exp(_ * d_ij)/[G P_i^" E P_j^ß exp(_ * d_ij)]

= P_j^ß exp(_ * d_ij)/[E P_j^ß exp(_ * d_ij)] , (2)

where N is the number of locations. In deriving equation (2), since origin i is fixed, P_i^" can be factored out of the summation expression appearing in the denominator, as can the constant of proportionality G. These two constant terms then appear in both the numerator and the denominator, and hence divide out.

The set of probabilities {P_ij|i fixed} constitute the MIF for origin i. This MIF will change as origin i changes. Thus, if there are N locations, then there are N MIFs.

These aforementioned features of a MIF are illustrated by the example presented here. Consider a geographic landscape consisting of the following configuration:

Inter_location distances

* A A B C D

* B A 0.0 10.44 9.49 14.14

B 10.44 0.0 9.22 7.0

* C C 9.49 9.22 0.0 7.07

* D D 14.14 7.0 7.07 0.0

Suppose that ß = 0 and * = 1/4. Then

P_ij = exp(_ d_ij/4)/[E (_ d_ij/4)] .

Consequently

Case I: for origin A fixed

P_AA = exp(_0/4)/[exp(_0/4) + exp(_10.44/4) + exp(_9.49/4) + exp(_14.14/4)]

= 1/(1 + 0.0735 + 0.0932 + 0.0292) = 1/1.1959 = 0.8362

P_AB = 0.0735/1.1959 = 0.0615

P_AC = 0.0932/1.1959 = 0.0779

P_AD = 0.0292/1.1959 = 0.0244

Case II: for origin B fixed

P_BA = exp(_10.44/4)/[exp(_10.44/4) + exp(_0/4) + exp(_9.22/4) + exp(_7/4)]

= 0.0735/(0.0735 + 1 + 0.0998 + 0.1738) = 0.0735/1.3471 = 0.0546

P_BB = 1/1.3471 = 0.7423

P_BC = 0.0998/1.3471 = 0.0741

P_BD = 0.1738/1.3471 = 0.1290

Case III: for origin C fixed

P_CA = exp(_9.49/4)/[exp(_9.49/4) + exp(_9.22/4) + exp(_0/4) + exp(_7.07/4)]

= 0.0932/(0.0932 + 0.0998 + 1 + 0.1708) = 0.0932/1.3638 = 0.0683

P_CB = 0.0998/1.3638 = 0.0732

P_CC = 1/1.3638 = 0.7333

P_CD = 0.1708/1.3638 = 0.1252

Case IV: for origin D fixed

P_DA = exp(_14.14/4)/[exp(_14.14/4) + exp(_7/4) + exp(_7.07/4) + exp(_0/4)]

= 0.0292/(0.0292 + 0.1738 + 0.1708 + 1) = 0.0292/1.3738 = 0.0213

P_DB = 0.1738/1.3738 = 0.1265

P_DC = 0.1708/1.3738 = 0.1243

P_DD = 1/1.3738 = 0.7279

In equation (2), if ß > 0 and * = 0, then the exponential term reduces to unity. Accordingly, equation (2) becomes

P_ij = P_j^ß/(E P_j^ß) , (3)

which is a purely hierarchical model. In contradistinction, as * 6 4, P_ii 6 1, implying that no interaction leaves origin i for any destination.

The random process mechanism used for simulation purposes takes each MIF in turn, and assigns integers from a given sequence [0, N_1] to the cells of this MIF. The number of integers assigned to a given MIF cell is proportional to the probability affiliated with the destination corresponding to that cell. Allocation of integers is done in such a way that the resulting groups are mutually exclusive and collectively exhaustive of the interval [0, N_1]. For the sake of simplicity, consecutive integers may be used.

This procedure may be exemplified with the foregoing Case I illustration. Suppose 1000 integers are to be used. Then the relevant interval is [0,999], or in terms of three_digit numbers [000,999]. Now the spatial distribution of probabilities in Case I is

* 0.8362

* 0.0615

* 0.0779

* 0.0244

One possible (and convenient) allocation of integers to this MIF for Location A would be

* 000_835

* 836_897

* 898_975

* 976_999

One should note that each integer is assigned to one and only one location, and that all integers in the interval [0,999] have been assigned to at least one of the locations.

row\column

1-4 5-8 9-12 13-16 17-20 21-24 25-28 29-32 33-36 37-40

1 08 54 86 21 23 60 37 45 16 15 84 73 40 10 01 00 90 02 71 55

2 29 01 44 01 92 09 31 97 13 17 56 67 91 47 44 77 76 87 31 00

3 97 91 86 46 69 16 60 24 78 91 45 10 20 84 39 63 95 07 94 65

4 20 19 25 90 95 63 45 83 42 17 53 07 47 22 64 96 34 77 44 40

5 65 88 31 95 35 47 90 75 60 27 74 82 66 09 15 55 21 90 25 22

6 49 39 81 12 05 72 67 07 29 54 80 42 78 15 71 59 14 21 64 90

7 77 49 83 21 89 66 29 97 57 19 30 62 42 16 44 42 96 41 97 98

8 21 06 20 53 58 01 02 15 67 65 89 76 87 37 17 06 26 36 44 49

9 30 67 97 44 90 44 43 02 49 59 33 29 50 68 30 09 98 11 57 66

10 95 08 12 94 20 85 30 10 08 37 03 83 48 04 54 30 87 15 43 67

11 82 06 07 84 99 79 68 19 88 10 59 33 26 45 77 96 28 13 24 14

12 38 66 52 15 59 60 28 15 50 98 72 40 89 64 33 12 88 14 44 22

13 73 11 94 24 87 96 00 79 03 68 55 47 47 10 52 97 78 70 10 34

14 75 55 71 46 19 85 68 64 62 15 60 18 74 56 27 76 50 86 36 89

15 92 24 20 68 63 09 41 88 65 64 10 70 35 09 60 55 67 87 41 21

16 22 82 29 99 85 47 97 63 36 12 97 74 09 57 96 94 88 28 36 52

17 02 18 79 01 88 35 02 94 99 01 69 40 40 01 58 23 49 47 19 96

18 20 75 00 61 54 04 01 19 55 86 92 48 28 79 88 22 34 72 85 96

19 67 54 06 02 87 33 39 33 64 80 86 53 49 61 37 04 46 86 47 73

20 94 45 07 19 03 91 58 86 83 50 49 62 09 97 80 99 23 07 44 61

21 33 01 19 11 78 59 02 79 95 42 35 98 34 77 50 34 33 83 96 41

22 43 17 01 26 64 36 77 30 12 31 67 59 97 99 57 39 24 86 61 55

23 86 96 06 73 79 41 79 27 21 96 16 76 00 32 91 75 62 61 94 08

24 76 13 53 48 92 62 51 14 36 97 44 21 34 09 40 00 62 03 63 67

25 81 04 49 46 32 93 91 48 20 23 20 60 64 60 44 76 17 33 54 54

In this case Location A is the origin. If one consults a random numbers table (see Table 1), then starting in row 11 and columns 18 through 20, and moving down through the table, for instance, three_digit numbers must be drawn. If 10 integers are selected, the resulting sequence here would be

810, 98, 368, 215, 564, 612, 901, 586, 480, 350 .

Therefore, Location A would interact with itself nine times (denoted by integers 810, 98, 368, 215, 564, 612, 586, 480, 350, which are contained in the interval [0,835]), Location A would not interact with Location B (since there are no integers from the interval [836,897]), Location A would interact with Location C one time (denoted by the integer 910, which is contained in the interval [898,975]), and Location A would not interact with Location D (since there are no integers from the interval [976,999]).

Cumulative diffusion over the landscape can be summarized with a logistic curve. This graph can be constructed by plotting the percentage of successful contacts against the number of transpired time periods. The resulting diagram would appear as

percent of 100 .

locations to .

which diffusion .

has occurred .

0 ................................... time

If P_i^" exp(_ * d_ii) = P_i^" dominates the probabilities defined by equation (2), then the slope of this S_shaped curve will be quite shallow. As * 6 0, the slope of this logistic curve should tend to become steeper.

INITIAL CONDITIONS

A landscape is provided that consists of a 13_by_11 set of Cartesian co_ordinates: (1,1), (1,2), ..., (10,11), (11,11). The spatial distribution of population over this punctiform surface conforms to a truncated K = 3 central place structure (see Figure 1), a clustered population distribution (MC = 0.76520; see Figure 2), or a random population distribution. In Figure 1, hexagons can be drawn by connecting the nearest neighbor locations having the same number. Euclidean distance between points is used in calculating the MIFs.

+ + + + + + + + + + + + +

1 + 2 1 3 1 2 1 2 1 3 1 2 +

2 + 1 2 1 2 1 3 1 2 1 2 1 +

3 + 2 1 5 1 2 1 2 1 5 1 2 +

4 + 1 2 1 2 1 3 1 2 1 2 1 +

5 + 2 1 3 1 2 1 2 1 3 1 2 +

6 + 1 2 1 2 1 9 1 2 1 2 1 +

7 + 2 1 3 1 2 1 2 1 3 1 2 +

8 + 1 2 1 2 1 3 1 2 1 2 1 +

9 + 2 1 5 1 2 1 2 1 5 1 2 +

10 + 1 2 1 2 1 3 1 2 1 2 1 +

11 + 2 1 3 1 2 1 2 1 3 1 2 +

12 + 2 1 3 1 2 1 2 1 3 1 2 +

13 + 1 2 1 2 1 5 1 2 1 2 1 +

14 + + + + + + + + + + + + +

1 2 3 4 5 6 7 8 9 1 1 1

0 1 2