Dear Dr. Excoffier,
I am currently conducting analyses using fastsimcoal28 but have encountered an issue with two of the four models I am testing. Specifically, I receive the error: "Simulations did not converge. MRCA not found!" in models 3 and 4.
I am using SFS data derived from RADseq, which I processed using easySFS. For each model, I am performing 50 replicate runs, with each run involving 100,000 simulations. I have attached an image detailing the models.
Could you review the parameters for models 3 and 4 to identify any potential errors? I have been unable to locate the issue myself.
Additionally, for models 1 and 2, I retrieve identical values for MaxObsLhood across replicates. Is this behavior expected?
Thank you for your time and assistance,
Ana
The files are the following:
MODEL 1
//Number of population samples (demes)
3
//Population effective sizes (number of genes)
N_POP0
N_POP1
N_POP2
//Samples sizes
12
72
34
//Growth rates: negative growth implies population expansion
0
0
0
//Number of migration matrices
3
//Migration matrix 0
0.000 0.000 0.000
0.000 0.000 0.000
0.000 0.000 0.000
//Migration matrix 1
0.000 0.000 0.000
0.000 0.000 0.000
0.000 0.000 0.000
//Migration matrix 2
0.000 0.000 0.000
0.000 0.000 0.000
0.000 0.000 0.000
//historical event: time, source, sink, migrants, new deme size, new growth rate, migration matrix index
2 historical event
TDIV01 0 1 1 RELANC01 0 0
TDIV2_ANC01 2 1 1 RELANC21 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 6.0277e-9 OUTEXP
// Priors and rules file
// *********************
[PARAMETERS]
//#isInt? #name #dist.#min #max
//all Ns are in number of haploid individuals
1 N_POP0 unif 10 2e6 output
1 N_POP1 unif 10 2e6 output
1 N_POP2 unif 10 2e6 output
1 N_ANC01 unif 10 2e6 output
1 N_ANCALL unif 10 2e6 output
1 TDIV01 unif 10 2e6 output
1 TPLUS01 unif 10 2e6 output
[COMPLEX PARAMETERS]
1 TDIV2_ANC01 = TDIV01+TPLUS01 output
0 RELANC21 = N_ANCALL/N_ANC01 hide
0 RELANC01 = N_ANC01/N_POP1 hide
MODEL 2
// Priors and rules file
// *********************
[PARAMETERS]
//#isInt? #name #dist.#min #max
//all Ns are in number of haploid individuals
1 N_POP0 unif 10 2e6 output
1 N_POP1 unif 10 2e6 output
1 N_POP2 unif 10 2e6 output
1 N_ANC01 unif 10 2e6 output
1 N_ANCALL unif 10 2e6 output
1 TDIV01 unif 10 2e6 output
1 TPLUS01 unif 10 2e6 output
[COMPLEX PARAMETERS]
1 TDIV2_ANC01 = TDIV01+TPLUS01 output
0 RELANC21 = N_ANCALL/N_ANC01 hide
0 RELANC01 = N_ANC01/N_POP1 hide
// Priors and rules file
// *********************
[PARAMETERS]
//#isInt? #name #dist.#min #max
//all Ns are in number of haploid individuals
1 N_POP0 unif 10 2e6 output
1 N_POP1 unif 10 2e6 output
1 N_POP2 unif 10 2e6 output
1 N_ANC01 unif 10 2e6 output
1 N_ANCALL unif 10 2e6 output
1 TDIV01 unif 10 2e6 output
1 TPLUS01 unif 10 2e6 output
0 NM2A logunif 1e-3 20 hide
0 NMA2 logunif 1e-3 20 hide
[COMPLEX PARAMETERS]
0 MIGA2 = NMA2/N_ANC01 output
0 MIG2A = NM2A/N_POP2 output
1 TDIV2_ANC01 = TDIV01+TPLUS01 output
0 RELANC21 = N_ANCALL/N_ANC01 hide
0 RELANC01 = N_ANC01/N_POP1 hide
MODEL 3
//Number of population samples (demes)
3
//Population effective sizes (number of genes)
N_POP0
N_POP1
N_POP2
//Samples sizes
12
72
34
//Growth rates: negative growth implies population expansion
0
0
0
//Number of migration matrices
3
//Migration matrix 0
0.000 MIG01 MIG02
MIG10 0.000 0.000
MIG20 0.000 0.000
//Migration matrix 1
0.000 0.000 0.000
0.000 0.000 0.000
0.000 0.000 0.000
//Migration matrix 2
0.000 0.000 0.000
0.000 0.000 0.000
0.000 0.000 0.000
//historical event: time, source, sink, migrants, new deme size, new growth rate, migration matrix index
2 historical event
TDIV01 0 1 1 RELANC01 0 0
TDIV2_ANC01 2 1 1 RELANC21 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 6.0277e-9 OUTEXP
// Priors and rules file
// *********************
[PARAMETERS]
//#isInt? #name #dist.#min #max
//all Ns are in number of haploid individuals
1 N_POP0 unif 10 2e6 output
1 N_POP1 unif 10 2e6 output
1 N_POP2 unif 10 2e6 output
1 N_ANC01 unif 10 2e6 output
1 N_ANCALL unif 10 2e6 output
1 TDIV01 unif 10 2e6 output
1 TPLUS01 unif 10 2e6 output
0 NM01 logunif 1e-3 20 hide
0 NM02 logunif 1e-3 20 hide
0 NM10 logunif 1e-3 20 hide
0 NM20 logunif 1e-3 20 hide
[COMPLEX PARAMETERS]
0 MIG01 = NM01/N_POP0 output
0 MIG02 = NM02/N_POP0 output
0 MIG10 = NM10/N_POP1 output
0 MIG20 = NM20/N_POP2 output
1 TDIV2_ANC01 = TDIV01+TPLUS01 output
0 RELANC21 = N_ANCALL/N_ANC01 hide
0 RELANC01 = N_ANC01/N_POP1 hide
MODEL 4
//Number of population samples (demes)
3
//Population effective sizes (number of genes)
N_POP0
N_POP1
N_POP2
//Samples sizes
12
72
34
//Growth rates: negative growth implies population expansion
0
0
0
//Number of migration matrices
3
//Migration matrix 0
0.000 MIG01 MIG02
MIG10 0.000 0.000
MIG20 0.000 0.000
//Migration matrix 1
0.000 0.000 0.000
0.000 0.000 MIGA2
0.000 MIG2A 0.000
//Migration matrix 2
0.000 0.000 0.000
0.000 0.000 0.000
0.000 0.000 0.000
//historical event: time, source, sink, migrants, new deme size, new growth rate, migration matrix index
2 historical event
TDIV01 0 1 1 RELANC01 0 0
TDIV2_ANC01 2 1 1 RELANC21 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 6.0277e-9 OUTEXP
// Priors and rules file
// *********************
[PARAMETERS]
//#isInt? #name #dist.#min #max
//all Ns are in number of haploid individuals
1 N_POP0 unif 10 2e6 output
1 N_POP1 unif 10 2e6 output
1 N_POP2 unif 10 2e6 output
1 N_ANC01 unif 10 2e6 output
1 N_ANCALL unif 10 2e6 output
1 TDIV01 unif 10 2e6 output
1 TPLUS01 unif 10 2e6 output
0 NM01 logunif 1e-3 20 hide
0 NM02 logunif 1e-3 20 hide
0 NM10 logunif 1e-3 20 hide
0 NM20 logunif 1e-3 20 hide
0 NM2A logunif 1e-3 20 hide
0 NMA2 logunif 1e-3 20 hide
[COMPLEX PARAMETERS]
0 MIG01 = NM01/N_POP0 output
0 MIG02 = NM02/N_POP0 output
0 MIG10 = NM10/N_POP1 output
0 MIG20 = NM20/N_POP2 output
0 MIGA2 = NMA2/N_ANC01 output
0 MIG2A = NM2A/N_POP2 output
1 TDIV2_ANC01 = TDIV01+TPLUS01 output
0 RELANC21 = N_ANCALL/N_ANC01 hide
0 RELANC01 = N_ANC01/N_POP1 hide