sim.metapopgen.dioecious.multilocus <- function(input.type,demographic.data,N1_M,N1_F,sigma_M,sigma_F,phi_M,phi_F,List_gene,List_mu,r=0.5,delta,recr.dd="settlers",kappa0,T_max,save.res=F,save.res.T=seq(1,T_max),verbose=F) {
##########################################################################
# Initial definitions
##########################################################################
if (input.type=="data.frame") {
print("Input type = data.frame")
a <- metapopgen.input.convert.dioecious(demographic.data)
N1_M <- a[[1]]
N1_F <- a[[2]]
sigma_M <- a[[3]]
sigma_F <- a[[4]]
phi_M <- a[[5]]
phi_F <- a[[6]]
rm(a)
} else {
if (input.type == "array") {
print("Input type = array")
} else {
stop("Unknown value for argument input.type. It must be either data.frame, array or txt")
}
}
# Define basic variables
m <- dim(N1_M)[1] # Number of genotypes
l <- length(Proba[,1]) # Number of alleles ### What is l in the multilocus version? Number of GAMETES? 07/09/18
n <- dim(N1_M)[2] # Number of demes
z <- dim(N1_M)[3] # Number of age-classes
# Here there should be the warning about one age-class and recr.dd=adults
# Define usefull matrix:
index<-genotype.index.multilocus(List_gene)
print ("creating probability matrix")
Proba<-create.probability.matrix(index,List_gene,r,List_mu)
print("probability matrix created")
##########################################################################
# Check if input data are time-dependent or not
##########################################################################
# Male Survival
if (is.na(dim(sigma_M)[4])) {
sigma_M <- array(rep(sigma_M,T_max),c(m,n,z,T_max))
}
# Female Survival
if (is.na(dim(sigma_F)[4])) {
sigma_F <- array(rep(sigma_F,T_max),c(m,n,z,T_max))
}
# Female fecundity
if (is.na(dim(phi_F)[4])) {
phi_F <- array(rep(phi_F,T_max),c(m,n,z,T_max))
}
# Male fecundity
if (is.na(dim(phi_M)[4])) {
phi_M <- array(rep(phi_M,T_max),c(m,n,z,T_max))
}
# Dispersal
if (is.na(dim(delta)[3])) {
delta <- array(rep(delta,T_max),c(n,n,T_max))
}
# Carrying capacity
if (is.vector(kappa0)) {
kappa0 <- array(rep(kappa0,T_max),c(n,T_max))
}
##########################################################################
# Initialize state variables
##########################################################################
print("Initializing variables...")
if (save.res){
N_M <- N1_M
N_F <- N1_F
rm(N1_M,N1_F)
} else {
N_M <- array(NA,dim=c(m,n,z,T_max))
N_M[,,,1] <- N1_M
N_F <- array(NA,dim=c(m,n,z,T_max))
N_F[,,,1] <- N1_F
rm(N1_M,N1_F)
L_M <- array(NA,dim=c(m,n,T_max))
L_F <- array(NA,dim=c(m,n,T_max))
S_M <- array(0,dim=c(m,n,T_max))
S_F <- array(0,dim=c(m,n,T_max))
}
##########################################################################
# Define functions
##########################################################################
# Survival
surv <-
function(sigma,N) {
rbinom(1,N,sigma)
}
#Type Gamete
Type_gamete<-function(fec,Proba){
l=length(Proba[,1])
m=length(Proba[1,])
Gamete<-array(0,dim=c(l,m))
for (k in 1:m){
err<- try (as.vector(rmultinom(1,fec[k],Proba[,k])),silent=T) #meiose
if (class(err)=="try-error") {
prob.mvr <- fec[k] * Proba[,k] # Vector of means of the multivariate normal distribution
var.mvr <- fec[k]* Proba[,k] * (1-Proba[,k]) # Vector of variances of the multivariate normal distribution
sigma.mvr <- diag(var.mvr, l) # Variance-covariance matrix of the multivariate normal distribution
for (i.mvr in 1 : l){
for (j.mvr in 1 : l) {
if (i.mvr == j.mvr) next
sigma.mvr[i.mvr,j.mvr] <- -fec[k] * Proba[i.mvr,k] * Proba[j.mvr,k]
}
}
err<-try(as.vector(round(loinorm(1,prob.mvr,sigma.mvr))),silent=T)
if (class(err)=="try-error"){
err<- fec[k]*Proba[,k]
}
}
Gamete[,k]<-err
}
type<-array(0,l)
for (j in 1:l){
type[j]<-sum(Gamete[j,])
}
return(type)
}
#Reproduction
repr <-
function(Nprime_M,Nprime_F,phi_F,phi_M,l,m,Proba) {
# Adjust dimensions if there is only one age-class and/or one deme
#if (length(dim(phi_F))==2) dim(phi_F)[3]<-1
#if (length(dim(phi_M))==2) dim(phi_M)[3]<-1
#dim(phi_F) <- c(m,n,z)
#dim(phi_M) <- c(m,n,z)
# n and z are not passed to the function? But it works...
#if (verbose) print("Calculates total number of female gametes")
fecx<- array(0,dim=c(m,z)) # Number of female gametes produced by all the individuals of each genotype in each age class
fec<- array(0,dim=m) # Number of female gametes produced by all the individuals of each genotype
for (k in 1 : m) {
for (x in 1 : z) {
fecx[k,x] <- sum(as.numeric(rpois(Nprime_F[k,x],phi_F[k,x]))) # This is the contribution of variation in reproductive success among individuals to genetic drift
}
fec[k] <- sum(fecx[k,])
}
#if (verbose) print("Calculate number of gametes for each allele")
G_F <-Type_gamete(fec,Proba)
#if (verbose) print("Calculates total number of male gametes")
fecx <- array(0,dim=c(m,z)) # Number of male gametes produced by all the individuals of each genotype in each age class
fec <- array(0,dim=m) # Number of male gametes produced by all the individuals of each genotype
for (k in 1 : m) {
for (x in 1 : z) {
fecx[k,x] <- sum(as.numeric(rpois(Nprime_M[k,x],phi_M[k,x]))) # This is the contribution of variation in reproductive success among individuals to genetic drift
}
fec[k] <- sum(fecx[k,])
}
#if (verbose) print("Calculate number of gametes for each allele")
G_M <- Type_gamete(fec,Proba)
#if (verbose) print("Union of gametes to form zygotes")
if (sum(G_F) <= sum(G_M)) {
mat_geno<- array(0,dim=c(l,l))
Gprime_M<- G_M
for (j in 1 : l) {
in_dist<- Gprime_M
odds<- array(1,dim=l)
ndraws<- G_F[j]
# If the mutlivariate hypergeometric does not work, use the multinomial
# If the multinomial does not work, use the multivariate normal
err1<- try(rMWNCHypergeo(1,in_dist,ndraws,odds),silent=T)
if (class(err1)=="try-error") {
err2<- try(as.numeric(rmultinom(1,ndraws,in_dist)),silent=T)
if (class(err2)=="try-error") {
# Use multivariate normal
prob <- in_dist/sum(in_dist)
mu.mvr <- ndraws * prob # Vector of means of the multivariate normal distribution
var.mvr <- ndraws * prob * (1-prob) # Vector of variances of the multivariate normal distribution
sigma.mvr <- diag(var.mvr, l) # Variance-covariance matrix of the multivariate normal distribution
for (i.mvr in 1 : l){
for (j.mvr in 1 : l) {
if (i.mvr == j.mvr) next
sigma.mvr[i.mvr,j.mvr] <- -ndraws * prob[i.mvr] * prob[j.mvr]
}
}
extr <- as.vector(round(loinorm(1,mu.mvr,sigma.mvr)))
} else {
extr<- err2 # Use multinomial
}
} else {
extr <- err1
}
mat_geno[j,]<- extr
Gprime_M<- Gprime_M - extr
}
mat_geno_l<- mat_geno
mat_geno_l[upper.tri(mat_geno_l)]<- 0
mat_geno_u<- mat_geno
mat_geno_u[lower.tri(mat_geno_u,diag=T)]<- 0
mat_geno_f<- mat_geno_l + t(mat_geno_u)
L <- mat_geno_f[lower.tri(mat_geno_f,diag=T)]
} else {
mat_geno <- array(0,dim=c(l,l))
Gprime_F <- G_F
for (j in 1 : l) {
in_dist <- Gprime_F
odds <- array(1,dim=l)
ndraws <- G_M[j]
# If the mutlivariate hypergeometric does not work, use the multinomial
# If the multinomial does not work, use the multivariate normal
err1 <- try(rMWNCHypergeo(1,in_dist,ndraws,odds),silent=T)
if (class(err1)=="try-error") {
err2 <- try(as.numeric(rmultinom(1,ndraws,in_dist)),silent=T)
if (class(err2)=="try-error") {
# Use multivariate normal
prob <- in_dist/sum(in_dist)
mu.mvr <- ndraws * prob # Vector of means of the multivariate normal distribution
var.mvr <- ndraws * prob * (1-prob) # Vector of variances of the multivariate normal distribution
sigma.mvr <- diag(var.mvr, l) # Variance-covariance matrix of the multivariate normal distribution
for (i.mvr in 1 : l){
for (j.mvr in 1 : l) {
if (i.mvr == j.mvr) next
sigma.mvr[i.mvr,j.mvr] <- -ndraws * prob[i.mvr] * prob[j.mvr]
}
}
extr <- as.vector(round(loinorm(1,mu.mvr,sigma.mvr)))
} else {
extr<- err2 # Use multinomial
}
} else {
extr <- err1 # Use multivariate hypergeometric
}
mat_geno[j,]<- extr
Gprime_F<- Gprime_F - extr
}
mat_geno_l<- mat_geno
mat_geno_l[upper.tri(mat_geno_l)]<- 0
mat_geno_u<- mat_geno
mat_geno_u[lower.tri(mat_geno_u,diag=T)]<- 0
mat_geno_f<- mat_geno_l + t(mat_geno_u)
L<- mat_geno_f[lower.tri(mat_geno_f,diag=T)]
}
#if (verbose) print("Determine sex")
# Sex differentiation is independent of genotype, with equal proportion of male and females produced
LL <- array(NA,dim=c(m,2))
for (i in 1 : m){
LL[i,1] <- rbinom(1,L[i],0.5)
LL[i,2] <- L[i] - LL[i,1]
}
return(LL)
}
# Dispersal
disp <-
function(L,delta){
delta_lost <- max(0,1 - sum(delta)) # The maximum function is needed to avoid errors due to precision
delta_add <- c(delta,delta_lost)
S <- rmultinom(1,L,delta_add)
return(S)
}
# Recruitment:
# S_M[,i],S_F[,i],Nprime_M[,i,],Nprime_F[,i,],m,kappa0[i,t]
# S, S_othersex. Number of settlers of all genotypes. Dimension: m
# N_M, N_F. Number of adults of all genotyps and age-classes. Dimensions: n*z
# m. Number of genotypes
# kappa0. Carrying capacity. Scalar.
recr <-
function(S,S_othersex,N_M,N_F,m,kappa0,recr.dd) {
switch(recr.dd,
# Dependence on settler density
settlers = {
Ntot <- sum(S+S_othersex)
sigma <- settler.survival(Ntot,kappa0)
},
# Dependence on adult density
adults = {
Ntot <- sum(N_M[,1:(z-1)] + N_F[,1:(z-1)])
Stot <- sum(S+S_othersex)
Recr <- kappa0 - Ntot
if (Recr <= 0){
sigma <- 0
} else {
sigma <- Recr / Stot
}
if (sigma > 1) sigma <- 1
},
# No-match: error
stop("Unknown value for argument recr.dd. Valid values: 'settlers', 'adults'")
)
# Use recruitment Probability to calculate the number of recruits
R <- array(0,dim=m)
for (k in 1 : m) {
R[k] <- rbinom(1,S[k],sigma)
}
return(R)
}
##########################################################################
# Simulate metapopulation genetics
##########################################################################
if (save.res){
dir.res.name <- paste(getwd(),format(Sys.time(), "%Y-%b-%d-%H.%M.%S"),sep="/")
dir.create(dir.res.name)
if (1 %in% save.res.T) {
file.name <- "N1.RData"
save(N_M,N_F,file=paste(dir.res.name,file.name,sep="/"))
}
}
print("Running simulation...")
for (t in 1 : (T_max-1)) {
print(t)
m <- dim(N_M)[1] # Number of genotypes
l <- length(Proba [,1]) # Nuber of alleles
n <- dim(N_M)[2] # Number of demes
z <- dim(N_M)[3] # Number of age-classes
# At each time-step, redefine variable Nprime
# If save.res, redefine also larval and settlers numbers
if (save.res) {
Nprime_M <- array(NA,dim=c(m,n,z))
Nprime_F <- array(NA,dim=c(m,n,z))
L_M <- array(NA,dim=c(m,n))
L_F <- array(NA,dim=c(m,n))
S_M <- array(0,dim=c(m,n))
S_F <- array(0,dim=c(m,n))
} else {
Nprime_M <- array(NA,dim=c(m,n,z))
Nprime_F <- array(NA,dim=c(m,n,z))
}
### Survival
# If there is only one age-class, we must force the third dimension. What if only one year?
#if (length(dim(sigma_M))==2) dim(sigma_M)[3] <- 1
# If there is only one age-class, we must force the third dimension. What if only one year?
#if (length(dim(sigma_F))==2) dim(sigma_F)[3] <- 1
if (verbose) print("Apply survival function")
for (i in 1 : n) {
for (x in 1 : z) {
for (k in 1 : m) {
if (save.res){
Nprime_M[k,i,x] = surv(sigma_M[k,i,x,t],N_M[k,i,x])
Nprime_F[k,i,x] = surv(sigma_F[k,i,x,t],N_F[k,i,x])
} else {
Nprime_M[k,i,x] = surv(sigma_M[k,i,x,t],N_M[k,i,x,t])
Nprime_F[k,i,x] = surv(sigma_F[k,i,x,t],N_F[k,i,x,t])
}
if(is.na(Nprime_M[k,i,x])){
Nprime_M[k,i,x]<-0
cat("NA adjustment in deme",i,"age",x,"genotype",k,"males \n")
}
if(is.na(Nprime_F[k,i,x])){
Nprime_F[k,i,x]<-0
cat("NA adjustment in deme",i,"age",x,"genotype",k,"females \n")
}
}
}
}
if (verbose) print("Apply reproduction function")
# If there is only one age-class, we must force the third dimension
if (length(dim(phi_F))==2) dim(phi_F)[3] <- 1
if (length(dim(phi_M))==2) dim(phi_M)[3] <- 1
for (i in 1 : n) {
if (save.res) {
if (sum(Nprime_M[,i,],Nprime_F[,i,])==0) { # To save computing time. In the older version it was: if (sum(Nprime_M[,i,])==0 | sum(Nprime_F[,i,])==0)
L_M[,i] = 0
L_F[,i] = 0
next
} else {
LL <- repr(Nprime_M[,i,],Nprime_F[,i,],phi_F[,i,,t],phi_M[,i,,t],l,m,Proba)
L_M[,i] <- LL[,1]
L_F[,i] <- LL[,2]
rm(LL)
}
} else {
if (sum(Nprime_M[,i,],Nprime_F[,i,])==0) { # To save computing time. In the older version it was: if (sum(Nprime_M[,i,])==0 | sum(Nprime_F[,i,])==0)
L_M[,i,t] = 0
L_F[,i,t] = 0
next
} else {
LL <- repr(Nprime_M[,i,],Nprime_F[,i,],phi_F[,i,,t],phi_M[,i,,t],l,m,Proba)
L_M[,i,t] <- LL[,1]
L_F[,i,t] <- LL[,2]
}
}
}
if (verbose) print("Apply dispersal function")
for (i in 1 : n) {
for (k in 1 : m) {
if (save.res) {
y = disp(L_M[k,i],delta[,i,t])
S_M[k,] <- S_M[k,] + y[1:n]
y = disp(L_F[k,i],delta[,i,t])
S_F[k,] <- S_F[k,] + y[1:n]
} else {
y = disp(L_M[k,i,t],delta[,i,t])
S_M[k,,t] <- S_M[k,,t] + y[1:n]
y = disp(L_F[k,i,t],delta[,i,t])
S_F[k,,t] <- S_F[k,,t] + y[1:n]
}
}
}
if (verbose) print("Apply recruitment function")
for (i in 1 : n) {
if (save.res) {
N_M[,i,1] <- recr(S_M[,i],S_F[,i],Nprime_M[,i,],Nprime_F[,i,],m,kappa0[i,t],recr.dd) # Pass the abundance of other classes too
N_F[,i,1] <- recr(S_F[,i],S_M[,i],Nprime_M[,i,],Nprime_F[,i,],m,kappa0[i,t],recr.dd)
} else {
N_M[,i,1,t+1] <- recr(S_M[,i,t],S_F[,i,t],Nprime_M[,i,],Nprime_F[,i,],m,kappa0[i,t],recr.dd)
N_F[,i,1,t+1] <- recr(S_F[,i,t],S_M[,i,t],Nprime_M[,i,],Nprime_F[,i,],m,kappa0[i,t],recr.dd)
}
}
if (verbose) print("Calculates N at t+1")
for (i in 1 : n){
for (x in 1 : z) {
if (x == 1) next
for (k in 1 : m) {
if (save.res) {
N_M[k,i,x] <- Nprime_M[k,i,x-1]
N_F[k,i,x] <- Nprime_F[k,i,x-1]
} else {
N_M[k,i,x,t+1] <- Nprime_M[k,i,x-1]
N_F[k,i,x,t+1] <- Nprime_F[k,i,x-1]
}
}
}
}
# Save results if save.res=T
if (save.res){
if ((t+1) %in% save.res.T) {
file.name <- paste("N",(t+1),".RData",sep="")
save(N_M,N_F,file=paste(dir.res.name,file.name,sep="/"))
}
}
}
print(T_max)
print("...done")
if (save.res==F) return(list(N_M=N_M,N_F=N_F))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.