Imagine we simulate a group of 5 demes connected by the following migration pattern (looking back in time, because our simulation use the coalescent):

     [,1]   [,2]   [,3]   [,4]   [,5]  
[1,] 0.0000 0.0000 0.0000 0.0000 0.0005
[2,] 0.0005 0.0000 0.0000 0.0000 0.0000
[3,] 0.0000 0.0005 0.0000 0.0000 0.0000
[4,] 0.0000 0.0000 0.0005 0.0000 0.0000
[5,] 0.0000 0.0000 0.0000 0.0005 0.0000

The interpretation of this migration pattern forward in time would be:

     [,1]   [,2]   [,3]   [,4]   [,5]  
[1,] 0.0000 0.0005 0.0000 0.0000 0.0000
[2,] 0.0000 0.0000 0.0005 0.0000 0.0000
[3,] 0.0000 0.0000 0.0000 0.0005 0.0000
[4,] 0.0000 0.0000 0.0000 0.0000 0.0005
[5,] 0.0005 0.0000 0.0000 0.0000 0.0000

In words, we see that forward in time, migration from \(deme 1 \rightarrow deme 2 = 0.0005\), \(deme 2 \rightarrow deme 3 = 0.0005\), \(deme 3 \rightarrow deme 4 = 0.0005\), \(deme 4 \rightarrow deme 5 = 0.0005\) and \(deme 5 \rightarrow deme 1 = 0.0005\). This is equivilent to all demes being connected in a circular migration pattern in the direction of deme 1 to 5 as below.

Simulation

Now lets see how divMigrate resolves this pattern. First generate a .par file for use in fastsimcoal2

# load the custom parMaker function
source("parMaker.R")
# create a migration matrix with backward in time migration pattern
x <- matrix(0, ncol = 5, nrow = 5)
x[2, 1] <- 0.0005
x[3, 2] <- 0.0005
x[4, 3] <- 0.0005
x[5, 4] <- 0.0005
x[1, 5] <- 0.0005
# run parMaker
parMaker(outfile = "circ_pops", npops = 5, popSizes = 2000,
         sampSizes = 500, numMigMat = 1, migMat = x,
         numHistEvents = 0,
         numIndependentChromo = c(50, 0),
         lociData = list(locType = "MICROSAT",
                         numLoc = 1,
                         recombRate = 0,
                         mu = 0.001,
                         addPar = 0))

The above code generates the following fastsimcoal2 parameter file:

//Number of population samples (demes)
5 samples to simulate :
//Population effective sizes (number of genes)
2000
2000
2000
2000
2000
//Samples sizes
500
500
500
500
500
//Growth rates  : negative growth implies population expansion
0
0
0
0
0
//Number of migration matrices : 0 implies no migration between demes
1
//migration matrix
0.0000 0.0000 0.0000 0.0000 0.0005
0.0005 0.0000 0.0000 0.0000 0.0000
0.0000 0.0005 0.0000 0.0000 0.0000
0.0000 0.0000 0.0005 0.0000 0.0000
0.0000 0.0000 0.0000 0.0005 0.0000
//historical event: time, source, sink, migrants, new size, new growth rate, migr. matrix 
0  historical event
//Number of independent loci [chromosome] 
50 0
//Per chromosome: Number of linkage blocks
1
//per Block: data type, num loci, rec. rate and mut rate + optional parameters
MICROSAT 1 0 0.001 0

Now we can run fastsimcoal2:

system("./fastsimcoal21 -i circ_pops.par -n 1 -g")

To analyse the results of this simulation using divMigrate we need to convert the arlequin file into genepop format:

library(diveRsity)
arp2gen("circ_pops/circ_pops_1_1.arp")

From here it is possible to either analyse the generated genepop file online with divMigrate online or it can be analysed using the function divMigrate in the package diveRsity as below:

res <- divMigrate("circ_pops/circ_pops_1_1.gen", plot_network = TRUE, 
                  boots = 1000, para = TRUE, stat = "d")

DivMigrate generates two figures one with only the significant asymmetric migration calculated from 1000 bootstrap iterations and one with all values. It is further possible to filter the result such that only values above a certain value is included in the figures.

res <- divMigrate("circ_pops/circ_pops_1_1.gen", plot_network = TRUE, 
                  boots = 1000, para = TRUE, stat = "d",filter_threshold = 0.5)

We can see that the method resolves the correct direction of the forward time migration specified above. Although some migration rates in unspecified directions are found to be significantly asymmetric, these can be filtered out using the filter_threshold option in divMigrate. These false positive are not too surprising given the complexity of the migration pattern simulated.