Supplementary Mathematica code for simulations and analytical plots of the meiotic-constraint model
Data files
Jul 20, 2026 version files 27.87 KB
-
README.md
12.30 KB
-
Supplementary_Code_S1_Mathematica_simulation.wl
7.55 KB
-
Supplementary_Code_S1_README.txt
1.15 KB
-
Supplementary_Code_S2_Mathematica_Fig4b_plot.wl
5.80 KB
-
Supplementary_Code_S2_README.txt
1.07 KB
Abstract
Most multicellular animals practice anisogamy (fertilization between eggs and sperm). When mothers produce sons and daughters at a 1:1 ratio, the “twofold cost of males” arises because males do not directly contribute to population growth. If thelytokous parthenogens producing only daughters invade a population, they should spread rapidly. Although thelytoky has repeatedly evolved across invertebrate and vertebrate taxa, it remains a minority. Why? The evolutionary transition from anisogamy to thelytoky requires eggs to initiate embryonic development without fertilization. However, in metazoan animals, meiotic metaphase (MM) arrest halts oogenesis midway and normally resumes only after stimulation by sperm penetration. Empirical and experimental evidence indicates that release of MM arrest without fertilization is extremely difficult, providing a strong mechanistic barrier against parthenogenesis. Even if MM arrest were released, oogenesis would proceed to produce either a haploid embryo or a diploid embryo through refusion with the second polar body (terminal fusion automixis). Outbred species typically accumulate more than one lethal equivalent of recessive deleterious alleles per genome as heterozygotes. Upon transition to haploid or automictic development, these recessive lethals normally masked in outbred diploids would be exposed simultaneously, causing embryonic death and creating the next barrier. Thus, thelytoky cannot be achieved simply by modification of the existing meiotic system; instead, other mechanisms, such as apomixis, that bypass meiosis are required. Mathematical models and simulations support this “meiotic constraint” hypothesis. Combined with recently proposed immediate benefits of anisogamy and traditional genetic benefits (e.g., Red Queen), it may largely explain the maintenance of costly anisogamy.
Dataset DOI: 10.5061/dryad.sbcc2frpg
Description of the data and file structure
Dataset title:
Supplementary Mathematica code for "The twofold cost of sex reconsidered: Meiotic mechanisms protect anisogamous populations from invasion by thelytoky"
Description:
This dataset contains Mathematica/Wolfram Language scripts and README files used to reproduce the analytical plot and simulation results shown in Fig. 4b-d of the manuscript. No new empirical raw data were generated in this study. The empirical parameter values used in Fig. 4 were taken from the published sources cited in the manuscript and are defined directly in the deposited code.
Files:
1. Supplementary_Code_S1_Mathematica_simulation.wl
Count-only Mathematica simulation corresponding to Fig. 4c,d. The code compares haploid development, terminal-fusion automixis, and anisogamy outcrossing.
2. Supplementary_Code_S1_README.txt
README file explaining the purpose, assumptions, parameters, and output of Supplementary Code S1.
3. Supplementary_Code_S2_Mathematica_Fig4b_plot.wl
Mathematica code for the analytical phase-diagram style plot corresponding to Fig. 4b.
4. Supplementary_Code_S2_README.txt
README file explaining the purpose, analytical boundaries, parameters, and output of Supplementary Code S2.
Software:
The scripts were written for Mathematica / Wolfram Language 14.1.
Reproducibility:
Supplementary Code S1 uses a fixed random seed, SeedRandom[12345], so the same stochastic output can be reproduced in the same Mathematica/Wolfram Language environment. Supplementary Code S2 is deterministic and does not use random sampling.
How to use:
Open each .wl file in Mathematica and evaluate the full script. For Supplementary Code S1, set parameterSet = "Fig4c" or parameterSet = "Fig4d" to reproduce the corresponding simulation panel. Supplementary Code S2 directly produces the analytical background plot for Fig. 4b.
License:
The files are original code and documentation produced by the author and are compatible with the CC0 license waiver required by Dryad.
Files and variables
Files
1. Supplementary_Code_S1_Mathematica_simulation.wl
This file contains Mathematica/Wolfram Language code for the count-only simulation corresponding to Fig. 4c,d of the manuscript. The code compares three reproductive/developmental routes: haploid development, terminal-fusion automixis, and anisogamy outcrossing. It generates survival/death counts and a paired histogram for either Fig. 4c or Fig. 4d, depending on the selected parameter set.
Main variables and abbreviations used in this file:
- H: haploid development.
- TF: terminal-fusion automixis.
- AN: anisogamy outcrossing.
- L: number of functional loci. In the code, L is represented by nTot.
- numchromosome: number of chromosomes used in the simplified count-only model.
- numlocus: number of loci per chromosome.
- nTot: total number of functional loci, calculated as numchromosome*numlocus.
- p: standing per-locus probability that a haploid genome carries a recessive lethal allele.
- r: recombination fraction. In the code, this is represented by recombination.
- recombination: recombination fraction used to divide the genome into recombined and unrecombined regions.
- qSurv: per-gamete probability of carrying allele 0 among gametes sampled from surviving diploid parents. Defined as qSurv = p/(1+p).
- allele 0: recessive lethal allele.
- allele 1: normal allele.
- iteration: number of simulation replicates.
- lethalThreshold: threshold number of expressed lethal alleles/loci at or above which an embryo is counted as dead. In this code, lethalThreshold = 1.
- tableH: sampled counts of exposed lethal alleles in haploid development.
- tableTF: sampled counts of lethal loci in the unrecombined region under terminal-fusion automixis.
- tableAN: sampled counts of homozygous lethal loci under anisogamy outcrossing.
- motherAlive, motherDead: diagnostic counts of viable/inviable diploid mothers.
- hAlive, hDead: survival/death counts for haploid development.
- tfAlive, tfDead: survival/death counts for terminal-fusion automixis.
- anAlive, anDead: survival/death counts for anisogamy outcrossing.
- summary: Mathematica Association containing the main parameter values and survival/death counts.
- histFig4cd: paired histogram corresponding to Fig. 4c or Fig. 4d.
Units:
All variables are dimensionless counts, probabilities, or proportions. L, numchromosome, numlocus, nTot, iteration, and alive/dead values are counts. p, r, recombination, qSurv, and survival probabilities are dimensionless probabilities.
Missing values:
No missing values are used in this file.
2. Supplementary_Code_S1_README.txt
This text file explains the purpose, assumptions, parameters, output, and reproducibility of Supplementary Code S1.
Main abbreviations defined in this file:
- H: haploid development.
- TF: terminal-fusion automixis.
- AN: anisogamy outcrossing.
- L: number of functional loci.
- p: standing per-locus probability that a haploid genome carries a recessive lethal allele.
- r: recombination fraction.
- qSurv: p/(1+p), the per-gamete probability of carrying allele 0 among gametes sampled from surviving diploid parents.
Units:
All variables are dimensionless counts, probabilities, or proportions.
Missing values:
No missing values are used in this file.
3. Supplementary_Code_S2_Mathematica_Fig4b_plot.wl
This file contains Mathematica/Wolfram Language code for the analytical phase-diagram style plot corresponding to Fig. 4b of the manuscript. The code plots analytical boundary curves separating parameter regions in which terminal-fusion automixis and anisogamy differ in viability. It also superimposes approximate empirical reference points for Drosophila, human, and mouse.
Main variables and abbreviations used in this file:
- L: number of functional loci. In the code, this is represented by loci or l.
- p: standing per-locus probability that a haploid genome carries a recessive lethal allele.
- r: recombination fraction. Fixed at r = 0.1.
- qSurv: per-gamete probability of carrying allele 0 among gametes sampled from surviving diploid parents. Defined as qSurv = p/(1+p).
- TF: terminal-fusion automixis.
- AN: anisogamy outcrossing.
- tfBoundary05: analytical boundary function for TF survival = 0.5.
- anBoundary001: analytical boundary function for AN survival = 0.01.
- lMin, lMax: minimum and maximum plotted values of L.
- pMin, pMax: minimum and maximum plotted values of p.
- humanL: approximate number of functional loci for human, represented by annotated protein-coding gene count.
- humanMhap: estimated number of recessive lethal mutations per haploid genome for human.
- humanP: humanMhap/humanL.
- flyL: approximate number of functional loci for Drosophila melanogaster.
- flyMhap: estimated number of recessive lethal mutations per haploid genome for Drosophila.
- flyP: flyMhap/flyL.
- mouseL: approximate number of functional loci for mouse.
- mouseMhapLow, mouseMhapHigh: lower and upper proxy values for the estimated number of recessive lethal mutations per haploid genome in mouse.
- mousePLow, mousePHigh: corresponding lower and upper values of p for mouse.
- mousePMid: geometric midpoint of mousePLow and mousePHigh.
- parameterSummary: Mathematica Association containing the empirical parameter values used in the plot.
- fig4bSquare: analytical phase-diagram style plot corresponding to Fig. 4b.
Units:
All variables are dimensionless counts, probabilities, or proportions. L, humanL, flyL, and mouseL are counts of loci. p, r, qSurv, humanP, flyP, mousePLow, mousePHigh, and mousePMid are dimensionless probabilities.
Missing values:
No missing values are used in this file.
4. Supplementary_Code_S2_README.txt
This text file explains the purpose, analytical boundaries, definitions, empirical reference values, output, and reproducibility of Supplementary Code S2.
Main abbreviations defined in this file:
- L: number of functional loci.
- p: standing per-locus probability that a haploid genome carries a recessive lethal allele.
- r: recombination fraction.
- qSurv: p/(1+p), the per-gamete probability of carrying allele 0 among gametes sampled from surviving diploid parents.
- TF: terminal-fusion automixis.
- AN: anisogamy outcrossing.
- Mhap: estimated number of recessive lethal mutations per haploid genome.
Units:
All variables are dimensionless counts, probabilities, or proportions.
Missing values:
No missing values are used in this file.
Compressed archives and folders
No compressed archives or folders are included in this dataset. All files are provided directly as individual files.
Code/software
The deposited .wl files are plain-text Mathematica/Wolfram Language scripts and can be viewed with any standard text editor. The README files are plain-text files and can also be viewed with any standard text editor.
The scripts were written and tested in Mathematica / Wolfram Language 14.1. No additional Mathematica packages were loaded or required.
Workflow:
1. Supplementary_Code_S1_Mathematica_simulation.wl
This script reproduces the count-only simulations corresponding to Fig. 4c,d. Open the script in Mathematica / Wolfram Language 14.1 and evaluate the full file. To reproduce Fig. 4c or Fig. 4d, set parameterSet = "Fig4c" or parameterSet = "Fig4d" near the beginning of the script before evaluation. The script prints the parameter values, survival/death counts, and survival probabilities, and returns a paired histogram.
2. Supplementary_Code_S2_Mathematica_Fig4b_plot.wl
This script reproduces the analytical phase-diagram style plot corresponding to Fig. 4b. Open the script in Mathematica / Wolfram Language 14.1 and evaluate the full file. The script prints the parameter summary and returns the analytical background plot with empirical reference points.
No open-source replacement for Mathematica / Wolfram Language was used or validated for executing these scripts. However, all code files are plain-text files and can be inspected without proprietary software.
Access information
No other public repository location is currently available for this dataset. The files deposited here are original Mathematica/Wolfram Language scripts and README documentation created by the author.
No third-party data files, copyrighted figures, copyrighted tables, screenshots, or article PDFs are redistributed in this Dryad dataset.
Some empirical parameter values used in Supplementary Code S2 were derived from published sources and public genome annotation pages cited in the manuscript. Specifically, approximate protein-coding gene counts for Drosophila melanogaster, Mus musculus, and Homo sapiens were taken from Ensembl annotation pages, and the estimated recessive-lethal burdens were taken from the published sources cited in the manuscript. These values are manually defined as scalar parameters in the code only; no external database file or copyrighted table is reproduced.
Sources used for empirical parameter values:
- Ensembl annotation page for Drosophila melanogaster protein-coding gene count.
- Ensembl annotation page for Mus musculus protein-coding gene count.
- Ensembl annotation page for Homo sapiens protein-coding gene count.
- Simmons & Crow (1977), cited in the manuscript, for the Drosophila recessive-lethal burden.
- Bittles & Neel (1994) and Wade et al. (2023), cited in the manuscript, for the human recessive-lethal burden.
- The mouse value is treated as a proxy range in the code, as described in the manuscript.
License information:
The deposited files themselves are original code and documentation produced by the author and are compatible with the CC0 license waiver required by Dryad. Ensembl integrated data are publicly accessible; only scalar summary values from the annotation pages are used in the code. Published literature sources are cited in the manuscript, but no copyrighted text, figures, or tables from those sources are included in this dataset.