1554.py

import math
import numpy as np
import msprime

seed = 1102
divergence_time = 50
####Variables####
# These are the DIPLOID counts
sequence_length = 100000
Ne_A = 10000
Ne_B = Ne_A / 2
peripheral_pops = 500
pi_1 = 0.01
pi_2 = 0.1
tol = 0.01
###################
demography = msprime.Demography()
descendant_pops = []
demography.add_population(name="ancestral_population", initial_size=Ne_A)
for i in range(0, peripheral_pops):
    name = "peripheral_%s" % i
    demography.add_population(name=name, initial_size=Ne_B)
    demography.add_mass_migration(
        time=divergence_time, source=name, dest="ancestral_population", proportion=1
    )
    descendant_pops.append("peripheral_%s" % i)

# Do a recombination map with a single recombination hotspot in the middle (ask Gil why!)
positions = range(0, sequence_length)
rate_map = np.zeros(sequence_length - 1)
# Replace 0 with hotspot in middle
rate_map[math.floor(sequence_length / 2)] = 0.9

tree = msprime.sim_ancestry(
    recombination_rate=msprime.RateMap(position=positions, rate=rate_map),
    sequence_length=sequence_length - 1,
    demography=demography,
    samples={i: 500 for i in descendant_pops},
    random_seed=seed,
)