MaxEstLhood zero when adding migration matrix

[Alex van Rensburg]

Dear Laurent, 

I'm running a simple three population model to estimate divergence time. I'm using derived SFS estimated from ~9000 independent RAD loci. The test models with migration set to 0 works with no problems, but as soon as I set a migration model MaxEstLhood is always zero. I assume fsc finds the correct obs file given that my previous test models worked. Do you have some suggestions of what I could change?

fsc26 -t BA-Model1.tpl -e BA-Model1.est -L50 -n100000 -d -M  -q -0 -C1 --multiSFS -c1 -B1

.tpl file: 

//Parameters for the coalescence simulation program : fsimcoal2.exe

3 samples to simulate :

//Population effective sizes (number of genes)

NNEW

NSOUTH

NHOD

//Samples sizes and samples age 

10

10

6

//Growth rates : negative growth implies population expansion

0

0

0

//Number of migration matrices : 0 implies no migration between demes

2

//migration matrix 0

0 0 MIG02

0 0 0

MIG20 0 0

//migration matrix 1

0 MIG01 0

MIG10 0 MIG12

0 MIG21 0

//historical event: time, source, sink, migrants, new deme size, new growth rate, migration matrix index

2

TDIV1 0 2 1 RESIZE 0 0

TDIV2 2 1 1 RESIZE 0 1

//Number of independent loci [chromosome] 

1 0

//Per chromosome: Number of contiguous linkage Block: a block is a set of contiguous loci

1

//per Block:data type, number of loci, per generation recombination and mutation rates and optional parameters

FREQ  1   0   2.5e-9 OUTEXP

.est file

// Priors and rules file

// *********************

[PARAMETERS]

//#isInt? #name #dist.#min #max

//all N are in number of haploid individuals

1 ANCSIZE unif 1000 1000000 output

1 NNEW unif 1000 1000000 output

1 NSOUTH unif 1000 1000000 output

1 NHOD unif 1000 1000000 output

1 TDIV1 unif 10 100 output

1 TPLUSDIV unif 10 10000 output

0 MIG01 logunif 1e-5 1e-2 output

0 MIG10 logunif 1e-5 1e-2 output

0 MIG02 logunif 1e-5 1e-2 output

0 MIG20 logunif 1e-5 1e-2 output

0 MIG12 logunif 1e-5 1e-2 output

0 MIG21 logunif 1e-5 1e-2 output

0 MUTRATE unif 2.9e-9 5.5e-9 output

[RULES]

[COMPLEX PARAMETERS]

0 RESIZE = ANCSIZE/NSOUTH output

0 TDIV2 = TDIV1+TPLUSDIV output

Thanks for the help

Alex

[Chrystelle Delord]

Dear Alex,

I would suggest to try rewriting your migration matrices and add a third matrix with no migration (i.e., full of zeros) to set at TDIV2.

Indeed, currently, at TDIV1, you had set migration matrix 0 which defines migrations rates between deme 0 and 2, but since these have actually merged together; migration rates between these two are irrelevant.

Backward in time you have:

Deme 0, 1 and 2 with migration between potentially all possible pairs, defined by matrix 0.

Then at TDIV1, deme 0 merges with deme 2, so the new migration matrix (matrix 1) should contain only migration rates between sink deme 2, and the remaining deme 1.

Finally at TDIV2, deme 2 merges with deme 1 which is the only deme left. There is no migration anymore because there is only one population. Thus the new migration matrix at TDIV2 (matrix 2) should be

0 0 0

0 0 0

0 0 0.

Try writing:

//Number of migration matrices : 0 implies no migration between demes

3
//migration matrix 0
0 MIG01 MIG02
MIG10 0 MIG12
MIG20 MIG21 0
//migration matrix 1
0 0 0

0 0 MIG12

0 MIG21 0

//migration matrix 2
0 0 0
0 0 0
0 0 0

//historical event: time, source, sink, migrants, new deme size, new growth rate, migration matrix index
2

TDIV1 0 2 1 RESIZE 0 1

TDIV2 2 1 1 RESIZE 0 2

I hope this will work,

Best regards,

Chrys

[Alex van Rensburg]

Hi Chrys,

Thank you - that's solved the problem! I realise I misunderstood that the matrix 0 gets used initially, so any changes from that would be specified by the other matrices - hence I need 3 matrices, not two like I had. This makes more sense now. 

Thank you

Alex