In angiosperm self-incompatibility systems, pollen with an allele matching the pollen recipient at the self-incompatibility locus is rejected. Extreme allelic polymorphism is maintained by frequency-dependent selection favoring rare alleles. However, two challenges result in a "chicken-egg"problem for the spread of a new allele (a tightly linked haplotype in this case) under the widespread "collaborative non-self recognition" mechanism. A novel pollen-function mutation alone would merely grant compatibility with a nonexistent style-function allele: a neutral change at best. A novel pistil-function mutation alone could only be fertilized by pollen with a nonexistent pollen-function allele: a deleterious change that would eliminate all seed set. However, a pistil-function mutation complementary to a previously neutral pollen mutation may spread if it restores self-incompatibility to a self-compatible intermediate. We show that novel haplotypes can also drive elimination of existing ones with fewer siring opportunities. We calculate relative probabilities of increase and collapse in haplotype number given the initial collection of incompatibility haplotypes and the population gene conversion rate. Expansion in haplotype number is possible when population gene conversion rate is large, but large contractions are likely otherwise. A Markov chain model derived from these expansion and collapse probabilities generates a stable haplotype number distribution in the realistic range of 10--40 under plausible parameters. However, smaller populations might lose many haplotypes beyond those lost by chance during bottlenecks.

The file equilibria.txt, which contains the equilibrium frequencies of self-compatible intermediates, was created by numerically solving genotype frequency recursions in Mathematica using the uploaded notebook sc-balance.nb. All other data were generated by simulation or deterministic iteration using custom R scripts. Data were used either to calculate derived values or as starting points for further simulations, as described in the readme.

The transition matrices contain NA values, which are the result of a reflecting lower boundary: the state at the lower boundary is never reached, so transitions out from it were never calculated. See the README.txt for the workflow for generating the data.