#---------------------------------------------- # Temperature effect #---------------------------------------------- # temperature effect in GR3D # from Rosso temperatureEffect = function(tempWater, Tmin, Topt, Tmax, effectFunction = Rosso){ return(effectFunction(tempWater, Tmin, Topt, Tmax)) } modifiedLehman = function(temperature, Tmin, Topt, Tmax){ # Svirezhev, Yu. M., Krysanova, V. P. & Voinov, A. A. (1984) Mathematical modelling of a fish pond ecosystem. # Ecological Modelling 21, 315–337. # rate(Tmin or Tmax) = exp(-4.6) = 1 % rate = rep(0, length(temperature)) rate[temperature < Topt] = exp(-4.6 * ((Topt - temperature[temperature < Topt]) / (Topt - Tmin))^4) rate[temperature >= Topt] = exp(-4.6 * ((temperature[temperature >= Topt] - Topt) / (Tmax - Topt))^4) return(rate) } Rosso = function(tempWater, Tmin, Topt, Tmax){ response = (tempWater - Tmin) * (tempWater - Tmax) / ((tempWater - Tmin) * (tempWater - Tmax) - (tempWater - Topt)^2) response[tempWater <= Tmin | tempWater >= Tmax] = 0 return(response) } thermalRange = function(pct = 0.8, Tmin, Topt, Tmax,...){ lower = uniroot(function(x) temperatureEffect(x, Tmin, Topt, Tmax,...) - pct, interval = c(Tmin, Topt))\$root upper = uniroot(function(x) temperatureEffect(x, Tmin, Topt, Tmax,...) - pct, interval = c(Topt, Tmax))\$root return(c(lower = lower,upper = upper)) } # thermalRange(Tmin = 3, Topt = 17, Tmax = 27, effectFunction = modifiedLehman ) #---------------------------------------------- # Growth simulation #--------------------------------------------- # von Bertalaffy : L = f(age) vonBertalanffyGrowth = function(age, L0, Linf, K){ t0 = log(1 - L0 / Linf) / K return(Linf * (1 - exp(-K * (age - t0)))) } # von Bertalanffy inverse : age = f(L) # with L0 rather the usual t0 parameter vonBertalanffyInverse = function(L, L0, Linf, K){ t0 = log(1 - L0/Linf)/K return(t0 - log(1 - L/Linf)/K) } # von Bertalanffy increment # pas cohérent avec la temperature effet sur le coeff de croissance mais ca marche vonBertalanffyIncrement = function(nStep, L0, Linf, K, timeStepDuration, sigma, withTempEffect=FALSE, TrefAtSea = c(9.876946, 13.489854, 15.891487, 11.554104) ){ tempEffect = temperatureEffect( TrefAtSea , 3, 17, 26) L = matrix(nrow = nStep + 1) L[1] = L0 for (i in 1:nStep) { mu = log((Linf - L[i]) * (1 - exp(-K * timeStepDuration))) - sigma * sigma / 2 increment = exp(rnorm(1, mu, sigma)) if (withTempEffect) { increment = increment * tempEffect[((i - 1) %% 4) + 1] } L[i + 1] = L[i] + increment } return(L) } vonBertalanffyWithNextIncrement = function(L, Linf, K, timeStepDuration, sigma, tempEffect ){ if (sigma == 0) { mu = log((Linf - L) * (1 - exp(-K * timeStepDuration))) increment = exp(mu) } else { mu = log((Linf - L) * (1 - exp(-K * timeStepDuration))) - (sigma * sigma) / 2 increment = exp(rnorm(length(mu), mu, sigma)) } L = L + increment * tempEffect return(L) } # generate a cohort of length trajectories computeMultipleLengthTrajectories = function(temperaturePattern, Nind = 10, growPar, ...){ ages = temperaturePattern\$age res <- expand_grid(ind = seq.int(Nind), age = ages) %>% inner_join(temperaturePattern, by = 'age') %>% arrange(age, ind) %>% mutate(temperatureEffect = temperatureEffect(temperature, growPar\$tempMinGrow, growPar\$tempOptGrow, growPar\$tempMaxGrow, ...), L = if_else(age == 0, growPar\$lengthAtHatching, 0)) for (i in 2:length(ages)) { previousAge = ages[i - 1] currentAge = ages[i] tempEffect = res %>% filter(age == currentAge) %>% select(temperatureEffect) %>% unlist(use.names = FALSE) previousL <- res %>% filter(age == previousAge) %>% select(L) %>% unlist(use.names = FALSE) currentL <- vonBertalanffyWithNextIncrement(L = previousL, Linf = growPar\$Linf, K = growPar\$kOpt, timeStepDuration = currentAge - previousAge, sigma = growPar\$sigmaDeltaLVonBert, tempEffect = tempEffect) res = res %>% mutate(L =replace(L, age == ages[i], currentL) ) } return(res) } computeMultipleLengthTrajectoriesWithRandomSeed = function(temperaturePattern, Nind = 10, growPar, RNGseed =1, ...){ set.seed(RNGseed) ages = temperaturePattern %>% distinct(age) %>% unlist(use.names = FALSE) res <- expand_grid(ind = seq.int(Nind), age = ages) %>% mutate(random = rnorm(Nind * length(ages))) %>% inner_join(temperaturePattern, by = 'age') %>% arrange(age, ind) %>% mutate(temperatureEffect = temperatureEffect(temperature, growPar\$tempMinGrow, growPar\$tempOptGrow, growPar\$tempMaxGrow, ...), L = if_else(age == 0, growPar\$lengthAtHatching, 0)) for (i in 2:length(ages)) { previousAge = ages[i - 1] currentAge = ages[i] tempEffect <- res %>% filter(age == currentAge) %>% select(temperatureEffect) %>% unlist(use.names = FALSE) previousL <- res %>% filter(age == previousAge) %>% select(L) %>% unlist(use.names = FALSE) rnd <- res %>% filter(age == currentAge) %>% select(random) %>% unlist(use.names = FALSE) currentL <- vonBertalanffyWithRandomVector(L = previousL, Linf = growPar\$Linf, K = growPar\$kOpt, timeStepDuration = currentAge - previousAge, randomVector = rnd, sigma = growPar\$sigmaDeltaLVonBert, tempEffect = tempEffect) res = res %>% mutate(L = replace(L, age == ages[i], currentL) ) } return(res) } #computeMultipleLengthTrajectoriesWithRandomSeed(temperaturePattern, Nind = 10, #growPar = growParUnisex, # RNGseed = 1) vonBertalanffyWithRandomVector = function(L, Linf, K, timeStepDuration, randomVector, sigma, tempEffect ){ if (sigma == 0) { mu = if_else(L < Linf, log((Linf - L) * (1 - exp(-K * tempEffect * timeStepDuration))),-Inf) #mu = log((Linf - L) * (1 - exp(-K * tempEffect * timeStepDuration))) increment = exp(mu) } else { mu = if_else(L < Linf, log((Linf - L) * (1 - exp(-K * tempEffect * timeStepDuration))) - (sigma * sigma) / 2, -Inf) # mu = log((Linf - L) * (1 - exp(-K * tempEffect * timeStepDuration))) - (sigma * sigma) / 2 increment = exp(randomVector * sigma + mu) } L = pmin(Linf, L + increment) return(L) } computeOgive = function(lengthTrajectories, lengthAtMaturity){ ogive <- lengthTrajectories %>% group_by(age) %>% summarise(nTotal = n()) %>% left_join(lengthTrajectories %>% group_by(age) %>% filter(L >= lengthAtMaturity) %>% summarise(taller = n(), .groups = 'drop'), by = 'age') %>% replace_na(list(taller = 0)) %>% mutate(mature = c(0,diff(taller)), immature = nTotal - cumsum(mature), p = if_else(mature + immature > 0 , mature / (mature + immature), 1)) %>% select(age, p) return(ogive) } # deprecated computeOgive3 = function(lengthTrajectories, lengthAtMaturity) { ogive <- lengthTrajectories %>% group_by(age) %>% summarise(nTotal=n()) %>% left_join(lengthTrajectories %>% group_by(age) %>% filter(L > lengthAtMaturity) %>% summarise(mature = n(), .groups = 'drop'), by = 'age') %>% replace_na(list(mature = 0)) %>% mutate(p = mature/nTotal) %>% select(age, p) return(ogive) } # ------------------------------------------------------- # Dispersal # see see (Chapman et al., 2007; Holloway et al., 2016) # ------------------------------------------------------- logitKernel = function(distance, alpha0, alpha1, meanInterDistance, standardDeviationInterDistance, basinSurface = 0, alpha2 = 0, meanBasinSurface = 0, standardDeviationBasinSurface = 1, fishLength = 0, alpha3 = 0, meanLengthAtRepro = 0, standardDeviationLengthAtRepro = 1) { logitW = alpha0 - alpha1 * (distance - meanInterDistance) / standardDeviationInterDistance + alpha2 * (basinSurface - meanBasinSurface) / standardDeviationBasinSurface + alpha3 * (fishLength - meanLengthAtRepro) / standardDeviationInterDistance; W = 1 / (1+ exp(-logitW)) return(W) } # extended negative exponential kernel eneKernel = function(distance, alpha_D, beta_D){ W = exp(-alpha_D * distance ^ beta_D) return(W) } strayerSurvival = function(distance, m ){ # m: extra mortality rate for strayers (L-1) return(exp(-m * distance)) } #------------------------------------------------------- # Spwaner survival before reproduction # Dome-shape curve with temperature effect #------------------------------------------------------- spawnerSurvivalPreReproduction <- function(Triver,Tmin, Topt, Tmax, parSurv){ #River survival survProbOptRiver = parSurv\$survProbOptRiver spRiver = survProbOptRiver * temperatureEffect(Triver, Tmin, Topt, Tmax) #StockAfterSurv = S * SpRiver return(spRiver) } #---------------------------------------------- # Stock Recruitement relationship #See Rougier et al, 2014; 2015 and Jatteau et al., 2017 #---------------------------------------------- #survival of larvae 14dph from Jatteau et al, 2017 stockRecruitementRelationship <- function(Triver,Tmin, Topt, Tmax, survivalStock, parRep, parSurv, surfaceWatershed = 84810) { #surfaceWatershed = 84810 Surface de la Garonne lambda = parRep\$lambda deltaT = parRep\$deltaT survOptRep = parRep\$survOptRep eta = parRep\$eta #simule l'effet allee ratioS95_S50 = parRep\$ratioS95_S50 alpha = parRep\$a #species fecundity survProbOptRiver = parSurv\$survProbOptRiver periodAtSea = 5 - deltaT ###--------------------- SR relationship ----------------- #parametre c de la RS de BH int?grant un effet du BV consid?r? cj = lambda/surfaceWatershed #parametre b repr?sentant la mortalit? densit? d?pendante de la RS de BH int?grant un effet de la temperature # bj = (-(1/deltaTrecruitement))* # log(survOptRep * temperatureEffect(temp, 9.8, 20.0, 26.0)) bj = -log(survOptRep * temperatureEffect(Triver, Tmin , Topt , Tmax)) / deltaT #parametre a (f?condit? de l'esp?ce) de la RS de BH int?grant un effet de la temperature alphaj = (bj * exp(-bj * deltaT)) / (cj * (1 - exp(-bj * deltaT))) alphaj[is.na(alphaj)] <- 0 #Bj = param?tre de la relation SR int?grant l'effet de la temp?rature betaj = bj/(alpha*cj*(1 - exp(-bj*deltaT))) #p = proportion de g?niteurs participant ? la reproduction en focntion de la quantit? de g?niteur total #p = 1/(1+exp(-log(19)*(S-n)/(Ratio*surfaceWatershed))) S95 = eta * surfaceWatershed S50 = S95/ratioS95_S50 p = 1/(1 + exp(-log(19)*(survivalStock - S50)/(S95 - S50))) #relation Stock Recrutement ie calcul le nombre de recrues en fonction du nombre de g?niteurs et de la T en int?grant l'effet Allee #R0 = aj * S * p alleeEffect = 1/(1 + exp(-log(19)*(survivalStock - eta/ratioS95_S50*surfaceWatershed)/(eta*surfaceWatershed - eta/ratioS95_S50*surfaceWatershed))) Rj = (alphaj * survivalStock * alleeEffect)/(betaj + survivalStock * alleeEffect) #Rj = ((aj * S) * p)/(Bj +S * p) stockRecruitement = as.vector(Rj) return(stockRecruitement) } # ---------------------------------------------- # Spawner survival after reproduction #---------------------------------------------- #Logit with temperature effect depending on a Triver spawnerSurvivalPostReproductionTempRiver <- function(Triver, coeffa, coeffb){ spRiverPostSpawn = (coeffb/ (1 + exp(coeffa*(Triver)))) return(spRiverPostSpawn) } #Logit with temperature effect depending on the difference between the river temperature and Tref spawnerSurvivalPostReproductionTempRef <- function(Triver,Tref, coeffa, coeffb){ spRiverPostSpawn = (coeffb/ (1 + exp(coeffa*(Triver-Tref)))) return(spRiverPostSpawn) } #Logit for survival after reproduction with log19 #similar to spawnerSurvivalPostReproductionWithTempRef which include a logit2 = function(Triver, Tref,minTempForIteroparity){ return( 1/ (1+exp((log(19)/(minTempForIteroparity-Tref))*(Triver-Tref)))) } #Dome-shape curve with temperature effect spawnerSurvivalPostReproductionWithBellCurve <- function(Triver, Tmin, Topt, Tmax, coeffb){ #P1: #SpRiverPostSpawn = probSurvAfterRepro/(probSurvAfterRepro + (1 - probSurvAfterRepro)*exp(-coeffLogit*(TemperatureEffect(Triver, Tmin, Topt, Tmax)))) #P2: #SpRiverPostSpawn = probSurvAfterRepro/(probSurvAfterRepro + (1 - probSurvAfterRepro)*exp(-coeffLogit)) #SpRiverPostSpawnTempEffect = (probSurvAfterRepro/(probSurvAfterRepro + (1 - probSurvAfterRepro)*exp(-coeffLogit))) * TemperatureEffect(Triver, Tmin, Topt, Tmax) #P3: spRiverPostSpawn = coeffb * temperatureEffect(Triver, Tmin, Topt, Tmax) #stockAfterSpawn = S * spRiverPostSpawn return(spRiverPostSpawn) }