Problem with SFS & estimation?

[Ivan Scotti]

Dear Laurent and the group,

we are stuck with a model that does not seem to want to run, but we're having a hard time figuring out where the error comes from - whether it is  a syntax problem, or a problem with the ranges of the priors, or a true hardware problem.

This is the TPL we're using: four populations that first sequentially merge into one, and then undergo a series of changes in growth rates:

//Number of population samples (demes)

4 samples to simulate

//Population effective sizes (number of genes)

POPSIZE0

POPSIZE1

POPSIZE2

POPSIZE3

//Sample sizes

40 0 -0.05253664

50 0 -0.04238771

50 0 -0.04314575

52 0 -0.04139063

//Growth rates : negative growth implies population expansion

GROWTHRATE0

GROWTHRATE0

GROWTHRATE0

GROWTHRATE0

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

1

//migration matrix

0.000 M01 M02 M03 

M12 0.000 M12 M13 

M02 M12 0.000 M23 

M03 M13 M23 0.000

//historical event: time, source, sink, migrants, new size, growth rate, migr. matrix NOTE: order of event 1 and 2 to be estimated

8 historical event

TIME1 0 1 1 1 GROWTHRATE0 0

TIME2 2 1 1 1 GROWTHRATE0 0

TIME3 3 1 1 1 GROWTHRATE0 0

TIMEA 1 1 1 1 0 0

TIMEB 1 1 1 1 GROWTHRATE1 0

TIMEC 1 1 1 1 0 0

TIMED 1 1 1 1 GROWTHRATE2 0

TIMEE 1 1 1 1 0 0

//Number of independent loci [chromosome]

12 0

//Per chromosome: Number of linkage blocks

1

//per Block: data type, num loci, rec. rate and mut rate + optional parameters

FREQ 1 0 1e-8 OUTEXP

//Per chromosome: Number of linkage blocks

1

//per Block: data type, num loci, rec. rate and mut rate + optional parameters

FREQ 1 0 1e-8 OUTEXP

//Per chromosome: Number of linkage blocks

1

//per Block: data type, num loci, rec. rate and mut rate + optional parameters

FREQ 1 0 1e-8 OUTEXP

//Per chromosome: Number of linkage blocks

1

//per Block: data type, num loci, rec. rate and mut rate + optional parameters

FREQ 1 0 1e-8 OUTEXP

//Per chromosome: Number of linkage blocks

1

//per Block: data type, num loci, rec. rate and mut rate + optional parameters

FREQ 1 0 1e-8 OUTEXP

//Per chromosome: Number of linkage blocks

1

//per Block: data type, num loci, rec. rate and mut rate + optional parameters

FREQ 1 0 1e-8 OUTEXP

//Per chromosome: Number of linkage blocks

1

//per Block: data type, num loci, rec. rate and mut rate + optional parameters

FREQ 1 0 1e-8 OUTEXP

//Per chromosome: Number of linkage blocks

1

//per Block: data type, num loci, rec. rate and mut rate + optional parameters

FREQ 1 0 1e-8 OUTEXP

//Per chromosome: Number of linkage blocks

1

//per Block: data type, num loci, rec. rate and mut rate + optional parameters

FREQ 1 0 1e-8 OUTEXP

//Per chromosome: Number of linkage blocks

1

//per Block: data type, num loci, rec. rate and mut rate + optional parameters

FREQ 1 0 1e-8 OUTEXP

//Per chromosome: Number of linkage blocks

1

//per Block: data type, num loci, rec. rate and mut rate + optional parameters

FREQ 1 0 1e-8 OUTEXP

//Per chromosome: Number of linkage blocks

1

//per Block: data type, num loci, rec. rate and mut rate + optional parameters

FREQ 1 0 1e-8 OUTEXP

 and this is the EST where the parameters are taken from:

// Priors and rules file

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

[PARAMETERS]

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

//all N are in number of haploid individuals

1 POPSIZE0 logunif 100 30000 output

1 POPSIZE1 logunif 100 30000 output

1 POPSIZE2 logunif 100 30000 output

1 POPSIZE3 logunif 50 30000 output

0 GROWTHRATE0 unif -0.5 -0.05 output

0 M01 logunif 0.001 0.1 output

0 M02 logunif 0.001 0.1 output

0 M03 logunif 0.001 0.1 output

0 M12 logunif 0.001 0.1 output

0 M13 logunif 0.001 0.1 output

0 M23 logunif 0.001 0.1 output

1 TIMEA unif 2 20 output

1 TIME3 logunif 1 TIMEA output paramInRange

1 TIME1 logunif 1 TIMEA output paramInRange

1 TIME2 logunif 1 TIMEA output paramInRange

0 MUTRATE logunif 4e–09 6.56e–09 output

1 TIMEB unif TIMEA 30 output paramInRange

1 TIMEC unif TIMEB 100 output paramInRange

1 TIMED unif TIMEC 300 output paramInRange

1 TIMEE unif TIMED 1000 output paramInRange

0 GROWTHRATE1 unif 1e-9 0.1 output

0 GROWTHRATE2 unif -0.01 -1e-9 output

[COMPLEX PARAMETERS]

And finally, the command line we're using to produce simulations:

fsc2702 -t inputTemplate_v9.tpl -n 1 -m -e inputTemplate_v9.est -E 1000

and the screen output we got from the last attempt at running it:

fastsimcoal was invoked with the following command line arguments:

fsc2702 -t inputTemplate_v9.tpl -n 1 -m -e inputTemplate_v9.est -E 1000 

Reading parameter settings file inputTemplate_v9.est...

done reading parameter settings file!

Random generator seed : 279725

Deme sizes

Deme 0 706

Deme 1 1320

Deme 2 200

Deme 3 69

Sample sizes

Deme 0 40

Deme 1 50

Deme 2 50

Deme 3 52

Sample ages

Deme 0 0

Deme 1 0

Deme 2 0

Deme 3 0

Growth rates

Deme 0 -0.0532047

Deme 1 -0.0532047

Deme 2 -0.0532047

Deme 3 -0.0532047

Migration matrix 0 (note that migration rates have been pre-multiplied by 1e+06n

 0.0000000 24116.6129606 38186.1167473 1081.5045148 

2856.0462804  0.0000000 2856.0462804 22620.5963647 

38186.1167473 2856.0462804  0.0000000 34041.1724768 

1081.5045148 22620.5963647 34041.1724768  0.0000000 

Historical events

Event 0

#Time              : 5

#Source            : 3

#Sink              : 1

#Migrants          : 1.0000000

#New relative size : 1.0000000

#New growth rate   : -0.0532047

#New migr. matrix  : 0

Event 1

#Time              : 6

#Source            : 2

#Sink              : 1

#Migrants          : 1.0000000

#New relative size : 1.0000000

#New growth rate   : -0.0532047

#New migr. matrix  : 0

Event 2

#Time              : 13

#Source            : 0

#Sink              : 1

#Migrants          : 1.0000000

#New relative size : 1.0000000

#New growth rate   : -0.0532047

#New migr. matrix  : 0

Event 3

#Time              : 16

#Source            : 1

#Sink              : 1

#Migrants          : 1.0000000

#New relative size : 1.0000000

#New growth rate   : 0.0000000

#New migr. matrix  : 0

Event 4

#Time              : 22

#Source            : 1

#Sink              : 1

#Migrants          : 1.0000000

#New relative size : 1.0000000

#New growth rate   : 0.0695284

#New migr. matrix  : 0

Event 5

#Time              : 87

#Source            : 1

#Sink              : 1

#Migrants          : 1.0000000

#New relative size : 1.0000000

#New growth rate   : 0.0000000

#New migr. matrix  : 0

Event 6

#Time              : 234

#Source            : 1

#Sink              : 1

#Migrants          : 1.0000000

#New relative size : 1.0000000

#New growth rate   : -0.0085402

#New migr. matrix  : 0

Event 7

#Time              : 988

#Source            : 1

#Sink              : 1

#Migrants          : 1.0000000

#New relative size : 1.0000000

#New growth rate   : 0.0000000

#New migr. matrix  : 0

Number of independent loci to simulate : 12

with the same chromosomal structure

Number of linkage blocks to simulate in structure 1: 1

Data type: FREQ   :  recombination and mutation rates

0.0000000000   0.0000000100   

   Output expected MAF spectrum

   Do not output expected DAF spectrum

Population growth detected in input file

fastsimcoal2 is building 12 genealogies ...

Computing expected SFS using 12 batches and 1 threads.

free(): invalid next size (fast)

Abandon (core dumped)

The example of one of the observed SFS is here:
(we did not find a way to attach text files to the message)

1 observation

d0_0 d0_1 d0_2 d0_3 d0_4 d0_5 d0_6 d0_7 d0_8 d0_9 d0_10 d0_11 d0_12 d0_13 d0_14 d0_15 d0_16 d0_17 d0_18 d0_19 d0_20

d1_0 4458 2406 617 217 95 24 11 6 0 0 1 0 0 0 0 0 0 0 0 0 0

d1_1 3556 1475 690 341 154 74 29 16 5 2 0 1 1 0 0 0 0 0 0 0 0

d1_2 1179 970 607 382 199 101 70 47 21 7 2 0 0 0 0 0 0 0 0 0 0

d1_3 452 564 529 403 255 166 92 50 24 11 8 3 1 0 0 0 0 0 0 0 0

d1_4 220 306 380 339 270 192 115 83 38 18 15 6 1 1 1 0 0 0 0 0 0

d1_5 86 194 251 275 269 198 135 100 52 31 24 9 7 4 0 2 0 0 0 0 0

d1_6 47 103 190 218 225 197 156 96 91 37 25 13 14 3 0 0 0 1 0 0 0

d1_7 20 84 110 142 183 160 155 116 95 77 40 26 14 4 3 0 2 0 0 0 0

d1_8 10 35 66 101 142 169 126 128 120 67 42 39 29 9 8 5 2 2 0 1 1

d1_9 6 18 41 47 87 124 126 130 121 83 93 49 37 32 10 5 2 0 0 3 0

d1_10 1 9 25 31 60 96 111 101 108 104 84 67 44 24 16 10 6 2 2 1 0

d1_11 0 5 11 26 52 78 77 113 114 109 93 66 53 35 34 19 8 5 5 1 1

d1_12 0 8 10 14 34 52 64 100 82 92 97 77 69 48 36 25 21 9 8 2 0

d1_13 0 4 2 16 21 23 68 60 77 85 87 92 77 68 66 49 33 21 5 11 7

d1_14 0 0 0 9 6 18 45 63 53 73 80 84 90 70 55 44 41 29 14 5 2

d1_15 0 0 1 2 4 17 16 34 54 52 72 85 61 59 68 54 53 48 26 10 4

d1_16 0 0 0 2 1 14 12 23 38 35 60 81 75 96 97 56 53 37 29 19 9

d1_17 0 0 0 0 3 10 8 10 21 26 48 49 62 88 60 69 69 47 43 44 19

d1_18 0 1 0 1 1 1 6 5 20 27 35 48 60 64 87 84 58 64 59 58 28

d1_19 0 0 0 1 1 0 3 5 14 22 29 43 53 66 69 62 82 62 63 57 24

d1_20 0 0 0 0 0 0 3 4 5 16 14 36 46 53 51 72 82 78 86 67 32

d1_21 0 0 0 0 0 0 5 5 6 8 12 20 21 48 53 66 74 96 105 100 56

d1_22 0 0 0 0 0 0 1 1 2 3 12 19 15 48 52 55 66 104 113 137 66

d1_23 0 0 0 0 1 0 0 0 1 1 5 11 15 21 49 61 88 90 111 169 75

d1_24 0 0 0 0 0 0 0 2 2 2 6 18 20 17 43 63 62 105 167 161 151

d1_25 0 0 0 0 0 0 0 0 1 2 6 3 6 15 15 28 43 51 59 99 125

Notice that we have tested a run replacing the FREQ data type with DNA (taking the line at page 12 of the manual to describe the crhomosomes) and the simulations pan out nicely. So the problem must be somewhere between the way we write the chromosomes in the TPL and the way we write the OBS files.

[Notice that it is the first time that we see the one-before-last line appearing on our screen! Previously, we only got the 'core dumped' message. So this perhaps adds a sudden turn to our story...]

We'll be happy with any help you might provide us. For the time being, we'll keep trying changing things here and there to see whether we can find the culprit empirically.

Many thanks, yours

Ivan Scotti, Andrea Modica
INRAE Avignon, France

[Andrea Modica]

Dear Laurent and the group,

Here are also the SFSs, if they can be of some use in understanding the nature of the issue.

Thanks for any help or suggestion, yours

Ivan Scotti, Andrea Modica

INRAE Avignon, France

[Laurent Excoffier]

Hi,

if you want to estimate parameters from the data, youdo not need to specify multiple FREQ data types

The last lines of your tpl file shoudl look like:

//Number of independent loci [chromosome]

1 0

//Per chromosome: Number of linkage blocks
1
//per Block: data type, num loci, rec. rate and mut rate + optional parameters
FREQ 1 0 1e-8 OUTEXP

and that's it!

Moreover, it seems as if you had unfolded SFSs instead of folded ones

Then the way you call the program is not appropriate

Indeed of 

fsc2702 -t inputTemplate_v9.tpl -n 1 -m -e inputTemplate_v9.est -E 1000

it shoudl look rather like

fsc2702 -t inputTemplate_v9.tpl -n 1000000 -d -e inputTemplate_v9.est -M -q

-M will launch parameter estimations, -E is not correct here, and I change -m to -d if you have unfolded SFS

Finally, in your observed SFS, you do not seem to have many monomoprhic sites, making me wonder if you have really analysed random DNA sequences. I looks more like a SNP array to me...

Hope it helps

laurent 

[Ivan Scotti]

Dear Laurent,

thank you for you reply and your patience. 

Indeed, our command line was wrong and, moreover, the SFSs missed all the monomorphic sites (indeed, we have >3Mb of sequence in this project, but after computing the frequencies in each pop for the variants, we sort of "forgot" (code-wise) to count the invariant sites again...).

We're going to run the analyses after the corrections and will report back here about how it goes.

Kind regards,

Ivan

[Andrea Modica]

Dear Laurent and the group,

thanks for your kind replies and patience.

As Ivan said, after having corrected the files and added the missing monomorphic loci to [0,0], we ran the analyses many times with the command line you suggested.

Mostly the software returned a core dumped, but sometimes it returned errors different from this one and sometimes it ran for a while before going core dumped.

We attach again all the datafiles (TPL, EST, OBS) we are using, including the CLI we write in the terminal, hoping they can help you in helping us by noticing some hidden bug that escaped our eyes. We wrote the screen errors in the file VariousErrors.txt, if they can testify something more about the issue that's going on in our analysis. 

Again, many thanks for your time and advice.

Kindest regards,

Andrea Modica, Ivan Scotti

[Laurent Excoffier]

Hi,

I had a look at your files and you have several problems

- negative inbreeding coefficients (they are not FIS)

- potentially asymmetric migration matrix (ok, but I symmetrized it)

- number of independent loci set to 12, but shoudl be 1 for FREQ data

- Inconsistent sample sizes e.g. pop0 has 50 in tpl file and 25 in obs file

- unused parameter in est file (MUTRATE) (okay but will slow things down)

However, the real trouble is that you have too negative growth rate, such that deme size goes to zero while still having lineages in it. If deme size goes to zero, probability of coalescences go to infinite, and therefore in fsc27, demes of size zero are killed

Maybe it is a bug on my side, but I think your groath rates are too large and you shoudl make sure that their range is within reasonable values.

When setting growth rates to zero and using the attached corrected tpl files, the program doe not crash, but still is running super slow given the high migration rates you are specifying. Lineages keep moving every generations and take ages before coalescing

Hope it helps

laurent

[Ivan Scotti]

Dear Laurent,

thank you once again for your feedback - this only goes to prove that the process of learning how to use FastSimCoal is a never-ending one (which adds to the fun!), and that the number of mistakes one can do is also rather large (which is rather humbling).
Indeed, there was a series of macroscopic errors with our SFS tables, which, being both multi-population and based on minor alleles, needed to be computed over the entire number of sampled chromosomes in each population, but for some reason I had my scripts all wrong. So we've fixed that, and also reviewed our TPL and EST files for potential errors. Yet our full model still would not run.

We then took a more systematic approach to find what was wrong, starting from the simplest model and progressively making it more complicated.

We indeed managed to successfully run the most basic model - four pops, three events corresponding to population mergers, growth rate set to zero, no other "frills" - based on the attached files and with the command line:

fsc2702 -t inputTemplate_v11.tpl -n 1000 -m -e inputTemplate_v11.est -M

We correctly obtained best likelihoods and parameter estimates, after about 250 iterations (and 4 hours of machine time, using 12 batches and 1 thread).

Then, based on this encouraging result, we started making the model slightly more complex: we started the simulation with a negative growth rate, that stayed the same throughout the population merger events, and then fell very quickly to zero after the mergers. As a consequence, the TPL was modified as follows (changes relative to the attached file are in bold):

//Number of population samples (demes)

4 samples to simulate

//Population effective sizes (number of genes)

POPSIZE0

POPSIZE1

POPSIZE2

POPSIZE3

//Sample sizes

40

50

50

52

//Growth rates : negative growth implies population expansion

GROWTHRATE0

GROWTHRATE0

GROWTHRATE0

GROWTHRATE0

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

1

//migration matrix

0.000 M01 M02 M03 

M01 0.000 M12 M13 

M02 M12 0.000 M23 

M03 M13 M23 0.000

//historical event: time, source, sink, migrants, new size, growth rate, migr. matrix NOTE: order of event 1 and 2 to be estimated

4 historical event

TIME1 0 1 1 1 GROWTHRATE0 0

TIME2 2 1 1 1 GROWTHRATE0 0

TIME3 3 1 1 1 GROWTHRATE0 0

TIMEA 1 1 0 1 0 0

//Number of independent loci [chromosome]

1 0

//Per chromosome: Number of linkage blocks

1

//per Block: data type, num loci, rec. rate and mut rate + optional parameters

FREQ 1 0 1e-8 OUTEXP

and the EST was changed as follows  (changes relative to the attached file are in bold):

// Priors and rules file

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

[PARAMETERS]

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

//all N are in number of haploid individuals

1 POPSIZE0 logunif 100 30000 output

1 POPSIZE1 logunif 100 30000 output

1 POPSIZE2 logunif 100 30000 output

1 POPSIZE3 logunif 50 30000 output

0 M01 logunif 0.001 0.1 output

0 M02 logunif 0.001 0.1 output

0 M03 logunif 0.001 0.1 output

0 M12 logunif 0.001 0.1 output

0 M13 logunif 0.001 0.1 output

0 M23 logunif 0.001 0.1 output

1 TIMEA unif 2 20 output

1 TIME3 logunif 1 TIMEA output paramInRange

1 TIME1 logunif 1 TIMEA output paramInRange

1 TIME2 logunif 1 TIMEA output paramInRange

0 GROWTHRATE0 unif -0.5 -0.05 output

[COMPLEX PARAMETERS]

This does not look like a major change, but clearly we must have introduced a mistake even with our infinitesimal change: with the same command line as above, the program immediately core-dumped. 

Notice that all events happen over a very short time, so it seems unlikely that, e.g., the new growth rate(s) have done anything weird to the populations (plus, they expand going backward, so they cannot disappear).

So back to square one. Or, should I say, to square 2, because we managed at least one successful model! 

If you still have some patience left for us, we'd be happy to check with you what may still be wrong.

Kind regards,

Ivan, Andrea

[Ivan Scotti]

...and of course, misunderstanding what the manual says about forward and backward things does not help. The growth rates must be positive if I want expansion backward. So I have modified the PAR, which now looks like this (I took the opportunity to reduce mig rates as well, just in case):

// Priors and rules file

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

[PARAMETERS]

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

//all N are in number of haploid individuals

1 POPSIZE0 logunif 100 30000 output

1 POPSIZE1 logunif 100 30000 output

1 POPSIZE2 logunif 100 30000 output

1 POPSIZE3 logunif 50 30000 output

0 M01 logunif 0.001 0.01 output

0 M02 logunif 0.001 0.01 output

0 M03 logunif 0.001 0.01 output

0 M12 logunif 0.001 0.01 output

0 M13 logunif 0.001 0.01 output

0 M23 logunif 0.001 0.01 output

1 TIMEA unif 2 20 output

1 TIME3 logunif 1 TIMEA output paramInRange

1 TIME1 logunif 1 TIMEA output paramInRange

1 TIME2 logunif 1 TIMEA output paramInRange

0 GROWTHRATE0 unif 0.01 0.05 output

[COMPLEX PARAMETERS]

So now we're out of "straight" core dumps, but still something's wrong - it seems that the program does not "see" the historical event setting the growth rate at zero, and populations (all of them, even though they should have merged) go to infinity and eventually we have our dear core dump, of course.

Here's an excerpt of the screen output:

fsc2702 -t inputTemplate_v11.tpl -n 1000 -m -e inputTemplate_v11.est -M

fastsimcoal was invoked with the following command line arguments:

fsc2702 -t inputTemplate_v11.tpl -n 1000 -m -e inputTemplate_v11.est -M 

Unspecified number of loops to perform. Setting it to 20

Reading parameter settings file inputTemplate_v11.est...

done reading parameter settings file!

Random generator seed : 675180

Deme sizes

Deme 0 5960

Deme 1 2307

Deme 2 916

Deme 3 9184

Sample sizes

Deme 0 40

Deme 1 50

Deme 2 50

Deme 3 52

Sample ages

Deme 0 0

Deme 1 0

Deme 2 0

Deme 3 0

Growth rates

Deme 0 0.0496717

Deme 1 0.0496717

Deme 2 0.0496717

Deme 3 0.0496717

Migration matrix 0 (note that migration rates have been pre-multiplied by 1e+06n

 0.0000000 4462.4938328 4312.7245585 2396.4115292 

4462.4938328  0.0000000 2834.9554236 6641.9203710 

4312.7245585 2834.9554236  0.0000000 1684.3880678 

2396.4115292 6641.9203710 1684.3880678  0.0000000 

Historical events

Event 0

#Time              : 3

#Source            : 0

#Sink              : 1

#Migrants          : 1.0000000

#New relative size : 1.0000000

#New growth rate   : 0.0496717

#New migr. matrix  : 0

Event 1

#Time              : 10

#Source            : 2

#Sink              : 1

#Migrants          : 1.0000000

#New relative size : 1.0000000

#New growth rate   : 0.0496717

#New migr. matrix  : 0

Event 2

#Time              : 12

#Source            : 3

#Sink              : 1

#Migrants          : 1.0000000

#New relative size : 1.0000000

#New growth rate   : 0.0496717

#New migr. matrix  : 0

Event 3

#Time              : 18

#Source            : 1

#Sink              : 1

#Migrants          : 0.0000000

#New relative size : 1.0000000

#New growth rate   : 0.0000000

#New migr. matrix  : 0

Number of independent loci to simulate : 1

with the same chromosomal structure

Number of linkage blocks to simulate in structure 1: 1

Data type: FREQ   :  recombination and mutation rates

0.0000000000   0.0000000100   

   Output expected MAF spectrum

   Do not output expected DAF spectrum

Population growth detected in input file

Estimating model parameters using 12 batches and 1 thread

TDemeCollection::adjustDemeSizesCT: Deme size set to MAXFLOAT in deme 1 at time 1880 as it was reaching INFINITY!

Check your par file...!!!

[...]

TDemeCollection::adjustDemeSizesCT: Deme size set to MAXFLOAT in deme 1 at time 13645 as it was reaching INFINITY!

Check your par file...!!!

TDemeCollection::adjustDemeSizesCT: Deme size set to MAXFLOAT in deme 2 at time 13645 as it was reaching INFINITY!

Check your par file...!!!

TDemeCollection::adjustDemeSizesCT: Deme size set to MAXFLOAT in deme 0 at time 13652 as it was reaching INFINITY!

Check your par file...!!!

TDemeCollection::adjustDemeSizesCT: Deme size set to MAXFLOAT in deme 1 at time 13652 as it was reaching INFINITY!

Check your par file...!!!

TDemeCollection::adjustDemeSizesCT: Deme size set to MAXFLOAT in deme 2 at time 13652 as it was reaching INFINITY!

Check your par file...!!!

Erreur de segmentation (core dumped)

Well, I guess we're making progress, but still...

Cheers

Ivan

[Ivan Scotti]

I think I've found the problem: it had escaped me that populations having merged with another one had to be actively killed (P. 29 of the manual) or else they stay around and continue their zombie life - and undergo all possible messy positive and negative growth leftover from previous steps, causing the program to crash. This version of the TPL ran smoothly and in less than two minutes:

7 historical event

TIME1 0 1 1 1 GROWTHRATE0 0

TIME1 0 0 0 0 GROWTHRATE0 0 // kills pop0

TIME2 2 1 1 1 GROWTHRATE0 0

TIME2 2 2 0 0 GROWTHRATE0 0 // kills pop2

TIME3 3 1 1 1 GROWTHRATE0 0

TIME3 3 3 0 0 GROWTHRATE0 0 // kills pop3

TIMEA 1 1 1 1 0 0 nomig

//Number of independent loci [chromosome]

1 0

//Per chromosome: Number of linkage blocks

1

//per Block: data type, num loci, rec. rate and mut rate + optional parameters

FREQ 1 0 1e-8 OUTEXP

Now I think we can proceed with the rest of the model...

Kind regards

Ivan