The issue

During our discussion yesterday Yaniv pointed out that the -em argument in ms was designed to set a hybridization from the specified time and extending rootward from that point. After some testing, this is indeed what is happening. In my opinion, this leaves two options to represent the kind of gene flow we are interested in within the ms framework: 1) Use two -em arguments: the first closer to the tips turns migration on; the second closer to the roots turns migration off. 2) Use a combination of -es and -ej arguments: as we move from a sampled tip we split the pipe into one that follows the backbone tree and another that joins with another pipe at some point further in the past.

I believe that these should have identical effects if 1) the es and ej statements are very close in time (this eliminates the opportunity for drift in the hybridizing lineage) and 2) the two -em arguments are very close in time (I think to truly be equal it would em would have to turn on migration for just a single generation. I did a simple test of this just comparing the frequency of sister relationships among tips which can be seen below.

Here is the backbone tree that we are working with

To evaluate the effect of our different ms runs we can just look at the frequency of sister pairings of tips 1,2 or 1,3 or 2,3.

We can compare 3 sets of ms commands and we will simulate 1000 gene trees under each command to get a feel for what is going on.

NO MIGRATION

ms 3 1000 -I 3 1 1 1 -ej 1.625 1 2 -ej 8.901 2 3 -T > no.mig.tree

MIGRATION START STOP

ms 3 1000 -I 3 1 1 1 -em 0.12 3 1 50 -em 0.13 3 1 0 -ej 1.625 1 2 -ej 8.901 2 3 -T   > em.mig.tree

ES MIGRATION

ms 3 1000 -I 3 1 1 1 -es 0.12 1 0.5 -ej 0.13 4 3 -ej 1.625 1 2 -ej 8.901 2 3 -T   > es.mig.tree

and this is what the results look like (again values are % of gene trees having the indicated sister pairings on the left):

no.migration -em.migration -es.migration
1 and 2 100 60.4 51.9
1 and 3 0 39.4 48.1
2 and 3 0 0.2 0.0

That result is actually after playing with the intensity argument for -em. I am still not a 100% sure what this is in units of. It looks to me that the more straightforward solution is to simply use the es ej option. Josh and I will work on implementing that solution and update everybody next week cheers.