Model of meiotic drive in ZW systems

Author

Darren Irwin

Published

May 25, 2025

Here I will develop a model for simulating how meiotic drive affects diversity within a population and sweeps between populations.

The code blocks below are written in the Julia programming language. If you’ve never used Julia, you can learn more and easily install for free here.

This website is a Quarto project and each page is a Quarto notebook, which can run and display the results of Julia (or other) code blocks, along with text narration, and output in html, pdf, Word, etc.To learn more about Quarto websites visit https://quarto.org/docs/websites.

Citation

The scripts, data, and figures shown in this website were used as the basis for the research article listed below, which should be cited as the source of information from this website:

Irwin, D. 2025. The driving W hypothesis for low within-population mitochondrial DNA diversity and between-population mitochondrial capture. BioRxiv.

Load required Julia packages

If you don’t already have these packages installed, you will need to load packages used in the script, so run what is in this section below. It will take some time to install and precompile the packages:

import Pkg
Pkg.add("StatsBase")
Pkg.add("Plots")

Now actually load those packages into the Julia session:

using StatsBase  # needed for "sample" function
using Plots  # for graphing results

Set a random number seed, so we get the same series of random numbers as in the paper (reproducible randomness! :)

using Random
Random.seed!(379)
TaskLocalRNG()

Building the model

In the code blocks below, I will create data types (using the struct command) that represent chromosomes, individuals, and populations. I also build functions that create and govern interactions between these entities.

Create a data structure to store info about a chromosome

This will store:

  • kind of chromosome (Z or W, or another name representing an autosome)
  • ability to connect to spindle to egg
  • survival fitness loss of having this chromosome (to be subtracted from an idealized chromosome of fitness = 1)
  • ability to suppress drivers

The cell below creates a data type (Chrom) that can store four fields with the names give (corresponding to the above):

struct Chrom
    kind
    spindleAttraction
    fitnessLoss
    suppressStrength
end
Now actually create some chromosome types

We now create some specific kinds of chromosomes, storing them as instances of the Chrom data type:

driveStrength = 1.5  # sets this to be the drive strength of each driving chrom type
boringW = Chrom("W", 1.0, 0.0, 0.0) # a simple nothing-special W chromosome
drivingW = Chrom("W", driveStrength, 0.0, 0.0) # a driving W
boringZ = Chrom("Z", 1.0, 0.0, 0.0) 
drivingZ = Chrom("Z", driveStrength, 0.0, 0.0)
boringAutosome1 = Chrom("1", 1.0, 0.0, 0.0)
drivingAutosome1 = Chrom("1", driveStrength, 0.0, 0.0)
driveSuppressorAutosome1 = Chrom("1", 1.0, 0.0, 0.5) # autosome causing partial suppression of drive (half suppression when homozygous)
driveCancellerAutosome1 = Chrom("1", 1.0, 0.0, 1.0)  # when homozygous, completely suppresses drive
driveStrengthStrong = 3.0 # sets this as the stronger driving strength for some sims
drivingWstrong = Chrom("W", driveStrengthStrong, 0.0, 0.0)
driveStrengthWeak = 1.1 # sets this as the weaker driving strength for some sims
drivingWweak = Chrom("W", driveStrengthWeak, 0.0, 0.0)
Chrom("W", 1.1, 0.0, 0.0)

In the above, “boring” means non-driving and non-suppressing chromosomes.

Store some of these kinds of chromosomes in a vector, for easy reference later:

chromKinds = [boringZ, 
            drivingZ,
            boringW,
            drivingW,
            drivingWstrong,
            drivingWweak,
            boringAutosome1,
            drivingAutosome1]
8-element Vector{Chrom}:
 Chrom("Z", 1.0, 0.0, 0.0)
 Chrom("Z", 1.5, 0.0, 0.0)
 Chrom("W", 1.0, 0.0, 0.0)
 Chrom("W", 1.5, 0.0, 0.0)
 Chrom("W", 3.0, 0.0, 0.0)
 Chrom("W", 1.1, 0.0, 0.0)
 Chrom("1", 1.0, 0.0, 0.0)
 Chrom("1", 1.5, 0.0, 0.0)

Make a struct for storing info about an individual

This will store:

  • Sex chromosome kind from egg.
  • Sex chromosome kind from sperm.
  • Autosome 1 kind from egg.
  • Autosome 1 kind from sperm.
struct Indiv
    sexChromFromEgg 
    sexChromFromSperm
    autosome1FromEgg
    autosome1FromSperm
end
Now create some kinds of individuals

Below, we construct individuals by specifying which kind of chromosome is in each of the four slots above:

begin
femaleBoring = Indiv(boringW, boringZ, boringAutosome1, boringAutosome1)
femaleDriveW = Indiv(drivingW, boringZ, boringAutosome1, boringAutosome1)
femaleDriveWstrong = Indiv(drivingWstrong, boringZ, boringAutosome1, boringAutosome1)
femaleDriveWweak = Indiv(drivingWweak, boringZ, boringAutosome1, boringAutosome1)
femaleDriveZ = Indiv(boringW, drivingZ, boringAutosome1, boringAutosome1)
femaleDriveWZ = Indiv(drivingW, drivingZ, boringAutosome1, boringAutosome1)
femaleFirstDriveA = Indiv(boringW, boringZ, drivingAutosome1, boringAutosome1)
femaleSecondDriveA = Indiv(boringW, boringZ, boringAutosome1, drivingAutosome1)
femaleTwoDriveA = Indiv(boringW, boringZ, drivingAutosome1, drivingAutosome1)
femaleTwoDriveSuppressorA = Indiv(boringW, boringZ, driveSuppressorAutosome1, driveSuppressorAutosome1)
femaleTwoDriveCancellerA = Indiv(boringW, boringZ, driveCancellerAutosome1, driveCancellerAutosome1)
maleBoring = Indiv(boringZ, boringZ, boringAutosome1, boringAutosome1)
maleFirstDriveZ = Indiv(drivingZ, boringZ, boringAutosome1, boringAutosome1)
maleSecondDriveZ = Indiv(boringZ, drivingZ, boringAutosome1, boringAutosome1)
maleTwoDriveZ = Indiv(drivingZ, drivingZ, boringAutosome1, boringAutosome1)
maleFirstDriveA = Indiv(boringZ, boringZ, drivingAutosome1, boringAutosome1)
maleSecondDriveA = Indiv(boringZ, boringZ, boringAutosome1, drivingAutosome1)
maleTwoDriveA = Indiv(boringZ, boringZ, drivingAutosome1, drivingAutosome1)
maleTwoDriveSuppressorA = Indiv(boringZ, boringZ, driveSuppressorAutosome1, driveSuppressorAutosome1)
maleTwoDriveCancellerA = Indiv(boringW, boringZ, driveCancellerAutosome1, driveCancellerAutosome1)
end
Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 1.0), Chrom("1", 1.0, 0.0, 1.0))

Write a mating function

Given one female (ZW) and one male (ZZ) parent, the function will produce one offspring, incorporating transmission advantage as calculated by the ratio of spindleAttraction of each chromosome compared to total spindleAttraction (modified by suppression).

function makeOneOffspring(female, male) 
    # determine sex chromosome of egg, weighted by spindleAttraction
    # but moderated by suppressors
    driveAdvantage = female.sexChromFromEgg.spindleAttraction / 
                    (female.sexChromFromEgg.spindleAttraction + female.sexChromFromSperm.spindleAttraction)
    # determine drive suppression strength by averaging two alleles
    suppressionStrength = (female.autosome1FromEgg.suppressStrength +   
                            female.autosome1FromSperm.suppressStrength) / 2
    if suppressionStrength > 1.0
        suppressionStrength = 1.0
    end
    transProbW = suppressionStrength * 0.5 + (1 - suppressionStrength) * 
                driveAdvantage 
    eggSexChrom = sample([female.sexChromFromEgg, female.sexChromFromSperm], 
        Weights([transProbW, 1 - transProbW]))
    # determine sex chromosome of sperm, with no drive (Mendelian)
    spermSexChrom = sample([male.sexChromFromEgg, male.sexChromFromSperm])
    # determine autosome1 in egg, weighted by spindleAttraction
    # but moderated by suppressors (I don't think I will use this much, but for completeness putting it here now)
    driveAdvantageAutosome1 = female.autosome1FromEgg.spindleAttraction / 
                    (female.autosome1FromEgg.spindleAttraction +
                    female.autosome1FromSperm.spindleAttraction)
    # use same suppression strength calculated above
    transProbChromFromEgg = suppressionStrength * 0.5 + (1 - suppressionStrength) * 
                driveAdvantageAutosome1
    eggAutosome1 = sample([female.autosome1FromEgg, female.autosome1FromSperm], 
                    Weights([transProbChromFromEgg, 1 - transProbChromFromEgg]))
    # determine autosome1 in sperm, Mendelian
    spermAutosome1 = sample([male.autosome1FromEgg, male.autosome1FromSperm])
    
    return Indiv(eggSexChrom, spermSexChrom, eggAutosome1, spermAutosome1)
end
makeOneOffspring (generic function with 1 method)

Make a function to test the function above:

function test_makeOneOffspring(female, male, reps)
    offspringSex = fill("N", reps)
    for i in 1:reps
        kid = makeOneOffspring(female, male)
        if kid.sexChromFromEgg.kind == "W"
            offspringSex[i] = "F"
        else offspringSex[i] = "M"
        end
    end
    proportionFemale = (sum(offspringSex .== "F") / reps)
    println("The proportion of offspring who are female is " * string(proportionFemale))
end
test_makeOneOffspring (generic function with 1 method)

Now run some tests:

test_makeOneOffspring(femaleDriveW, maleBoring, 1_000_000)
The proportion of offspring who are female is 0.599241
test_makeOneOffspring(femaleDriveWZ, maleBoring, 1_000_000)
The proportion of offspring who are female is 0.499607
femaleDriveWandSuppressorTwoCopies = Indiv(drivingW, boringZ, driveSuppressorAutosome1, driveSuppressorAutosome1)
test_makeOneOffspring(femaleDriveWandSuppressorTwoCopies, maleBoring, 1_000_000)
The proportion of offspring who are female is 0.550221
femaleDriveWandCancellerTwoCopies = Indiv(drivingW, boringZ, driveCancellerAutosome1, driveCancellerAutosome1)
test_makeOneOffspring(femaleDriveWandCancellerTwoCopies, maleBoring, 1_000_000)
The proportion of offspring who are female is 0.499685

Make a Population type for storing info about whole populations

This will contain two vectors (females, and males), in which the cells contain individuals (i.e., of type Indiv), each of which contains chromosomes (i.e., of type Chrom).

mutable struct Population
    popFemales
    popMales
end
Now construct some populations:

pop1 has 5% driving W:

pop1setPropDrivingW = 0.05
pop1SizeFemale = 10_000
pop1SizeMale = 10_000
pop1Females = vcat(fill(femaleDriveW, Int(round(pop1setPropDrivingW * pop1SizeFemale))), fill(femaleBoring, Int(round((1-pop1setPropDrivingW) * pop1SizeFemale))))
pop1Males = vcat(fill(maleBoring, pop1SizeMale))
pop1 = Population(pop1Females, pop1Males)
Population(Indiv[Indiv(Chrom("W", 1.5, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.5, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.5, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.5, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.5, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.5, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.5, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.5, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.5, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.5, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0))  …  Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0))], Indiv[Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0))  …  Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0))])

pop1B has 5% driving W (of the weak variety):

pop1BsetPropDrivingWweak = 0.05
pop1BSizeFemale = 10_000
pop1BSizeMale = 10_000
pop1BFemales = vcat(fill(femaleDriveWweak, Int(round(pop1BsetPropDrivingWweak * pop1BSizeFemale))), fill(femaleBoring, Int(round((1-pop1BsetPropDrivingWweak) * pop1BSizeFemale))))
pop1BMales = vcat(fill(maleBoring, pop1BSizeMale))
pop1B = Population(pop1BFemales, pop1BMales)
Population(Indiv[Indiv(Chrom("W", 1.1, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.1, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.1, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.1, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.1, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.1, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.1, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.1, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.1, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.1, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0))  …  Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0))], Indiv[Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0))  …  Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0))])

pop2 has 5% driving Z. Here we set the proportions of males with different kinds of Z by assuming Hardy-Weinberg equilibrium:

pop2setPropDrivingZ = 0.05
pop2SizeFemale = 10_000
pop2SizeMale = 10_000
pop2Females = vcat(fill(femaleDriveZ, Int(round(pop2setPropDrivingZ * pop2SizeFemale))), 
                fill(femaleBoring, Int(round((1-pop2setPropDrivingZ) * pop2SizeMale))))
pop2numMaleTwoDriveZ = Int(round(pop2setPropDrivingZ^2 * pop2SizeMale))
pop2numMaleFirstDriveZ = Int(round(pop2setPropDrivingZ * (1-pop2setPropDrivingZ) * pop2SizeMale))
pop2numMaleSecondDriveZ = Int(round((1-pop2setPropDrivingZ) * pop2setPropDrivingZ * pop2SizeMale))
pop2numMaleBoring = Int(round((1-pop2setPropDrivingZ)^2 * pop2SizeMale))
pop2Males = vcat(fill(maleTwoDriveZ, pop2numMaleTwoDriveZ),
                fill(maleFirstDriveZ, pop2numMaleFirstDriveZ),
                fill(maleSecondDriveZ, pop2numMaleSecondDriveZ),
                fill(maleBoring, pop2numMaleBoring))
pop2 = Population(pop2Females, pop2Males)
Population(Indiv[Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.5, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.5, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.5, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.5, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.5, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.5, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.5, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.5, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.5, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.5, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0))  …  Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0))], Indiv[Indiv(Chrom("Z", 1.5, 0.0, 0.0), Chrom("Z", 1.5, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.5, 0.0, 0.0), Chrom("Z", 1.5, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.5, 0.0, 0.0), Chrom("Z", 1.5, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.5, 0.0, 0.0), Chrom("Z", 1.5, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.5, 0.0, 0.0), Chrom("Z", 1.5, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.5, 0.0, 0.0), Chrom("Z", 1.5, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.5, 0.0, 0.0), Chrom("Z", 1.5, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.5, 0.0, 0.0), Chrom("Z", 1.5, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.5, 0.0, 0.0), Chrom("Z", 1.5, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.5, 0.0, 0.0), Chrom("Z", 1.5, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0))  …  Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0))])

pop3 has 5% driving autosome (A). We again use HW equilibrium for the initial setup of the frequencies of different individual types:

pop3setPropDrivingA = 0.05
pop3SizeFemale = 10_000
pop3SizeMale = 10_000

pop3numFemaleTwoDriveA = Int(round(pop3setPropDrivingA^2 * pop3SizeFemale))
pop3numFemaleFirstDriveA = Int(round(pop3setPropDrivingA * (1-pop3setPropDrivingA) * pop3SizeFemale))
pop3numFemaleSecondDriveA = Int(round((1-pop3setPropDrivingA) * pop3setPropDrivingA * pop3SizeFemale))
pop3numFemaleBoring = Int(round((1-pop3setPropDrivingA)^2 * pop3SizeFemale))

pop3Females = vcat(fill(femaleTwoDriveA, pop3numFemaleTwoDriveA),
                fill(femaleFirstDriveA, pop3numFemaleFirstDriveA),
                fill(femaleSecondDriveA, pop3numFemaleSecondDriveA),
                fill(femaleBoring, pop3numFemaleBoring))

pop3numMaleTwoDriveA = Int(round(pop3setPropDrivingA^2 * pop3SizeMale))
pop3numMaleFirstDriveA = Int(round(pop3setPropDrivingA * (1-pop3setPropDrivingA) * pop3SizeMale))
pop3numMaleSecondDriveA = Int(round((1-pop3setPropDrivingA) * pop3setPropDrivingA * pop3SizeMale))
pop3numMaleBoring = Int(round((1-pop3setPropDrivingA)^2 * pop3SizeMale))

pop3Males = vcat(fill(maleTwoDriveA, pop3numMaleTwoDriveA),
                fill(maleFirstDriveA, pop3numMaleFirstDriveA),
                fill(maleSecondDriveA, pop3numMaleSecondDriveA),
                fill(maleBoring, pop3numMaleBoring))
pop3 = Population(pop3Females, pop3Males)
Population(Indiv[Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.5, 0.0, 0.0), Chrom("1", 1.5, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.5, 0.0, 0.0), Chrom("1", 1.5, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.5, 0.0, 0.0), Chrom("1", 1.5, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.5, 0.0, 0.0), Chrom("1", 1.5, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.5, 0.0, 0.0), Chrom("1", 1.5, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.5, 0.0, 0.0), Chrom("1", 1.5, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.5, 0.0, 0.0), Chrom("1", 1.5, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.5, 0.0, 0.0), Chrom("1", 1.5, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.5, 0.0, 0.0), Chrom("1", 1.5, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.5, 0.0, 0.0), Chrom("1", 1.5, 0.0, 0.0))  …  Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0))], Indiv[Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.5, 0.0, 0.0), Chrom("1", 1.5, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.5, 0.0, 0.0), Chrom("1", 1.5, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.5, 0.0, 0.0), Chrom("1", 1.5, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.5, 0.0, 0.0), Chrom("1", 1.5, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.5, 0.0, 0.0), Chrom("1", 1.5, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.5, 0.0, 0.0), Chrom("1", 1.5, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.5, 0.0, 0.0), Chrom("1", 1.5, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.5, 0.0, 0.0), Chrom("1", 1.5, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.5, 0.0, 0.0), Chrom("1", 1.5, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.5, 0.0, 0.0), Chrom("1", 1.5, 0.0, 0.0))  …  Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0))])

pop4 has 5% driving W and 5% driving Z:

pop4setPropDrivingW = 0.05
pop4setPropDrivingZ = 0.05
pop4SizeFemale = 10_000
pop4SizeMale = 10_000

pop4numFemaleDriveWZ = Int(round(pop4setPropDrivingW * pop4setPropDrivingZ * pop4SizeFemale))
pop4numFemaleDriveW = Int(round(pop4setPropDrivingW * (1-pop4setPropDrivingZ) * pop4SizeFemale))
pop4numFemaleDriveZ = Int(round((1-pop4setPropDrivingW) * pop4setPropDrivingZ * pop4SizeFemale))
pop4numFemaleBoring = Int(round((1-pop4setPropDrivingW) * (1-pop4setPropDrivingZ) * pop4SizeFemale))

pop4Females = vcat(fill(femaleDriveWZ, pop4numFemaleDriveWZ),
                fill(femaleDriveW, pop4numFemaleDriveW),
                fill(femaleDriveZ, pop4numFemaleDriveZ),
                fill(femaleBoring, pop4numFemaleBoring))

pop4numMaleTwoDriveZ = Int(round(pop4setPropDrivingZ^2 * pop4SizeMale))
pop4numMaleFirstDriveZ = Int(round(pop4setPropDrivingZ * (1-pop4setPropDrivingZ) * pop4SizeMale))
pop4numMaleSecondDriveZ = Int(round((1-pop4setPropDrivingZ) * pop4setPropDrivingZ * pop4SizeMale))
pop4numMaleBoring = Int(round((1-pop4setPropDrivingZ)^2 * pop4SizeMale))
pop4Males = vcat(fill(maleTwoDriveZ, pop4numMaleTwoDriveZ),
                fill(maleFirstDriveZ, pop4numMaleFirstDriveZ),
                fill(maleSecondDriveZ, pop4numMaleSecondDriveZ),
                fill(maleBoring, pop4numMaleBoring))
pop4 = Population(pop4Females, pop4Males)
Population(Indiv[Indiv(Chrom("W", 1.5, 0.0, 0.0), Chrom("Z", 1.5, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.5, 0.0, 0.0), Chrom("Z", 1.5, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.5, 0.0, 0.0), Chrom("Z", 1.5, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.5, 0.0, 0.0), Chrom("Z", 1.5, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.5, 0.0, 0.0), Chrom("Z", 1.5, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.5, 0.0, 0.0), Chrom("Z", 1.5, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.5, 0.0, 0.0), Chrom("Z", 1.5, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.5, 0.0, 0.0), Chrom("Z", 1.5, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.5, 0.0, 0.0), Chrom("Z", 1.5, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.5, 0.0, 0.0), Chrom("Z", 1.5, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0))  …  Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0))], Indiv[Indiv(Chrom("Z", 1.5, 0.0, 0.0), Chrom("Z", 1.5, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.5, 0.0, 0.0), Chrom("Z", 1.5, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.5, 0.0, 0.0), Chrom("Z", 1.5, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.5, 0.0, 0.0), Chrom("Z", 1.5, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.5, 0.0, 0.0), Chrom("Z", 1.5, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.5, 0.0, 0.0), Chrom("Z", 1.5, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.5, 0.0, 0.0), Chrom("Z", 1.5, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.5, 0.0, 0.0), Chrom("Z", 1.5, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.5, 0.0, 0.0), Chrom("Z", 1.5, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.5, 0.0, 0.0), Chrom("Z", 1.5, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0))  …  Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0))])

pop5 has 5% driving W (the stronger one) and 5% drive-suppressing autosome. For simplicity in setting this up, 5% of both females and males will start with 2 copies of drive-suppressing autosome, and otherwise boring chromosomes.

pop5setPropDrivingWstrong = 0.05
pop5setPropSuppressorA = 0.05
pop5SizeFemale = 10_000
pop5SizeMale = 10_000

pop5numFemaleTwoDriveSuppressorA = Int(round(pop5setPropSuppressorA * pop5SizeFemale))
pop5numFemaleDriveWstrong = Int(round(pop5setPropDrivingWstrong * pop5SizeFemale))
pop5numFemaleBoring = pop5SizeFemale - pop5numFemaleTwoDriveSuppressorA - pop5numFemaleDriveWstrong

pop5Females = vcat(fill(femaleTwoDriveSuppressorA, pop5numFemaleTwoDriveSuppressorA), 
                fill(femaleDriveWstrong, pop5numFemaleDriveWstrong), 
                fill(femaleBoring, pop5numFemaleBoring))

pop5numMaleTwoDriveSuppressorA = Int(round(pop5setPropSuppressorA * pop5SizeMale))
pop5numMaleBoring = pop5SizeMale - pop5numMaleTwoDriveSuppressorA

pop5Males = vcat(fill(maleTwoDriveSuppressorA, pop5numMaleTwoDriveSuppressorA),
                fill(maleBoring, pop5numMaleBoring))
pop5 = Population(pop5Females, pop5Males)
Population(Indiv[Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.5), Chrom("1", 1.0, 0.0, 0.5)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.5), Chrom("1", 1.0, 0.0, 0.5)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.5), Chrom("1", 1.0, 0.0, 0.5)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.5), Chrom("1", 1.0, 0.0, 0.5)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.5), Chrom("1", 1.0, 0.0, 0.5)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.5), Chrom("1", 1.0, 0.0, 0.5)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.5), Chrom("1", 1.0, 0.0, 0.5)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.5), Chrom("1", 1.0, 0.0, 0.5)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.5), Chrom("1", 1.0, 0.0, 0.5)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.5), Chrom("1", 1.0, 0.0, 0.5))  …  Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0))], Indiv[Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.5), Chrom("1", 1.0, 0.0, 0.5)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.5), Chrom("1", 1.0, 0.0, 0.5)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.5), Chrom("1", 1.0, 0.0, 0.5)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.5), Chrom("1", 1.0, 0.0, 0.5)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.5), Chrom("1", 1.0, 0.0, 0.5)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.5), Chrom("1", 1.0, 0.0, 0.5)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.5), Chrom("1", 1.0, 0.0, 0.5)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.5), Chrom("1", 1.0, 0.0, 0.5)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.5), Chrom("1", 1.0, 0.0, 0.5)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.5), Chrom("1", 1.0, 0.0, 0.5))  …  Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0))])

pop6 has 5% driving W (the stronger one) and 5% drive-cancelling autosome. For simplicity in setting this up, 5% of both females and males will start with 2 copies of drive-cancelling autosome, and otherwise boring chromosomes.

pop6setPropDrivingWstrong = 0.05
pop6setPropCancellerA = 0.05
pop6SizeFemale = 10_000
pop6SizeMale = 10_000

pop6numFemaleTwoDriveCancellerA = Int(round(pop6setPropCancellerA * pop6SizeFemale))
pop6numFemaleDriveWstrong = Int(round(pop6setPropDrivingWstrong * pop6SizeFemale))
pop6numFemaleBoring = pop6SizeFemale - pop6numFemaleTwoDriveCancellerA - pop6numFemaleDriveWstrong

pop6Females = vcat(fill(femaleTwoDriveCancellerA, pop6numFemaleTwoDriveCancellerA), 
                fill(femaleDriveWstrong, pop6numFemaleDriveWstrong), 
                fill(femaleBoring, pop6numFemaleBoring))

pop6numMaleTwoDriveCancellerA = Int(round(pop6setPropCancellerA * pop6SizeMale))
pop6numMaleBoring = pop6SizeMale - pop6numMaleTwoDriveCancellerA

pop6Males = vcat(fill(maleTwoDriveCancellerA, pop6numMaleTwoDriveCancellerA),
                fill(maleBoring, pop6numMaleBoring))
pop6 = Population(pop6Females, pop6Males)
Population(Indiv[Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 1.0), Chrom("1", 1.0, 0.0, 1.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 1.0), Chrom("1", 1.0, 0.0, 1.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 1.0), Chrom("1", 1.0, 0.0, 1.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 1.0), Chrom("1", 1.0, 0.0, 1.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 1.0), Chrom("1", 1.0, 0.0, 1.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 1.0), Chrom("1", 1.0, 0.0, 1.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 1.0), Chrom("1", 1.0, 0.0, 1.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 1.0), Chrom("1", 1.0, 0.0, 1.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 1.0), Chrom("1", 1.0, 0.0, 1.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 1.0), Chrom("1", 1.0, 0.0, 1.0))  …  Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0))], Indiv[Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 1.0), Chrom("1", 1.0, 0.0, 1.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 1.0), Chrom("1", 1.0, 0.0, 1.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 1.0), Chrom("1", 1.0, 0.0, 1.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 1.0), Chrom("1", 1.0, 0.0, 1.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 1.0), Chrom("1", 1.0, 0.0, 1.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 1.0), Chrom("1", 1.0, 0.0, 1.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 1.0), Chrom("1", 1.0, 0.0, 1.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 1.0), Chrom("1", 1.0, 0.0, 1.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 1.0), Chrom("1", 1.0, 0.0, 1.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 1.0), Chrom("1", 1.0, 0.0, 1.0))  …  Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0))])

Make a function for breeding for one generation

We will assume random mating, constant population size, and non-overlapping generations.

function breedOneGeneration(population::Population; k = length(population.popFemales) + length(population.popMales))
    # can set the number of offspring (k) in function call, or leave out and it will set k according to the current population size
    # cycle through offspring, choosing random parents
    offspringFemales = []
    offspringMales = []
    for i in 1:k
        mother = sample(population.popFemales)
        father = sample(population.popMales)
        offspring = makeOneOffspring(mother, father)
        if offspring.sexChromFromEgg.kind == "W" # so offspring is female
            push!(offspringFemales, offspring)
        elseif offspring.sexChromFromEgg.kind == "Z" # so offspring is male
            push!(offspringMales, offspring)
        end
    end
    return Population(offspringFemales, offspringMales)
end
breedOneGeneration (generic function with 2 methods)

Test the above function with one generation of breeding:

popNext = breedOneGeneration(pop1)
println("The offspring generation has ", string(length(popNext.popFemales)), " females and ", string(length(popNext.popMales)), " males.")
The offspring generation has 10005 females and 9995 males.

Create some functions for tracking population proportions

Function for calculating population proportion of driving W
function getPropDrivingW(population::Population)
    driving = []
    for i in population.popFemales # all females have a single W, making next line OK
        if i.sexChromFromEgg.spindleAttraction > 1.0
            push!(driving, true)
        else push!(driving, false)
        end
    end
    return sum(driving) / length(driving)
end
getPropDrivingW (generic function with 2 methods)

Test the function above:

getPropDrivingW(pop1)
0.05
Function for calculating population proportion of driving Z
function getPropDrivingZ(population::Population)
    driving = []
    # have to cycle through females and males to count up all the driving Z
    for i in population.popFemales
        if i.sexChromFromSperm.spindleAttraction > 1.0
            push!(driving, true)
        else push!(driving, false)
        end
    end
    for i in population.popMales
        if i.sexChromFromEgg.spindleAttraction > 1.0
            push!(driving, true)
        else push!(driving, false)
        end
        if i.sexChromFromSperm.spindleAttraction > 1.0
            push!(driving, true)
        else push!(driving, false)
        end     
    end
    return sum(driving) / length(driving)       
end
getPropDrivingZ (generic function with 2 methods)

Test the function above:

getPropDrivingZ(pop2)
0.05
Function for calculating population proportion of driving autosome
function getPropDrivingA(population::Population)
    driving = []
    # cycle through females and males to count up all the driving autosomes
    for i in population.popFemales
        if i.autosome1FromEgg.spindleAttraction > 1.0
            push!(driving, true)
        else push!(driving, false)
        end
        if i.autosome1FromSperm.spindleAttraction > 1.0
            push!(driving, true)
        else push!(driving, false)
        end     
    end
    for i in population.popMales
        if i.autosome1FromEgg.spindleAttraction > 1.0
            push!(driving, true)
        else push!(driving, false)
        end
        if i.autosome1FromSperm.spindleAttraction > 1.0
            push!(driving, true)
        else push!(driving, false)
        end     
    end
    return sum(driving) / length(driving)
end
getPropDrivingA (generic function with 1 method)

Test the function above:

getPropDrivingA(pop3)
0.05
Function for calculating population proportion of suppressing autosome
function getPropSuppressingA(population::Population)
    suppressing = []
    # have to cycle through females and males to count up all the autosomes of this type
    for i in population.popFemales
        if i.autosome1FromEgg.suppressStrength > 0.0
            push!(suppressing, true)
        else push!(suppressing, false)
        end
        if i.autosome1FromSperm.suppressStrength > 0.0
            push!(suppressing, true)
        else push!(suppressing, false)
        end     
    end
    for i in population.popMales
        if i.autosome1FromEgg.suppressStrength > 0.0
            push!(suppressing, true)
        else push!(suppressing, false)
        end
        if i.autosome1FromSperm.suppressStrength > 0.0
            push!(suppressing, true)
        else push!(suppressing, false)
        end     
    end
    return sum(suppressing) / length(suppressing)
end
getPropSuppressingA (generic function with 1 method)

Function for survival to adulthood

This is not discussed in the paper, but I’ve built possible survival costs into the model for expansion in the future. There are two types:

  • First, that specified by the fitnessLoss field of the Chrom struct.

  • Also, in the below we can specify a fitness cost (in survival to adulthood) for high spindle attraction during mitosis. This takes the total spindle attraction above 1 and multiplies it by a parameter driveCost to get the fitness cost. Here’s the function to calculate this spindle attraction cost:

function getSpindleAttractionFitnessCost(indiv::Indiv; driveCost = 0.0)
    excessSpindleAttraction = (indiv.sexChromFromEgg.spindleAttraction - 1) +
                            (indiv.sexChromFromSperm.spindleAttraction - 1) +
                            (indiv.autosome1FromEgg.spindleAttraction - 1) +
                            (indiv.autosome1FromSperm.spindleAttraction - 1)
    spindleAttractionFitnessCost = driveCost * excessSpindleAttraction
    if spindleAttractionFitnessCost > 1.0
        spindleAttractionFitnessCost = 1.0
    end
    return spindleAttractionFitnessCost
end
getSpindleAttractionFitnessCost (generic function with 1 method)

Here’s the complete function for deterimining the population of individuals that survived to adulthood:

function surviveToAdulthood(population::Population; driveCost = 0.0)
    femalePopSize = length(population.popFemales) 
    malePopSize = length(population.popMales)
    # cycle through individuals, determining survival
    survivingFemales = []
    survivingMales = []
    for i in population.popFemales  
        probSurvival = (1.0 - i.sexChromFromEgg.fitnessLoss) * (1.0 - i.sexChromFromSperm.fitnessLoss) * (1.0 - getSpindleAttractionFitnessCost(i; driveCost = driveCost))
        if probSurvival > rand()  # the individual survives
            push!(survivingFemales, i)
        end
    end
    for i in population.popMales
        probSurvival = (1.0 - i.sexChromFromEgg.fitnessLoss) * (1.0 - i.sexChromFromSperm.fitnessLoss) * (1.0 - getSpindleAttractionFitnessCost(i; driveCost = driveCost))
        if probSurvival > rand()  # the individual survives
            push!(survivingMales, i)        
        end
    end
    return Population(survivingFemales, survivingMales)
end
surviveToAdulthood (generic function with 2 methods)

In the paper, there are no survival fitness costs. (So all fitnessLoss fields of chromosomes and driveCost in the above function were set to 0.0)

Function for running for many generations

Keep track of proportion of females, driving W, driving Z, driving A, and suppressors.

function runManyGens(population::Population, gens::Int; k = length(population.popFemales) + length(population.popMales), driveCost = 0)
    # The above sets k (the number of offspring in each generation) as the starting population size
    propFemale = fill(-0.9999, gens+1)  # start with obviously wrong values
    propDrivingZ = fill(-0.9999, gens+1)
    propDrivingW = fill(-0.9999, gens+1)
    propDrivingA = fill(-0.9999, gens+1)
    propSuppressingA = fill(-0.9999, gens+1)
    # record proportions in starting population (generation 0; index 1)
    propFemale[1] = length(population.popFemales) / 
                        (length(population.popFemales) + length(population.popMales))
    propDrivingW[1] = getPropDrivingW(population)
    propDrivingZ[1] = getPropDrivingZ(population)
    propDrivingA[1] = getPropDrivingA(population)
    propSuppressingA[1] = getPropSuppressingA(population)
    for i in 1:gens
        # produce offspring
        offspringPopulation = breedOneGeneration(population; k)
        # determine survival to adulthood
        population = surviveToAdulthood(offspringPopulation; driveCost = driveCost)
        # record proportions in the offspring adult generation
        propFemale[i+1] = length(population.popFemales) / 
                        (length(population.popFemales) + length(population.popMales))
        propDrivingW[i+1] = getPropDrivingW(population)
        propDrivingZ[i+1] = getPropDrivingZ(population)
        propDrivingA[i+1] = getPropDrivingA(population)
        propSuppressingA[i+1] = getPropSuppressingA(population)
    end
    return (propFemale = propFemale, 
            propDrivingW = propDrivingW, 
            propDrivingZ = propDrivingZ,
            propDrivingA = propDrivingA, 
            propSuppressingA = propSuppressingA, 
            population = population)
end
runManyGens (generic function with 2 methods)

Run many generations, and plot proportions over time of:

  • females
  • driving W
  • driving Z
  • suppressing Autosome 1
runLength = 100
pop1run = runManyGens(pop1, runLength, driveCost = 0.0)
pop1plot = plot(0:runLength, pop1run.propFemale, color = :purple, label = "Females", xlabel = "Generations", ylabel = "Proportion", title = "Driving W, Drive Strength = " * string(driveStrength))
xlims!(0, 100)
ylims!(0, 1.05)
plot!(0:runLength, pop1run.propDrivingW, color = :red, label = "Driving W")
pop1Brun = runManyGens(pop1B, 200, driveCost = 0.0)
pop1Bplot = plot(0:200, pop1Brun.propFemale, color = :purple, label = "Females", xlabel = "Generations", ylabel = "Proportion", title = "Driving W, Drive Strength = " * string(driveStrengthWeak))
xlims!(0, 200)
ylims!(0, 1.05)
plot!(0:200, pop1Brun.propDrivingW, color = :red, label = "Driving W")
pop2run = runManyGens(pop2, runLength, driveCost = 0.0)
pop2plot = plot(0:runLength, pop2run.propFemale, color = :purple, label = "Female", xlabel = "Generations", ylabel = "Proportion", title = "Driving Z, Drive Strength = " * string(driveStrength))
xlims!(0,100)
ylims!(0, 1.05)
plot!(0:runLength, pop2run.propDrivingZ, color = :blue, label = "Driving Z")
pop3run = runManyGens(pop3, runLength, driveCost = 0.0)
pop3plot = plot(0:runLength, pop3run.propFemale, color = :purple, label = "Female", xlabel = "Generations", ylabel = "Proportion", title = "Driving Autosome, Drive Strength = " * string(driveStrength))
xlims!(0, 100)
ylims!(0, 1.05)
plot!(0:runLength, pop3run.propDrivingA, color = :green, label = "Driving autosome")
pop4run = runManyGens(pop4, runLength, driveCost = 0.0)
pop4plot = plot(0:runLength, pop4run.propFemale, color = :purple, label = "Female", xlabel = "Generations", ylabel = "Proportion", title = "Driving W and Z, Drive Strength = " * string(driveStrength))
xlims!(0, 100)
ylims!(0, 1.05)
plot!(0:runLength, pop4run.propDrivingW, color = :red, label = "Driving W")
plot!(0:runLength, pop4run.propDrivingZ, color = :blue, label = "Driving Z")

Combine some of above for Fig. 2 of paper

Fig2A = plot(0:runLength, pop1run.propFemale, color = :purple, linewidth = 2, label = "Females", xlabel = "Generations", ylabel = "Proportion", title = "A. Driving W", titlefontsize = 14, titlelocation = :left)
xlims!(0, 100)
ylims!(0, 1.05)
plot!(0:runLength, pop1run.propDrivingW, color = :red, linewidth = 4, label = "Driving W")
annotate!(35, 0.85, ("Driving W", :left, 8, :red))
annotate!(80, 0.5, ("females", :left, 8, :purple))

Fig2B = plot(0:runLength, pop2run.propFemale, color = :purple, linewidth = 2, label = "Female", xlabel = "Generations", ylabel = "Proportion", title = "B. Driving Z", titlefontsize = 14, titlelocation = :left)
xlims!(0, 100)
ylims!(0, 1.05)
plot!(0:runLength, pop2run.propDrivingZ, color = :blue, linewidth = 4, label = "Driving Z")
annotate!(70, 0.67, ("Driving Z", :left, 8, :blue))
annotate!(80, 0.3, ("females", :left, 8, :purple))

Fig2C = plot(0:runLength, pop3run.propFemale, color = :purple, linewidth = 2, label = "Female", xlabel = "Generations", ylabel = "Proportion", title = "C. Driving autosome")
xlims!(0, 100)
ylims!(0, 1.05)
plot!(0:runLength, pop3run.propDrivingA, color = :green, linewidth = 4, label = "Driving autosome", titlefontsize = 14, titlelocation = :left)
annotate!(55, 0.82, ("Driving autosome", :left, 8, :green))
annotate!(80, 0.4, ("females", :left, 8, :purple))

Figure2 = plot(Fig2A, Fig2B, Fig2C, layout=(3,1), legend=false, size = (400, 600))

To save the figure, run the following (after changing the false to true):

if false  # set to true to save plot
    savefig(Figure2, "Figure2.svg")
end 

Run some simulations with suppressors

First, with a 50% autosomal suppressor (when homozygous, half that when heterozygous) of any drivers:

pop5run = runManyGens(pop5, 400, driveCost = 0.0)
pop5plot = plot(0:400, pop5run.propFemale, color = :purple, label = "Female", xlabel = "Generations", ylabel = "Proportion", title = "Driving W and 50% Suppressor, Drive Strength = " * string(driveStrengthStrong))
xlims!(0, 400)
ylims!(0, 1.05)
plot!(0:400, pop5run.propDrivingW, color = :red, label = "Driving W")
plot!(0:400, pop5run.propSuppressingA, color = :grey, label = "Autosomal drive suppressor")

Now with a full supressor when homozygous (i.e., a “canceller”):

pop6run = runManyGens(pop6, 400, driveCost = 0.0)
pop6plot = plot(0:400, pop6run.propFemale, color = :purple, label = "Female", xlabel = "Generations", ylabel = "Proportion", title = "Driving W and Canceller, Drive Strength = " * string(driveStrengthStrong))
xlims!(0, 400)
ylims!(0, 1)
plot!(0:400, pop6run.propDrivingW, color = :red, label = "Driving W")
plot!(0:400, pop6run.propSuppressingA, color = :grey, label = "Autosomal drive canceller")

Build a deterministic version of the model

The above simulations were individual-based, with finite population size, hence had stochasticity (i.e., slightly different results each time, although running with large population size—20,000—makes that effect small.)

Below, we can develop a deterministic version that has no randomness. We do that by assuming infinite population sizes (not realistic of course, but handy if interested in a non-stochastic outcome). There are also efficiency benefits in terms of computation (i.e., faster and needing less memory).

Rather than have individuals in our simulations, we will keep track of proportions of individuals (of infinite population) that are of these four female types:

  • femaleBoring
  • femaleDriveW
  • femaleDriveZ
  • femaleDriveWZ

And these four males types:

  • maleBoring
  • maleFirstDriveZ
  • maleSecondDriveZ
  • maleTwoDriveZ

We will allow sex ratio to diverge from 50:50, by having all of the above add up to frequency of 1.

Using this approach, we will not (yet) incorporate suppressors of drive.

Will make a new data type called PopInfinite, to store data describing individual type frequencies in an infinite population.

mutable struct PopInfinite
    indivKinds # a vector of `Indiv` kinds
    freqs # a vector of frequencies of those `indKinds` in the population
end

Write a function for setting up a PopInfinite population, assuming HW equilibrium:

function startPopInfinite(startFreqfemale, startFreqDriveW, startFreqDriveZ)
    
    indivKindsStart = [femaleBoring,
            femaleDriveW,
            femaleDriveZ,
            femaleDriveWZ,
            maleBoring,
            maleFirstDriveZ,
            maleSecondDriveZ,
            maleTwoDriveZ]

    startKindFreqs = [startFreqfemale * (1-startFreqDriveW) * (1-startFreqDriveZ),
                startFreqfemale * startFreqDriveW * (1-startFreqDriveZ),
                startFreqfemale * (1-startFreqDriveW) * startFreqDriveZ,
                startFreqfemale * startFreqDriveW * startFreqDriveZ,
                (1-startFreqfemale) * (1-startFreqDriveZ) * (1-startFreqDriveZ),
                (1-startFreqfemale) * startFreqDriveZ * (1-startFreqDriveZ),
                (1-startFreqfemale) * (1-startFreqDriveZ) * startFreqDriveZ,
                (1-startFreqfemale) * startFreqDriveZ * startFreqDriveZ]

    return PopInfinite(indivKindsStart, startKindFreqs)
end
startPopInfinite (generic function with 1 method)

Set up some infinite populations (to correspond to some of the finite ones above):

popInf1 = startPopInfinite(0.5, 0.05, 0)  # driving W
popInf2 = startPopInfinite(0.5, 0, 0.05)  # driving Z
popInf4 = startPopInfinite(0.5, 0.05, 0.05)  # driving W & Z
popInf5 = startPopInfinite(0.5, 0.05, 1.0)  # driving W & already fixed driving Z
PopInfinite(Indiv[Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.5, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.0, 0.0, 0.0), Chrom("Z", 1.5, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("W", 1.5, 0.0, 0.0), Chrom("Z", 1.5, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.5, 0.0, 0.0), Chrom("Z", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.0, 0.0, 0.0), Chrom("Z", 1.5, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0)), Indiv(Chrom("Z", 1.5, 0.0, 0.0), Chrom("Z", 1.5, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0), Chrom("1", 1.0, 0.0, 0.0))], [0.0, 0.0, 0.475, 0.025, 0.0, 0.0, 0.0, 0.5])

Function for breeding one generation using a PopInfinite

This adds a method to the existing function breedOneGeneration()

Will build this functions below using the idea of producing gametes in their proportions, and then combining them randomly. This will get the probability of W transmission, based on spindleAttraction of W and Z:

function produceEggProbW(femaleIndiv::Indiv)
transProbW = femaleIndiv.sexChromFromEgg.spindleAttraction /
            (femaleIndiv.sexChromFromEgg.spindleAttraction + femaleIndiv.sexChromFromSperm.spindleAttraction)
end
produceEggProbW (generic function with 1 method)

Now make a function to carry out one generation of breeding:

function breedOneGeneration(population::PopInfinite, chromKinds)
    # generate egg and sperm genotype frequencies
    eggFreqs = zeros(length(chromKinds))  # these are frequencies as proportions of all gametes
    spermFreqs = zeros(length(chromKinds))  # these are frequencies as proportions of all gametes
    for i in 1:length(population.freqs)  # cycle through kinds of individuals
        if population.indivKinds[i].sexChromFromEgg.kind == "W"  # female individual
            eggProbW = produceEggProbW(population.indivKinds[i]) # get prob of W egg
            chromKindRow = findfirst(isequal(population.indivKinds[i].sexChromFromEgg), chromKinds) # get row for recording the egg type
            eggFreqs[chromKindRow] = eggFreqs[chromKindRow] + (population.freqs[i] * eggProbW) # add prob to proper row of egg frequency table 
            # now consider prob of Z egg
            eggProbZ = 1 - eggProbW
            chromKindRow = findfirst(isequal(population.indivKinds[i].sexChromFromSperm), chromKinds) # get row for recording the egg type, as carrying mother's chromosome that came from sperm (in her parent)
            eggFreqs[chromKindRow] = eggFreqs[chromKindRow] + (population.freqs[i] * eggProbZ)
        
        elseif population.indivKinds[i].sexChromFromEgg.kind == "Z"  # male individual, so no centromere drive
            # first consider transferring chromosome that came from individual's mother
            chromKindRow = findfirst(isequal(population.indivKinds[i].sexChromFromEgg), chromKinds) # get row for recording the sperm type
            spermFreqs[chromKindRow] = spermFreqs[chromKindRow] + (population.freqs[i] * 0.5)
            # now consider transferring chromosome that came from individual's father
            chromKindRow = findfirst(isequal(population.indivKinds[i].sexChromFromSperm), chromKinds) # get row for recording sperm type
            spermFreqs[chromKindRow] = spermFreqs[chromKindRow] + (population.freqs[i] * 0.5)
        end
    end

    # calculate egg freqs as a proportion of eggs
    eggProps = eggFreqs ./ sum(eggFreqs)

    # calculate sperm freqs as a proportion of sperm
    spermProps = spermFreqs ./ sum(spermFreqs)
    
    # now combine eggs and sperm to make new population
    newIndivKindFreqs = zeros(length(population.freqs))
    for i in 1:length(eggProps)  # loop through egg types
        for j in 1:length(spermProps)  # loop through sperm types
            # find row for the resulting kind of individual when combining these gametes
            focalIndivKind = Indiv(chromKinds[i], chromKinds[j], boringAutosome1, boringAutosome1)
            indivKindRow = findfirst(isequal(focalIndivKind), population.indivKinds)
            if eggProps[i] > 0.0 && spermProps[j] > 0.0
                # add frequency of this egg/sperm combo to table of new individual freqs
                newIndivKindFreqs[indivKindRow] = newIndivKindFreqs[indivKindRow] + (eggProps[i] * spermProps[j])
                #println(focalIndivKind)
                #println(indivKindRow)
            end
        end
    end

    return PopInfinite(population.indivKinds, newIndivKindFreqs)
end
breedOneGeneration (generic function with 2 methods)

Add method to surviveToAdulthood() function, for PopInfinite object

function surviveToAdulthood(population::PopInfinite)
    probSurvivals = fill(-9.0, length(population.indivKinds))
    for i in 1:length(population.indivKinds) # cycle through kinds of individuals
        indivKind = population.indivKinds[i]
        probSurvivals[i] = (1.0 - indivKind.sexChromFromEgg.fitnessLoss) * (1.0 - indivKind.sexChromFromSperm.fitnessLoss)
    end
    return PopInfinite(population.indivKinds, (population.freqs .* probSurvivals))
end
surviveToAdulthood (generic function with 2 methods)

(Note that in all the simulations presented in the paper, there is 100% survival to adulthood—because fitnessLoss is zero. Hence prior to using fitnessLoss of other values, I recommend careful testing.)

Add methods to functions for calculating proportions

The below with add methods to existing functions, for PopInfinite objects:

function getPropDrivingW(population::PopInfinite)
    freqDriveW = 0.0
    freqBoringW = 0.0
    for i in 1:length(population.freqs)
        if population.indivKinds[i].sexChromFromEgg.kind == "W"  # females
            if population.indivKinds[i].sexChromFromEgg.spindleAttraction > 1.0
                freqDriveW = freqDriveW + population.freqs[i]
            else 
                freqBoringW = freqBoringW + population.freqs[i]
            end
        end
    end
    return freqDriveW / (freqDriveW + freqBoringW)
end
getPropDrivingW (generic function with 2 methods)
function getPropDrivingZ(population::PopInfinite)
    freqDriveZ = 0.0
    freqBoringZ = 0.0
    for i in 1:length(population.freqs)
        if population.indivKinds[i].sexChromFromEgg.kind == "W"  # females
            if population.indivKinds[i].sexChromFromSperm.spindleAttraction > 1.0
                freqDriveZ = freqDriveZ + population.freqs[i]
            else 
                freqBoringZ = freqBoringZ + population.freqs[i]
            end
        end

        if population.indivKinds[i].sexChromFromEgg.kind == "Z"  # males
            # check first Z (from egg)
            if population.indivKinds[i].sexChromFromEgg.spindleAttraction > 1.0
                freqDriveZ = freqDriveZ + population.freqs[i]
            else 
                freqBoringZ = freqBoringZ + population.freqs[i]
            end
            # check second Z (from sperm)
            if population.indivKinds[i].sexChromFromSperm.spindleAttraction > 1.0
                freqDriveZ = freqDriveZ + population.freqs[i]
            else 
                freqBoringZ = freqBoringZ + population.freqs[i]
            end
        end
    end
    return freqDriveZ / (freqDriveZ + freqBoringZ)
end
getPropDrivingZ (generic function with 2 methods)
function getPropFemale(population::PopInfinite)
    freqFemale = 0.0
    for i in 1:length(population.freqs)
        if population.indivKinds[i].sexChromFromEgg.kind == "W"  # females
            freqFemale = freqFemale + population.freqs[i]
        end
    end
    return freqFemale
end
getPropFemale (generic function with 1 method)

Add method to runManyGens() function, for PopInfinite object

function runManyGens(population::PopInfinite, chromKinds, gens::Int)
    propFemale = fill(-0.9999, gens+1)  # start with obviously wrong values
    propDrivingZ = fill(-0.9999, gens+1)
    propDrivingW = fill(-0.9999, gens+1)
    # record proportions in starting population (generation 0; index 1)
    propFemale[1] = getPropFemale(population)
    propDrivingW[1] = getPropDrivingW(population)
    propDrivingZ[1] = getPropDrivingZ(population)
    for i in 1:gens
        # breed and survive to adulthood
        populationOffspring = breedOneGeneration(population, chromKinds)
        population = surviveToAdulthood(populationOffspring)
        # record proportions in the offspring adult generation
        propFemale[i+1] = getPropFemale(population)
        propDrivingW[i+1] = getPropDrivingW(population)
        propDrivingZ[i+1] = getPropDrivingZ(population)
    end
    return (propFemale = propFemale, 
            propDrivingW = propDrivingW, 
            propDrivingZ = propDrivingZ,
            population = population)
end
runManyGens (generic function with 2 methods)

Run some deterministic simulations

popInf1Run = runManyGens(popInf1, chromKinds, runLength)
popInf1plot = plot(0:runLength, popInf1Run.propFemale, color = :purple, label = "Female", xlabel = "Generations", ylabel = "Proportion", title = "Driving W, Drive Strength = " * string(driveStrength))
xlims!(0,100)
ylims!(0, 1)
plot!(0:runLength, popInf1Run.propDrivingW, color = :red, label = "Driving W")
popInf2Run = runManyGens(popInf2, chromKinds, runLength)
popInf2plot = plot(0:runLength, popInf2Run.propFemale, color = :purple, label = "Female", xlabel = "Generations", ylabel = "Proportion", title = "Driving Z, Drive Strength = " * string(driveStrength))
xlims!(0,100)
ylims!(0, 1)
plot!(0:runLength, popInf2Run.propDrivingZ, color = :blue, label = "Driving Z")
popInf4Run = runManyGens(popInf4, chromKinds, runLength)
popInf4plot = plot(0:runLength, popInf4Run.propFemale, color = :purple, label = "Female", xlabel = "Generations", ylabel = "Proportion", title = "Driving W & Z, Drive Strength = " * string(driveStrength))
xlims!(0,100)
ylims!(0, 1)
plot!(0:runLength, popInf4Run.propDrivingW, color = :red, label = "Driving W")
plot!(0:runLength, popInf4Run.propDrivingZ, color = :blue, label = "Driving Z")
popInf5Run = runManyGens(popInf5, chromKinds, runLength)
popInf5plot = plot(0:runLength, popInf5Run.propFemale, color = :purple, label = "Female", xlabel = "Generations", ylabel = "Proportion", title = "Driving W and fixed driving Z, Drive Strength = " * string(driveStrength))
xlims!(0,100)
ylims!(0, 1)
plot!(0:runLength, popInf5Run.propDrivingW, color = :red, label = "Driving W")
plot!(0:runLength, popInf5Run.propDrivingZ, color = :blue, label = "Driving Z")

Put some of the above together as Fig. 3 of the paper

Fig3A = plot(0:runLength, popInf1Run.propFemale, color = :purple, linewidth = 2, label = "Female", xlabel = "Generations", ylabel = "Proportion", title = "A. Driving W", titlefontsize = 13, titlelocation = :left)
xlims!(0,100)
ylims!(0, 1.05)
plot!(0:runLength, popInf1Run.propDrivingW, color = :red, linewidth = 4, linealpha = 0.75, label = "Driving W")

Fig3B = plot(0:runLength, popInf2Run.propFemale, color = :purple, linewidth = 2, label = "Female", xlabel = "Generations", ylabel = "Proportion", title = "B. Driving Z", titlefontsize = 13, titlelocation = :left)
xlims!(0,100)
ylims!(0, 1.05)
plot!(0:runLength, popInf2Run.propDrivingZ, color = :blue, linewidth = 4, linealpha = 0.75, label = "Driving Z")

Fig3C = plot(0:runLength, popInf4Run.propFemale, color = :purple, linewidth = 2, label = "Female", xlabel = "Generations", ylabel = "Proportion", title = "C. Driving W & Z", titlefontsize = 13, titlelocation = :left)
xlims!(0,100)
ylims!(0, 1.05)
plot!(0:runLength, popInf4Run.propDrivingW, color = :red, linewidth = 4, linealpha = 0.75, label = "Driving W")
plot!(0:runLength, popInf4Run.propDrivingZ, color = :blue, linewidth = 4, linealpha = 0.75, label = "Driving Z")

Fig3D = plot(0:runLength, popInf5Run.propFemale, color = :purple, linewidth = 2, label = "Female", xlabel = "Generations", ylabel = "Proportion", title = "D. Driving W; fixed driving Z", titlefontsize = 13, titlelocation = :left)
xlims!(0,100)
ylims!(0, 1.05)
plot!(0:runLength, popInf5Run.propDrivingZ, color = :blue, linewidth = 4, linealpha = 0.75, label = "Driving Z")
plot!(0:runLength, popInf5Run.propDrivingW, color = :red, linewidth = 4, linealpha = 0.75, label = "Driving W")

Fig3E = plot(0:400, pop5run.propFemale, color = :purple, linewidth = 2, label = "Female", xlabel = "Generations", ylabel = "Proportion", title = "E. Driving W; partial suppressor", titlefontsize = 13, titlelocation = :left)
xlims!(0, 400)
ylims!(0, 1.05)
plot!(0:400, pop5run.propDrivingW, color = :red, linealpha = 0.75, linewidth = 4, label = "Driving W")
plot!(0:400, pop5run.propSuppressingA, color = :grey, linealpha = 0.75, linewidth = 4, label = "Autosomal drive suppressor")
annotate!(25, 0.92, ("Driving W", :left, 8, :red))
annotate!(84, 0.3, ("Autosomal drive suppressor (partial)", :left, 8, :grey))
annotate!(310, 0.58, ("females", :left, 8, :purple))

Fig3F = plot(pop6run.propFemale, color = :purple, linewidth = 2, label = "Female", xlabel = "Generations", ylabel = "Proportion", title = "F. Driving W; complete suppressor", titlefontsize = 13, titlelocation = :left)
xlims!(0, 400)
ylims!(0, 1.05)
plot!(0:400, pop6run.propDrivingW, color = :red, linealpha = 0.75, linewidth = 4, label = "Driving W")
plot!(0:400, pop6run.propSuppressingA, color = :grey, linealpha = 0.75, linewidth = 4, label = "Autosomal drive negator")
annotate!(25, 0.92, ("Driving W", :left, 8, :red))
annotate!(55, 0.3, ("Autosomal drive suppressor (complete)", :left, 8, :grey))
annotate!(310, 0.48, ("females", :left, 8, :purple))

l = @layout [[a{0.9w};b{0.9w};c{0.9w};d{0.9w}] [e;f]]
Figure3 = plot(Fig3A, Fig3B, Fig3C, Fig3D, Fig3E, Fig3F, layout=l, legend=false, size = (800, 800))

To save the figure, run the following (after changing the false to true):

if false  # set to true to save plot
    savefig(Figure3, "Figure3.svg")
end 

This is a Quarto website.

To learn more about Quarto websites visit https://quarto.org/docs/websites.