#### MODÈLE TK-BLM ------------------------------------------------------------------------------ library(deSolve) # pour utiliser ode pour calculer Ci library(ggplot2) library(dplyr) library(tidyr) #************************************************************************************************ ### Fixer tous les paramètres ------------------------------------------------------------------- #************************************************************************************************ # Temps t_debut <- 0 t_fin <- 21 dt <- 0.1 # avec 1 ça fait des angles temps <- seq(t_debut, t_fin, by = dt) # Concentrations eau (mg/L) Ca_mgL <- 80 Mg_mgL <- 26 Na_mgL <- 6.5 pH <- (7.2 * 3.5 / 24) + (8.6 * 20.5 / 24) # Conversion en mol/L Ca_mol <- Ca_mgL / 40.078 / 1000 Mg_mol <- Mg_mgL / 24.305 / 1000 Na_mol <- Na_mgL / 22.990 / 1000 H_mol <- 10^(-pH) #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ # Paramètres pour les métaux # Paramètres BLM pour le Nickel params_Ni <- list( Jmax= 1240.478, # µg/g/j capacité max de transport ke = 0.212, # j⁻¹ constante d'élimination K_MBL = 10^4.7, # affinité Ni pour ligand (mol/L)⁻¹ K_CaBL = 10^3.29, # affinité Ca pour ligand K_MgBL = 10^2.91, # affinité Mg pour ligand K_NaBL = 0, #pas de logK_NaBL pour le Ni2+ K_HBL = 10^5.75, Ceau_mol = 0.00378 / 58.693 / 1000) # concentration Ni eau en mol/L # concentration Ni eau / masse molaire / 1000 # Paramètres BLM pour le Cadmium params_Cd <- list( Jmax = 90.485, ke = 0.044, K_MBL = 10^6.5, K_CaBL = 10^3.41, K_MgBL = 10^2.91, K_NaBL = 10^2.05, K_HBL = 10^5.75, Ceau_mol = 0.0111 / 112.41 / 1000) # Paramètres BLM pour le Plomb params_Pb <- list( Jmax = 788.028, ke = 0.0548, K_MBL = 10^6.9, K_CaBL = 10^3.4, K_MgBL = 0, #10^2.8, K_NaBL = 0, #10^2.2, K_HBL = 10^5.8, Ceau_mol = 0.01 / 207.2 / 1000) #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ # Paramètres métal + paramètres eau pour modèle TK-BLM params_env_Ni <- list( params = params_Ni, Ca_mol = Ca_mol, Mg_mol = Mg_mol, Na_mol = Na_mol, H_mol = H_mol) params_env_Cd <- list( params = params_Cd, Ca_mol = Ca_mol, Mg_mol = Mg_mol, Na_mol = Na_mol, H_mol = H_mol) params_env_Pb <- list( params = params_Pb, Ca_mol = Ca_mol, Mg_mol = Mg_mol, Na_mol = Na_mol, H_mol = H_mol) #************************************************************************************************ #### Faire les modèles TK-BLM et TK ------------------------------------------------------------- #************************************************************************************************ # Condition initiale: gammare non contaminé au départ donc Ci = 0 state_init_Ni <- c(Ci = 0) state_init_Cd <- c(Ci = 0) state_init_Pb <- c(Ci = 0) # Calcul de J(t) pour le modèle TK-BLM calculer_J <- function(params, Ca_mol, Mg_mol, Na_mol, H_mol) { # calculer le dénominateur avant le gros calcul denominateur <- 1 + params$K_MBL * params$Ceau_mol + params$K_CaBL * Ca_mol + params$K_MgBL * Mg_mol + params$K_NaBL * Na_mol + params$K_HBL * H_mol J <- (params$Jmax * params$K_MBL * params$Ceau_mol) / denominateur return(J) } # Modèle TK-BLM modele_TKBLM <- function(t, state, params_env) { Ci <- state["Ci"] J <- calculer_J(params_env$params, params_env$Ca_mol, params_env$Mg_mol, params_env$Na_mol, params_env$H_mol) dCi <- J - params_env$params$ke * Ci return(list(dCi)) } #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ # Résolution numérique avec ode() pour avoir les Ci (résolution d'une EDO) ---- sol_Ni <- ode(y = state_init_Ni, times = temps, func = modele_TKBLM, parms = params_env_Ni) sol_Cd <- ode(y = state_init_Cd, times = temps, func = modele_TKBLM, parms = params_env_Cd) sol_Pb <- ode(y = state_init_Pb, times = temps, func = modele_TKBLM, parms = params_env_Pb) #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ # Tout convertir en tableau ---- tab_Ni <- as.data.frame(sol_Ni) %>% mutate(metal = "Nickel", calcium_mgL = Ca_mgL) tab_Cd <- as.data.frame(sol_Cd) %>% mutate(metal = "Cadmium", calcium_mgL = Ca_mgL) tab_Pb <- as.data.frame(sol_Pb) %>% mutate(metal = "Plomb", calcium_mgL = Ca_mgL) #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ # Graphiques ---- bind_rows(tab_Ni, tab_Cd, tab_Pb) %>% ggplot(aes(x = time, y = Ci, color = metal)) + geom_line(linewidth = 1.2) + scale_color_manual(values = c("Nickel" = "darkmagenta", "Cadmium" = "sandybrown", "Plomb" = "aquamarine2")) + labs(x = "Temps (jours)", y = "Concentration interne Ci (µg/g PS)", color = "Métal") + ggtitle("Bioaccumulation des métaux") + scale_y_log10( breaks = 10^(0:5), labels = c("1", "10", "100", "1 000", "10 000", "100 000")) + theme_minimal() + theme(plot.title = element_text(face = "bold", hjust = 0.5)) #************************************************************************************************ #### Ci en fonction de Ca2+ pour voir les effets protecteurs ------------------------------------ #************************************************************************************************ # Sur le modèle TK de base, le calcium n'a aucun effet vu qu'il est pas dans les paramètres # On peut pas voir les effets protecteurs du calcium sur l'Arsenic et le Chrome # Faire une échelle de Ca Ca_gamme_mgL <- seq(1, 200, by = 1) Ca_gamme_mol <- Ca_gamme_mgL / 40.078 / 1000 calculer_Ci_eq <- function(params, Ca_gamme_mol, Mg_mol, Na_mol, H_mol) { Ci_eq_vec <- sapply(Ca_gamme_mol, function(ca) { J_eq <- calculer_J(params, ca, Mg_mol, Na_mol, H_mol) Ci_eq <- J_eq / params$ke return(Ci_eq) }) return(Ci_eq_vec) } #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ # Appel pour chaque métal ---- tab_eq_Cd <- data.frame( calcium_mgL = Ca_gamme_mgL, Ci_eq = calculer_Ci_eq(params_Cd, Ca_gamme_mol, Mg_mol, Na_mol, H_mol), metal = "Cadmium") tab_eq_Ni <- data.frame( calcium_mgL = Ca_gamme_mgL, Ci_eq = calculer_Ci_eq(params_Ni, Ca_gamme_mol, Mg_mol, Na_mol, H_mol), metal = "Nickel") tab_eq_Pb <- data.frame( calcium_mgL = Ca_gamme_mgL, Ci_eq = calculer_Ci_eq(params_Pb, Ca_gamme_mol, Mg_mol, Na_mol, H_mol), metal = "Plomb") # Coller tous les métaux ensemble tab_eq_tous <- bind_rows(tab_eq_Cd, tab_eq_Ni, tab_eq_Pb) #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ # Graphique ggplot(tab_eq_tous, aes(x = calcium_mgL, y = Ci_eq, color = metal, group = metal)) + geom_line(linewidth = 1.2) + scale_y_log10( breaks = 10^(0:5), labels = c("1", "10", "100", "1 000", "10 000", "100 000")) + labs(title = "Effet protecteur du calcium sur l'accumulation des métaux", x = "Concentration en calcium (mg/L)", y = "Ci à l'équilibre (µg/g PS)", color = "Métal") + scale_color_manual(values = c("Nickel" = "darkmagenta", "Cadmium" = "sandybrown", "Plomb" = "aquamarine2")) + theme_minimal() + theme(plot.title = element_text(face = "bold", hjust = 0.5)) #************************************************************************************************ #### J en fonction de Ca2+, Mg2+, Na+ (Effet de la compétition) --------------------------------- #************************************************************************************************ calculer_J_gamme_ion <- function(params, ion_gamme_mol, ion_variable, Ca_mol_fixe, Mg_mol_fixe, Na_mol_fixe, H_mol) { J_vec <- sapply(ion_gamme_mol, function(val) { ca <- if (ion_variable == "Ca") val else Ca_mol_fixe mg <- if (ion_variable == "Mg") val else Mg_mol_fixe na <- if (ion_variable == "Na") val else Na_mol_fixe calculer_J(params, ca, mg, na, H_mol) }) return(J_vec) } # Thème commun pour les 3 graphiques J = f(ion) theme_J <- theme_minimal() + theme( axis.title = element_text(size = 16), axis.text = element_text(size = 16), legend.title = element_text(size = 16, face = "bold"), legend.text = element_text(size = 16)) couleurs_metaux <- c("Nickel" = "darkmagenta", "Cadmium" = "sandybrown", "Plomb" = "aquamarine2") #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ # J = f([Ca2+]) # Gamme de calcium: 1 à 200 mg/L Ca_gamme_mgL <- seq(1, 200, by = 1) Ca_gamme_mol <- Ca_gamme_mgL / 40.078 / 1000 tab_J_Ca <- bind_rows( data.frame(ion_mgL = Ca_gamme_mgL, J = calculer_J_gamme_ion(params_Cd, Ca_gamme_mol, "Ca", Ca_mol, Mg_mol, Na_mol, H_mol), metal = "Cadmium"), data.frame(ion_mgL = Ca_gamme_mgL, J = calculer_J_gamme_ion(params_Ni, Ca_gamme_mol, "Ca", Ca_mol, Mg_mol, Na_mol, H_mol), metal = "Nickel"), data.frame(ion_mgL = Ca_gamme_mgL, J = calculer_J_gamme_ion(params_Pb, Ca_gamme_mol, "Ca", Ca_mol, Mg_mol, Na_mol, H_mol), metal = "Plomb")) fig_J_Ca <- ggplot(tab_J_Ca, aes(x = ion_mgL, y = J, color = metal)) + geom_line(linewidth = 1.2) + scale_y_log10(labels = scales::label_comma()) + scale_color_manual(values = couleurs_metaux) + labs(x = "Concentration en calcium (mg/L)", y = "Coefficient J (µg/g/j) (échelle log)", color = "Métal") + theme_J print(fig_J_Ca) #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ # J = f([Mg2+]) # Gamme de magnésium: 1 à 100 mg/L Mg_gamme_mgL <- seq(1, 100, by = 0.5) Mg_gamme_mol <- Mg_gamme_mgL / 24.305 / 1000 tab_J_Mg <- bind_rows( data.frame(ion_mgL = Mg_gamme_mgL, J = calculer_J_gamme_ion(params_Cd, Mg_gamme_mol, "Mg", Ca_mol, Mg_mol, Na_mol, H_mol), metal = "Cadmium"), data.frame(ion_mgL = Mg_gamme_mgL, J = calculer_J_gamme_ion(params_Ni, Mg_gamme_mol, "Mg", Ca_mol, Mg_mol, Na_mol, H_mol), metal = "Nickel"), data.frame(ion_mgL = Mg_gamme_mgL, J = calculer_J_gamme_ion(params_Pb, Mg_gamme_mol, "Mg", Ca_mol, Mg_mol, Na_mol, H_mol), metal = "Plomb")) fig_J_Mg <- ggplot(tab_J_Mg, aes(x = ion_mgL, y = J, color = metal)) + geom_line(linewidth = 1.2) + scale_y_log10(labels = scales::label_comma()) + scale_color_manual(values = couleurs_metaux) + labs(x = "Concentration en magnésium (mg/L)", y = "Coefficient J (µg/g/j) (échelle log)", color = "Métal") + theme_J print(fig_J_Mg) #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ # J = f([Na+]) # Gamme de sodium: 1 à 150 mg/L Na_gamme_mgL <- seq(1, 150, by = 0.5) Na_gamme_mol <- Na_gamme_mgL / 22.990 / 1000 tab_J_Na <- bind_rows( data.frame(ion_mgL = Na_gamme_mgL, J = calculer_J_gamme_ion(params_Cd, Na_gamme_mol, "Na", Ca_mol, Mg_mol, Na_mol, H_mol), metal = "Cadmium"), data.frame(ion_mgL = Na_gamme_mgL, J = calculer_J_gamme_ion(params_Ni, Na_gamme_mol, "Na", Ca_mol, Mg_mol, Na_mol, H_mol), metal = "Nickel"), data.frame(ion_mgL = Na_gamme_mgL, J = calculer_J_gamme_ion(params_Pb, Na_gamme_mol, "Na", Ca_mol, Mg_mol, Na_mol, H_mol), metal = "Plomb")) fig_J_Na <- ggplot(tab_J_Na, aes(x = ion_mgL, y = J, color = metal)) + geom_line(linewidth = 1.2) + scale_y_log10(labels = scales::label_comma()) + scale_color_manual(values = couleurs_metaux) + labs(x = "Concentration en sodium (mg/L)", y = "Coefficient J (µg/g/j) (échelle log)", color = "Métal",) + theme_J print(fig_J_Na) #************************************************************************************************ #### Inverser le modèle pour avoir Ceau = f(Ci) ------------------------------------------------- #************************************************************************************************ # même logique que juste avant, c'est juste pas la même fonction calculer_Ceau <- function(Ci_t, params, Ca_mol, Mg_mol, Na_mol, H_mol, t) { somme <- 1 + params$K_CaBL * Ca_mol + params$K_MgBL * Mg_mol + params$K_NaBL * Na_mol + params$K_HBL * H_mol Ceau_mol <- ( Ci_t * params$ke * somme ) / ( params$K_MBL * (params$Jmax * (1 - exp(-params$ke * t)) - Ci_t * params$ke)) return(Ceau_mol) } t_expo <- 7 # Gamme de Ci observées Ci_gamme <- seq(0.1, 50, by = 0.1) # µg/g PS Ceau_predite_Cd <- sapply(Ci_gamme, function(ci) { calculer_Ceau(ci, params_Cd, Ca_mol, Mg_mol, Na_mol, H_mol, t_expo) }) Ceau_predite_Ni <- sapply(Ci_gamme, function(ci) { calculer_Ceau(ci, params_Ni, Ca_mol, Mg_mol, Na_mol, H_mol, t_expo) }) Ceau_predite_Pb <- sapply(Ci_gamme, function(ci) { calculer_Ceau(ci, params_Pb, Ca_mol, Mg_mol, Na_mol, H_mol, t_expo) }) # Conversion mol/L -> µg/L Ceau_predite_ugL_Cd <- Ceau_predite_Cd * 112.41 * 1e6 Ceau_predite_ugL_Ni <- Ceau_predite_Ni * 58.693 * 1e6 Ceau_predite_ugL_Pb <- Ceau_predite_Pb * 207.2 * 1e6 tab_inverse_Cd <- data.frame( Ci_ugkg = Ci_gamme, Ceau_ugL = Ceau_predite_ugL_Cd, metal = "Cadmium") tab_inverse_Ni <- data.frame( Ci_ugkg = Ci_gamme, Ceau_ugL = Ceau_predite_ugL_Ni, metal = "Nickel") tab_inverse_Pb <- data.frame( Ci_ugkg = Ci_gamme, Ceau_ugL = Ceau_predite_ugL_Pb, metal ="Plomb") tab_inverse_tous <- bind_rows(tab_inverse_Cd, tab_inverse_Ni, tab_inverse_Pb) ggplot(tab_inverse_tous, aes(x = Ci_ugkg, y = Ceau_ugL, color = metal, group = metal)) + geom_line(linewidth = 1.2) + scale_color_manual( values = c("Cadmium" = "sandybrown", "Nickel" = "darkmagenta", "Plomb" = "aquamarine2")) + labs(title = "Concentration dans l'eau en fonction de la concentration interne", x = "Concentration interne Ci (µg/g PS)", y = "Concentration dans l'eau Ceau (µg/L)", color = "Métal") + theme_minimal() + theme(plot.title = element_text(face = "bold", hjust = 0.5))