Supplemental Materials for: To accept or reject heterospecific mates
Daizaburo Shizuka and Emily J. Hudson
This supplemental report is generated using R Markdown. All code required to reprodue Figures 1, 2, and 4 appear in code chunks.
Basic setup
First, generate normal distribution functions for conspecifics and heterospecifics. Conspecific distribution has mean = 0 and sd = 0.1. Heterospecifics have mean = 0.5 and sd = 0.15. This is just to generate overlapping distributions.
g = function(x) dnorm(x, mean = 0, sd = 0.2)
b = function(x) dnorm(x, mean = 0.5, sd = 0.2)We can visualize the trait distribution curves.
plot(g, ylim = c(0, 2), xlim = c(-0.5, 1), lwd = 2, ylab = "G(t) or B(t)", xlab = "Trait values (t)")
par(new = T)
plot(b, ylim = c(0, 2), xlim = c(-0.5, 1), lwd = 2, col = "red", ylab = "",
xlab = "")
Trait distributions of two species. Black=‘conspecifics’, red=‘heterospecifics’
We now calculate what proportion of conspecifics and heterospecifics are accepted under varying values of t. We implement this with a brute-force method with 100 values of t ranging from 0 to 1 (along the x-axis on above figure).
t = seq(0, 1, length = 1001)
Gt = list(0)
Bt = list(0)
# proportion of conspecifics accepted for each t
for (i in t) Gt[[which(t == i)]] = integrate(g, lower = -1, upper = i)$value
Gt.value = unlist(Gt)
# proportion of heterospecifics accepted for each t
for (i in t) Bt[[which(t == i)]] = integrate(b, lower = -1, upper = i)$value
Bt.value = unlist(Bt)As illustration, we can choose 4 threshold values–0.1, 0.25, 0.5, and 0.75–and show the proportion of conspecifics and heterospecifics that are accepted based on each value.
dat = data.frame(t = c(0.1, 0.25, 0.5, 0.75), prop_conspecific = Gt.value[c(101,
251, 501, 751)], prop_heterospecific = Bt.value[c(101, 251, 501, 751)],
line_on_figure = c("(1)", "(2)", "(3)", "(4)"))
dat## t prop_conspecific prop_heterospecific line_on_figure
## 1 0.10 0.6914622 0.02275013 (1)
## 2 0.25 0.8943499 0.10564977 (2)
## 3 0.50 0.9937900 0.50000000 (3)
## 4 0.75 0.9999113 0.89435023 (4)Figure 2 visualizes the threshold values above:
Locations of thresholds as shown in table above
Figure 1
Figure 1 shows how the optimal acceptance threshold (along the x-axis) shifts with the relative fitness consequence of accepting a heterospecific mate (i.e., producing hybrid offspring) vs. accepting a conspecific mate.
The following code will calculate which of the 1000 threshold values we have defined above generates the optimal fitness for the actor, and plot that value along the x-axis.
Note that here and throughout, we use “V” in place of parameter F* and “v” in place of parameter f because “F” means FALSE in R.
# set up we will plot 20 points
times = 20
best.t = vector(length = times)
# constant parameter values. 'V' stands in for 'F' in main text because 'F'
# means FALSE in R.
V = 10
P = 0.7
Cs = 1
# now, we vary the fitness consequence of accepting heterospecific mates (f,
# represented here by v), ranging from 0 (hybrids are inviable) to V (there
# is no fitness reduction for accepting heterospecific mates)
v = seq(0, V, length = times)
# now, calculate fitness for each threshold given a specific value of f, and
# estimate what threshold value gives the optimal fitness.
for (i in 1:times) {
Wt = (V * P * Gt.value + v[i] * (1 - P) * Bt.value - Cs)/(P * Gt.value +
(1 - P) * Bt.value)
best.t[i] = t[which.max(Wt)]
}
# generate plot
plot(best.t, v/V, type = "p", xlim = c(0, 1), lwd = 2, las = 1, pch = 19, xlab = "Optimal Acceptance Threshold",
ylab = "Hybrid Fitness / Pure-bred Fitness")
Figure 2
Figure 2 shows how the rates of acceptance error can vary based on search costs or skews in abundance of conspecifics and heterospecifics. Here, we use parameter values that match Reeve (1989) Figure 5.
g = function(x) dnorm(x, mean = 0, sd = 0.1) #conspecifics
b = function(x) dnorm(x, mean = 0.5, sd = 0.15) #heterospecifics
t = seq(0, 1, length = 1001) #set up 1000 possible thresholds between 0 and 1
Gt = list(0)
Bt = list(0)
for (i in t) Gt[[which(t == i)]] = integrate(g, lower = -1, upper = i)$value
Gt.value = unlist(Gt)
for (i in t) Bt[[which(t == i)]] = integrate(b, lower = -1, upper = i)$value
Bt.value = unlist(Bt)# set up parameter values for Figure 2a
V = 10
v = 8
P = 0.7
# vary costs between 0 and 1.5 (i.e., 0 to 0.15 the fitness of mating with
# conspecific)
Cs = seq(0, 1.5, length = 100)
accept.error = vector(length = 100) #set up vector
# acceptance error is proportion of heterospecifics accepted.
for (i in 1:100) {
Wt = (V * P * Gt.value + v * (1 - P) * Bt.value - Cs[i])/(P * Gt.value +
(1 - P) * Bt.value)
accept.error[i] = Bt.value[which.max(Wt)]
}
plot(Cs/V, accept.error, xlim = c(0, 0.15), ylim = c(0, 1), type = "l", lwd = 2,
ylab = "Acceptance Error", xlab = "Search cost as proportion of maximum fitness (c/F)",
las = 1)
## Figure 2b
V = 10
v = 8
Cs = 0.2
# vary proportion population that are conspecifics
P = seq(0.01, 0.99, 0.01)
accept.error = vector(length = length(P))
for (i in 1:length(P)) {
Wt = (V * P[i] * Gt.value + v * (1 - P[i]) * Bt.value - Cs)/(P[i] * Gt.value +
(1 - P[i]) * Bt.value)
accept.error[i] = Bt.value[which.max(Wt)]
}
plot(P, accept.error, xlim = c(0, 0.99), ylim = c(0, 1), type = "l", ylab = "Acceptance Error",
xlab = "Proportion of encounters with conspecific", las = 1, lwd = 2, main = "",
xaxp = c(0, 1, 5))
Figure 4
Figure 4 shows the results of a sample run of a simple extension of the Reeve (1989) acceptance threshold model where we introduce a third type of recipient–the F1 hybrid. For simplicity, we assume that the hybrids have intermediate phenotype, and backcrosses have intermediate fitness. We then plot the expected acceptance error along threshold trait values under different hybrid frequencies.
#### new model: three-category model to assess effect of F1 hybrids
# There are three trait distributions
g2 = function(x) dnorm(x, mean = 0, sd = 0.2)
b2 = function(x) dnorm(x, mean = 0.5, sd = 0.2)
# hybrid has intermediate phenotype distribution
h2 = function(x) dnorm(x, mean = 0.25, sd = 0.2)
# set up results list. Here, we will evaluate results for 5 different
# frequencies of hybrids
results.list = list()
for (j in 1:5) {
t2 = seq(0, 1, length = 101)
Gt2 = list(0)
Bt2 = list(0)
Ht2 = list(0)
for (i in t2) Gt2[[which(t2 == i)]] = integrate(g2, lower = -1, upper = i)$value
Gt.value2 = unlist(Gt2)
for (i in t2) Bt2[[which(t2 == i)]] = integrate(b2, lower = -1, upper = i)$value
Bt.value2 = unlist(Bt2)
for (i in t2) Ht2[[which(t2 == i)]] = integrate(h2, lower = -1, upper = i)$value
Ht.value2 = unlist(Ht2)
# plot 20 points along the trait axis
times = 20
accept.error2 = vector(length = times)
# set up new parameters fitness of mating with conspecific (F)
Vd = 10
# proportion of population that are conspecifics ('desired recipient) ranges
# from 0.1 to 0.5
Pd = j/10
# proportion of population that are heterospecifics ('undesired recipient)
# ranges from 0.1 to 0.5
Pu = j/10
# proportion of population that are hybrids ranges from 0.8 to 0
Ph = 1 - (Pd + Pu)
Cs2 = 0.1
# fitness of mating with heterospecific (f) varies from 0 (hybrids are
# inviable or sterile) to Vd (no fitness reduction compared to mating with
# conspecific)
Vu = seq(0, Vd, length = times)
# assume backcrosses have intermediate fitness
Vh = (Vd + Vu)/2
for (i in 1:times) {
Wt2 = (Vd * Pd * Gt.value2 + Vu[i] * Pu * Bt.value2 + Vh[i] * Ph * Ht.value2 -
Cs2)/(Pd * Gt.value2 + Pu * Bt.value2 + Ph * Ht.value2)
accept.error2[i] = Bt.value2[which.max(Wt2)]
}
results.list[[j]] = accept.error2
}Now, we can set up the plot
#### Figure 4
require(RColorBrewer)
colors = brewer.pal(6, "BuPu") #set up color gradient
plot(Vu/Vd, results.list[[1]], type = "l", xlim = c(0, 1), col = colors[2],
lwd = 2, las = 1, pch = 19, ylab = "Acceptance Error", xlab = "Relative fitness of F1 Hybrids",
cex.lab = 1.4, cex.axis = 1.2)
for (k in 2:5) {
points(Vu/Vd, results.list[[k]], type = "l", xlim = c(0, 1), col = colors[k +
1], lwd = 2, las = 1, pch = 19)
}
# legend
prophybrids = c(0.8, 0.6, 0.4, 0.2, 0)
legend(x = 0, y = 0.9, lty = 1, lwd = 2, col = colors[2:6], legend = prophybrids,
title = "Proportion F1 Hybrids \n in Population", bty = "n")