# Advanced Exercises in Dispersal - Simulating an Epidemic

Throughout the exercises in dispersal, the focus has been on examining and comparing the mechanisms and assumptions that would lead to approximate exponential or power law distribution, or related distributions of propagules or disease as a function of distance from source. We can also study dispersal by asking questions that pertain to the source of inoculum. For instance, when would the starting distribution of infections be random? The distribution of initial infection is a function of the type of inoculum that initiates an epidemic. If infested seed is the source of initial infection, this would result in a random distribution to the extent that infested seed is randomly mixed with clean seed. If aerial dispersal of propagules into a field initiates the epidemic, the resulting distribution would tend to be approximately random if propagules have moved a substantial distance from the original source. By contrast, if the source of propagules is an adjacent field, propagules will tend to be deposited more heavily near that field. The distribution of infection will also differ from random, even if the distribution of propagules is random, when the biotic and abiotic environment is not homogeneous. For example, if some areas of a field are wetter, these patches will favor infection for many types of pathogens. Also, if different genotypes or species are planted across the field, such that there are different levels of resistance present, infection may be patchy.

Such concepts about dispersal can be studied via simulation. In the R code that follows, dispersal within a field will be explored, beginning with the case of random positioning of initial infections. First, a simple model of disease spread will be considered. For these exercises, we recommend that you work through this section-by-section, especially spending time to understand what the R code does. The information presented in this section on dispersal also provides an introduction to how exploration of disease processes may be handled through spatial statistics (Sparks et al. (2008).

To begin, suppose you would like to generate a simulated map of individual plants and their infection status, where infections are randomly assigned to individuals in the initial map and the number of infections per individual is tallied. One way to do this in R is as follows. Check out the contents of the different objects along the way. Use the `help(sample)` command to understand the sample function better.

```# First generate a matrix of dimension 10x10 containing
#   zeroes, where zero indicates no infection
map1 <- matrix(0,10,10)

# Randomly select row and column coordinates for n1 initial
#   infections, where n1 can be assigned any value
n1 <- 3
row1 <- sample(1:10,size=n1, replace=T)
col1 <- sample(1:10, size=n1, replace=T)

# Assign the infections to the plants corresponding to the
#   randomly drawn coordinates, and look at the new map
for(j in 1:n1) {  map1[row1[j],col1[j]] <- map1[row1[j],col1[j]] + 1}
map1
```

For example, here is one representation of three initial randomly infected plants:

```map1

[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,]    0    0    0    0    0    0    0    0    1     0
[2,]    0    0    0    0    0    0    0    1    0     0
[3,]    0    0    0    0    0    0    1    0    0     0
[4,]    0    0    0    0    0    0    0    0    0     0
[5,]    0    0    0    0    0    0    0    0    0     0
[6,]    0    0    0    0    0    0    0    0    0     0
[7,]    0    0    0    0    0    0    0    0    0     0
[8,]    0    0    0    0    0    0    0    0    0     0
[9,]    0    0    0    0    0    0    0    0    0     0
[10,]    0    0    0    0    0    0    0    0    0     0
```

## Suggested exercise

Repeat this step a few times each for different values of `n1`, the initial number of infections, to get a feel for the appearance of a random distribution of infection. If you pick very large numbers of `n1`, the results may start to become unrealistic for many host-pathogen systems because host tissue would become limiting.

As a first step to modeling dispersal of pathogens, you might consider even simpler functions to describe the dispersal gradient than the exponential and power law models. Suppose we start with a scenario much simpler than dispersal following some probability distribution. For most dispersal modeling, nearer dispersal is more likely and distant dispersal is less likely. Instead, suppose that pathogens always disperse in equal numbers to the plants directly adjacent to an infected plant on a square grid. This is sometimes described as one "rook's move", with reference to the moves possible for the chess piece.

To consider this type of movement, start by generating a map of a small number of initial infections using the method above to generate `map1`.

Next simulate one generation of pathogen reproduction and dispersal. The following simulation has five successful propagules produced for each existing infection and one of these disperses one step to the north, one disperses one step to the south, one disperses one step to the east, one disperses one step to the west, and one infects the same host. Note that following this method, some pathogen propagules will effectively be lost off the edge of the map.

```map2 <- map1
for(xrow in 1:10){
for(xcol in 1:10){
numprop <- map1[xrow,xcol]
if(xrow > 1){      map2[xrow-1,xcol] <- map2[xrow-1,xcol] + numprop}
if(xrow < 10){      map2[xrow+1,xcol] <- map2[xrow+1,xcol] + numprop}
if(xcol < 10){      map2[xrow,xcol+1] <- map2[xrow,xcol+1] + numprop}
if(xcol > 1){      map2[xrow,xcol-1] <- map2[xrow,xcol-1] + numprop}
map2[xrow,xcol] <- map2[xrow,xcol] + numprop
}
}
#See what map2 looks like now

map2
```

Here is an example of `map2` for one set of randomly generated initial infection locations.

```      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,]    0    0    0    0    0    0    0    2    2     1
[2,]    0    0    0    0    0    0    2    2    2     0
[3,]    0    0    0    0    0    1    2    2    0     0
[4,]    0    0    0    0    0    0    1    0    0     0
[5,]    0    0    0    0    0    0    0    0    0     0
[6,]    0    0    0    0    0    0    0    0    0     0
[7,]    0    0    0    0    0    0    0    0    0     0
[8,]    0    0    0    0    0    0    0    0    0     0
[9,]    0    0    0    0    0    0    0    0    0     0
[10,]    0    0    0    0    0    0    0    0    0     0
```

To simulate `n2` generations of dispersal, try the following. This also relies on having map1 generated as above with a small number of initial infections. Check out the objects created in each step (for each step that is not within a loop). Note that this is essentially the same set of commands as the example for the creation of `map2`, but executed `n2 - 1` times.

```#Assign the number of generations of#   dispersal to be studied, n2
n2 <- 5

#Set up mapn as a list that includes n2 maps
mapn <- as.list(1:n2)

#Initialize each of the n2 maps to contain zeroes
for(j in 1:n2){mapnj <- matrix(0,10,10)}

#For each of the next maps
mapn[] <- map1
for(j in 2:n2){
#Temporarily make the map the same as one generation back
mapnj <- mapn[[j-1]]
#Then add the new infections following the rooks' moves
for(xrow in 1:10){
for(xcol in 1:10){
numprop <- mapn[[j-1]][xrow,xcol]
if(xrow > 1){mapnj[xrow-1,xcol] <-
mapnj[xrow-1,xcol] + numprop}
if(xrow < 10){mapnj [xrow+1,xcol] <-
mapnj[xrow+1,xcol] + numprop}
if(xcol < 10){mapnj [xrow,xcol+1] <-
mapnj[xrow,xcol+1] + numprop}
if(xcol > 1){mapnj [xrow,xcol-1] <-
mapnj[xrow,xcol-1] + numprop}
mapnj[xrow,xcol] <- mapnj[xrow,xcol] + numprop
}
}
}

#See what mapn looks like
mapn
```

Here is one example of `mapn`:

```1
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,]    0    0    0    0    0    0    0    0    1     0
[2,]    0    0    0    0    0    0    0    1    0     0
[3,]    0    0    0    0    0    0    1    0    0     0
[4,]    0    0    0    0    0    0    0    0    0     0
[5,]    0    0    0    0    0    0    0    0    0     0
[6,]    0    0    0    0    0    0    0    0    0     0
[7,]    0    0    0    0    0    0    0    0    0     0
[8,]    0    0    0    0    0    0    0    0    0     0
[9,]    0    0    0    0    0    0    0    0    0     0
[10,]    0    0    0    0    0    0    0    0    0     0

2
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,]    0    0    0    0    0    0    0    2    2     1
[2,]    0    0    0    0    0    0    2    2    2     0
[3,]    0    0    0    0    0    1    2    2    0     0
[4,]    0    0    0    0    0    0    1    0    0     0
[5,]    0    0    0    0    0    0    0    0    0     0
[6,]    0    0    0    0    0    0    0    0    0     0
[7,]    0    0    0    0    0    0    0    0    0     0
[8,]    0    0    0    0    0    0    0    0    0     0
[9,]    0    0    0    0    0    0    0    0    0     0
[10,]    0    0    0    0    0    0    0    0    0     0

3
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,]    0    0    0    0    0    0    4    8    9     4
[2,]    0    0    0    0    0    3    8   12    8     3
[3,]    0    0    0    0    1    4   10    8    4     0
[4,]    0    0    0    0    0    2    4    3    0     0
[5,]    0    0    0    0    0    0    1    0    0     0
[6,]    0    0    0    0    0    0    0    0    0     0
[7,]    0    0    0    0    0    0    0    0    0     0
[8,]    0    0    0    0    0    0    0    0    0     0
[9,]    0    0    0    0    0    0    0    0    0     0
[10,]    0    0    0    0    0    0    0    0    0     0

4
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,]    0    0    0    0    0    7   24   41   38    20
[2,]    0    0    0    0    4   18   45   56   44    18
[3,]    0    0    0    1    6   24   44   45   24     7
[4,]    0    0    0    0    3   12   24   18    7     0
[5,]    0    0    0    0    0    3    6    4    0     0
[6,]    0    0    0    0    0    0    1    0    0     0
[7,]    0    0    0    0    0    0    0    0    0     0
[8,]    0    0    0    0    0    0    0    0    0     0
[9,]    0    0    0    0    0    0    0    0    0     0
[10,]    0    0    0    0    0    0    0    0    0     0

5
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,]    0    0    0    0   11   56  141  200  181    96
[2,]    0    0    0    5   32  116  232  287  224   107
[3,]    0    0    1    8   44  128  226  232  151    56
[4,]    0    0    0    4   24   78  128  116   56    14
[5,]    0    0    0    0    6   24   44   32   11     0
[6,]    0    0    0    0    0    4    8    5    0     0
[7,]    0    0    0    0    0    0    1    0    0     0
[8,]    0    0    0    0    0    0    0    0    0     0
[9,]    0    0    0    0    0    0    0    0    0     0
[10,]    0    0    0    0    0    0    0    0    0     0
```

The output `mapn` is a list in which each matrix is a map after a new generation of dispersal.

## Suggested exercise

Try this for a series of different initial maps generated as `map1`. To do this, you could cut and paste the material from above when needed. Alternatively, you could create a function that performs these steps, as below.

```disperse <- function(map1,n2){
mapn <- as.list(1:n2)
for(j in 1:n2){mapnj <- matrix(0,10,10)}
mapn[] <- map1
for(j in 2:n2){
mapnj <- mapn[[j-1]]
for(xrow in 1:10){
for(xcol in 1:10){
numprop <- mapn[[j-1]][xrow,xcol]
if(xrow > 1){          mapnj[xrow-1,xcol] <-              mapnj[xrow-1,xcol] + numprop}
if(xrow < 10){          mapnj [xrow+1,xcol] <-              mapnj[xrow+1,xcol] + numprop}
if(xcol < 10){          mapnj [xrow,xcol+1] <-              mapnj[xrow,xcol+1] + numprop}
if(xcol > 1){          mapnj [xrow,xcol-1] <-              mapnj[xrow,xcol-1] + numprop}
mapnj[xrow,xcol] <-            mapnj[xrow,xcol] + numprop
}
}
}
mapn
}
```

You could use this function as in the following examples (which assume you have a matrix `map1` already created). You could try this command after generating different versions of `map1` to compare the patterns.

```disperse(map1,n2=3)
disperse(map1,n2=5)
```

In the above examples, locations were randomly selected individual hosts where the number of infections was tallied. It might be more realistic to consider a situation where the initial infections are randomly applied to a continuous map. That is, `nstart` infections are assigned to continuous x and y coordinates in a 1000 m by 1000 m map. The random assignment is performed using the `runif` function. This function generates random numbers with a uniform distribution across the range requested. Try `help(runif)` for more information.

```#Assign the number of initial infections, nstart
nstart <- 20

#Randomly draw x and y coordinates for the initial infections
row1 <- runif(n=nstart, min=0, max=1000)
col1 <- runif(n=nstart, min=0, max=1000)

#Plot the initial infections
plot(row1, col1, xlab="East-West",     ylab="North-South", pch=19, col="red")
```

Next, suppose that each infection produces nnew new infections in the next pathogen generation. These new infections are equally likely to occur in any direction from the source infection. (When would this assumption not be reasonable?) The coordinates for a randomly-drawn direction can be generated as follows, again using the `runif` function, but ranging from -? to ?. Here the random variables are transformed with functions `sin` and `cos` to generate points in all possible directions from a source in two dimensions.

```#Assign the new number of infections from#   each initial infection, nnew
nnew <- 50

#Generate nnew random directions from a source
temp <- runif(n=nnew, min=-pi, max=pi)

#Plot the points in the different directions,#   for wno equidistant from source
plot(sin(temp),cos(temp))
```

Suppose the distance from the source at which a new infection occurs is described by an exponential model. R supplies a function to generate random numbers generated from an exponential distribution, as well as the uniform distribution used above and many others (see Garrett et al. In Press). The `rexp` function generates random variables drawn from an exponential distribution with parameter rate. You can get a feel for how the rate parameter affects the shape of the distribution by trying the following. The `hist` command produces a histogram of the variables generated. If n=500 or more variables are generated, this gives a pretty clear idea of the shape of the distribution. If a small number of variables are generated, the shape is not so clear. But it is also informative to try drawing just a small number n of samples from the distribution to get a feel for what it means to "draw from a distribution".

```hist(rexp(n=500, rate=1),col="blue")

```

Now, combine the initial random map of infections, the random direction of new infections, and the exponential model describing the distance at which new infections occur. The following function is the first step for generating maps of infection for ngen generations.

```# nstart is the number of initial infections
# nnew is the number of new infections produced by each
#      pre-existing infection each generation
# ngen is the number of generations
# exppar is the rate parameter for#    the exponential distribution
disperse2 <- function(nstart, nnew, ngen, exppar){
coord <- matrix(runif(n=2*nstart, min=0, max=1000),
nrow=nstart, ncol=2)
ntot <- nstart
for(i in 1:ngen){
for(j in 1:ntot){
tempx <- coord[j,1]
tempy <- coord[j,2]
newdir <- runif(n=nnew, min=-pi, max=pi)
newdist <- rexp(n=nnew, rate=exppar)
newx <- tempx + newdist * sin(newdir)
newy <- tempy + newdist * cos(newdir)
newcoord <- cbind(newx,newy)
coord <- rbind(coord,newcoord)
}
ntot <- length(coord[,1])
}
coord
}
```

Try the following commands, varying the values of the input parameters. Especially try different values of the exponential paramater, `exppar`. More fun than a video game?

```tmat <- disperse2(nstart=20, nnew=10, ngen=3, exppar=0.1)
plot(tmat, xlim=c(0,1000), ylim=c(0,1000), xlab="East-West",     ylab="North-South", col="Red")
```

This leads to the following graph:

Output

As a more advanced exercise, try generating different maps using R functions that generate random numbers from the related gamma and Weibull distributions, too: `rgamma` and `rweibull`. Try different values of the shape and scale parameters for both these distributions too. Here are examples to remind you of the shape of those distributions; you can change the shape parameters to plot histograms representing the different parameter combinations.

```hist(rgamma(n=500, shape=1, scale=1), col="blue")
```
```hist(rweibull(n=500, shape=1, scale=1), col="blue")
```

Note that there are several different ways to get similar-looking spatial patterns. In a practical application, such as trying to model a real epidemic, it may not matter greatly which of these models you use. This exercise is intended to give you an idea of some of the tools available.

Next, concluding remarks