---
title: SIT resistance dynamics with proportional release 
# author: "Dick Gomulkiewicz and Jim Bull"
# date: "10 March 2022"
output:
  pdf_document: default
  html_document: default
  word_document: default
---

In the haploid sexual model, the population dynamics of the discriminating genotype is described by 
\begin{equation*}
    N^\prime_A = R_A \cdot N_A \label{disc.change}
\end{equation*}
where
\begin{equation*}
    R_A = \frac{b}{2} \left[ \frac{\pi}{\pi + d} (1+p) + \frac{1}{1 + d} (1-p)\right] .
\end{equation*}

```{r, echo=F}
b <- 2
pi <-2.5
```

Plot with $\pi = `r pi`$, $b = `r b`$:

```{r, echo=F}
#Discriminator finite growth rate
R_A <- function(p,b,pi,d){
  return(b/2*(pi/(pi + d)*(1+p) + 1/(1 + d)*(1-p)))
}

#Wild type finite growth rate
R_a <- function(q,b,pi,d){
  return(b/2*(1/(1 + d)*(1+q) + pi/(pi + d)*(1-q)))
}

ps <- seq(0, 1, by = 0.01)
d.1 <- 1.02
d.2 <- 2
d.3 <- 3

growth.rates.1 <- sapply(ps, function(p) R_A(p,b=b, pi=pi, d=d.1))
growth.rates.2 <- sapply(ps, function(p) R_A(p,b=b, pi=pi, d=d.2))
growth.rates.3 <- sapply(ps, function(p) R_A(p,b=b, pi=pi, d=d.3))
y.low <- min(growth.rates.1,growth.rates.2,growth.rates.3)
y.hi <- max(growth.rates.1,growth.rates.2,growth.rates.3)

# plot(ps,growth.rates.1,type="l",xlab="p",ylab=expression(R[A]), lty = 1, ylim = c(y.low,y.hi))
# lines(ps,growth.rates.2, col=2)
# lines(ps,growth.rates.3, col=3, type = "l")
# abline(h=1, lty = 3)
# legend("bottomright",c(paste(d.1),paste(d.2),paste(d.3)),title = "d", lty=1, col=1:3,bty="n",cex=0.75)
```


```{r, echo=F}

#Trajectories

last.gen <- 25 

#Discriminator frequencies for the 3 cases
p.1 <- matrix(NA,nrow=last.gen,ncol=1)
p.2 <- matrix(NA,nrow=last.gen,ncol=1)
p.3 <- matrix(NA,nrow=last.gen,ncol=1)

#initial densities
NA.0 <- 1 #discriminator
Na.0 <- 100 #wild type
N.A <- matrix(NA.0,nrow=last.gen+1,ncol=1)
N.a <- matrix(Na.0,nrow=last.gen+1,ncol=1)

####### CASE 1 ########
for(t in 1:last.gen){
  p.1[t] <- N.A[t]/(N.A[t]+N.a[t])
  q <- 1-p.1[t]
  N.A[t+1] <- R_A(p=p.1[t],b=b,pi=pi,d=d.1)*N.A[t]
  N.a[t+1] <- R_a(q=q,b=b,pi=pi,d=d.1)*N.a[t]
  }

####### CASE 2 ########
N.A[1] <- NA.0
N.a[1] <- Na.0
for(t in 1:last.gen){
  p.2[t] <- N.A[t]/(N.A[t]+N.a[t])
  q <- 1-p.2[t]
  N.A[t+1] <- R_A(p=p.2[t],b=b,pi=pi,d=d.2)*N.A[t]
  N.a[t+1] <- R_a(q=q,b=b,pi=pi,d=d.2)*N.a[t]
  }

####### CASE 3 ########
N.A[1] <- NA.0
N.a[1] <- Na.0
for(t in 1:last.gen){
  p.3[t] <- N.A[t]/(N.A[t]+N.a[t])
  q <- 1-p.3[t]
  N.A[t+1] <- R_A(p=p.3[t],b=b,pi=pi,d=d.3)*N.A[t]
  N.a[t+1] <- R_a(q=q,b=b,pi=pi,d=d.3)*N.a[t]
  }

####### DISCRIMINATOR REALIZED FINITE GROWTH RATES
RA.1 <- sapply(p.1, function(p) R_A(p,b=b, pi=pi, d=d.1))
RA.2 <- sapply(p.2, function(p) R_A(p,b=b, pi=pi, d=d.2))
RA.3 <- sapply(p.3, function(p) R_A(p,b=b, pi=pi, d=d.3))

##### Constant Release #######

#Numerically stable version of Discriminator finite growth rate
R_A.C <- function(N,p,b,pi,S){
  return(b/2*(pi*N/(pi*N + S)*(1+p) + N/(N + S)*(1-p)))
}

#Numerically stable version of Wild type finite growth rate
R_a.C <- function(N,q,b,pi,S){
  return(b/2*(N/(N + S)*(1+q) + pi*N/(pi*N + S)*(1-q)))
}

d.C <- 1.02
p.C <- matrix(NA,nrow=last.gen,ncol=1)
RA.C <-matrix(NA,nrow=last.gen,ncol=1)
S <- d.C*(NA.0 + Na.0)
N.A[1] <- NA.0
N.a[1] <- Na.0

for(t in 1:last.gen){
  N <- N.A[t]+N.a[t]
  p.C[t] <- N.A[t]/N
  q <- 1-p.C[t]
  RA.C[t]<- R_A.C(N=N,p=p.C[t],b=b,pi=pi,S=S)
  N.A[t+1] <- RA.C[t]*N.A[t]
  N.a[t+1] <- R_a.C(N=N, q=q,b=b,pi=pi,S=S)*N.a[t]
  }

####### PLOT  ######
y.low <- 0
plot(ps,growth.rates.1,type="l",xlab=expression(italic(p)),ylab=expression(italic(R[A])), lty = 1, ylim = c(y.low,y.hi),xlim = c(0,1),lwd=2)
lines(ps,growth.rates.2, col=2,lwd=2)
lines(ps,growth.rates.3, col=3,lwd=2)
abline(h=1, lty = 3)
points(p.1,RA.1,pch=16)
points(p.2,RA.2,pch=16, col=2)
points(p.3,RA.3,pch=16, col=3)
lines(p.C,RA.C)
points(p.C,RA.C)
legend(x=0.78,y=0.75,c(paste(d.1),paste(d.2),paste(d.3),"constant S"),title = expression(italic(d)), lty=1, col=1:3, pch = c(16, 16, 16, 1),bty="y",cex=0.9)

#save figure as jpeg file
# jpeg('discriminator.growth.rates.jpg')
# plot(ps,growth.rates.1,type="l",xlab="p",ylab=expression(R[A]), lty = 1, ylim = c(y.low,y.hi),lwd=2)
# lines(ps,growth.rates.2, col=2,lwd=2)
# lines(ps,growth.rates.3, col=3,lwd=2)
# abline(h=1, lty = 3)
# points(p.1,RA.1,pch=16)
# points(p.2,RA.2,pch=16, col=2)
# points(p.3,RA.3,pch=16, col=3)
# # points(p.C,RA.C,type = "l")
# legend("bottomright",c(paste(d.1),paste(d.2),paste(d.3)),title = "d", lty=1, col=1:3,bty="n",cex=0.75)
# dev.off()

#save figure as pdf file
pdf('prop.v.const.pdf')
y.low <- 0
plot(ps,growth.rates.1,type="l",xlab="p",ylab=expression(italic(R[A])), lty = 1, ylim = c(y.low,y.hi),xlim = c(0,1),lwd=2)
lines(ps,growth.rates.2, col=2,lwd=2)
lines(ps,growth.rates.3, col=3,lwd=2)
abline(h=1, lty = 3)
points(p.1,RA.1,pch=16)
points(p.2,RA.2,pch=16, col=2)
points(p.3,RA.3,pch=16, col=3)
lines(p.C,RA.C)
points(p.C,RA.C)
legend(x=0.78,y=0.75,c(paste(d.1),paste(d.2),paste(d.3),"constant S"),title = expression(italic(d)), lty=1, col=1:3, pch = c(16, 16, 16, 1),bty="y",cex=0.9)
dev.off()

```

