Split and bottleneck model

[Panos Zerv]

Hi Ryan and all wonderful dadi community,

I have 1D models on two populations, one of them follows a two-epoch model and the other a three-epoch. I am trying to integrate these into a 2D split model to calculate the time of the split between the populations. So here are my questions:

1. Do I need to integrate the 1D models in the first place? Or Just a simple split model would do the trick? I am not sure if recent demographic events would affect the way dadi calculates past events.

2. Assuming that I have to do it, how would I have to specify the time of each event? For example, let's assume that at a time T1 in the past they were split. Then population 1 started going through a bottleneck. Following, population 2 started going through a bottleneck which lasts until today. Finally, population 1 recovers.

def split_three_epoch(params, ns, pts):
        """
Population size change after split with two periods of migration after size change.
        nu1: Size of population 1 after split
        nu1b: bottleneck of pop1

        nu1a: Final size of population 1.
        nu2: Size of population 2 after split.
        nu2a: Final size of poulation size 2.
        F1: Inbreeding in population 1.
        F2: Inbreeding in population 2.
        T1: Time of split (in units of 2*Na generations)

        T2: Start of bottleneck (population 1) 

        T3: Start of bottleneck (population 2)

        T4: Recovery (population 1)

        """
nu1, nu2, nu1a, nu1b, nu2a, F1, F2, T1, T2, T3, T4= params
        xx = dadi.Numerics.default_grid(pts)
        phi = dadi.PhiManip.phi_1D(xx)
        phi = dadi.PhiManip.phi_1D_to_2D(xx, phi)
        phi = dadi.Integration.two_pops(phi, xx, T1, nu1=nu1, nu2=nu2) # sizes after split
        phi = dadi.Integration.two_pops(phi, xx, T2, nu1=nu1b, nu2=nu2) # start of bottleneck pop1
        phi = dadi.Integration.two_pops(phi, xx, T3, nu1=nu1b, nu2=nu2a) # start of bottleneck pop2
        phi = dadi.Integration.two_pops(phi, xx, T4, nu1=nu1a, nu2=nu2a) # recovery pop1
        fs = dadi.Spectrum.from_phi_inbreeding(phi, ns, (xx, xx), (F1, F2), (2, 2))
        return fs

So to calculate the total time since the split I would have to add up T1+T2+T3+T4? Or is time calculated separately for the two branches after the split? My tests so far show that this is not the case. Finally, would I need to use an if... statement and modify the model accordingly in case e.g. the bottleneck of population 2 came after the recovery of population 1?

Thank you,

Panagiotis

[Ryan Gutenkunst]

Hello Panagiotis,

1. Mis-modeled recent demographic events can bias the inference of older events.

2. In answer to your question, Integration operations define time intervals. So your model below would run for time T1+T2+T3+T4 since the split. It’s actually a bit complex to write the model you propose, since you have multiple events that can switch orders, while you also need to constrain parameters to make sense. (For example, ensuring that the bottleneck lengths fit in the time since the split.) It’s also worth thinking about whether you need that complexity. If you include a three epoch model after the split, you’re really modeling 4 size epochs for that population. You might consider simplifying the model so that there’s a size change at the split for each population, then just a single additional size change for population 1.

Best,
Ryan

> --
> You received this message because you are subscribed to the Google Groups "dadi-user" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to dadi-user+...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/dadi-user/425c0be0-9011-4224-8912-f29f9043b727n%40googlegroups.com.