Trouble generating expected PSMC output from constant population model

3 views
Skip to first unread message

Jared Grummer

unread,
Feb 18, 2025, 8:23:37 PMFeb 18
to msprime-users
Hello,

Before testing an empirical model in PSMC, we are trying to simulate a simple demographic model to make sure we are getting out of PSMC what we believe should match the simulated model. One of the main problems we're having now is getting the x-axis time value inferred from the psmc to match what we simulate in msprime. For example, we have this model:

demographic_model = msprime.Demography()
demographic_events=[
demographic_model.add_population(initial_size = 5000, growth_rate = 0), #Size at t=0, then shrinks to ~11k individuals by gen. 40k
demographic_model.add_population_parameters_change(time=600000, growth_rate = 0, population = 0), #Stop growth at 200k gens ago
demographic_model.add_population_parameters_change(time=1000000, initial_size=10000, population = 0)
]

Output from the demography debugger:

Screenshot 2025-02-18 at 5.53.47 PM.png

I assume this means the population's origin should go back to 1e6 or 6e5 generations ago. We use these values:

mut_rate = 8e-8
r_chrom = 1e-7
chr_length = int(1e7)

But the psmc output looks like this:

Screenshot 2025-02-18 at 5.57.07 PM.png

It's a little curious that the Ne is ~10x less than that simulated, and the x-axis is ~10x less than what it should be. But I'm not quite sure why. This model used multiple (independently recombining) chromosomes, but we get similar results with a model simulating a single chromosome. I wonder if there aren't enough recombination events to infer the time that long ago? I set the chromosome length and rec. rate to be the inverse so there is ~1 recom/chr. expected.

What am I missing? Any insight would be greatly appreciated!

Thanks in advance,
Jared

Reply all
Reply to author
Forward
0 new messages