Ve = 0.66
Vm = 0.12

library("deSolve")
library("reshape2")
library("ggplot2"); theme_set(theme_bw())

## Loading required package: methods

alpha <- c(seq(from = 2, to = 7, by = 1))
n.alpha <- length(alpha)
n.stage <- 10
n.I <- 2 + n.stage
v3 = 10^alpha

parameters <- c(
)

c = 1.25
lam1 = 4
lam3 = 4/3
Bmax = 0.317
K = 13938
a = 1.02
Dmax = 25.4
h = 0.41
D50 = 3058
Beta1unadj = 2.76
Beta3unadj = 0.76

##this allows you to change initial value
inimean = 3.5
inisd = 0.05
iniI = 0.001

yini <- c(
    S = 1 - iniI,
    i1 = iniI * dnorm(alpha, mean = inimean, sd = inisd)/sum(dnorm(alpha, mean = inimean, sd = inisd)),
    i2 = matrix(c(rep(0, n.alpha * n.stage)),
                ncol = n.alpha,
                nrow = n.stage,
                byrow = TRUE),
    i3 = rep(0, n.alpha)
)



norm2 <- function(x){dnorm(alpha, mean = x, sd = Vm)}

p0 <- matrix(c(sapply(alpha, norm2)),
             nrow = n.alpha,
             byrow = FALSE)

p <- p0/rowSums(p0)

pheno <- c(seq(from = 0, to = 9, by = 0.5))
n.pheno = length(pheno)

norm3 <- function(x){dnorm(pheno, mean = x, sd = Ve)}

p.gm0 <- matrix(c(sapply(alpha, norm3)),
               nrow = n.pheno,
               ncol = n.alpha)
p.gm <- sweep(p.gm0, 2, colSums(p.gm0), "/")

v3 = 10^pheno

Beta1 = Beta1unadj * c/ (Beta1unadj + c + lam1)
lam2 = (D50^h + v3^h)/(Dmax * D50^h)
Beta2unadj = Bmax * v3^a / (v3^a + K^a)
Beta2 = Beta2unadj * c/ (Beta2unadj + c + lam2)

Beta3 = Beta3unadj * c/ (Beta3unadj + c + lam3)

g2 <- function(t,yini,parameters) {
  with(as.list(c(yini,parameters)), {

      Imat <- matrix(yini[-1],ncol=n.I)

      infrate <- (S*(Beta1 * Imat[,1] +
                        Beta2 %*% sweep(p.gm, 2, rowSums(Imat[,2:(n.I-1)]), "*") +Beta3 * Imat[,n.I])) %*% p
      dS = sum(lam3 * Imat[,n.I]) - sum(infrate)

      di1 = infrate - lam1 * Imat[,1]

      dI2mat = matrix(NA,nrow=n.alpha,ncol=n.stage)

      dI2mat[,1] = lam1 * Imat[,1] - n.stage * lam2 %*% sweep(p.gm, 2, Imat[,2], "*")

      for(i in 2:n.stage){
        dI2mat[,i] = n.stage * lam2 %*% sweep(p.gm , 2 , (Imat[,i]-Imat[,(i+1)]) , "*")
      }

      di3 = n.stage * lam2 %*% sweep(p.gm, 2, Imat[,(n.I - 1)], "*") - lam3 * Imat[,n.I]

      list(c(dS, di1, dI2mat, di3))
  })
} 

tvec <- seq(from = 0, to = 400, by = 0.01)
r2 <- ode(y= yini,times=tvec,func=g2,parms=parameters)
mm2 <- melt(cbind(as.data.frame(r2)), id.var = "time")
ggplot(mm2, aes(time, value, colour = variable)) + geom_line()

Here's the evolution! Mean alpha:

sum_aI2 <- rowSums(sweep(r2[,-(1:2)],2,rep(alpha,n.I),"*"))
tot_I2 <- rowSums(r2[,-(1:2)])
plot(sum_aI2/tot_I2,type="l")

save.image("HIV_Mutation_long.RData")