Reference

HybridZones

A Julia implementation of the multilocus hybrid zone simulation and inference framework originally developed by Mallet and Barton in Turbo Pascal during the late 1980s, with subsequent extensions in PowerBASIC and R.

This package provides forward simulation of allele and genotype frequency clines under selection and migration, along with maximum-likelihood cline fitting and disequilibrium estimation. Validated against historical reference implementations including the original Pascal simulators and the R reimplementation in Rosser, Dasmahapatra & Mallet (2014).

References

• Mallet, J. and Barton, N. H. (1989). Inference from clines stabilized by frequency-dependent selection. Genetics 122: 967-976.
• Mallet, J., Barton, N., Lamas, G., Santisteban, J., Muedas, M. and Eeley, H. (1990). Estimates of selection and gene flow from measures of cline width and linkage disequilibrium in Heliconius hybrid zones. Genetics 124: 921-936.
• Rosser, N., Dasmahapatra, K. K. and Mallet, J. (2014). Stable Heliconius butterfly hybrid zones are correlated with a local rainfall peak at the edge of the Amazon basin. Evolution 68: 3470-3484.

package_version()

Return the version of the HybridZones package as a VersionNumber.

This function exists primarily so the package has a callable export from its first commit, allowing CI workflows to verify that the package loads and exports a function correctly. Real simulation functionality will be added in subsequent commits.

Examples

julia> using HybridZones

julia> v = package_version();

julia> v isa VersionNumber
true

GeneticArchitecture

Abstract supertype for genetic architecture specifications. Concrete subtypes carry locus count, ploidy, allele identities, and dominance relationships used by selection and migration models to dispatch the correct update rules.

OneLocusDiploid <: GeneticArchitecture

A single autosomal diploid locus with two alleles.

dominance specifies how alleles combine in heterozygotes:

• :codominant: both alleles contribute equally (additive)
• :dominant_first: the first allele is fully expressed; the second is masked
• :recessive_first: the first allele is masked; the second is fully expressed

allele_names provides human-readable labels for the two alleles, used in display and output. Defaults to ("A1", "A2").

Examples

julia> using HybridZones

julia> arch = OneLocusDiploid();

julia> arch isa GeneticArchitecture
true

julia> dominance(arch)
:codominant

julia> allele_names(arch)
("A1", "A2")

dominance(arch::OneLocusDiploid) -> Symbol

Return the dominance mode of arch: one of :codominant, :dominant_first, or :recessive_first.

n_loci(arch::OneLocusDiploid) -> Int

Return the number of loci. Always 1 for OneLocusDiploid.

n_alleles(arch::OneLocusDiploid) -> Int

Return the number of alleles per locus. Always 2 for OneLocusDiploid.

allele_names(arch::OneLocusDiploid) -> NTuple{2,String}

Return the two allele labels.

n_genotypes(arch::OneLocusDiploid) -> Int

Return the number of genotype classes. Always 3 for OneLocusDiploid (A1A1, A1A2, A2A2).

Examples

julia> using HybridZones

julia> arch = OneLocusDiploid();

julia> n_genotypes(arch)
3

SelectionModel

Abstract supertype for selection model specifications. Concrete subtypes carry selection parameters and implement select! to compute post-selection genotype frequencies from pre-selection genotype frequencies.

FrequencyDependentSelection <: SelectionModel

Frequency-dependent selection model following Mallet & Barton (1989, eq. 1–3).

Fitness of each phenotype is:

w_i = 1 - s × (1 - p_i)

where p_i is the frequency of phenotype i in the local deme and s is the selection coefficient. Rare phenotypes have lower fitness because predators do not recognise them as aposematic warning patterns. Common phenotypes are favoured, creating positive frequency-dependent selection that stabilises Müllerian mimicry clines in hybrid zones.

Fields

• s::Float64: selection coefficient in [0, 1]

Examples

julia> using HybridZones

julia> model = FrequencyDependentSelection(s=0.3)
FrequencyDependentSelection(s=0.3)

julia> model isa SelectionModel
true

julia> selection_coefficient(model)
0.3

select!(freq_next, freq_current, model::FrequencyDependentSelection, arch::OneLocusDiploid)

Apply one generation of frequency-dependent selection in-place.

Computes post-selection genotype frequencies from freq_current and writes them into freq_next. Both matrices must have the same size.

Matrix layout convention

State matrices use the convention rows = genotypes, columns = demes. For a OneLocusDiploid architecture, this means a 3×n_demes matrix:

This layout is chosen for cache efficiency: Julia's column-major storage means each deme's three genotype frequencies are contiguous in memory. Users coming from R or Python may need to transpose their input data, since those languages typically use rows-as-observations conventions.

Genotype-to-phenotype mapping

The mapping from genotypes to phenotypes is determined by arch.dominance:

• :codominant: each genotype is its own phenotype (three phenotypes total)
• :dominant_first: A1A1 and A1A2 share the A1 phenotype; A2A2 is distinct
• :recessive_first: A1A1 is distinct; A1A2 and A2A2 share the A2 phenotype

Contract

• freq_next is overwritten; freq_current is read-only.
• Each column of freq_current must sum to 1.0 (valid genotype frequencies).
• After return, each column of freq_next sums to 1.0 within floating-point precision.

Examples

julia> using HybridZones

julia> model = FrequencyDependentSelection(s=0.0);

julia> arch = OneLocusDiploid();

julia> freq = reshape([0.25, 0.50, 0.25], 3, 1);

julia> freq_next = similar(freq);

julia> select!(freq_next, freq, model, arch);

julia> freq_next ≈ freq  # s=0 → selection has no effect
true

select!(g_next, g_current, model::SemiDominantFrequencyDependentSelection, arch::OneLocusDiploid)

Apply one generation of semi-dominant frequency-dependent selection in-place.

Computes post-selection genotype frequencies from g_current using Pascal's semi-dominant phenotype mixing (WC1SEDO.PAS lines 720–735) and writes them into g_next. Both matrices must have the same size. Only :codominant dominance is supported; other modes raise an ArgumentError.

Semi-dominant mixing

The perceived phenotype frequency for each genotype includes contributions from related phenotypes:

sim_A1A1 = p_A1A1 + 0.5 * p_A1A2
sim_A1A2 = p_A1A2 + 0.5 * (p_A1A1 + p_A2A2)
sim_A2A2 = p_A2A2 + 0.5 * p_A1A2

Fitness for each genotype is w_G = 1 - s * (1 - sim_G).

See FrequencyDependentSelection for the strict codominant alternative where each genotype is its own distinct phenotype.

Examples

julia> using HybridZones

julia> model = SemiDominantFrequencyDependentSelection(0.0);

julia> arch = OneLocusDiploid(dominance = :codominant);

julia> freq = reshape([0.25, 0.50, 0.25], 3, 1);

julia> freq_next = similar(freq);

julia> select!(freq_next, freq, model, arch);

julia> freq_next ≈ freq  # s=0 → selection has no effect
true

selection_coefficient(model::FrequencyDependentSelection) -> Float64

Return the selection coefficient s of model.

selection_coefficient(model::SemiDominantFrequencyDependentSelection) -> Float64

Return the selection coefficient s of model.

SemiDominantFrequencyDependentSelection(s)

Frequency-dependent selection with semi-dominant phenotype mixing.

This model implements the selection regime used in Mallet's WC1SEDO Pascal program. Unlike FrequencyDependentSelection, which treats each genotype as its own distinct phenotype under :codominant dominance, this model treats heterozygotes as having partial similarity to both homozygous parental phenotypes.

For each genotype, a "perceived phenotype frequency" is computed that includes contributions from related phenotypes via semi-dominance mixing:

• Each homozygote receives a contribution of half the heterozygote frequency (the heterozygote looks half-like each parent)
• The heterozygote receives a contribution of half of each homozygote frequency (symmetric, since heterozygote looks half-like both parents)

The fitness of each genotype is then computed from its perceived phenotype frequency:

w_G = 1 - s * (1 - sim_G)

Selection acts proportionally on each genotype, with normalization by mean fitness.

Arguments

• s: selection coefficient, in [0, 1]

Biological interpretation

Semi-dominant frequency-dependent selection is biologically appropriate for warning-color systems where heterozygotes have visually intermediate phenotypes that share partial similarity with both parental forms, and predators have partial recognition across similar patterns. This is the model used by Mallet & Barton (1989) and the broader Pascal simulator tradition for Heliconius hybrid zones where the heterozygote wing pattern is intermediate between the parental races.

For systems where heterozygotes have a distinctly different phenotype recognized as its own pattern class, use FrequencyDependentSelection with :codominant dominance instead.

Dispatch

This model dispatches only on OneLocusDiploid with :codominant dominance. Combining it with :dominant_first or :recessive_first will raise a MethodError. The semi-dominance is intrinsic to the selection model itself.

Examples

julia> using HybridZones

julia> arch = OneLocusDiploid(dominance = :codominant);

julia> sel = SemiDominantFrequencyDependentSelection(0.3);

MigrationModel

Abstract supertype for migration model specifications. Concrete subtypes carry the migration kernel and parameters needed to redistribute frequencies across demes via migrate!.

BinomialStepping <: MigrationModel

Binomial stepping-stone migration model as used by Mallet & Barton (1989).

The migration kernel is a Binomial(2n, α) PMF evaluated at displacements -n:n from each source deme, where n is the maximum migration distance derived from the migration variance σ² and symmetry parameter α:

n = ceil(Int, σ² / (2α(1-α)))

The kernel is stored in full: kernel[n+1+i] holds the probability of displacement i for i ∈ -n:n, so the kernel has length 2n+1 and sums to 1.0.

Boundary demes receive the frequency of the nearest in-bounds deme to fill the kernel's reach beyond the transect edge (reflective boundary, matching the Pascal source's Reflection procedure).

Fields

• σ²::Float64: migration variance per generation
• α::Float64: symmetry parameter (0.5 = symmetric)
• n::Int: maximum migration distance (half-width of kernel)
• kernel::Vector{Float64}: migration weights of length 2n+1

Examples

julia> using HybridZones

julia> m = BinomialStepping(10.0)
BinomialStepping(σ²=10.0, α=0.5, n=20)

julia> m isa MigrationModel
true

julia> migration_variance(m)
10.0

julia> max_distance(m)
20

julia> isapprox(sum(HybridZones.MigrationModels.kernel(m)), 1.0; atol=1e-12)
true

migrate!(freq_next, freq_current, m::BinomialStepping)

Apply one generation of binomial stepping-stone migration in-place.

Writes post-migration frequencies into freq_next from freq_current. Both vectors must have the same length (number of demes). Demes beyond the transect boundaries are treated as having the frequency of the nearest edge deme (reflective boundary).

Arguments

• freq_next::AbstractVector: output buffer; overwritten on return
• freq_current::AbstractVector: input frequencies, one value per deme
• m::BinomialStepping: migration model supplying kernel and range

Examples

julia> using HybridZones

julia> m = BinomialStepping(10.0);

julia> freq = fill(0.3, 100);

julia> freq_next = similar(freq);

julia> migrate!(freq_next, freq, m);

julia> all(isapprox.(freq_next, 0.3; atol=1e-10))
true

migration_variance(m::BinomialStepping) -> Float64

Return the migration variance per generation σ².

max_distance(m::BinomialStepping) -> Int

Return n, the maximum displacement (in demes) that the kernel reaches.

kernel(m::BinomialStepping) -> Vector{Float64}

Return the full migration weight vector of length 2n+1. Element n+1+i is the probability of displacement i for i ∈ -n:n.

MatingModel

Abstract supertype for mating model specifications. Concrete subtypes implement mate! to compute post-mating genotype frequencies from pre-mating genotype frequencies, typically restoring Hardy-Weinberg proportions within each deme.

RandomMating <: MatingModel

Random mating within each deme, restoring Hardy-Weinberg (HW) proportions each generation.

This is a zero-field struct: random mating is parameter-free. The type exists to dispatch mate!, not to carry data.

Under random mating, the three diploid genotype frequencies at any deme are fully determined by the allele frequency p (frequency of A1):

A1A1 = p²
A1A2 = 2p(1-p)
A2A2 = (1-p)²

Migration mixing across demes introduces Wahlund-effect heterozygote deficits relative to HW expectations. Applying mate! each generation restores HW within demes, matching the biological convention used in Mallet's Pascal simulators and the broader population genetics literature.

Examples

julia> using HybridZones

julia> mat = RandomMating()
RandomMating()

julia> mat isa MatingModel
true

mate!(g_next, g_current, model::RandomMating, arch::OneLocusDiploid)

Apply one generation of random mating in-place, restoring Hardy-Weinberg proportions within each deme.

For each deme j, computes the A1 allele frequency:

p = g_current[1, j] + 0.5 × g_current[2, j]

then writes HW genotype frequencies into g_next:

g_next[1, j] = p²
g_next[2, j] = 2p(1-p)
g_next[3, j] = (1-p)²

Random mating is dominance-independent: the HW restoration formula is the same regardless of whether the locus is :codominant, :dominant_first, or :recessive_first. Dispatching on OneLocusDiploid rather than the abstract GeneticArchitecture mirrors the select! pattern and leaves room for architecture-specific mating logic in future subtypes.

Contract

• g_next is overwritten; g_current is read-only.
• Both matrices must have the same size.
• The first dimension of g_current must be 3 (one row per genotype).
• After return, each column of g_next sums to 1.0 within floating-point precision.

Examples

julia> using HybridZones

julia> mat = RandomMating();

julia> arch = OneLocusDiploid();

julia> g = reshape([0.0, 1.0, 0.0], 3, 1);  # all heterozygotes, p=0.5

julia> g_next = similar(g);

julia> mate!(g_next, g, mat, arch);

julia> g_next
3×1 Matrix{Float64}:
 0.25
 0.5
 0.25

simulate(arch, mat, sel, mig, initial_state; n_generations, lifecycle_order=:mate_migrate_select, save_trajectory=true)

Run a forward simulation of a hybrid zone for n_generations generations.

Each generation applies mating, selection, and migration in the order specified by lifecycle_order. The default :mate_migrate_select ordering — mating, then migration, then selection — matches Mallet's Pascal reference implementation and is the recommended default for warning-color hybrid zone models.

Lifecycle orderings

• :mate_migrate_select (default): mating restores Hardy-Weinberg within demes, then migration redistributes frequencies across demes, then selection acts on local frequencies. This matches Jim Mallet's standard Pascal ordering.
• :mate_select_migrate: mating, then selection, then migration. An alternative ordering for organisms where dispersal occurs after selection.
• :select_migrate_mate: selection first, then migration, then mating at the end of the generation. Biologically valid for organisms where dispersal occurs late in the life cycle after selection has acted, with mating completing the generation.

Other values raise an ArgumentError.

Breaking change from previous API

The previous signature was simulate(arch, sel, mig, initial_state; ...). The new signature inserts mat::MatingModel as the second positional argument: simulate(arch, mat, sel, mig, initial_state; ...). No backward-compatible default method is provided — the behavior change (HW restoration each generation) is significant enough to require explicit user choice.

Buffer reuse

The inner loop preallocates working buffers and updates them in place, avoiding per-generation heap allocation. For save_trajectory = true, two intermediate scratch buffers are allocated alongside the trajectory array. For save_trajectory = false, three buffers are rotated per operation each generation.

Migration convention

Migration is applied row-by-row (one genotype at a time) by calling migrate! for each genotype row. This keeps migrate! operating on AbstractVector as documented, matching the per-genotype stepping-stone logic in the Pascal reference implementation.

Arguments

• arch::OneLocusDiploid: genetic architecture
• mat::MatingModel: mating model (e.g. RandomMating())
• sel::SelectionModel: selection model
• mig::MigrationModel: migration model
• initial_state::AbstractMatrix: n_genotypes × n_demes genotype frequency matrix; each column must sum to approximately 1.0

Keyword arguments

• n_generations::Int: number of generations to simulate; must be positive
• lifecycle_order::Symbol: :mate_migrate_select (default), :mate_select_migrate, or :select_migrate_mate
• save_trajectory::Bool: if true (default), return all generations; if false, return only the final state (more memory-efficient for long runs)

Returns

• save_trajectory = true: Array{Float64, 3} of shape (n_genotypes, n_demes, n_generations + 1). Slice [:, :, 1] is initial_state; slice [:, :, g + 1] is the state after generation g.
• save_trajectory = false: Matrix{Float64} of shape (n_genotypes, n_demes) with the state after n_generations generations.

Examples

julia> using HybridZones

julia> arch = OneLocusDiploid(dominance = :codominant);

julia> mat = RandomMating();

julia> sel = FrequencyDependentSelection(s = 0.3);

julia> mig = BinomialStepping(30.0);

julia> initial = secondary_contact(arch);

julia> trajectory = simulate(arch, mat, sel, mig, initial; n_generations = 10);

julia> size(trajectory)
(3, 101, 11)

julia> trajectory[:, :, 1] == initial
true

julia> final_cline = allele_frequencies(trajectory[:, :, end], arch);

julia> length(final_cline)
101

secondary_contact(arch::OneLocusDiploid; n_demes=101) -> Matrix{Float64}

Return an initial genotype frequency matrix for the secondary contact scenario.

The transect is split into three zones:

• Left half (demes 1 to (n_demes - 1) ÷ 2): fixed A1A1 (freq[1, j] = 1)
• Contact deme ((n_demes + 1) ÷ 2): all heterozygotes (freq[2, j] = 1)
• Right half (remaining demes): fixed A2A2 (freq[3, j] = 1)

With the default n_demes = 101, the contact zone is at deme 51, the geometric centre of the transect. An odd number of demes is recommended so the centre is an integer index and the two flanking regions are equal in length.

The contact deme is set to all heterozygotes rather than a 50/50 mix of the two homozygotes, following the Pascal reference implementation convention. A heterozygote-only contact is the natural starting point for studying introgression between two differentiated populations meeting for the first time.

Examples

julia> using HybridZones

julia> arch = OneLocusDiploid();

julia> state = secondary_contact(arch);

julia> size(state)
(3, 101)

julia> state[:, 1]
3-element Vector{Float64}:
 1.0
 0.0
 0.0

julia> state[:, 51]
3-element Vector{Float64}:
 0.0
 1.0
 0.0

julia> state[:, 101]
3-element Vector{Float64}:
 0.0
 0.0
 1.0

allele_frequencies(state::AbstractMatrix, arch::OneLocusDiploid) -> Vector{Float64}

Compute the frequency of the A2 allele in each deme from a genotype frequency matrix.

For OneLocusDiploid, the A2 frequency at deme j is:

p_A2[j] = freq[3, j] + 0.5 × freq[2, j]

which counts A2A2 homozygotes fully and heterozygotes at half weight.

Arguments

• state: 3×n_demes genotype frequency matrix (rows: A1A1, A1A2, A2A2; columns: demes)
• arch: genetic architecture; currently only OneLocusDiploid is supported

Examples

julia> using HybridZones

julia> arch = OneLocusDiploid();

julia> state = reshape([0.25, 0.50, 0.25], 3, 1);

julia> allele_frequencies(state, arch)
1-element Vector{Float64}:
 0.5

allele_frequencies(trajectory, arch::OneLocusDiploid) -> Matrix{Float64}

Compute A2 allele frequencies across all time points in a simulation trajectory.

Applies allele_frequencies(state, arch) slice-by-slice to a 3D trajectory array of shape (n_genotypes, n_demes, n_time).

Returns

A matrix of shape (n_demes, n_time) where entry [j, t] is the frequency of A2 at deme j at time point t. For a trajectory returned by simulate with n_generations generations, t = 1 is the initial state and t = n_generations + 1 is the final state.

Arguments

• trajectory: 3D array of shape (n_genotypes, n_demes, n_time) from simulate
• arch: genetic architecture