Beak Of The Finch

We recently visited the Galapagos islands and to prepare I read The Beak of the Finch by Jonathan Weiner: a Pulitzer prize winning examination of evolution following the work of the Grants on the finches on Daphne Major. You can see the island as you land in the airport at Baltra: remarkable that they lived on it for months at a time.

The surprising bit of the book is how fast evolution works in the Galapagos. My intuition from my education was that evolution and speciation were very slow processes. Apparently not necessarily! The book showcases that the strong selection effects on the islands can substantially change the phenotypic makeup of the finches in just a few (quite harsh) selection episodes. They even have a unit for it: “the darwin”.

When going through the book I was mostly following the logic but it struck me that this is perfectly suited for a simulation to really drive home the dynamics. And thus…

The first dynamic that we’re introduce to in the book is character displacement. As Weiner notes:

Darwinian competition is not only the clas of stag horns, the gore on the jaws of lions, nature red in tooth and claw. Competition also can be a silent race, side by side, for the last food on a desert island, where te competititors never fight one another, and the only sound of battle is the occasional crack of a Tribulus seed.

and to quote Peter Thiel: “competition is for losers”. Thus what happens in nature is that similar species evolve to differ from each other in order to minimize competition for limited resources.

Let’s set up a simulation: we’ll imagine that we have supplies of resources \(X_{i} \sim Uniform(0,1)\). Individual birds utilize the resources from the range \([a_{i} - \delta, a_{i} + \delta]\) . When they reproduce they pass along their \(a_{i}\) plus some additional noise \(\epsilon \sim N(0, 0.05)\)

using DataFrames
using Distributions
using Random
using AlgebraOfGraphics
using CairoMakie

CARRYING_CAPACITY = 1000

struct Bird
    a::Any
    delta::Any
    species::Any
end

function reproduce(bird::Bird)
    return Bird(bird.a + rand(Normal(0, 0.05)), bird.delta, bird.species)
end

function step(birds)
    resources = rand(CARRYING_CAPACITY)

    next_generation = Bird[]
    for bird in shuffle(birds)
        food = findfirst(x -> abs(bird.a - x) <= bird.delta, resources)
        isnothing(food) && continue # bird starves
        deleteat!(resources, food)
        push!(next_generation, reproduce(bird))
        push!(next_generation, reproduce(bird))
    end

    return next_generation
end

function simulate(n_iters)
    subpopulation_left = Bird.(rand(Normal(0.45, 0.1), 500), 0.05, 1)
    subpopulation_right = Bird.(rand(Normal(0.55, 0.1), 500), 0.05, 2)

    population = [subpopulation_left; subpopulation_right]

    dfs = []

    for iter in 1:n_iters
        push!(dfs,
              DataFrame(
                  :iter => iter,
                  :targets => [bird.a for bird in population],
                  :species => [bird.species for bird in population],
              ),
              )
        population = step(population)
    end

    return vcat(dfs...)
end

And we see starrting from a fairly close distribution we see the species start to seprate into two separate clusters. We can make this even starker by introducing a non-uniform resource distribution

function step(birds)
    resources = rand(MixtureModel([Normal(0.2, 0.05), Normal(0.8, 0.05)], [0.5, 0.5]), CARRYING_CAPACITY)

    next_generation = Bird[]
    for bird in shuffle(birds)
        food = findfirst(x -> abs(bird.a - x) <= bird.delta, resources)
        isnothing(food) && continue # bird starves
        deleteat!(resources, food)
        push!(next_generation, reproduce(bird))
        push!(next_generation, reproduce(bird))
    end

    return next_generation
end

There is something unsatisfying with this simulation however: we’re not exactly seeing the creation of two species. Indeed you see a buch of Species 1 overlapping with Species 2; species only comes in as a historical accident but otherwise has no effect. So let’s fix that.

Apparently the definition of a species is a little more complicated then I thought: originally I was under the impression that species can’t interbreed. Turns out it’s more like they usually don’t interbreed. You can get this through sexual selection: if hybrids are less fit then the birds that prefer to, um, flock with birds of the same feather will have an advantage.

struct Bird
    a::Any
    delta::Any
    species_preference::Any
    species::Any
end

function reproduce(bird_a::Bird, bird_b::Bird)
    return Bird((bird_a.a + bird_b.a) / 2.0 + rand(Normal(0, 0.05)),
                bird_a.delta,
                bird_a.species_preference + rand(Normal(0, 0.05)),
                bird_a.species)
end

function step(birds)
    resources = rand(MixtureModel([Normal(0.2, 0.05), Normal(0.8, 0.05)]), CARRYING_CAPACITY)

    survivors = []
    for bird in shuffle(birds)
        food = findfirst(x -> abs(bird.a - x) <= bird.delta, resources)
        isnothing(food) && continue # bird starves
        push!(survivors, bird)
        deleteat!(resources, food)
    end

    next_generation = []
    for bird in survivors
        mate_species = rand() <= bird.species_preference ? bird.species : 3 - bird.species
        others = shuffle(survivors)
        mate = findfirst(x -> x.species == mate_species, others)
        isnothing(mate) && continue
        push!(next_generation, reproduce(bird, others[mate]))
        push!(next_generation, reproduce(bird, others[mate]))
    end

    return next_generation
end

function simulate(n_iters)
    subpopulation_left = Bird.(rand(Normal(0.45, 0.1), 500), 0.05, rand(Normal(0.9, 0.1), 500), 1)
    subpopulation_right = Bird.(rand(Normal(0.55, 0.1), 500), 0.05, rand(Normal(0.9, 0.1), 500), 2)

    population = [subpopulation_left; subpopulation_right]

    dfs = []

    for iter = 1:n_iters
        push!(
            dfs,
            DataFrame(
                :iter => iter,
                :targets => [bird.a for bird in population],
                :species => [bird.species for bird in population],
            ),
        )
        population = step(population)
    end

    return vcat(dfs...)
end

Of course we’re being a little prescriptive with our definition of species: it’s not really an enumeration but rather a cluster. Let’s give each bird a “trait affinity” for how they select a partner: higher affinity means they strongly prefer mates with similar beak size. Mates are selected with probability proportional to the product of each bird’s preference for the other.

struct Bird
    a::Float64
    delta::Float64
    affinity::Float64  # mate preference strength: higher → prefers more similar-looking mates
end

function mate_weight(bird_a::Bird, bird_b::Bird)
    diff_sq = (bird_a.a - bird_b.a)^2
    return exp(-(bird_a.affinity + bird_b.affinity) * diff_sq)
end

function reproduce(bird_a::Bird, bird_b::Bird)
    child_a = (bird_a.a + bird_b.a) / 2.0 + rand(Normal(0, 0.05))
    child_affinity = (bird_a.affinity + bird_b.affinity) / 2.0 + rand(Normal(0, 0.3))
    return Bird(child_a, bird_a.delta, max(0.0, child_affinity))
end

function step(birds; bimodal_resources = true)
    resources = bimodal_resources ?
        rand(MixtureModel([Normal(0.2, 0.05), Normal(0.8, 0.05)], [0.5, 0.5]), CARRYING_CAPACITY) :
        rand(Uniform(0, 1), CARRYING_CAPACITY)

    survivors = Bird[]
    for bird in shuffle(birds)
        food = findfirst(x -> abs(bird.a - x) <= bird.delta, resources)
        isnothing(food) && continue
        push!(survivors, bird)
        deleteat!(resources, food)
    end

    isempty(survivors) && return Bird[], Bird[]

    next_generation = Bird[]
    for bird in survivors
        weights = [mate_weight(bird, other) for other in survivors]
        total = sum(weights)
        total == 0.0 && continue
        mate = survivors[rand(Distributions.Categorical(weights ./ total))]
        push!(next_generation, reproduce(bird, mate))
        push!(next_generation, reproduce(bird, mate))
    end

    return survivors, next_generation
end

function simulate(n_iters; bimodal_resources = true)
    subpopulation_left  = Bird.(rand(Normal(0.2, 0.05), 500), 0.05, clamp.(rand(Normal(2.0, 0.5), 500), 0.0, Inf))
    subpopulation_right = Bird.(rand(Normal(0.8, 0.05), 500), 0.05, clamp.(rand(Normal(2.0, 0.5), 500), 0.0, Inf))

    population = [subpopulation_left; subpopulation_right]
    dfs = []

    for iter = 1:n_iters
        survivors, population = step(population; bimodal_resources)
        isempty(survivors) && break

        push!(dfs, DataFrame(
            :iter     => iter,
            :a        => [bird.a for bird in survivors],
            :affinity => [bird.affinity for bird in survivors],
        ))
    end

    return vcat(dfs...)
end

We can verify that affinity is being selected for by comparing against a control where resources are distributed uniformly (no penalty for hybridization). Without the penalty, being picky about mates is costly as less picky birds are more likely to mate so affinity should drift downward. With the penalty, avoiding costly crosses is advantageous, so affinity should rise.

It’s cool to see how straightforward it is to express these ideas in simulation. Reading the book was fascinating but actually implementing this helped me understand the dynamics much better. And of course seeing the outcomes in person was wild: the naturalists showed us the subtle differences in all of lizards as we went from island to island. Truly enchanted.