mmc2.r
In lemna: Lemna Ecotox Effect Model

##
## A mechanistic combined toxicokinetic–toxicodynamic (TK/TD) and growth
## model for the aquatic macrophytes Lemna spp. as published by
## Schmitt et al. (2013), DOI: 10.1016/j.ecolmodel.2013.01.017
##
## Parts of the implementation were commented out for reasons of brevity.
##

#
# options(warn=-1)
#
# #********** Control of calculation *********************************************
# # Simulation name  (will be used for naming output files !!!)
# Simname = 'Test'
#
# # Working directory
# workdir = 'C:/Users/zfwsm/Documents/4_R Projects/Lemna Model/Growth model_official'
#
# # F i l e s to be loaded
# # Name ='' for unused files
# Conc_file = 'Peak_2d@t=0.txt'
# Temperature_file = ''
# Irradiation_file = ''
# ##Data_file = 'Data.txt',
#
# # Save results after execution?
# fileout = F       # T=true/F=false
#
# #---- Settings of simulation ---------------------------------------------------
# # Define method to be used
# callfun = 'CS'   # CS = single simulation, DR = Dose response, MC = MC Simulation
# # Set  timepoints for output
# outtimes = seq(0,100,1)
# # DOSE RESPONSE: Define exposure steps/factors (only used for "Dose Response")
# expfact = 10^seq(0,1.5,0.25)
# # MC SIMULATION: Define list of parameters to vary and No of replicates
# MCrep = 100
# parlist = list(
#     parnames=c(
#       'k_G_max',
#       'k_mort',
# #      'Emax',
#       'EC50'),
# #      'b',
# #      'BM50',
# #      'P_up'),
#     pardis=list(
#     c(0.47,0.2,1),
#     c(0.05,0.1,1),
# #    c(1,0,1),
#     c(.48,.1,2))
# #    c(6.4,.1,1),
# #    c(176,.1,1),
# #    c(.3,.2,2))
# )


# Set initial values
vars <- c(
BM       = 50,     # [g_dw/m?]   Dry Biomass dryweight per m2
E        = 1,      # (0-1)      (Toxic) Effect = Factor on growth rate  (Range: 0 - 1, 1=no effect)
M_int    = 0       # [?g]       Amount of toxicant in biomass
)

#---- Parameter Setting --------------------------------------------------------
param <- list(
#
#      - Effect -
Conc     = 1,   #  [any]      Concentration of toxicant (may also be a table)
Emax     = .784,    #   maximum Effect
EC50     = 0.3,     #  [same as conc. data]      Midpoint of effect curve
b        = 4.16,    #  [-]         Slope of effect curve
#
#       - Toxicokinetics -
P_up     = .0054,   # [cm/d]      Permeability for uptake
AperBM   = 1000,   # [cm?/g_dw]   A_leaf / d_leaf = 1/d_leaf (for circular disc, d=0.05 cm) [Ref. HARLAN-022]
Kbm      = .75,   # []          Biomass(fw):water partition coefficient
P_Temp = F,       # Switch for temperature dependence of cuticle permeability
MolWeight = 381,  # Molmass of molecule (determines Q10_permeability)



#       - Fate of biomass -
k_phot_fix  = F,   #  T/F          If True k_G_max is not changed by environmental factors
k_phot_max  = 0.47, #  [1/d]       Maximum growth rate of biomass + kmort     [Ref. F 0191, Harlan-011]
k_resp   = 0.05, #  [1/d]       Rate of mortality                 [Ref. Harlan-011, rough estimate]
k_loss   = 0.0, #  [1/d]       Some rate of loss (e.g. Flow rate)
#
#      - Temperatur dependence -
# k_phot
Temp     = 12,   #  [?C]        Current temperature (may also be a table)
Tmin     = 8.0  , #  [?C]        Minimum growth temperature      [Ref. F 0191, data re-evaluated  incl. kmort(T)]
Tmax     = 40.5 , #  [?C]        Maximum growth temperature      [Ref. F 0191, data re-evaluated  incl. kmort(T)]
Topt     = 26.7 , #  [?C]        Optimum growth temperature      [Ref. F 0191, data re-evaluated  incl. kmort(T)]
# k_resp
t_ref       = 25,   # temperature at which t_mort is effective
Q10         = 2,
#
#      - Light dependence (linear dependence on global radiation (see Hodgeson 1969)
Rad 	    = 15000 , #  [kJ/m?/d]  Radiation  (may also be given as table)
k_0     	= 3 ,     #  [1/d]      Intercept of linear part
a_k     	= 5E-5 ,  #  [(1/d)/(kJ/m?/d)]        Slope of linear part

#      - Phosphorus dependence (Hill like dependence) -
C_P      = 0.3,          #   [mg/L]       Phosporus concentration in water
CP50     = 0.0043,        #   [mg/L]       P-conc. where growth rate is halfened   [Data from L??nd, 1983 evaluated with monod model]
a_P      = 1,        #   []            Hill coefficient
KiP      = 101,         #   [mg/L]       P-inhibition constant for very high P-conc.   [Ref. F 0191]

#      - Nitrogen dependence (Hill like dependence) -
C_N      = 0.6,         #   [mg/L]       Nitrogen concentration in water
CN50     = 0.034,          #   [mg/L]       N-conc. where growth rate is halfened   [Data from L??nd, 1983 evaluated with monod model]
a_N      = 1,        #   []            Hill coefficient
KiN      = 604,         #   [mg/L]       n-inhibition constant for very high P-conc.   [Ref. F 0191]

#      - Density dependence -
BM50     = 176,    #  [g_dw/m?]             Cut off BM   [Ref. F 0191]

#      - Others -
mass_per_frond     = 0.0001,   #  [g_dw/frond]  Dryweight per frond [Ref. HARLAN-022]
BMw2BMd     = 16.7   #  [g_fresh/g_dry]  Fresh- / dryweight [Ref. F 0191]

)



# **************** Call of scripts ********************************************
# +++++++  Do not change the following !!!! +++++++++++++++++++++++++++++++++++

# if (Conc_file!=''){Conc_file = paste(workdir,'/',Conc_file,sep="")}
# if (Temperature_file!=''){Temperature_file = paste(workdir,'/',Temperature_file,sep="")}
# if (Irradiation_file!=''){Irradiation_file = paste(workdir,'/',Irradiation_file,sep="")}
# #Data_file = 'Data.txt',
# filenames = c(
# Conc_file=Conc_file,
# Temperature_file=Temperature_file,
# Irradiation_file=Irradiation_file
# )
#
# # Function to read source files (needs to be defined before first call)
# sourceDir <- function(path, trace = TRUE, ...) {
#    for (nm in list.files(path, pattern = "\\.[RrSsQq]$")) {
#       if(trace) cat(nm,":")
#       source(file.path(path, nm), ...)
#       if(trace) cat("\n")
#     }
#  }
#
# # Read source files
# sourceDir(paste(workdir,'/SourceFun',sep=""))
#
# if (callfun == 'CS'){
# # Call: calcgrowth(timepoints, vars, param, filenames,graph=T)
#   res <- calcgrowth(outtimes, vars, param, filenames,graph=T)
#   if(fileout){write.table(res,paste(Simname,'_out.dat',sep=''))}
# }
#
# if (callfun=='DR'){
# # Call:  dos_resp('growthmodel',doses,eff_frac,t_recov,timepoints,vars,param,filenames,graph=T)
#    dosres <- dos_resp('growthmodel',expfact,0.01,5,outtimes,vars,param,filenames,graph=F)
#    res <- dosres$Popcurves
#
#    if(fileout){
#       write.table(dosres$Popcurves,paste(outname,'_data.dat',sep=""))
#       write.table(dosres$Concentrations,paste(outname,'_conc.dat',sep=""))
#       write.table(dosres$Effects,paste(outname,'_eval.dat',sep=""))
#       }
# }
#
# if(callfun=='MC'){
# # Call:    MCsim(model,parlist,MCrep,simperiod,vars,param,filenames,graph,eff_frac,t_recov)
#     MCres <- MCsim('growthmodel',parlist,MCrep,outtimes,vars,param,filenames,graph=F,0.01,28)
# }