#Czech Republic transition zone
#Barton-Turelli wave equation: hyperbolic tangent function (Barton, 1979; Turelli & Hoffmann, 1991; Schofield, 2002)
f <- function(x, sigma, t){0.5 * (1 - tanh((x - sigma * (1 - 2 * 0) * sqrt(0.98) * t / 2)/(2 * sigma / sqrt(0.98))))}
#our infection frequencies
y = c(1,1,1,1,0.85,0.93,0.92,0.71,0.45,0.38,0.17,0.22,0.06,0.1,0,0,0,0.04,0.04)
#distance from location 1 (calculated from online distance calculaors based on coordinates)
x = c(0,3.965,6.019,6.624,11.405,11.938,15.641,16.835,24.671,24.896,26.812,26.873,28.389,31.036,33.566,35.583,40.507,44.793,46.648)
#confidence interval calculation see bellow *after CI calculation, infection frequency is subtracted by the low and high CI because this function just adds or substracts from the dot (for graphing purposes)
CIl = c(0.0739728,0.185302,0.1322746,0.4096164,0.3044711,0.2686845,0.3047962,0.1506348,0.1515412,0.2414207,0.14913747,0.19185503,0.04692076,0.06530188,0,0,0,0.034913532,0.034913532)
CIh = c(0,0,0,0,0.1307933,0.0681932,0.0778924,0.1204692,0.1632682,0.3042224,0.31413775,0.38009357,0.11195993,0.12657766,0.07397279,0.1088812,0.2316358,0.102540973,0.102540973)
#curve fitting
# a = dispersal (sigma), b = generations (t)
plot(x,y)
curve(f(x,5,7), add = TRUE) #test the function plot with arbitrary values
#find a best fitting sigma and t
sigma.t.fit <- nls(y~f(x,sigma,t), start = list(sigma=1,t=1))
sigma.t.fit
#plot the best fitting sigma and t
plot(x,y, main="fit sigma & t", ylab="I")
curve(f(x,3.759,11.906), add = TRUE, lwd =2, col = "red", lty = 2 )
#install.packages("plotrix")
library(plotrix)
#creates plot and saves it, run it without the first (34) and last (40) lines just for plot creation
tiff("CZ_Transition_Zone3.tiff", width = 8, height = 6, units = "in", res = 300) #this is for saving higher quality plots only
plotCI(x,y, main="Transition Zone CZ", ylab="Infection Frequency", xlab="Distance(km)", err = "y", uiw = CIh, liw = CIl, col = "grey", xlim = c(0,80), slty = 2, cex.axis=1.5,cex.lab=1.5, cex.main=1.5)
points(x,y, lwd = "2")
curve(f(x,3.759,11.906), add = TRUE, lwd =2, col = "red", lty = 2 )
text(x= 60, y= 0.8, "R^2=0.9821045", cex=1.75)
text(x= 60, y= 0.7, "sigma of best fit=3.759", cex=1.5)
dev.off()
summary(sigma.t.fit)
#residual sum of squares and R value
RSS <- sum(residuals(sigma.t.fit)^2)
TSS <- sum((y - mean(y))^2)
R.square <- 1 - (RSS/TSS)
R.square
#Hungary transition zone
f <- function(x, sigma, t){0.5 * (1 - tanh((x - sigma * (1 - 2 * 0) * sqrt(0.98) * t / 2)/(2 * sigma / sqrt(0.98))))}
#our infection frequencies
y = c(1,1,1,1,1,1,0.91,0.2,0.11,0.02,0.02,0.2)
#distance from location 1 (calculated from online distance calculaors based on coordinates)
x = c(0,15.829,22.618,23.804,28.809,34.178,36.949,43.85,47.878,55.076,65.199,71.926)
#confidence interval calculation see bellow
CIl = c(0.3084971,0.3084971,0.3694166,0.4592581,0.185302,0.1277029,0.1,0.1117939,0.07375163,0.019424761,0.019449764,0.10642427)
CIh = c(0,0,0,0,0,0,0.09,0.1486655,0.1256967,0.100241582,0.095271826,0.1391476)
#curve fitting
# a = dispersal (sigma), b = generations (t)
plot(x,y)
curve(f(x,5,7), add = TRUE) #test the function plot
#find a best fitting sigma and t
sigma.t.fit <- nls(y~f(x,sigma,t), start = list(sigma=2,t=1))
sigma.t.fit
#plot the best fitting sigma and t
plot(x,y, main="fit sigma & t", ylab="I")
curve(f(x,1.961,42.615), add = TRUE, lwd =2, col = "red", lty = 2 )
#creates plot and saves it, run it without the first (83) and last (89) lines just for plot creation
tiff("HUN_Transition_Zone4.tiff", width = 8, height = 6, units = "in", res = 300) #this is for saving higher quality plots only
plotCI(x,y, main="Transition Zone HUN", ylab="Infection Frequency", xlab="Distance(km)", err = "y", uiw = CIh, liw = CIl, col = "grey", xlim = c(0,80), slty = 2, cex.axis=1.5,cex.lab=1.5, cex.main=1.5)
points(x,y,lwd=2)
curve(f(x,1.961,42.615), add = TRUE, lwd =2, col = "red", lty = 2 )
text(x= 60, y= 0.8, "R^2=0.9792222", cex=1.75)
text(x= 60, y= 0.7, "sigma of best fit=1.9612", cex=1.5)
dev.off()
summary(sigma.t.fit)
#residual sum of square and R value
RSS <- sum(residuals(sigma.t.fit)^2)
TSS <- sum((y - mean(y))^2)
R.square <- 1 - (RSS/TSS)
R.square
#==================================================================================================================
##confidence intervals were calculated with the following code
#Clopper-Pearson exact CI
#install.packages("PropCIs")
library(PropCIs)
#example usage
exactci(34, 48, 0.95)
#34 infected flies
#48 total flies
#95% confidence interval
#==================================================================================================================
#model
#p(x,t) = 0.5 * (1 - tanh((x - sigma * (1 - 2 * p) * sqrt(Sh) * t / 2)/(2 * sigma / sqrt(Sh))))
x <- 0:80 #one dimensional distance (km) from point 1 (x axis in graph)
sigma <- 6 #dispersal/migration (best fit is calculated from non-linear best fit)
Sh <- 0.98 #level of CI determined from Boller et al., 1976
Sf <- 0.00
p <- Sf/Sh #Sf/Sh ;we assume Sf to be zero
t <- 4 #generations; this parameters switches the graph left and right only; one should find a best fit (not important parameter as we are looking at only one generation)
plot(0.5 * (1 - tanh((x - sigma * (1 - 2 * p) * sqrt(Sh) * t / 2)/(2 * sigma / sqrt(Sh)))), type = "l", ylim = c(0,1), xlim = c(0,80), lwd = 3, ylab = "I", xlab = "X" )
#==================================================================================================================
#References
#
#Barton NH. 1979 The dynamics of hybrid zones. Heredity, 43, 341-359. (doi: 10.1038/hdy.1979.87)
#Turelli M, Hoffmann AA. 1991 Rapid spread of an inherited incompatibility factor in California Drosophila. Nature 353, 440–442. (doi:10.1038/353440a0)
#Schofield P, 2002. Spatially explicit models of Turelli-Hoffmann Wolbachia invasive wave fronts. J. Theor. Biol., 215, 121-131. (doi:10.1006/jtbi.2001.2493)