up

src/polbarton.jl

@@ -39,11 +39,11 @@ end
 """
     MixedPloidyDeme{A,T}
 
-A single mixed-ploidy deme, most of the 'population
+Single mixed-ploidy deme, most of the 'population
 genetic' environment is implemented at this level (drift, selection,
 mutation). OV: Viability matrix, a symmetric matrix that contains the viability
 of offspring for each possible combination of gametes. Ug: Unreduced gamete
-formation matrix, a matrix that contains the probability of unreduced gametes
+formation rate matrix, a matrix that contains the probability of unreduced gamete formation
 for each level of ploidy in the population.
 
 - `K` : Carrying capacity
@@ -120,10 +120,6 @@ end
     Habitat2D{D}
 
 Two-dimensional habitat (matrix of connected demes).
-- `σ` : Variance of dispersal
-- `b` : Steepness of linear gradient
-- `θ` : Phenotypic optimum in the central deme of the habitat
-- `Dm` : Number of demes
 """
 @with_kw struct Habitat2D{D,T}
     demes::Matrix{D}
@@ -194,13 +190,6 @@ expected_heterozygosity(H₀, t, N) = ((1.0-1.0/N)^t)*H₀
 
 #Mating:
 
-"""
-	viability(a::Agent, b::Agent, d::deme)
-"""
-function viability(a::Agent, b::Agent, d::AbstractDeme)
-	return d.OV[ploidy(a), ploidy(b)]
-end
-
 """
 	recombination(a::Agent)
 
@@ -251,10 +240,19 @@ function gametogenesis(a::Agent, d::AbstractDeme)
     return Agent(a.loci[idx, :], a.d)
 end	
 
+"""
+	viability(a::Agent, b::Agent, d::deme)
+After fusing of gametes, calculate the viability of resulting individual. For example the viability of
+a triploid individual is by default set to 0, diploid and tetraploid individuals set to 1.
+"""
+function viability(a::Agent, b::Agent, d::AbstractDeme)
+	return d.OV[ploidy(a), ploidy(b)]
+end
+
 """
 	meiosis(a::Agent,c)
 
-Meiosis in a diploid.
+UNDER CONSTRUCTION. Meiosis in a diploid.
 - `c` : Recombination rate
 """
 function meiosis(a::Agent,c)

tutorial.jl

@@ -0,0 +1,284 @@
+### A Pluto.jl notebook ###
+# v0.14.7
+
+using Markdown
+using InteractiveUtils
+
+# ╔═╡ 1e540d75-0a9d-4bb2-85aa-02ea78f08766
+using Parameters, Random, Distributions, Plots, StatsBase, PlutoUI, ColorSchemes, DataFrames, StatsPlots
+
+# ╔═╡ 42a769f0-1bf8-11ec-2a76-e7cea7f2e77a
+import PolyStab:randagent
+
+# ╔═╡ 9fe6dd08-89c2-49c7-99ed-bf27e4a90242
+md"##### Generating agents"
+
+# ╔═╡ 111bd54b-d480-4c5e-9a92-aa6f38a2d074
+md""" Generating a single agent:"""
+
+# ╔═╡ ee0a9506-823f-449d-8bb9-7f864c8e73ae
+"""
+    randagent(p, α, n, k; d)
+
+-n: number of loci
+-k: ploidy l
+-α: allelic effect size
+-p: starting allele frequency
+"""
+a = randagent(0.5, 0.1, 100, 2, d=1.)
+
+# ╔═╡ 54ea2331-7430-48c6-b0b9-f82c4b80f3ff
+a.loci
+
+# ╔═╡ 1c163c4d-d8dd-49a4-a07a-372abf74ef21
+import PolyStab:trait
+
+# ╔═╡ 3f95bac8-5c98-472f-b53d-f5f2ebe2e113
+trait(a)
+
+# ╔═╡ 2e17fd8e-f360-421e-ba35-8921c9807c47
+md""" It is also possible to generate an array of agents:"""
+
+# ╔═╡ 4d4f91c4-116e-45bd-aca5-a4e728026159
+import PolyStab:randagent_p
+
+# ╔═╡ de22de55-116d-41fd-8865-a9a0ae5be0c9
+"""
+    randagent_P(p, α, n, k, N; d)
+
+-n: number of loci
+-k: ploidy l
+-α: allelic effect size
+-p: starting allele frequency
+-N: number of individuals
+"""
+b = randagent_p(0.5, 0.1, 100, [0.,1.,0.,0.], 10, d=1.)
+
+# ╔═╡ 7c3af315-31b4-40bb-af86-8a105f427b7c
+b
+
+# ╔═╡ 5e0aaf09-4453-4aea-b0f9-5f14b16200f5
+md"""#### Generating a single mixed-ploidy deme"""
+
+# ╔═╡ 870aab66-3e7a-4b47-917d-190a72f8e2de
+import PolyStab:MixedPloidyDeme
+
+# ╔═╡ 5d5393dd-bebe-4275-992b-7477a5c7834b
+"""
+    MixedPloidyDeme{A,T}
+
+Single mixed-ploidy deme, most of the 'population
+genetic' environment is implemented at this level (drift, selection,
+mutation). OV: Viability matrix, a symmetric matrix that contains the viability
+of offspring for each possible combination of gametes. Ug: Unreduced gamete
+formation rate matrix, a matrix that contains the probability of unreduced gamete formation for each level of ploidy in the population.
+
+- `K` : Carrying capacity
+- `θ` : Environmental optimum
+- `rm`: Growth rate per capita
+- `Vs`: Variance of stabilizing selection
+- `α` : Allelic effect size
+- `μ` : Mutation rate
+- `OV` : Offspring viability
+- `UG` : Unreduced gamete probabilities
+"""
+d_p1 = MixedPloidyDeme(agents = randagent_p(0.5, 0.2, 200, [1., 0., 0., 0.],10), OV = [1. 0. 0. 0. ; 0. 1. 0. 0. ; 0. 0. 0. 0. ; 0. 0. 0. 0.], UG = [1. 0. 0. 0. ; 1. 0. 0. 0. ; 0. 0. 0. 0. ; 0. 1. 0. 0.], θ = 20., Vs = 0.1)
+
+# ╔═╡ 6b866b78-80d7-43b1-a4c5-33f6c7908a03
+d_p1
+
+# ╔═╡ 1dc54ec2-f4d3-4237-8fb2-b331e0bf887f
+d_p1.agents
+
+# ╔═╡ 6e07be7d-8af6-4bf1-9591-f42ab3ea9db0
+trait.(d_p1.agents)
+
+# ╔═╡ a24f94e4-c012-4c82-ba4c-93484d0a7263
+md"""#### Simulating a single mixed-ploidy deme"""
+
+# ╔═╡ e2ced766-0032-4cb2-b73c-fd1a086e0950
+import PolyStab:evolving_selectiondeme
+
+# ╔═╡ 7ac2f20c-db6a-4927-a76b-d7ded656adbc
+"""
+	evolving_selectiondeme(d::AbstractDeme, ngen)
+
+Simulate a single deme with mixed ploidy, malthusian fitness and unreduced
+gamete formation.
+
+function evolving_selectiondeme(d::MixedPloidyDeme, ngen; 
+	heterozygosities_p=heterozygosities_p, fit=malthusian_fitness, trait_mean = trait_mean, allelefreqs_p = 
+	allelefreqs_p, pf = ploidy_freq, fta = f_trait_agents)
+	het = [heterozygosities_p(d)]
+	pop = [length(d)]
+	tm = [trait_mean(d)]
+	af = [allelefreqs_p(d)]
+	p2 = [ploidy_freq(d)[2]]
+	p3 = [ploidy_freq(d)[3]]
+	p4 = [ploidy_freq(d)[4]]
+	fta = [f_trait_agents(d)]
+	
+	for n=1:ngen
+		d = mating_PnB_x(d)
+		d = mutate(d) 
+		push!(het, heterozygosities_p(d))
+		push!(pop, length(d))
+		push!(tm, trait_mean(d))
+		push!(af, allelefreqs_p(d))
+		push!(p2, ploidy_freq(d)[2])
+		push!(p3, ploidy_freq(d)[3])
+		push!(p4, ploidy_freq(d)[4])
+		push!(fta, f_trait_agents(d))
+		
+	end
+	(pop=pop, deme=d, p2=p2, p3=p3, p4=p4, ngen=ngen, het=het,tm=tm, af=af, fta=fta) 
+end
+"""
+
+
+# ╔═╡ bb6aa892-6ef1-480a-8a2e-be1c131c0d5b
+evolving_selectiondeme(d_p1, 500)
+
+# ╔═╡ f91908fb-8202-4087-9ee0-ed8207a41378
+md""" ###### Components of the mating function"""
+
+# ╔═╡ a50b4da8-39d5-4977-a1de-a8a172c0e6ac
+"""
+	mating_PnB_x(d::MixedPloidyDeme{A})
+
+Mating in a mixed ploidy deme with unreduced gamete formation and partner
+choice weighted by malthusian fitness.
+
+function mating_PnB_x(d::AbstractDeme{A}) where A
+	new_agents = A[]
+	fitnesses = exp.(malthusian_fitness(d))
+	for i=1:length(d)
+		B1 = d.agents[i]
+		noff = number_of_offspring(d, B1)
+        Bs = sample(d.agents, weights(fitnesses), noff) 
+        offspring = filter(!ismock, map(B2->mate_p(B1, B2, d), Bs))
+        new_agents = vcat(new_agents, offspring)
+    end
+    d(new_agents)
+end
+"""
+
+# ╔═╡ d6eeb3b0-fbce-4ed8-b32b-0adbae74019c
+"""
+	malthusian_fitness(d::AbstractDeme)
+Return the Malthusian fitness (density dependence and stabilizing selection) of each agent in a deme.
+
+function malthusian_fitness(d::MixedPloidyDeme)
+    N = length(d)
+    fitnesses = Float64[]
+	for agent in d.agents
+		z = trait(agent)
+    	f = d.rm*(1-(N/d.K))-((z-d.θ)^2)/(2*d.Vs)
+		push!(fitnesses, f)
+	end
+	fitnesses
+end
+"""
+
+# ╔═╡ 3db37b1c-693b-4c34-9e18-b9126b0fe235
+"""
+	number_of_offspring(d::AbstractDeme, a::Agent)
+
+function number_of_offspring(d::MixedPloidyDeme, a::Agent)
+    logw = malthusian_fitness(d, a)
+    rand(Poisson(exp(logw)))
+end
+"""
+
+# ╔═╡ 776a2d51-ca39-4368-9431-5a3ad8afe14d
+"""
+	mate_p(a::Agent, b::Agent, d::AbstractDeme)
+
+Mating in a mixed ploidy deme. Assumes that different cytotypes can be
+compatible with a decrease in viability (cfr. OffspringViability matrix) (i.e.
+when #individuals with different ploidy hybridize, they generate have a
+probability p to generate no viable offspring).  Selfing is allowed without
+cost. This influences the dynamics by including hybrid offspring that can
+compete for space (and affects the malthusian fitness/density dependence of
+selection).
+	
+function mate_p(a::Agent, b::Agent, d::AbstractDeme)
+    #gamete formation
+    # XXX we should have recombinatoin *before* gamete formation ?!
+    ag = gametogenesis(recombination(a), d)
+    bg = gametogenesis(recombination(b), d)
+    #combine gametes and assign viability
+    via = viability(ag, bg, d) 
+    # not sure if returning 0 is the best idea in terms of type stability
+    # perhaps you should return a mock-agent (for instance an agent with length
+    # zero genome) when there is no viable offspring
+	# if d is defined at agent level, how do we cope with it here? 
+    return rand() < via ? Agent(loci=[ag.loci ; bg.loci], d=a.d) : 0
+end
+"""
+
+# ╔═╡ 82e7cf93-ba18-483a-9880-4dd396416510
+"""
+	recombination(a::Agent)
+
+Free recombination between loci in a mixed ploidy population. 
+
+function recombination(a::Agent)
+    genome = similar(a.loci)
+    for i in 1:nloci(a)
+        genome[:,i] = shuffle(a.loci[:,i])
+	end
+    Agent(loci=genome, d=a.d) 
+end
+
+
+	gametogenesis(a::Agent, d::AbstractDeme)
+
+Gamete formation in a mixed ploidy population (1n to 4n).
+
+function gametogenesis(a::Agent, d::AbstractDeme)
+	num = sample(1:4, weights(d.UG[ploidy(a),:]))
+    idx = sample(1:ploidy(a), num, replace=false)
+    return Agent(a.loci[idx, :], a.d)
+end	
+
+
+	viability(a::Agent, b::Agent, d::deme)
+After fusing of gametes, calculate the viability of resulting individual. For example the viability of
+a triploid individual is by default set to 0, diploid and tetraploid individuals set to 1.
+
+function viability(a::Agent, b::Agent, d::AbstractDeme)
+	return d.OV[ploidy(a), ploidy(b)]
+end
+"""
+
+
+# ╔═╡ Cell order:
+# ╠═1e540d75-0a9d-4bb2-85aa-02ea78f08766
+# ╠═42a769f0-1bf8-11ec-2a76-e7cea7f2e77a
+# ╟─9fe6dd08-89c2-49c7-99ed-bf27e4a90242
+# ╟─111bd54b-d480-4c5e-9a92-aa6f38a2d074
+# ╠═ee0a9506-823f-449d-8bb9-7f864c8e73ae
+# ╠═54ea2331-7430-48c6-b0b9-f82c4b80f3ff
+# ╠═1c163c4d-d8dd-49a4-a07a-372abf74ef21
+# ╠═3f95bac8-5c98-472f-b53d-f5f2ebe2e113
+# ╟─2e17fd8e-f360-421e-ba35-8921c9807c47
+# ╠═4d4f91c4-116e-45bd-aca5-a4e728026159
+# ╠═de22de55-116d-41fd-8865-a9a0ae5be0c9
+# ╠═7c3af315-31b4-40bb-af86-8a105f427b7c
+# ╟─5e0aaf09-4453-4aea-b0f9-5f14b16200f5
+# ╠═870aab66-3e7a-4b47-917d-190a72f8e2de
+# ╠═5d5393dd-bebe-4275-992b-7477a5c7834b
+# ╠═6b866b78-80d7-43b1-a4c5-33f6c7908a03
+# ╠═1dc54ec2-f4d3-4237-8fb2-b331e0bf887f
+# ╠═6e07be7d-8af6-4bf1-9591-f42ab3ea9db0
+# ╟─a24f94e4-c012-4c82-ba4c-93484d0a7263
+# ╠═e2ced766-0032-4cb2-b73c-fd1a086e0950
+# ╠═7ac2f20c-db6a-4927-a76b-d7ded656adbc
+# ╠═bb6aa892-6ef1-480a-8a2e-be1c131c0d5b
+# ╟─f91908fb-8202-4087-9ee0-ed8207a41378
+# ╠═a50b4da8-39d5-4977-a1de-a8a172c0e6ac
+# ╠═d6eeb3b0-fbce-4ed8-b32b-0adbae74019c
+# ╠═3db37b1c-693b-4c34-9e18-b9126b0fe235
+# ╠═776a2d51-ca39-4368-9431-5a3ad8afe14d
+# ╠═82e7cf93-ba18-483a-9880-4dd396416510