This code runs just one simulation

**For Fig 6.7 (10-patch neutral model)**: `J`

= 100, `init.1`

= 0.5*`J`

, `num.patch`

= 10, `num.years`

= 50, `fit.ratio.avg[]`

= 1, `freq.dep[]`

= 0, `m`

varies among panels

**For Fig 6.8 (spatially variable selection - symmetrical)**: as in Fig. 6.7 except `J`

= 1000, `num.patch`

= 2, `fit.ratio.avg[] <- c(1.2,(1/1.2))`

**For Fig 6.9 (spatially variable selection - asymmetrical)**: as in Fig. 6.8 except `num.years`

= 100, `fit.ratio.avg[] <- c(1.5,(1/1.1))`

```
## specify parameters, initial conditions, and output matrix
num.years <- 50
num.patch <- 10
freq.1.mat <- matrix(nrow = num.years, ncol = num.patch)
J <- 100 # number of individuals PER PATCH
init.1 <- 0.5*J
COM <- matrix(nrow=J, ncol=num.patch)
COM[1:init.1,] <- 1; COM[(init.1+1):J,] <- 2
year <- 2
m <- 0
fit.ratio.avg <- vector(length=num.patch)
fit.ratio.avg[] <- 1
freq.dep <- vector(length=num.patch)
freq.dep[] <- 0
## record data (frequency of species 1) for year 1
freq.1.mat[1,] <- init.1/J
## run simulation
for (i in 1:(J*num.patch*(num.years-1))) {
## choose a patch where a death even will occur
patch <- sample(1:num.patch,1)
## calculate Pr.1 if dispersal occurs
if (runif(1) < m) {
Pr.1 <- sum(COM==1)/(J*num.patch)
} else {
## calculate Pr.1 if local reproduction (not dispersal)
freq.1 <- sum(COM[,patch]==1)/J; freq.2 <- 1 - freq.1
fit.ratio <- exp(freq.dep[patch]*(freq.1-0.5) + log(fit.ratio.avg[patch]))
Pr.1 <- fit.ratio*freq.1/(fit.ratio*freq.1 + freq.2)
}
COM[ceiling(J*runif(1)),patch] <- sample(c(1,2), 1, prob=c(Pr.1,1-Pr.1))
## record data
if (i %% (J*num.patch) == 0) {
freq.1.mat[year,] <- colSums(COM==1)/J
year <- year + 1
}
}
## graph the results
plot(1:num.years, freq.1.mat[,1], type="l", xlab="Time",
ylab="Frequency of species 1", ylim=c(0,1))
for (i in 2:(num.patch)) {
lines(1:num.years,freq.1.mat[,i], type="l", lty=2, ylim=c(0,1))
}
```

- Unlike Boxes 1-4, we can now have more than two patches, with the number defined by
`num.patch`

.`freq.1.mat`

is a matrix for output, as before, but with columns that are patches rather than separate simulations.

```
num.years <- 50
num.patch <- 10
freq.1.mat <- matrix(nrow = num.years, ncol = num.patch)
```

`m`

is the dispersal parameter (`m`

for migration, borrowed from the population genetics tradition). It defines the probability that the individual chosen for reproduction comes from the entire metacommunity rather than the local community.

`m <- 0`

- With multiple patches, each can have its own value of the selection parameters, such that
`fit.ratio.avg`

and`freq.dep`

are now vectors with`num.patch`

elements. For example,`fit.ratio.avg[] <- 1`

gives all patches the same average fitness ratio (in this case 1). For two patches,`fit.ratio.avg[] <- c(1.2,(1/1.2))`

assigns the average fitness ratios of 1.2 and 1/1.2 to the two patches. For more patches, we add more values, separated by commas.

```
fit.ratio.avg <- vector(length=num.patch)
fit.ratio.avg[] <- 1
freq.dep <- vector(length=num.patch)
freq.dep[] <- 0
```

- The metacommunity contains
`J*num.patch`

individuals, so each year must involve`J*num.patch`

deaths. To choose a random individual for death, we first choose a patch by sampling an integer between 1 and`num.patch`

at random.

```
for (i in 1:(J*num.patch*(num.years-1))) {
## choose a patch where a death even will occur
patch <- sample(1:num.patch,1)
```

- Here is where dispersal happens. With probability
`m`

, we will define`Pr.1`

based only on the frequency of species 1 in the entire metacommunity, effectively making the reproducing individual a random selection from the metacommunity. The`else`

statement says that we will otherwise (i.e., with probability 1 -`m`

) choose the reproducing individual from the local community as in previous Boxes.

```
if (runif(1) < m) {
Pr.1 <- sum(COM==1)/(J*num.patch)
} else {
```

- The graphs in Figures 6.7-6.9 have slightly different formatting. After calling the
`lines`

function,`lty=2`

makes all the lines (except the first) dashed, as in Figures 6.8 and 6.9. For all the lines to be solid (as in Figure 6.7), we can set`lty=1`

instead.

```
plot(1:num.years, freq.1.mat[,1], type="l", xlab="Time",
ylab="Frequency of species 1", ylim=c(0,1))
for (i in 2:(num.patch)) {
lines(1:num.years,freq.1.mat[,i], type="l", lty=2, ylim=c(0,1))
```