diff --git a/ModelF/Changed/AFIRT_calculation.R b/ModelF/Changed/AFIRT_calculation.R new file mode 100644 index 0000000..2e698e5 --- /dev/null +++ b/ModelF/Changed/AFIRT_calculation.R @@ -0,0 +1,381 @@ +#---------------------------------------------------------------------------------- +# Helper function that returns a range of variable when performing +# sensitvity analysis +#---------------------------------------------------------------------------------- + +read.param.file = function(filename) { + d = read_excel(filename, 1) + param.as.double = as.numeric(d$Value) + names(param.as.double) = d$Parameter + param.as.double = param.as.double[model$pin] #keep only parameters used in ODE +} + +lseq = function(from, to, length.out){ + # Arguments: + # from : initial value of the variable + # to : teminal value of the variable + # length.out : fold number of + # Return : + # A vector of length + + sequence = seq(log(from), log(to), length.out=length.out) + sequence = exp(sequence) + return(sequence) +} + +#---------------------------------------------------------------------------------- +# Function computing lumped parameters from theory +#---------------------------------------------------------------------------------- + +lumped.parameters.theory = function(param.as.double = param.as.double, + dose.nmol = dose.nmol, + tau = tau, + soluble = FALSE){ + # Arguments: + # params_file_path: full path of the parameters file. + # dose.nmol: dosing amout in nmol + # tau: dosing interval in days + # soluble: flag saying whether or not drug is soluble. Default to FALSE. + # Return: + # A data frame of lumped parameters calculated from theory + + pars = as.data.frame(t(param.as.double)) + + # Calculate M3tot.ss, M30, S3tot.ss and S30 + numerator.DM3 = with(pars, k13DM*(VD1/VD3)* ksynM1 + (keDM1 + kshedDM1 + k13DM)* ksynM3) + numerator.M3 = with(pars, k13M *(VD1/VD3)* ksynM1 + (keM1 + kshedM1 + k13M )* ksynM3) + numerator.DM1 = with(pars, k31DM*(VD3/VD1)* ksynM3 + (keDM3 + kshedDM3 + k31DM)* ksynM1) + numerator.M1 = with(pars, k31M *(VD3/VD1)* ksynM3 + (keM3 + kshedM3 + k31M )* ksynM1) + + denomenator.DM_1or3 = with(pars, (keDM1 + kshedDM1 + k13DM)*(keDM3 + kshedDM3 + k31DM)-k31DM*k13DM) + denomenator.M_1or3 = with(pars, (keM1 + kshedM1 + k13M )*(keM3 + kshedM3 + k31M )-k31M *k13M ) + + M1tot.ss = numerator.DM1 / denomenator.DM_1or3 + M3tot.ss = numerator.DM3 / denomenator.DM_1or3 + M10 = numerator.M1 / denomenator.M_1or3 + M30 = numerator.M3 / denomenator.M_1or3 + + #note that this aligns with the numerator columns above and can be copied and pasted for comparison + numerator.DS3 = with(pars, k13DS*(VD1/VD3)*(ksynS1 + kshedDM1*M1tot.ss) + (keDS1 + k13DS)*(ksynS3 + kshedDM3*M3tot.ss)) + numerator.S3 = with(pars, k13S *(VD1/VD3)*(ksynS1 + kshedM1 *M10) + (keS1 + k13S) *(ksynS3 + kshedM3 *M30)) + + denomenator.DS3= with(pars, (keDS1 + k13DS)*(keDS3 + k31DS)-k31DS*k13DS) + denomenator.S3 = with(pars, (keS1 + k13S )*(keS3 + k31S )-k31S *k13S ) + + S3tot.ss = numerator.DS3 / denomenator.DS3 + S30 = numerator.S3 / denomenator.S3 + + + if (!soluble){ + Kssd = with(pars, (koff3 + keDM3 + kshedDM3 + k31DM)/kon3) + Kss = with(pars, (koff3 + keDM3 + kshedDM3 )/kon3) + Kd = with(pars, koff3 /kon3) + Tacc.tum = M3tot.ss / M30 + + } else { + Kssd = with(pars, (koff3 + keDS3 + + k31DS)/kon3) + Kss = with(pars, (koff3 + keDS3 )/kon3) + Kd = with(pars, koff3 /kon3) + Tacc.tum = S3tot.ss / S30 + } + + # Biodistribution coefficient (reference: ModelF_Appendix) + B = with(pars, (k13D/(keD3 + k31D) * (VD1/VD3))) + + # Clearance + CL = with(pars, (keD1*VD1)) + + # Average drug concentration in the central compartment + Cavg1 = dose.nmol/(CL*tau) + + # Compute various AFIRTs + AFIRT.Kssd = Kssd*Tacc.tum/(B*Cavg1) + AFIRT.Kss = Kss *Tacc.tum/(B*Cavg1) + AFIRT.Kd = Kd *Tacc.tum/(B*Cavg1) + + a0 = with(pars, keD1*k21D*k31D) + a1 = with(pars, keD1*k31D + k21D*k31D + k21D*k13D + keD1*k21D + k31D*k12D) + a2 = with(pars, keD1 + k12D + k13D + k21D + k31D) + + p = a1 - (a2^2)/3 + q = 2*(a2^3)/27 - a1*a2/3 + a0 + r1 = (-(p^3)/27)^0.5 + r2 = 2*(r1^(1/3)) + + phi = acos(-q/(2*r1))/3 + + + alpha = -(cos(phi) *r2 - a2/3) + beta = -(cos(phi + 2*pi/3)*r2 - a2/3) + gamma = -(cos(phi + 4*pi/3)*r2 - a2/3) + + V = with(pars, VD1) + A = with(pars, (1/V) * ((k21D - alpha)/(alpha - beta)) * ((k31D - alpha)/(alpha - gamma))) + B = with(pars, (1/V) * ((k21D - beta )/(beta - alpha)) * ((k31D - beta )/(beta - gamma ))) + C = with(pars, (1/V) * ((k21D - gamma)/(gamma - beta)) * ((k31D - gamma)/(gamma - alpha))) + + D = dose.nmol + Cmin1 = D*((A*exp(-alpha*tau))/(1 - exp(-alpha*tau)) + + (B*exp(-beta *tau))/(1 - exp(-beta *tau)) + + (C*exp(-gamma*tau))/(1 - exp(-gamma*tau))) + + if(!soluble){ + Tfold = M3tot.ss/M30 + }else{ + Tfold = S3tot.ss/S30 + } + + TFIRT.Kssd = Kssd*Tfold/(B*Cmin1) + TFIRT.Kss = Kss *Tfold/(B*Cmin1) + TFIRT.Kd = Kd *Tfold/(B*Cmin1) + + lumped_parameters_theory = data.frame(type = "theory", + M30 = M30, + S30 = S30, + M3tot.ss = M3tot.ss, + S3tot.ss = S3tot.ss, + Tacc.tum = Tacc.tum, + B = B, + Cavg1 = Cavg1, + Cavg3 = B*Cavg1, + Cmin = Cmin1, + Cmin1 = Cmin1, + AFIRT.Kssd = AFIRT.Kssd, + AFIRT.Kss = AFIRT.Kss, + AFIRT.Kd = AFIRT.Kd, + AFIRT = AFIRT.Kssd, + TFIRT.Kssd = TFIRT.Kssd, + TFIRT.Kss = TFIRT.Kss, + TFIRT.Kd = TFIRT.Kd, + TFIRT = TFIRT.Kssd, + stringsAsFactors = FALSE) + return(lumped_parameters_theory) + } + +#---------------------------------------------------------------------------------- +# Function simulates the lumped parameters +#---------------------------------------------------------------------------------- + +lumped.parameters.simulation = function(model = model, + param.as.double = param.as.double, + dose.nmol = dose.nmol, + tmax = tmax, + tau = tau, + compartment, + soluble = FALSE){ + + # Arguments: + # model_name: name of the model + # params_file_path: full path of the parameters file. + # dose.nmol: dosing amount in nmol + # tmax: maximum doing period in days + # tau: dosing interval in days + # compartment: compartment to which dosing is applied + # (in model F case, compartment=2) + # soluble: flag saying whether or not drug is soluble. Default to FALSE. + # Return: + # A data frame of lumped parameters calculated from simulation + + # Run simulation + #d <- xlsx::read.xlsx(params_file_path, 1) + #param.as.double = d$Value + #names(param.as.double) = d$Parameter + ev = eventTable(amount.units="nmol", time.units="days") + sample.points = c(seq(-7, tmax, 0.1), 10^(-3:0)) # sample time, increment by 0.1 + sample.points = sort(sample.points) + sample.points = unique(sample.points) + ev$add.sampling(sample.points) + ev$add.dosing(dose=dose.nmol, nbr.doses=floor(tmax/tau)+1, dosing.interval=tau, + dosing.to=compartment) + + #model = do.call(model_name, list()) # model file can contain only one model + init = model$init(param.as.double) + out = model$rxode$solve(param.as.double, ev, init) + out = model$rxout(out) + out = out %>% + mutate(Sfree.pct = S3/init["S3"], + Mfree.pct = M3/init["M3"], + dose.nmol = dose.nmol) + + # Calculate lumped parameters + initial_state = out %>% + filter(time==0) + M30 = initial_state$M3 + S30 = initial_state$S3 + + ## Assume the system reaches steady state during the last dosing period + steady_state = out %>% + filter(time > (floor(tmax/tau)-1)*tau & time % + filter(time > 0) + + # Average drug concentration in central compartment + Cavg1 = mean(steady_state$D1) + # Minimum drug concentration in central compartment + Cmin1 = min(steady_state$D1) + + # Average drug concentration in tumor compartment + Cavg3 = mean(steady_state$D3) + # Minimum drug concentration in tumor compartment + Cmin3 = min(steady_state$D3) + + # AFIRT and target accumulation + if (soluble) { + AFIRT = mean(steady_state$Sfree.pct) + Tacc.tum = S3tot.ss / S30 + } else { + AFIRT = mean(steady_state$Mfree.pct) + Tacc.tum = M3tot.ss / M30 + } + + # Simulation of TFIRT + if (soluble) { + TFIRT = max(steady_state$Sfree.pct) + } else { + TFIRT = max(steady_state$Mfree.pct) + } + + + lumped_parameters_sim = data.frame(type = "simulation", + M30 = M30, + M3tot.ss = M3tot.ss, + S30 = S30, + S3tot.ss = S3tot.ss, + Tacc.tum = Tacc.tum, + Cavg1 = Cavg1, + Cavg3 = Cavg3, + Cmin = Cmin1, + Cmin1 = Cmin1, + Cmin3 = Cmin3, + B = Cavg3/Cavg1, + AFIRT = AFIRT, + AFIRT.sim = AFIRT, + TFIRT = TFIRT, + TFIRT.sim = TFIRT, + stringsAsFactors = FALSE) #having one named sim will be helpful later on in Task01, Task02, etc. + + return(lumped_parameters_sim) +} + +#---------------------------------------------------------------------------------- +# Function runs simulation +#---------------------------------------------------------------------------------- + + +simulation = function(model = model, + param.as.double = param.as.double, + dose.nmol = dose.nmol, + tmax = tmax, + tau = tau){ + #d <- xlsx::read.xlsx(params_file_path, 1) + #param.as.double = d$Value + #names(param.as.double) = d$Parameter + ev = eventTable(amount.units="nmol", time.units="days") + sample.points = c(seq(-7, tmax, 0.1), 10^(-3:0)) # sample time, increment by 0.1 + sample.points = sort(sample.points) + sample.points = unique(sample.points) + ev$add.sampling(sample.points) + ev$add.dosing(dose=dose.nmol, nbr.doses=floor(tmax/tau)+1, dosing.interval=tau, + dosing.to=2) + + init = model$init(param.as.double) + out = model$rxode$solve(param.as.double, ev, init) + out = model$rxout(out) + out = out %>% + mutate(Sfree.pct = S1/init["S1"], + Mfree.pct = M3/init["M3"], + dose.nmol = dose.nmol) + return(out) +} + +#---------------------------------------------------------------------------------- +# This function does the sensitivity analysis on the user inputted parameter +# and compares the theoretical result to the simulated result. +#---------------------------------------------------------------------------------- + +# Input: +# model - model system of ODE's solved with RxODE. In this project, it is 'ivsc_4cmtct_shedct'. +# param.as.double - read parameters from Excel file. read.param.file("file directory"). +# dose.nmol - dose in nmol. +# tmax - time of treatment in days +# tau - frequency of administering dose in days +# compartment - compartment where drug is administered +# param.to.change - parameter on which to do SA. This must be a string. +# param.to.change.range - range of parameter on which to do SA. The range must be symmetric in fold change. This must be a vector of odd length. +# soluble - boolean that is true/false if the drug is soluble/insoluble. Need this since soluble and insoluble are treated differently. + +# Output: +# Data frame of AFIRT vs parameter value + +compare.thy.sim = function(model = model, + param.as.double = param.as.double, + dose.nmol = dose.nmol, + tmax = tmax, + tau = tau, + compartment = compartment, + param.to.change = param.to.change, + param.to.change.range = param.to.change.range, + soluble = FALSE) { + + # Save initial parameter vector + param0 = param.as.double + param.to.change.name = param.to.change #character of parameter to + param.to.change.0 = ifelse(param.to.change.name=="dose", + dose.nmol, + param0[param.to.change.name]) + + + + # Simulation - iterate through values in range + df_sim = data.frame() + if (param.to.change == 'dose'){ + for (param.iter in param.to.change.range){ + row = lumped.parameters.simulation(model, param.as.double, param.iter, tmax, tau, compartment, soluble) + df_sim = rbind(df_sim, row) + } + } else { + for (param.iter in param.to.change.range){ + param.as.double[param.to.change] = param.iter + row = lumped.parameters.simulation(model, param.as.double, dose.nmol, tmax, tau, compartment, soluble) + df_sim = rbind(df_sim, row) + } + } + + #save parameters that were changed in the simulation output + df_sim = df_sim %>% + mutate(param.to.change = param.to.change.range, + fold.change = param.to.change.range/param.to.change.0) + + # Theory - Iterate through values in range. + df_thy = data.frame() + if (param.to.change == 'dose'){ + for (param.iter in param.to.change.range){ + row = lumped.parameters.theory(param.as.double, param.iter, tau, soluble) + df_thy = rbind(df_thy, row) + } + } else{ + for (param.iter in param.to.change.range){ + param.as.double[param.to.change] = param.iter + row = lumped.parameters.theory(param.as.double, dose.nmol, tau, soluble) + df_thy = rbind(df_thy, row) + } + } + + #save value that were change + df_sim = df_sim %>% + mutate(param.to.change = param.to.change.range, + fold.change = param.to.change.range/param.to.change.0) + + # Arrange theory and simulation in single data frame. + df_compare = bind_rows(df_thy,df_sim) + df_compare = df_compare %>% + mutate(param = param.to.change.name) %>% + arrange(param.to.change,type) %>% + mutate_if(is.numeric,signif,2) + + return(df_compare) +} \ No newline at end of file diff --git a/ModelF/Changed/Task05_Sensitivity_Analysis_pt2.Rmd b/ModelF/Changed/Task05_Sensitivity_Analysis_pt2.Rmd new file mode 100644 index 0000000..aeb6af6 --- /dev/null +++ b/ModelF/Changed/Task05_Sensitivity_Analysis_pt2.Rmd @@ -0,0 +1,283 @@ +--- +title: "Task05_Sensitivity_Analysis_pt2" +author: "Miandra Ellis" +date: "`r format(Sys.time(), '%d %B, %Y')`" +output: html_document +# output: +# html_document: +# toc: true +# toc_float: true +# code_folding: hide +--- + +# Preliminary stuff + +```{r, warning=FALSE} +# Need this to get rxode working. I don't know why I keep having to run this. +Sys.setenv(PATH = paste(Sys.getenv("PATH"), "C:/RBuildTools/3.4/bin/", + "C:/RBuildTools/3.4/mingw_64/bin", sep = ";")) +Sys.setenv(BINPREF = "C:/RBuildTools/3.4/mingw_64/bin/") + +# To be called at the top of every Rmd file. Initialization code and some useful constants. +suppressMessages(source("ams_initialize_script.R")) +``` + +# Create data frame for multiple drugs and parameters. + +```{r, include=FALSE} + +# Note: Insoluble and soluble drugs are treated differently. +# For insoluble, Tacc = M3totss/M30. For soluble, Tacc = S3totss/S30. +# We are also interested in different parameters in the SA for either case. +# See writeup for details. + +# -------------------------------------------------------------------------------- +# Initialize. +# -------------------------------------------------------------------------------- + +# Load model. +model = ivsc_4cmtct_shedct(target = TRUE) +# Drugs to explore. +drugs = c("Atezolizumab", "Herceptin", "Pembrolizumab", "Bevacizumab") +# List of parameters of interest. +parameters = c("dose", "ksynM3", "k13D", "kshedM3", "kshedDM3", "VD3") +# names(parameters) = parameters +# Create vector of soluble drugs. +soluble = "Bevacizumab" + +# Create paths to data files for each drug. +paths = NULL +for (i in 1:length(drugs)) { + paths[i] = paste0("../data/ModelF_", drugs[i],"_Params.xlsx") +} + +# Dose time, frequency, compartment, nominal dose +tmax = 26*7 #days +tau = 21 #days +compartment = 2 +dose.nmol = scale.mpk2nmol +joined = NULL +isSoluble = FALSE + +# -------------------------------------------------------------------------------- +# Iterate over all the drugs. +# -------------------------------------------------------------------------------- +i.dfs = 0 +for (i in 1:length(drugs)){ + # Load parameters. + param.as.double = read.param.file(paths[i]) + df_param = as.data.frame(t(param.as.double)) + + # Check if drug has soluble target. + if (drugs[i] %in% soluble) { + # Change any non-shed parameters including M to S. + parameters = sapply(parameters, + function(x) { + if (!grepl("shed|dose",x)) {return(gsub("M","S", x))} + return(x) + }) + isSoluble = TRUE + } + + # Set range for parameters of interest in SA. + # Check which parameters are nonzero, not including dose which isn't in df_param. + nnzero = df_param[parameters[which(parameters != "dose")]] != 0 + nnzero = colnames(nnzero)[which(nnzero)] + params.to.iterate = data.frame(dose = lseq(scale.mpk2nmol*.1,scale.mpk2nmol*10,7), + lapply(df_param[nnzero], function(x) lseq(x*0.1, x*10, 7))) + + dfs = list() + # Iterate all of the parameters for a single drug. + for (j in 1:ncol(params.to.iterate)){ + i.dfs = i.dfs+1 + param.to.change = names(params.to.iterate)[j] + param.to.change.range = params.to.iterate[[j]] + soluble = isSoluble + compare_i = compare.thy.sim(model = model, + param.as.double = param.as.double, + dose.nmol = dose.nmol, + tmax = tmax, + tau = tau, + compartment = compartment, + param.to.change = param.to.change, + param.to.change.range = param.to.change.range, + soluble = soluble) %>% + mutate(drug = drugs[i], isSol = isSoluble) + + dfs[[i.dfs]] = compare_i + } + + # Reset isSoluble to false since the default in compare.thy.sim is false. + isSoluble = FALSE +} +joined = bind_rows(dfs) + +``` + +# Plot the output AFIRT. + +```{r} +# I only wanted a subset of the parameters for gather. +joined.A = joined[c("fold.change","AFIRT.Kssd", "AFIRT.Kss", "AFIRT.Kd", "AFIRT.sim", "param","drug", "isSol")] + +# Use gather to make the long data frame for ggplot. +plots.A = joined.A %>% gather(key, AFIRT.value, -param, -fold.change,-drug, -isSol) %>% + filter(!is.na(AFIRT.value)) %>% + # Group by isSoluble to create to separate plots. + group_by(isSol)%>% + do(plot = ggplot(., aes(fold.change, AFIRT.value, color=key)) + + scale.x.log10() + + scale.y.log10() + + geom_line() + + geom_point() + + facet_grid(drug ~ param)) +plots.A$plot +``` + +# Plot the output TFIRT. + +```{r} +# I only wanted a subset of the parameters for gather. +joined.T = joined[c("fold.change","TFIRT.Kssd", "TFIRT.Kss", "TFIRT.Kd", "TFIRT.sim", "param","drug", "isSol")] + +# Use gather to make the long data frame for ggplot. +plots.T = joined.T %>% gather(key, TFIRT.value, -param, -fold.change,-drug, -isSol) %>% + filter(!is.na(TFIRT.value)) %>% + # Group by isSoluble to create to separate plots. + group_by(isSol)%>% + do(plot = ggplot(., aes(fold.change, TFIRT.value, color=key)) + + scale.x.log10() + + scale.y.log10() + + geom_line() + + geom_point() + + facet_grid(drug ~ param)) +plots.T$plot +``` + +# Create data frame for multiple drugs and parameters. +# This is redone for no target in order to compare Cmin. +# The only difference is 'target = FALSE' is the input for the model. + +```{r, include=FALSE} + +# Note: Insoluble and soluble drugs are treated differently. +# For insoluble, Tacc = M3totss/M30. For soluble, Tacc = S3totss/S30. +# We are also interested in different parameters in the SA for either case. +# See writeup for details. + +# -------------------------------------------------------------------------------- +# Initialize. +# -------------------------------------------------------------------------------- + +# Load model. +model = ivsc_4cmtct_shedct(target = FALSE) +# Drugs to explore. +drugs = c("Atezolizumab", "Herceptin", "Pembrolizumab", "Bevacizumab") +# List of parameters of interest. +parameters = c("dose", "ksynM3", "k13D", "kshedM3", "kshedDM3", "VD3") +# names(parameters) = parameters +# Create vector of soluble drugs. +soluble = "Bevacizumab" + +# Create paths to data files for each drug. +paths = NULL +for (i in 1:length(drugs)) { + paths[i] = paste0("../data/ModelF_", drugs[i],"_Params.xlsx") +} + +# Dose time, frequency, compartment, nominal dose +tmax = 26*7 #days +tau = 21 #days +compartment = 2 +dose.nmol = scale.mpk2nmol +joined = NULL +isSoluble = FALSE + +# -------------------------------------------------------------------------------- +# Iterate over all the drugs. +# -------------------------------------------------------------------------------- + +for (i in 1:length(drugs)){ + # Load parameters. + param.as.double = read.param.file(paths[i]) + df_param = as.data.frame(t(param.as.double)) + + # Check if drug has soluble target. + if (drugs[i] %in% soluble) { + # Change any non-shed parameters including M to S. + parameters = sapply(parameters, + function(x) { + if (!grepl("shed|dose",x)) {return(gsub("M","S", x))} + return(x) + }) + isSoluble = TRUE + } + + # Set range for parameters of interest in SA. + # Check which parameters are nonzero, not including dose which isn't in df_param. + nnzero = df_param[parameters[which(parameters != "dose")]] != 0 + nnzero = colnames(nnzero)[which(nnzero)] + params.to.iterate = data.frame(dose = lseq(scale.mpk2nmol*.1,scale.mpk2nmol*10,7), + lapply(df_param[nnzero], function(x) lseq(x*0.1, x*10, 7))) + + dfs = list() + # Iterate all of the parameters for a single drug. + for (j in 1:ncol(params.to.iterate)){ + i.dfs = i.dfs + 1 + dfs[[i.dfs]] = compare.thy.sim(model = model, + param.as.double = param.as.double, + dose.nmol = dose.nmol, + tmax = tmax, + tau = tau, + compartment = compartment, + param.to.change = names(params.to.iterate)[j], + param.to.change.range = params.to.iterate[[j]], + soluble = isSoluble) %>% + mutate(drug = drugs[i], isSol = isSoluble) + } + # Reset isSoluble to false since the default in compare.thy.sim is false. + isSoluble = FALSE +} +joined = bind_rows(dfs) + +``` + +# Plot the output Cmin. + +```{r} + +# I only wanted a subset of the parameters for gather. +joined.Cmin = joined[c("fold.change","Cmin", "param","drug", "isSol")] + +# Use gather to make the long data frame for ggplot. +plots.Cmin = joined.Cmin %>% gather(key, Cmin.value, -param, -fold.change,-drug, -isSol) %>% + filter(!is.na(Cmin.value)) %>% + # Group by isSoluble to create to separate plots. + group_by(isSol)%>% + do(plot = ggplot(., aes(fold.change, Cmin.value, color=key)) + + scale.x.log10() + + scale.y.log10() + + geom_line() + + geom_point() + + facet_grid(drug ~ param)) +plots.Cmin$plot + +``` + + + + + + + + + + + + + + + + + + diff --git a/ModelF/Task03_Explore_Bevacizumab.html b/ModelF/Task03_Explore_Bevacizumab.html new file mode 100644 index 0000000..63ff2a8 --- /dev/null +++ b/ModelF/Task03_Explore_Bevacizumab.html @@ -0,0 +1,487 @@ + + + + + + + + + + + + + + +Explore The Model When Soluble Target Is Present + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ + + + + + + + + + + + + + + + +
+
+
+
+
+ +
+ + + + + + + +
suppressMessages(source("ams_initialize_script.R"))
+suppressMessages(source("AFIRT_calculation.R"))
+suppressMessages(source("ivsc_4cmtct_shedct.R"))
+
+

Set initial parameters and functions

+

Hongshan says: I just assume (.1, 100) x scale.mpk2nmol is a reasonable dose Bevacizumab

+
dose.nmol.range = lseq(.1,100,7)*scale.mpk2nmol
+
+model = ivsc_4cmtct_shedct()
+
+tmax = 26*7 #days
+tau  = 21   #days
+compartment = 2
+
+
+

Import parameters for Bevacizumab

+
param.as.double =  read.param.file("../data/ModelF_Bevacizumab_Params.xlsx", model=model)
+
## Warning in read.param.file("../data/ModelF_Bevacizumab_Params.xlsx", model
+## = model): NAs introduced by coercion
+
pdb.bevacizumab = param.as.double
+p.bevacizumab = as.data.frame(t(param.as.double))
+p.bevacizumab$drug = "Bevacizumab"
+
+df_compare = compare.theory.sim(model = model,
+                                param.as.double = param.as.double,
+                                dose.nmol.range = dose.nmol.range,
+                                tmax = tmax,
+                                tau  = tau,
+                                compartment = compartment)
+
##         type        Dose M30 Mtot3.ss        S30  Stot3.ss Tacc.tum.M
+## 1 simulation    46.66667   0        0 0.01053504 0.1092132        NaN
+## 2 simulation   147.57296   0        0 0.01053504 0.1319732        NaN
+## 3 simulation   466.66667   0        0 0.01053504 0.1418579        NaN
+## 4 simulation  1475.72957   0        0 0.01053504 0.1453674        NaN
+## 5 simulation  4666.66667   0        0 0.01053504 0.1465208        NaN
+## 6 simulation 14757.29575   0        0 0.01053504 0.1468901        NaN
+## 7 simulation 46666.66667   0        0 0.01053504 0.1470073        NaN
+##   Tacc.tum.S       Cavg1       Cavg3         B AFIRT.M      AFIRT.S
+## 1   10.36667    10.56622    3.432662 0.3248715     NaN 0.1944385835
+## 2   12.52707    33.69863   11.064274 0.3283301     NaN 0.0727659277
+## 3   13.46534   106.99189   35.273053 0.3296797     NaN 0.0244848180
+## 4   13.79847   338.83469  111.862785 0.3301397     NaN 0.0079051795
+## 5   13.90795  1072.01157  354.073506 0.3302889     NaN 0.0025166066
+## 6   13.94300  3390.52880 1120.015163 0.3303364     NaN 0.0007975159
+## 7   13.95413 10722.32719 3542.136713 0.3303515     NaN 0.0002523666
+
## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character
+
## Warning in bind_rows_(x, .id): binding character and factor vector,
+## coercing into character vector
+
+## Warning in bind_rows_(x, .id): binding character and factor vector,
+## coercing into character vector
+
##          type        Dose M30 Mtot3.ss        S30  Stot3.ss Tacc.tum.S
+## 1      theory    46.66667   0        0 0.01053504 0.1470626   13.95938
+## 2      theory   147.57296   0        0 0.01053504 0.1470626   13.95938
+## 3      theory   466.66667   0        0 0.01053504 0.1470626   13.95938
+## 4      theory  1475.72957   0        0 0.01053504 0.1470626   13.95938
+## 5      theory  4666.66667   0        0 0.01053504 0.1470626   13.95938
+## 6      theory 14757.29575   0        0 0.01053504 0.1470626   13.95938
+## 7      theory 46666.66667   0        0 0.01053504 0.1470626   13.95938
+## 8  simulation    46.66667   0        0 0.01053504 0.1092132   10.36667
+## 9  simulation   147.57296   0        0 0.01053504 0.1319732   12.52707
+## 10 simulation   466.66667   0        0 0.01053504 0.1418579   13.46534
+## 11 simulation  1475.72957   0        0 0.01053504 0.1453674   13.79847
+## 12 simulation  4666.66667   0        0 0.01053504 0.1465208   13.90795
+## 13 simulation 14757.29575   0        0 0.01053504 0.1468901   13.94300
+## 14 simulation 46666.66667   0        0 0.01053504 0.1470073   13.95413
+##          Cavg1       Cavg3         B      AFIRT.S
+## 1     12.34568    4.115226 0.3333333 0.3208844543
+## 2     39.04046   13.013488 0.3333333 0.1014725741
+## 3    123.45679   41.152263 0.3333333 0.0320884454
+## 4    390.40465  130.134883 0.3333333 0.0101472574
+## 5   1234.56790  411.522634 0.3333333 0.0032088445
+## 6   3904.04649 1301.348831 0.3333333 0.0010147257
+## 7  12345.67901 4115.226337 0.3333333 0.0003208845
+## 8     10.56622    3.432662 0.3248715 0.1944385835
+## 9     33.69863   11.064274 0.3283301 0.0727659277
+## 10   106.99189   35.273053 0.3296797 0.0244848180
+## 11   338.83469  111.862785 0.3301397 0.0079051795
+## 12  1072.01157  354.073506 0.3302889 0.0025166066
+## 13  3390.52880 1120.015163 0.3303364 0.0007975159
+## 14 10722.32719 3542.136713 0.3303515 0.0002523666
+
## Warning: Transformation introduced infinite values in continuous y-axis
+

+
##           type  Dose      param   value
+## 1   simulation    47        M30 0.0e+00
+## 2       theory    47        M30 0.0e+00
+## 3   simulation   150        M30 0.0e+00
+## 4       theory   150        M30 0.0e+00
+## 5   simulation   470        M30 0.0e+00
+## 6       theory   470        M30 0.0e+00
+## 7   simulation  1500        M30 0.0e+00
+## 8       theory  1500        M30 0.0e+00
+## 9   simulation  4700        M30 0.0e+00
+## 10      theory  4700        M30 0.0e+00
+## 11  simulation 15000        M30 0.0e+00
+## 12      theory 15000        M30 0.0e+00
+## 13  simulation 47000        M30 0.0e+00
+## 14      theory 47000        M30 0.0e+00
+## 15  simulation    47   Mtot3.ss 0.0e+00
+## 16      theory    47   Mtot3.ss 0.0e+00
+## 17  simulation   150   Mtot3.ss 0.0e+00
+## 18      theory   150   Mtot3.ss 0.0e+00
+## 19  simulation   470   Mtot3.ss 0.0e+00
+## 20      theory   470   Mtot3.ss 0.0e+00
+## 21  simulation  1500   Mtot3.ss 0.0e+00
+## 22      theory  1500   Mtot3.ss 0.0e+00
+## 23  simulation  4700   Mtot3.ss 0.0e+00
+## 24      theory  4700   Mtot3.ss 0.0e+00
+## 25  simulation 15000   Mtot3.ss 0.0e+00
+## 26      theory 15000   Mtot3.ss 0.0e+00
+## 27  simulation 47000   Mtot3.ss 0.0e+00
+## 28      theory 47000   Mtot3.ss 0.0e+00
+## 29  simulation    47        S30 1.1e-02
+## 30      theory    47        S30 1.1e-02
+## 31  simulation   150        S30 1.1e-02
+## 32      theory   150        S30 1.1e-02
+## 33  simulation   470        S30 1.1e-02
+## 34      theory   470        S30 1.1e-02
+## 35  simulation  1500        S30 1.1e-02
+## 36      theory  1500        S30 1.1e-02
+## 37  simulation  4700        S30 1.1e-02
+## 38      theory  4700        S30 1.1e-02
+## 39  simulation 15000        S30 1.1e-02
+## 40      theory 15000        S30 1.1e-02
+## 41  simulation 47000        S30 1.1e-02
+## 42      theory 47000        S30 1.1e-02
+## 43  simulation    47   Stot3.ss 1.1e-01
+## 44      theory    47   Stot3.ss 1.5e-01
+## 45  simulation   150   Stot3.ss 1.3e-01
+## 46      theory   150   Stot3.ss 1.5e-01
+## 47  simulation   470   Stot3.ss 1.4e-01
+## 48      theory   470   Stot3.ss 1.5e-01
+## 49  simulation  1500   Stot3.ss 1.5e-01
+## 50      theory  1500   Stot3.ss 1.5e-01
+## 51  simulation  4700   Stot3.ss 1.5e-01
+## 52      theory  4700   Stot3.ss 1.5e-01
+## 53  simulation 15000   Stot3.ss 1.5e-01
+## 54      theory 15000   Stot3.ss 1.5e-01
+## 55  simulation 47000   Stot3.ss 1.5e-01
+## 56      theory 47000   Stot3.ss 1.5e-01
+## 57  simulation    47 Tacc.tum.S 1.0e+01
+## 58      theory    47 Tacc.tum.S 1.4e+01
+## 59  simulation   150 Tacc.tum.S 1.3e+01
+## 60      theory   150 Tacc.tum.S 1.4e+01
+## 61  simulation   470 Tacc.tum.S 1.3e+01
+## 62      theory   470 Tacc.tum.S 1.4e+01
+## 63  simulation  1500 Tacc.tum.S 1.4e+01
+## 64      theory  1500 Tacc.tum.S 1.4e+01
+## 65  simulation  4700 Tacc.tum.S 1.4e+01
+## 66      theory  4700 Tacc.tum.S 1.4e+01
+## 67  simulation 15000 Tacc.tum.S 1.4e+01
+## 68      theory 15000 Tacc.tum.S 1.4e+01
+## 69  simulation 47000 Tacc.tum.S 1.4e+01
+## 70      theory 47000 Tacc.tum.S 1.4e+01
+## 71  simulation    47      Cavg1 1.1e+01
+## 72      theory    47      Cavg1 1.2e+01
+## 73  simulation   150      Cavg1 3.4e+01
+## 74      theory   150      Cavg1 3.9e+01
+## 75  simulation   470      Cavg1 1.1e+02
+## 76      theory   470      Cavg1 1.2e+02
+## 77  simulation  1500      Cavg1 3.4e+02
+## 78      theory  1500      Cavg1 3.9e+02
+## 79  simulation  4700      Cavg1 1.1e+03
+## 80      theory  4700      Cavg1 1.2e+03
+## 81  simulation 15000      Cavg1 3.4e+03
+## 82      theory 15000      Cavg1 3.9e+03
+## 83  simulation 47000      Cavg1 1.1e+04
+## 84      theory 47000      Cavg1 1.2e+04
+## 85  simulation    47      Cavg3 3.4e+00
+## 86      theory    47      Cavg3 4.1e+00
+## 87  simulation   150      Cavg3 1.1e+01
+## 88      theory   150      Cavg3 1.3e+01
+## 89  simulation   470      Cavg3 3.5e+01
+## 90      theory   470      Cavg3 4.1e+01
+## 91  simulation  1500      Cavg3 1.1e+02
+## 92      theory  1500      Cavg3 1.3e+02
+## 93  simulation  4700      Cavg3 3.5e+02
+## 94      theory  4700      Cavg3 4.1e+02
+## 95  simulation 15000      Cavg3 1.1e+03
+## 96      theory 15000      Cavg3 1.3e+03
+## 97  simulation 47000      Cavg3 3.5e+03
+## 98      theory 47000      Cavg3 4.1e+03
+## 99  simulation    47          B 3.2e-01
+## 100     theory    47          B 3.3e-01
+## 101 simulation   150          B 3.3e-01
+## 102     theory   150          B 3.3e-01
+## 103 simulation   470          B 3.3e-01
+## 104     theory   470          B 3.3e-01
+## 105 simulation  1500          B 3.3e-01
+## 106     theory  1500          B 3.3e-01
+## 107 simulation  4700          B 3.3e-01
+## 108     theory  4700          B 3.3e-01
+## 109 simulation 15000          B 3.3e-01
+## 110     theory 15000          B 3.3e-01
+## 111 simulation 47000          B 3.3e-01
+## 112     theory 47000          B 3.3e-01
+## 113 simulation    47    AFIRT.S 1.9e-01
+## 114     theory    47    AFIRT.S 3.2e-01
+## 115 simulation   150    AFIRT.S 7.3e-02
+## 116     theory   150    AFIRT.S 1.0e-01
+## 117 simulation   470    AFIRT.S 2.4e-02
+## 118     theory   470    AFIRT.S 3.2e-02
+## 119 simulation  1500    AFIRT.S 7.9e-03
+## 120     theory  1500    AFIRT.S 1.0e-02
+## 121 simulation  4700    AFIRT.S 2.5e-03
+## 122     theory  4700    AFIRT.S 3.2e-03
+## 123 simulation 15000    AFIRT.S 8.0e-04
+## 124     theory 15000    AFIRT.S 1.0e-03
+## 125 simulation 47000    AFIRT.S 2.5e-04
+## 126     theory 47000    AFIRT.S 3.2e-04
+
## Warning: Transformation introduced infinite values in continuous y-axis
+
## Warning: Transformation introduced infinite values in continuous y-axis
+

+
+ + + +
+
+ +
+ + + + + + + + diff --git a/ModelF/Task04_Sensitivity_Analysis.html b/ModelF/Task04_Sensitivity_Analysis.html new file mode 100644 index 0000000..f8c81b0 --- /dev/null +++ b/ModelF/Task04_Sensitivity_Analysis.html @@ -0,0 +1,302 @@ + + + + + + + + + + + + + + +Task04_Sensitivity_Analysis + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ + + + + + + + + + + + + + +
+

Preliminary stuff

+
# Need this to get rxode working. I don't know why I keep having to run this.
+Sys.setenv(PATH = paste(Sys.getenv("PATH"), "C:/RBuildTools/3.4/bin/",
+            "C:/RBuildTools/3.4/mingw_64/bin", sep = ";"))
+Sys.setenv(BINPREF = "C:/RBuildTools/3.4/mingw_64/bin/")
+
+# To be called at the top of every Rmd file. Initialization code and some useful constants.
+suppressMessages(source("ams_initialize_script.R"))
+
+
+

Function to compute sensitivity analysis

+
# This function does the sensitivity analysis on the user inputted parameter and compares the theoretical result to the simulated result.
+
+# Input: 
+# model - model system of ODE's solved with RxODE. In this project, it is 'ivsc_4cmtct_shedct'.
+# param.as.double - read parameters from Excel file. read.param.file("file directory").
+# dose.nmol - dose in nmol.
+# tmax - time of treatment in days
+# tau - frequency of administering dose in days
+# compartment - compartment where drug is administered
+# param.to.change - parameter on which to do SA. This must be a string.
+# param.to.change.range - range of parameter on which to do SA. The range must be symmetric in fold change. This must be a vector of odd length.
+
+# Output:
+# Plot of theoretical and simulated AFIRT vs fold change in parameter
+# Data frame of AFIRT vs parameter value
+
+compare.thy.sim = function(model = model, 
+                           param.as.double = param.as.double,
+                           dose.nmol = dose.nmol,
+                           tmax = tmax,
+                           tau = tau,
+                           compartment = compartment,
+                           param.to.change = param.to.change,
+                           param.to.change.range = param.to.change.range) {
+  
+  # Evaluate model ------------------------------------------------------------------------------------------
+
+  # Put all simulations for different dose into one data frame.
+  df_sim = data.frame()
+  
+  # Iterate through values in simulation.
+  if (param.to.change == 'dose'){
+    for (param.iter in param.to.change.range){
+      row = lumped.parameters.simulation(model, param.as.double, param.iter, tmax, tau, compartment)
+      df_sim = rbind(df_sim, row)
+    }  
+  } else {
+    for (param.iter in param.to.change.range){
+    param.as.double[param.to.change] = param.iter
+    row = lumped.parameters.simulation(model, param.as.double, dose.nmol, tmax, tau, compartment)
+    df_sim = rbind(df_sim, row)
+    }
+  }
+  df_sim$param.to.change = param.to.change.range
+  
+  # Put all theoretical calculations for different dose into one data frame.
+  df_thy = data.frame()
+  
+  # Iterate through values in theory.
+  if (param.to.change == 'dose'){
+    for (param.iter in param.to.change.range){
+      row = lumped.parameters.theory(param.as.double, param.iter, tau)
+      df_thy = rbind(df_thy, row)
+    }
+  } else{
+    for (param.iter in param.to.change.range){
+      param.as.double[param.to.change] = param.iter
+      row = lumped.parameters.theory(param.as.double, dose.nmol, tau)
+      df_thy = rbind(df_thy, row)
+    }
+  }
+  df_thy = df_thy %>% mutate(param.to.change = param.to.change.range)
+  
+  # Arrange theory and simulation into single data frame.
+  df_compare = bind_rows(df_thy,df_sim)
+  df_compare = df_compare %>%
+    arrange(param.to.change,type) %>%
+    mutate_if(is.numeric,signif,2)
+  
+  # Plot results --------------------------------------------------------------------------------------------
+  
+  # Arrange data in a way suitable for plotting.
+  Param  = data.frame(Param = param.to.change.range/median(param.to.change.range))
+  AFIRT  = data.frame(AFIRT.sim  = df_sim$AFIRT.sim,
+                      AFIRT.Kssd = df_thy$AFIRT.Kssd,
+                      AFIRT.Kss  = df_thy$AFIRT.Kss,
+                      AFIRT.Kd   = df_thy$AFIRT.Kd)
+  Data  = data.frame(Param, AFIRT)
+  names = names(Data)
+  Data  = Data %>% gather(key, AFIRT.value, -c(get(names[1])))
+
+  # Plot all AFIRT's on same axes.
+  g = ggplot(Data, aes(Param, AFIRT.value, color=key)) +
+      # ylim(10^(-4), 1) + 
+      scale.x.log10() +
+      scale.y.log10(limits=c(10^(-3), 10^(1))) + 
+      geom_line() +
+      geom_point() + 
+      labs(x = param.to.change)
+  print(g)
+  
+  return(df_compare)
+}
+
+
+

Example of calling this function

+
# Load model and parameters.
+model           = ivsc_4cmtct_shedct()
+param.as.double = read.param.file("../data/ModelF_Atezolizumab_Params.xlsx")
+df_param        = as.data.frame(t(param.as.double))
+dose.nmol       = scale.mpk2nmol
+
+# Set range for parameters of interest in SA.
+dose.nmol.range = lseq(.1,10,7)*scale.mpk2nmol
+ksynM3.range    = lseq(df_param$ksynM3  *.1, df_param$ksynM3  *10, 7)
+kshedM3.range   = lseq(df_param$kshedM3 *.1, df_param$kshedM3 *10, 7)
+kshedDM3.range  = lseq(df_param$kshedDM3*.1, df_param$kshedDM3*10, 7)
+k13D.range      = lseq(df_param$k13D    *.1, df_param$k13D    *10, 7)
+# k13DM.range     = lseq(df_param$k13DM   *.1, df_param$k13DM   *10, 7)
+# k31DM.range     = lseq(df_param$k31DM   *.1, df_param$k31DM   *10, 7)
+VD3.range       = lseq(df_param$VD3     *.1, df_param$VD3     *10, 7)
+
+# Dose time, frequency, and compartment
+tmax = 26*7 #days
+tau  = 21   #days
+compartment = 2
+
+df_compare = compare.thy.sim(model = model,
+                             param.as.double = param.as.double,
+                             dose.nmol = dose.nmol,
+                             tmax = tmax,
+                             tau  = tau,
+                             compartment = compartment,
+                             param.to.change = 'kshedDM3',
+                             param.to.change.range = kshedDM3.range)
+
## Warning: package 'bindrcpp' was built under R version 3.2.5
+
## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character
+
## Warning in bind_rows_(x, .id): binding character and factor vector,
+## coercing into character vector
+
+## Warning in bind_rows_(x, .id): binding character and factor vector,
+## coercing into character vector
+

+
+ + + + +
+ + + + + + + + diff --git a/ModelF/Task04_Sensitivity_Analysis_Individual.html b/ModelF/Task04_Sensitivity_Analysis_Individual.html new file mode 100644 index 0000000..954c2f1 --- /dev/null +++ b/ModelF/Task04_Sensitivity_Analysis_Individual.html @@ -0,0 +1,224 @@ + + + + + + + + + + + + + + +Task04_Sensitivity_Analysis + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ + + + + + + + + + + + + + +
+

Preliminary stuff

+
# Need this to get rxode working. I don't know why I keep having to run this.
+Sys.setenv(PATH = paste(Sys.getenv("PATH"), "C:/RBuildTools/3.4/bin/",
+            "C:/RBuildTools/3.4/mingw_64/bin", sep = ";"))
+Sys.setenv(BINPREF = "C:/RBuildTools/3.4/mingw_64/bin/")
+
+# To be called at the top of every Rmd file. Initialization code and some useful constants.
+suppressMessages(source("ams_initialize_script.R"))
+
+
+

Individual sensitivity analysis

+
# Load model and parameters.
+model           = ivsc_4cmtct_shedct()
+drugs           = c("Atezolizumab", "Trastuzumab", "Pembrolizumab", "Bevacizumab"); paths = NULL; i = 4
+paths[i]        = paste0("../data/ModelF_", drugs[i],"_Params.xlsx")
+param.as.double = read.param.file(paths[i])
+param.as.double["keS3"] = 1
+df_param        = as.data.frame(t(param.as.double))
+dose.nmol       = scale.mpk2nmol
+
+# Soluble parameters require different treatment.
+soluble = "Bevacizumab"
+if (drugs[i] %in% soluble) {
+  isSoluble = TRUE
+  } else {
+  isSoluble = FALSE
+}
+
+# Set parameter of interest and range for SA.
+dose.nmol.range = lseq(.1,10,7)*scale.mpk2nmol
+param.to.change = 'keS3'
+param.to.change.range = lseq(as.numeric(df_param[param.to.change])*.001,   
+                             as.numeric(df_param[param.to.change])*1000, 
+                             7)
+
+# Dose time, frequency, compartment
+tmax = 26*7 #days
+tau  = 21   #days
+compartment = 2
+
+# Create data frame containing variables of interest for chaning parameter of interest.
+df_compare = compare.thy.sim(model = model,
+                             param.as.double = param.as.double,
+                             dose.nmol = dose.nmol,
+                             tmax = tmax,
+                             tau  = tau,
+                             compartment = compartment,
+                             param.to.change = param.to.change,
+                             param.to.change.range = param.to.change.range,
+                             soluble = isSoluble)
+df_compare = df_compare %>% mutate(drug = drugs[i], isSol = isSoluble)
+
+# Retrieve simulation and theory only for individual plotting.
+df_sim = subset(df_compare, type=='simulation')
+df_thy = subset(df_compare, type=='theory')
+
+# Plot.
+ggplot(df_sim, aes(x=Tacc.tum, y=AFIRT.sim)) + 
+  geom_line() + 
+  scale.x.log10() +
+  scale.y.log10() +
+  # scale.y.log10(limits=c(10^-4, 10^4)) +
+  # labs(x=param.to.change, y='AFIRT')
+  labs(x='Target Accumulation (Tfold)', y='Average Free Tissue target\nto Inital target Ratio (AFTIR)\nby simulation ')
+

+
ggsave("../results/Task04_Taccum.pdf",width = 4,height=3)
+
+ + + + +
+ + + + + + + + diff --git a/ModelF/Task05_Sensitivity_Analysis_Aggregate.Rmd b/ModelF/Task05_Sensitivity_Analysis_Aggregate.Rmd index b87038f..5f21c48 100644 --- a/ModelF/Task05_Sensitivity_Analysis_Aggregate.Rmd +++ b/ModelF/Task05_Sensitivity_Analysis_Aggregate.Rmd @@ -149,7 +149,7 @@ plots.A.Kssd = plots.A.Kssd %>% param = str_replace(param,"S","T"), key = str_replace(key,"_",""), key = str_replace(key,"AFIRT\\.",""), - key = str_replace(key,"Kssd","Theory (AFTIR*)"), + key = str_replace(key,"Kssd","Theory (AFTIR)"), key = str_replace(key,"sim","Simulation"), Target = ifelse(isSol,"Soluble","Membrane-Bd"), Target = factor(Target,levels=c("Soluble","Membrane-Bd"))) @@ -158,7 +158,7 @@ plots.T.Kssd = plots.T.Kssd %>% param = str_replace(param,"S","T"), key = str_replace(key,"_",""), key = str_replace(key,"TFIRT\\.",""), - key = str_replace(key,"Kssd","Theory (TFTIR*)"), + key = str_replace(key,"Kssd","Theory (TFTIR)"), key = str_replace(key,"sim","Simulation"), Target = ifelse(isSol,"Soluble","Membrane-Bd"), Target = factor(Target,levels=c("Soluble","Membrane-Bd"))) @@ -168,7 +168,7 @@ plots.A.Kss = plots.A.Kss %>% param = str_replace(param,"S","T"), key = str_replace(key,"_",""), key = str_replace(key,"AFIRT\\.",""), - key = str_replace(key,"Kss","Theory (AFTIR*)"), + key = str_replace(key,"Kss","Theory (AFTIR)"), key = str_replace(key,"sim","Simulation"), Target = ifelse(isSol,"Soluble","Membrane-Bd"), Target = factor(Target,levels=c("Soluble","Membrane-Bd"))) @@ -177,7 +177,7 @@ plots.T.Kss = plots.T.Kss %>% param = str_replace(param,"S","T"), key = str_replace(key,"_",""), key = str_replace(key,"TFIRT\\.",""), - key = str_replace(key,"Kss","Theory (TFTIR*)"), + key = str_replace(key,"Kss","Theory (TFTIR)"), key = str_replace(key,"sim","Simulation"), Target = ifelse(isSol,"Soluble","Membrane-Bd"), Target = factor(Target,levels=c("Soluble","Membrane-Bd"))) @@ -187,7 +187,7 @@ plots.A.Kd = plots.A.Kd %>% param = str_replace(param,"S","T"), key = str_replace(key,"_",""), key = str_replace(key,"AFIRT\\.",""), - key = str_replace(key,"Kd","Theory (AFTIR*)"), + key = str_replace(key,"Kd","Theory (AFTIR)"), key = str_replace(key,"sim","Simulation"), Target = ifelse(isSol,"Soluble","Membrane-Bd"), Target = factor(Target,levels=c("Soluble","Membrane-Bd"))) @@ -196,7 +196,7 @@ plots.T.Kd = plots.T.Kd %>% param = str_replace(param,"S","T"), key = str_replace(key,"_",""), key = str_replace(key,"TFIRT\\.",""), - key = str_replace(key,"Kd","Theory (TFTIR*)"), + key = str_replace(key,"Kd","Theory (TFTIR)"), key = str_replace(key,"sim","Simulation"), Target = ifelse(isSol,"Soluble","Membrane-Bd"), Target = factor(Target,levels=c("Soluble","Membrane-Bd"))) @@ -208,10 +208,9 @@ g = g + geom_line() g = g + facet_grid(drug + Target ~ param, switch = "y") g = g + scale_linetype_manual(values = c("solid","longdash")) g = g + labs(x="Fold Change in Parameter", - y="Average Free Tissue target to\nInitial target Ratio (AFTIR)", + y="Average Free Tissue target to\nInitial target Ratio (AFTIR)") # color = "Estimate", - # linetype = "Estimate", - caption = "Kssd") + # linetype = "Estimate") g = g + theme(legend.position="bottom", text=element_text(size=15), axis.title.x=element_text(size=20), @@ -224,30 +223,25 @@ height = 7 plot(g) ggsave("../results/Task05_AFTIR.Kssd.pdf",width = width,height=height) g = g %+% plots.T.Kssd -g = g + labs(y="Trough Free Tissue target to\ninitial target Ratio (TFTIR)", - caption = "Kssd") +g = g + labs(y="Trough Free Tissue target to\ninitial target Ratio (TFTIR)") plot(g) ggsave("../results/Task05_TFTIR.Kssd.pdf",width = width,height=height) g = g %+% plots.A.Kss -g = g + labs(y="Average Free Tissue target to\nInitial target Ratio (AFTIR)", - caption = "Kss") +g = g + labs(y="Average Free Tissue target to\nInitial target Ratio (AFTIR)") plot(g) ggsave("../results/Task05_AFTIR.Kss.pdf",width = width,height=height) g = g %+% plots.T.Kss -g = g + labs(y="Trough Free Tissue target to\ninitial target Ratio (TFTIR)", - caption = "Kss") +g = g + labs(y="Trough Free Tissue target to\ninitial target Ratio (TFTIR)") plot(g) ggsave("../results/Task05_TFTIR.Kss.pdf",width = width,height=height) g = g %+% plots.A.Kd -g = g + labs(y="Average Free Tissue target to\nInitial target Ratio (AFTIR)", - caption = "Kd") +g = g + labs(y="Average Free Tissue target to\nInitial target Ratio (AFTIR)") plot(g) ggsave("../results/Task05_AFTIR.Kd.pdf",width = width,height=height) g = g %+% plots.T.Kd -g = g + labs(y="Trough Free Tissue target to\ninitial target Ratio (TFTIR)", - caption = "Kd") +g = g + labs(y="Trough Free Tissue target to\ninitial target Ratio (TFTIR)") plot(g) ggsave("../results/Task05_TFTIR.Kd.pdf",width = width,height=height) ``` diff --git a/ModelF/Task05_Sensitivity_Analysis_Aggregate.html b/ModelF/Task05_Sensitivity_Analysis_Aggregate.html index b2dd28c..4f8107c 100644 --- a/ModelF/Task05_Sensitivity_Analysis_Aggregate.html +++ b/ModelF/Task05_Sensitivity_Analysis_Aggregate.html @@ -233,7 +233,7 @@

Task05_Sensitivity_Analysis_pt2

Miandra Ellis

-

29 June, 2018

+

01 August, 2018

@@ -254,55 +254,144 @@

Create data frame for multiple drugs and parameters.

Plot the output variable of interest

# I only wanted a subset of the parameters for gather.
-joined.A = joined[c("fold.change.param","AFIRT.Kssd_", "AFIRT.Kssd", "AFIRT.sim", "param","drug", "isSol")]
-joined.T = joined[c("fold.change.param","TFIRT.Kssd_", "TFIRT.Kssd", "TFIRT.sim", "param","drug", "isSol")]
+joined.A.Kssd = joined[c("fold.change.param", "AFIRT.Kssd_", "AFIRT.sim", "param", "drug", "isSol", "M30", "S30")]
+joined.T.Kssd = joined[c("fold.change.param", "TFIRT.Kssd_", "TFIRT.sim", "param", "drug", "isSol", "M30", "S30")]
+
+joined.A.Kss = joined[c("fold.change.param", "AFIRT.Kss_", "AFIRT.sim", "param", "drug", "isSol", "M30", "S30")]
+joined.T.Kss = joined[c("fold.change.param", "TFIRT.Kss_", "TFIRT.sim", "param", "drug", "isSol", "M30", "S30")]
+
+joined.A.Kd = joined[c("fold.change.param", "AFIRT.Kd_", "AFIRT.sim", "param", "drug", "isSol", "M30", "S30")]
+joined.T.Kd = joined[c("fold.change.param", "TFIRT.Kd_", "TFIRT.sim", "param", "drug", "isSol", "M30", "S30")]
 
 # Use gather to make the long data frame for ggplot. Group by isSoluble to create to separate plots.
-plots.A = joined.A %>% gather(key, Value, -param, -fold.change.param, -drug, -isSol) %>%
+plots.A.Kssd = joined.A.Kssd %>% gather(key, Value, -param, -fold.change.param, -drug, -isSol, -M30, -S30) %>%
   filter(!is.na(Value)) %>%
-  group_by(isSol) 
+  group_by(isSol)
+plots.T.Kssd = joined.T.Kssd %>% gather(key, Value, -param, -fold.change.param, -drug, -isSol, -M30, -S30) %>%
+  filter(!is.na(Value)) %>%
+  group_by(isSol)
 
-plots.T = joined.T %>% gather(key, Value, -param, -fold.change.param, -drug, -isSol) %>%
+plots.A.Kss = joined.A.Kss %>% gather(key, Value, -param, -fold.change.param, -drug, -isSol, -M30, -S30) %>%
+  filter(!is.na(Value)) %>%
+  group_by(isSol)
+plots.T.Kss = joined.T.Kss %>% gather(key, Value, -param, -fold.change.param, -drug, -isSol, -M30, -S30) %>%
+  filter(!is.na(Value)) %>%
+  group_by(isSol)
+
+plots.A.Kd = joined.A.Kd %>% gather(key, Value, -param, -fold.change.param, -drug, -isSol, -M30, -S30) %>%
+  filter(!is.na(Value)) %>%
+  group_by(isSol)
+plots.T.Kd = joined.T.Kd %>% gather(key, Value, -param, -fold.change.param, -drug, -isSol, -M30, -S30) %>%
   filter(!is.na(Value)) %>%
-  group_by(isSol) 
+  group_by(isSol)
 
 #Rename for plotting - map M->T and S->T (for target)
-plots.A = plots.A %>%
-  mutate(param = str_replace(param,"M","T"),
-         param = str_replace(param,"S","T"),
-         key   = str_replace(key,"_","*"),
-         key   = str_replace(key,"AFIRT\\.",""),
-         key   = str_replace(key,"Kssd","theory.Kssd"),
-         key   = str_replace(key,"sim","simulation"))
-plots.T = plots.T %>%
-  mutate(param = str_replace(param,"M","T"),
-         param = str_replace(param,"S","T"),
-         key   = str_replace(key,"_","*"),
-         key   = str_replace(key,"TFIRT\\.",""),
-         key   = str_replace(key,"Kssd","theory.Kssd"),
-         key   = str_replace(key,"sim","simulation"))
-
-g = ggplot(plots.A, aes(fold.change.param, Value, color=key,linetype=key))
-g = g + scale.x.log10(limits=c(.003,30))
-g = g + scale.y.log10()
+plots.A.Kssd = plots.A.Kssd %>%
+  mutate(param  = str_replace(param,"M","T"),
+         param  = str_replace(param,"S","T"),
+         key    = str_replace(key,"_",""),
+         key    = str_replace(key,"AFIRT\\.",""),
+         key    = str_replace(key,"Kssd","Theory (AFTIR)"),
+         key    = str_replace(key,"sim","Simulation"),
+         Target = ifelse(isSol,"Soluble","Membrane-Bd"),
+         Target = factor(Target,levels=c("Soluble","Membrane-Bd")))
+plots.T.Kssd = plots.T.Kssd %>%
+  mutate(param  = str_replace(param,"M","T"),
+         param  = str_replace(param,"S","T"),
+         key    = str_replace(key,"_",""),
+         key    = str_replace(key,"TFIRT\\.",""),
+         key    = str_replace(key,"Kssd","Theory (TFTIR)"),
+         key    = str_replace(key,"sim","Simulation"),
+         Target = ifelse(isSol,"Soluble","Membrane-Bd"),
+         Target = factor(Target,levels=c("Soluble","Membrane-Bd")))
+
+plots.A.Kss = plots.A.Kss %>%
+  mutate(param  = str_replace(param,"M","T"),
+         param  = str_replace(param,"S","T"),
+         key    = str_replace(key,"_",""),
+         key    = str_replace(key,"AFIRT\\.",""),
+         key    = str_replace(key,"Kss","Theory (AFTIR)"),
+         key    = str_replace(key,"sim","Simulation"),
+         Target = ifelse(isSol,"Soluble","Membrane-Bd"),
+         Target = factor(Target,levels=c("Soluble","Membrane-Bd")))
+plots.T.Kss = plots.T.Kss %>%
+  mutate(param  = str_replace(param,"M","T"),
+         param  = str_replace(param,"S","T"),
+         key    = str_replace(key,"_",""),
+         key    = str_replace(key,"TFIRT\\.",""),
+         key    = str_replace(key,"Kss","Theory (TFTIR)"),
+         key    = str_replace(key,"sim","Simulation"),
+         Target = ifelse(isSol,"Soluble","Membrane-Bd"),
+         Target = factor(Target,levels=c("Soluble","Membrane-Bd")))
+
+plots.A.Kd = plots.A.Kd %>%
+  mutate(param  = str_replace(param,"M","T"),
+         param  = str_replace(param,"S","T"),
+         key    = str_replace(key,"_",""),
+         key    = str_replace(key,"AFIRT\\.",""),
+         key    = str_replace(key,"Kd","Theory (AFTIR)"),
+         key    = str_replace(key,"sim","Simulation"),
+         Target = ifelse(isSol,"Soluble","Membrane-Bd"),
+         Target = factor(Target,levels=c("Soluble","Membrane-Bd")))
+plots.T.Kd = plots.T.Kd %>%
+  mutate(param  = str_replace(param,"M","T"),
+         param  = str_replace(param,"S","T"),
+         key    = str_replace(key,"_",""),
+         key    = str_replace(key,"TFIRT\\.",""),
+         key    = str_replace(key,"Kd","Theory (TFTIR)"),
+         key    = str_replace(key,"sim","Simulation"),
+         Target = ifelse(isSol,"Soluble","Membrane-Bd"),
+         Target = factor(Target,levels=c("Soluble","Membrane-Bd")))
+
+g = ggplot(plots.A.Kssd, aes(fold.change.param, Value, color=key,linetype=key))
+g = g + scale.x.log10(limits=c(.009,11))
+g = g + scale.y.log10(limits=c(.001,1))
 g = g + geom_line()
-g = g + facet_grid(drug ~ param,switch = "y")
-g = g + scale_linetype_manual(values  = c("solid","longdash","dashed"))
+g = g + facet_grid(drug + Target ~ param, switch = "y")
+g = g + scale_linetype_manual(values  = c("solid","longdash"))
 g = g + labs(x="Fold Change in Parameter",
-             y="AFTIR",
-             color = "Estimate",
-             linetype = "Estimate")
-width = 8
-height = 5
+             y="Average Free Tissue target to\nInitial target Ratio (AFTIR)")
+             # color = "Estimate",
+             # linetype = "Estimate")
+g = g + theme(legend.position="bottom", 
+              text=element_text(size=15), 
+              axis.title.x=element_text(size=20), 
+              axis.title.y=element_text(size=20))
+g = g + scale_color_manual(values = c("blue","black"))
+
+width  = 10
+height = 7
+
+plot(g)
+

+
ggsave("../results/Task05_AFTIR.Kssd.pdf",width = width,height=height)
+g = g %+% plots.T.Kssd
+g = g + labs(y="Trough Free Tissue target to\ninitial target Ratio (TFTIR)")
+plot(g)
+

+
ggsave("../results/Task05_TFTIR.Kssd.pdf",width = width,height=height)
+
+g = g %+% plots.A.Kss
+g = g + labs(y="Average Free Tissue target to\nInitial target Ratio (AFTIR)")
+plot(g)
+

+
ggsave("../results/Task05_AFTIR.Kss.pdf",width = width,height=height)
+g = g %+% plots.T.Kss
+g = g + labs(y="Trough Free Tissue target to\ninitial target Ratio (TFTIR)")
 plot(g)
-

-
ggsave("../results/Task05_AFTIR.pdf",width = width,height=height)
+

+
ggsave("../results/Task05_TFTIR.Kss.pdf",width = width,height=height)
 
-g = g %+% plots.T
-g = g + labs(y="TFTIR")
+g = g %+% plots.A.Kd
+g = g + labs(y="Average Free Tissue target to\nInitial target Ratio (AFTIR)")
+plot(g)

+

+
ggsave("../results/Task05_AFTIR.Kd.pdf",width = width,height=height)
+g = g %+% plots.T.Kd
+g = g + labs(y="Trough Free Tissue target to\ninitial target Ratio (TFTIR)")
 plot(g)

-

-
ggsave("../results/Task05_TFTIR.pdf",width = width,height=height)

+

+
ggsave("../results/Task05_TFTIR.Kd.pdf",width = width,height=height)

Create data frame for multiple drugs and parameters.

diff --git a/ModelF/Task06_ThieleModulus_Explore.Rmd b/ModelF/Task06_ThieleModulus_Explore.Rmd new file mode 100644 index 0000000..ff6095c --- /dev/null +++ b/ModelF/Task06_ThieleModulus_Explore.Rmd @@ -0,0 +1,44 @@ +--- +title: "Calculate Thiele Modulus" +author: "Andy Stein" +date: "`r format(Sys.time(), '%d %B, %Y')`" +output: + html_document: + toc: true + toc_float: true + code_folding: hide +--- + +```{r, warning=FALSE} +suppressMessages(source("ams_initialize_script.R")) +``` + +## Load the data +```{r} + +#function for reading the parameter file +read.param.file.full = function(filename) { + d = read_excel(filename, 1) + p = as.numeric(d$Value) + names(p) = d$Parameter + return(p) +} + +#read in a parameter + drugs = c("Herceptin","Atezolizumab","Pembrolizumab") + P = list() + for (i in 1:length(drugs)) { + drugi = drugs[i] + p = read.param.file.full(paste0("../data/ModelF_",drugi,"_Params.xlsx")) + P[[i]] = p + } + P = bind_rows(P) + P$drug = drugs + P = select(P,drug,everything()) + +#compute the Thiele modulus + Thiele = + + +``` + diff --git a/ModelF/Task07_CompareEigenvalues.html b/ModelF/Task07_CompareEigenvalues.html new file mode 100644 index 0000000..9e2b6ee --- /dev/null +++ b/ModelF/Task07_CompareEigenvalues.html @@ -0,0 +1,198 @@ + + + + + + + + + + + + + +CompareEigenvalues + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ + + + + + + + + + + + + + +
suppressMessages(source("ams_initialize_script.R"))
+
## Warning: package 'ggplot2' was built under R version 3.2.5
+
## Warning: package 'rmarkdown' was built under R version 3.2.5
+
## Warning: package 'knitr' was built under R version 3.2.5
+
## Warning: package 'stringr' was built under R version 3.2.5
+
## Warning: package 'dplyr' was built under R version 3.2.5
+
## Warning: package 'ggrepel' was built under R version 3.2.5
+
source("CheckFunctions.R")
+

For Trastuzumab

+
model = ivsc_4cmtct_shedct()
+
+param.as.double = read.param.file("../data/ModelF_Trastuzumab_Params.xlsx")
+
## Warning in read_xlsx_(path, sheet, col_names = col_names, col_types =
+## col_types, : unknown timezone 'zone/tz/2018c.1.0/zoneinfo/America/New_York'
+
#make sure keD3 is 0
+print(param.as.double["keD3"])
+
## keD3 
+##    0
+
compare_eigenvalues(param.as.double=param.as.double)
+
## [1] "Eigenvalues of the matrix A are"
+## [1] -0.85350748 -0.48538545 -0.05817667
+## [1] "alpha, beta, gamma from the formula"
+## [1] 0.05817667 0.85350748 0.48538545
+

For Atezo

+
param.as.double = read.param.file("../data/ModelF_Atezolizumab_Params.xlsx")
+print(param.as.double["keD3"])
+
## keD3 
+##    0
+
compare_eigenvalues(param.as.double = param.as.double)
+
## [1] "Eigenvalues of the matrix A are"
+## [1] -0.46475731 -0.34295555 -0.02596482
+## [1] "alpha, beta, gamma from the formula"
+## [1] 0.02596482 0.46475731 0.34295555
+

For Pembrolizumab

+
param.as.double = read.param.file("../data/ModelF_Pembrolizumab_Params.xlsx")
+
## Warning in read.param.file("../data/ModelF_Pembrolizumab_Params.xlsx"): NAs
+## introduced by coercion
+
print(param.as.double["keD3"])
+
## keD3 
+##    0
+
compare_eigenvalues(param.as.double = param.as.double)
+
## [1] "Eigenvalues of the matrix A are"
+## [1] -0.60418701 -0.44958686 -0.02677083
+## [1] "alpha, beta, gamma from the formula"
+## [1] 0.02677083 0.60418701 0.44958686
+ + + + +
+ + + + + + + + diff --git a/ModelF/Task09b_Compare_AFIRT_AFIRT_.Rmd b/ModelF/Task09b_Compare_AFIRT_AFIRT_.Rmd index c590bab..f4da2f7 100644 --- a/ModelF/Task09b_Compare_AFIRT_AFIRT_.Rmd +++ b/ModelF/Task09b_Compare_AFIRT_AFIRT_.Rmd @@ -136,32 +136,32 @@ g = g + scale_color_manual (values = c(simulation = "black", theory.Kssd = "red", theory.Kssd_ = "green4"), labels = c(simulation = "Simulation", - theory.Kssd = "Theory Kssd", - theory.Kssd_ = "Theory Kssd*")) + theory.Kssd = "Theory Kssd simple", + theory.Kssd_ = "Theory Kssd")) g = g + scale_shape_manual (values = c(simulation = 46, theory.Kssd = 15, theory.Kssd_ = 17), labels = c(simulation = "Simulation", - theory.Kssd = "Theory Kssd", - theory.Kssd_ = "Theory Kssd*")) + theory.Kssd = "Theory Kssd simple", + theory.Kssd_ = "Theory Kssd")) g = g + scale_size_manual (values = c(simulation = 1, theory.Kssd = 1, theory.Kssd_ = 1), labels = c(simulation = "Simulation", - theory.Kssd = "Theory Kssd", - theory.Kssd_ = "Theory Kssd*")) + theory.Kssd = "Theory Kssd simple", + theory.Kssd_ = "Theory Kssd")) g = g + scale_linetype_manual(values = c(simulation = "solid", theory.Kssd = "dotdash", theory.Kssd_ = "dashed"), labels = c(simulation = "Simulation", - theory.Kssd = "Theory Kssd", - theory.Kssd_ = "Theory Kssd*")) + theory.Kssd = "Theory Kssd simple", + theory.Kssd_ = "Theory Kssd")) g = g + scale_alpha_manual (values = c(simulation = .5, theory.Kssd = 1, theory.Kssd_ = 1), labels = c(simulation = "Simulation", - theory.Kssd = "Theory Kssd", - theory.Kssd_ = "Theory Kssd*")) + theory.Kssd = "Theory Kssd simple", + theory.Kssd_ = "Theory Kssd")) print(g) width = 10 diff --git a/ModelF/Task09b_Compare_AFIRT_AFIRT_.html b/ModelF/Task09b_Compare_AFIRT_AFIRT_.html new file mode 100644 index 0000000..150c7b5 --- /dev/null +++ b/ModelF/Task09b_Compare_AFIRT_AFIRT_.html @@ -0,0 +1,246 @@ + + + + + + + + + + + + + + +Task09b_Compare_AFIRT_AFIRT* + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ + + + + + + + + + + + + + +
# Need this to get rxode working. I don't know why I keep having to run this.
+Sys.setenv(PATH = paste(Sys.getenv("PATH"), "C:/RBuildTools/3.4/bin/",
+                        "C:/RBuildTools/3.4/mingw_64/bin", sep = ";"))
+Sys.setenv(BINPREF = "C:/RBuildTools/3.4/mingw_64/bin/")
+
+# To be called at the top of every Rmd file. Initialization code and some useful constants.
+suppressMessages(source("ams_initialize_script.R"))
+
+

Create data frame for multiple drugs and parameters.

+
+
+

Plot the output AFIRT —-

+
# I only wanted a subset of the parameters for gather.
+joined.AT = joined[c("fold.change.param", "AFIRT.Kssd", "AFIRT.Kssd_", "AFIRT.sim", "TFIRT.Kssd", "TFIRT.Kssd_", "TFIRT.sim", "param", "drug", "isSol", "M30", "S30")]
+
+# Use gather to make the long data frame for ggplot.
+plots.AT = joined.AT %>% 
+  gather(key, AFIRT.value, -param, -fold.change.param, -drug, -isSol, -M30, -S30) %>%
+  filter(!is.na(AFIRT.value)) %>%
+  mutate(Target = ifelse(isSol,"Soluble","Membrane-Bd"),
+         Target = factor(Target,levels=c("Soluble","Membrane-Bd")),
+         AT     = NA,
+         AT     = ifelse(str_detect(key,"AFIRT"),"AFIRT",AT),
+         AT     = ifelse(str_detect(key,"TFIRT"),"TFIRT",AT),
+         key    = str_replace(key,"\\wFIRT.",""),
+         key    = str_replace(key,"sim","simulation"),
+         key    = str_replace(key,"K","theory.K"),
+         # key    = str_replace(key,"_","*"),
+         key    = factor(key,levels=c("simulation","theory.Kssd","theory.Kssd_")),
+         T30    = ifelse(isSol,S30,M30),
+         T30    = paste("T0 =",signif(T30,1),"nM")) %>%
+  mutate(drug   = str_replace(drug,"Herceptin","Trastuzumab"),
+         drug   = factor(drug,levels=c("Bevacizumab","Pembrolizumab","Atezolizumab","Trastuzumab")))
+
+  plots.A = filter(plots.AT,AT!="TFIRT")
+  plots.T = filter(plots.AT,AT!="AFIRT")
+
+g = ggplot(plots.A, aes(fold.change.param, AFIRT.value, color=key, shape=key, linetype=key))
+g = g + scale.x.log10()
+g = g + scale.y.log10()
+g = g + geom_line(mapping=aes(size=key))
+#g = g + geom_point()
+g = g + facet_grid(.~drug + Target + T30)
+g = g + labs(x        = "Dose (mg/kg) every 3 weeks",
+             y        = "Average Free Tissue target to\nInitial target Ratio (AFTIR)",
+             color    = "",
+             shape    = "",
+             size     = "",
+             linetype = "")
+g = g + scale_color_manual   (values = c(simulation   = "black",
+                                         theory.Kssd  = "red",
+                                         theory.Kssd_ = "green4"),
+                              labels = c(simulation   = "Simulation", 
+                                         theory.Kssd  = "Theory Kssd simple", 
+                                         theory.Kssd_ = "Theory Kssd"))
+g = g + scale_shape_manual   (values = c(simulation   = 46,
+                                         theory.Kssd  = 15,
+                                         theory.Kssd_ = 17),
+                              labels = c(simulation   = "Simulation", 
+                                         theory.Kssd  = "Theory Kssd simple", 
+                                         theory.Kssd_ = "Theory Kssd"))
+g = g + scale_size_manual    (values = c(simulation   = 1,
+                                         theory.Kssd  = 1,
+                                         theory.Kssd_ = 1),
+                              labels = c(simulation   = "Simulation", 
+                                         theory.Kssd  = "Theory Kssd simple", 
+                                         theory.Kssd_ = "Theory Kssd"))
+g = g + scale_linetype_manual(values = c(simulation   = "solid",
+                                         theory.Kssd  = "dotdash",
+                                         theory.Kssd_ = "dashed"),
+                              labels = c(simulation   = "Simulation", 
+                                         theory.Kssd  = "Theory Kssd simple", 
+                                         theory.Kssd_ = "Theory Kssd"))
+g = g + scale_alpha_manual   (values = c(simulation   = .5,
+                                         theory.Kssd  = 1,
+                                         theory.Kssd_ = 1),
+                              labels = c(simulation   = "Simulation", 
+                                         theory.Kssd  = "Theory Kssd simple", 
+                                         theory.Kssd_ = "Theory Kssd"))
+
+print(g)
+

+
width = 10
+height = 2.7
+ggsave("../results/Task09b_Compare_AFTIR_AFTIRstar.pdf",width=width,height=height)
+
+g = g %+% plots.T
+g = g + labs(y="Trough Free Tissue target to\ninitial target Ratio (TFTIR)")
+print(g)
+

+
ggsave("../results/Task09b_Compare_TFTIR_TFTIRstar.pdf",width=width,height=height)
+
+ + + + +
+ + + + + + + + diff --git a/ModelF/Task09b_PAGE2018_Poster_AFIRTstar.R b/ModelF/Task09b_PAGE2018_Poster_AFIRTstar.R new file mode 100644 index 0000000..74a7296 --- /dev/null +++ b/ModelF/Task09b_PAGE2018_Poster_AFIRTstar.R @@ -0,0 +1,153 @@ +#HEADER +# Need this to get rxode working. I don't know why I keep having to run this. +Sys.setenv(PATH = paste(Sys.getenv("PATH"), "C:/RBuildTools/3.4/bin/", + "C:/RBuildTools/3.4/mingw_64/bin", sep = ";")) +Sys.setenv(BINPREF = "C:/RBuildTools/3.4/mingw_64/bin/") + +# To be called at the top of every Rmd file. Initialization code and some useful constants. +suppressMessages(source("ams_initialize_script.R")) + +# Note: Insoluble and soluble drugs are treated differently. +# For insoluble, Tacc = M3totss/M30. For soluble, Tacc = S3totss/S30. +# We are also interested in different parameters in the SA for either case. +# See writeup for details. + +# Initialize ---- + +# Load model. +model = ivsc_4cmtct_shedct(target = TRUE) +# Drugs to explore. +drugs = c("Atezolizumab", "Trastuzumab", "Pembrolizumab", "Bevacizumab") +# List of parameters of interest. +parameters = c("dose") +# names(parameters) = parameters +# Create vector of soluble drugs. +soluble = "Bevacizumab" + +# Create paths to data files for each drug. +paths = NULL +for (i in 1:length(drugs)) { + paths[i] = paste0("../data/ModelF_", drugs[i],"_Params.xlsx") +} + +# Dose time, frequency, compartment, nominal dose +tmax = 26*7 #days +tau = 21 #days +compartment = 2 +dose.nmol = scale.mpk2nmol +joined = NULL +isSoluble = FALSE + +# Iterate over all the drugs. ----------------- +for (i in 1:length(drugs)){ + # Load parameters. + param.as.double = read.param.file(paths[i]) + df_param = as.data.frame(t(param.as.double)) + + # Check if drug has soluble target. + if (drugs[i] %in% soluble) { + # Change any non-shed parameters including M to S. + parameters = sapply(parameters, + function(x) { + if (!grepl("shed|dose",x)) {return(gsub("M","S", x))} + return(x) + }) + isSoluble = TRUE + } + + # Set range for parameters of interest in SA. + # Check which parameters are nonzero, not including dose which isn't in df_param. + n.iter = 13 + params.to.iterate = data.frame(dose = lseq(scale.mpk2nmol*.1,scale.mpk2nmol*100,length.out = n.iter)) + + dfs = list() + # Iterate all of the parameters for a single drug. + for (j in 1:ncol(params.to.iterate)){ + dfs[[j]] = compare.thy.sim(model = model, + param.as.double = param.as.double, + dose.nmol = dose.nmol, + tmax = tmax, + tau = tau, + compartment = compartment, + param.to.change = names(params.to.iterate)[j], + param.to.change.range = params.to.iterate[[j]], + soluble = isSoluble) + dfs[[j]] = dfs[[j]] %>% mutate(drug = drugs[i], isSol = isSoluble) + } + joined = bind_rows(joined,dfs) + + # Reset isSoluble to false since the default in compare.thy.sim is false. + isSoluble = FALSE +} + +# Plot the output AFIRT ---- + +# I only wanted a subset of the parameters for gather. +joined.AT = joined[c("fold.change.param","AFIRT.Kssd_", "AFIRT.Kss_", "AFIRT.Kd_", "AFIRT.sim", "TFIRT.Kssd_", "TFIRT.Kss_", "TFIRT.Kd_", "TFIRT.sim","param","drug", "isSol","M30","S30")] +names(joined.AT) = names(joined.AT) %>% + str_replace("_","") + +# Use gather to make the long data frame for ggplot. +plots.AT = joined.AT %>% + gather(key, AFIRT.value, -param, -fold.change.param,-drug, -isSol, -M30, -S30) %>% + filter(!is.na(AFIRT.value)) %>% + mutate(Target = ifelse(isSol,"Soluble","Membrane-Bd"), + Target = factor(Target,levels=c("Soluble","Membrane-Bd")), + AT = NA, + AT = ifelse(str_detect(key,"AFIRT"),"AFIRT",AT), + AT = ifelse(str_detect(key,"TFIRT"),"TFIRT",AT), + key = str_replace(key,"\\wFIRT.",""), + key = str_replace(key,"sim","simulation"), + key = str_replace(key,"K","theory.K"), + key = factor(key,levels=c("simulation","theory.Kssd","theory.Kss","theory.Kd")), + T30 = ifelse(isSol,S30,M30), + T30 = paste("T0 =",signif(T30,1),"nM")) %>% + mutate(drug = str_replace(drug,"Herceptin","Trastuzumab"), + drug = factor(drug,levels=c("Bevacizumab","Pembrolizumab","Atezolizumab","Trastuzumab"))) + + plots.A = filter(plots.AT,AT!="TFIRT") + plots.T = filter(plots.AT,AT!="AFIRT") + +g = ggplot(plots.A, aes(fold.change.param, AFIRT.value, color=key, shape=key, linetype=key)) +g = g + scale.x.log10(limits=c(.05,30)) +g = g + scale.y.log10(limits = c(.001,1)) +g = g + geom_line(mapping=aes(size=key)) +#g = g + geom_point() +g = g + facet_grid(.~drug + Target + T30) +g = g + labs(x = "Dose (mg/kg) every 3 weeks", + y = "Average Free Tissue target to\nInitial target Ratio* (AFTIR*)", + color = "", + shape = "", + size = "", + linetype="") +g = g + scale_color_manual(values = c(simulation = "black", + theory.Kd = "red", + theory.Kss = "green4", + theory.Kssd= "blue")) +g = g + scale_shape_manual(values = c(simulation = 46, + theory.Kd = 15, + theory.Kss = 17, + theory.Kssd= 16)) +g = g + scale_size_manual(values = c(simulation = 1, + theory.Kd = 1, + theory.Kss = 1, + theory.Kssd= 1 )) +g = g + scale_linetype_manual(values = c(simulation = "solid", + theory.Kd = "dotdash", + theory.Kss = "dashed", + theory.Kssd = "longdash")) +g = g + scale_alpha_manual(values = c(simulation = .5, + theory.Kd = 1, + theory.Kss = 1, + theory.Kssd = 1)) + +print(g) +width = 7 +height = 2.7 +ggsave("../results/Task09b_PAGE_Figure_AFTIRstar.pdf",width = width,height=height) + +g = g %+% plots.T +g = g + labs(y="Trough Free Tissue target to\ninitial target Ratio* (TFTIR*)") +print(g) +ggsave("../results/Task09b_PAGE_Figure_TFTIRstar.pdf",width = width,height=height) + diff --git a/ModelF/Task09c_PAGE2018_Poster_AFIRTstar.R b/ModelF/Task09c_PAGE2018_Poster_AFIRTstar.R index 70b6a39..57a54d0 100644 --- a/ModelF/Task09c_PAGE2018_Poster_AFIRTstar.R +++ b/ModelF/Task09c_PAGE2018_Poster_AFIRTstar.R @@ -115,7 +115,7 @@ g = g + geom_line(mapping=aes(size=key)) #g = g + geom_point() g = g + facet_grid(.~drug + Target + T30) g = g + labs(x = "Dose (mg/kg) every 3 weeks", - y = "Average Free Tissue target to\nInitial target Ratio* (AFTIR*)", + y = "Average Free Tissue target to\nInitial target Ratio (AFTIR)", color = "", shape = "", size = "", diff --git a/ModelF/Task11_Jumps_In_Sensitivity_Analysis_v2.html b/ModelF/Task11_Jumps_In_Sensitivity_Analysis_v2.html new file mode 100644 index 0000000..ae33d7b --- /dev/null +++ b/ModelF/Task11_Jumps_In_Sensitivity_Analysis_v2.html @@ -0,0 +1,361 @@ + + + + + + + + + + + + + + +Task11_Jumps_In_Sensitivity_Analysis + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ + + + + + + + + + + + + + + + +
+
+
+
+
+ +
+ + + + + + + +

In the first plot of “Task05_Sensistivity_Analysis_pt5.Rmd”, it is strange that AFIRT.sim jumps as ksynM3 changes from the 4th value to the 5th value (Take pemb as an example), because as indicated in the AFTIR derivation part of the writeup, Keq does not depend on ksynM3, i.e. AFTIR does not depend on ksynM3. This is illustrated in the plots of AFTIR.Kssd, AFTIR.Kss, and AFTIR.Kd, where the dependences on ksynM3 are barely noticeable.

+

This indicates that our approximation of Keq does not apply when ksynM3 is large (or when target concentration is large). By definition \[ + K_{ssd, mem} = \frac{D_3\cdot M_3}{DM_3} +\] so it is reasonable that as ksynM3 gets big, \(K_{ssd, mem}\) gets big. Then, we will not be able assume (2.2.20 of the writeup) \[ + \frac{T_3}{T_{3tot}} = \frac{K_{eq}}{K_{eq} + D_3} \sim + \frac{K_{eq}}{D_3} +\]

+

Therefore, in order for the above approximation to be true, we will have to increase the dose for large ksynM3, so that \(Keq << D_3\).

+

I will run a sensitivity analysis on ksynM3 with the same range of value as in Task05, but I will significantly increase the dose. Then, we should be able to
+see better agreement between theory and simulation

+
suppressMessages(source("ams_initialize_script.R"))
+suppressMessages(source("ams_initialize_script.R"))
+model = ivsc_4cmtct_shedct(target=TRUE)
+
+# test on Pemb, and run senstivity analysis on ksynM3
+drug = "Pembrolizumab"
+drug.path = paste0("../data/ModelF_", drug, "_Params.xlsx")
+parameter = "ksynM3"
+
+# Set global variables
+tmax = 26*7
+tau = 21
+compartment = 2
+
+# Load model parameters for pemb
+param.as.double = read.param.file(drug.path)
+df_param = as.data.frame(t(param.as.double))
+
+# grep ksynM3
+ksynM3.init = df_param["ksynM3"]
+
+# set the range for ksynM3
+ksynM3.range = data.frame(lapply(ksynM3.init, function(x) lseq(x*0.01, x*100, 7)))
+

With low dose

+
dose.nmol = scale.mpk2nmol
+# run the sensitivity analysis on the above range
+AFIRT = compare.thy.sim(model=model,
+                    param.as.double=param.as.double,
+                    dose.nmol = dose.nmol,
+                    tmax = tmax,
+                    tau = tau,
+                    compartment = compartment,
+                    param.to.change = 'ksynM3',
+                    param.to.change.range = ksynM3.range[[1]],
+                    soluble = FALSE)
+
+# get the relevant data from the above frame
+relevant = AFIRT[c("fold.change", "AFIRT.Kssd", "AFIRT.Kss","AFIRT.Kd", "AFIRT.sim", "param")]
+
+relevant = relevant %>% 
+  gather(key, AFIRT.value, -param, -fold.change) %>%
+  filter(!is.na(AFIRT.value))
+
+g = ggplot(relevant, aes(fold.change, AFIRT.value, color=key, linetype=key)) +
+  scale.x.log10() + 
+  scale.y.log10() + 
+  geom_line() + 
+  geom_point()
+
+print(g)
+

+

With high dose (100 times higher)

+
dose.nmol = 1000*scale.mpk2nmol
+# run the sensitivity analysis on the above range
+AFIRT = compare.thy.sim(model=model,
+                    param.as.double=param.as.double,
+                    dose.nmol = dose.nmol,
+                    tmax = tmax,
+                    tau = tau,
+                    compartment = compartment,
+                    param.to.change = 'ksynM3',
+                    param.to.change.range = ksynM3.range[[1]],
+                    soluble = FALSE)
+
+# get the relevant data from the above frame
+relevant = AFIRT[c("fold.change", "AFIRT.Kssd", "AFIRT.Kss","AFIRT.Kd", "AFIRT.sim", "param")]
+
+relevant = relevant %>% 
+  gather(key, AFIRT.value, -param, -fold.change) %>%
+  filter(!is.na(AFIRT.value))
+
+g = ggplot(relevant, aes(fold.change, AFIRT.value, color=key, linetype=key)) +
+  scale.x.log10() + 
+  scale.y.log10() + 
+  geom_line() + 
+  geom_point()
+
+print(g)
+

+ + + +
+
+ +
+ + + + + + + + diff --git a/ModelF/Task12_Compare_AFIRTKssd_AFIRTStar_Simulation.html b/ModelF/Task12_Compare_AFIRTKssd_AFIRTStar_Simulation.html new file mode 100644 index 0000000..505f98f --- /dev/null +++ b/ModelF/Task12_Compare_AFIRTKssd_AFIRTStar_Simulation.html @@ -0,0 +1,201 @@ + + + + + + + + + + + + + + + +Compare AFIRT AFIRT* and Simulation + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ + + + + + + + + + + + + + +

The goal of this task is to compare AFIRT, AFIRT* and simulation over a range of dose

+
suppressMessages(source("ams_initialize_script.R"))
+
+dose.nmol.range = lseq(0.1, 100000, 7)
+model = ivsc_4cmtct_shedct()
+tmax = 26*7 #days
+tau = 21 #days
+compartment = 2
+

Load parameters for Atezo

+
param.as.double = read.param.file("../data/ModelF_Atezolizumab_Params.xlsx")
+p.atezo = as.data.frame(t(param.as.double))
+
df_sim = data.frame() # put all simulations for different dose into one data frame
+  for (dose.nmol in dose.nmol.range){
+    row = lumped.parameters.simulation(model,
+      param.as.double, dose.nmol, tmax, tau, compartment)
+    df_sim = rbind(df_sim, row)
+  }
+
+df_thy = data.frame() # put all theoretical calculations of lumped parameters at different dose together
+for (dose.nmol in dose.nmol.range){
+  row = lumped.parameters.theory(param.as.double=param.as.double, 
+                                dose.nmol=dose.nmol, 
+                                tau=tau)
+  df_thy = rbind(df_thy, row)
+}
+
+df_thy$Dose = dose.nmol.range
+

Make a data frame with 4 columns: Dose, AFIRT.Kssd, AFIRT_, AFIRT.sim

+
df_compare = df_thy %>%
+  select(Dose, AFIRT.Kssd, AFIRT.Kssd_)
+
+df_compare$AFIRT.sim = df_sim$AFIRT.sim
+

Plot the data

+
df_plot = df_compare %>%
+  gather(AFIRT, value, -Dose)
+
+g = ggplot(df_plot, aes(Dose, value, group=AFIRT, color=AFIRT)) + 
+  geom_line() + 
+  geom_point()+
+  scale.x.log10() +
+  scale.y.log10()
+
+print(g)
+

+ + + + +
+ + + + + + + + diff --git a/ModelF/Task13_Example_Individual_Run_Plot.Rmd b/ModelF/Task13_Example_Individual_Run_Plot.Rmd new file mode 100644 index 0000000..d6a277f --- /dev/null +++ b/ModelF/Task13_Example_Individual_Run_Plot.Rmd @@ -0,0 +1,105 @@ +--- +title: "Task13_Example_Individual_Run_Plot.R" +author: "Andy Stein" +date: "`r format(Sys.time(), '%d %B, %Y')`" +output: + html_document: + toc: true + toc_float: true + code_folding: hide +--- + +# Preliminary stuff + +```{r, warning=FALSE} +# Sameed need this to get rxode working. I don't know why I keep having to run this. +Sys.setenv(PATH = paste(Sys.getenv("PATH"), "C:/RBuildTools/3.4/bin/", + "C:/RBuildTools/3.4/mingw_64/bin", sep = ";")) +Sys.setenv(BINPREF = "C:/RBuildTools/3.4/mingw_64/bin/") + +# To be called at the top of every Rmd file. Initialization code and some useful constants. +source("ams_initialize_script.R") +``` + +# Individual simulation examlpe with rxode + +```{r} +# Load model and parameters. +model = ivsc_4cmtct_shedct() +drugs = c("Pembrolizumab") +filename = paste0("../data/ModelF_", drugs,"_Params.xlsx") +param.as.double = read.param.file(filename) +df_param = as.data.frame(t(param.as.double)) + +#Set dosing parameters +dose.nmol = scale.mpk2nmol #1 mg/kg (mpk) +tmax = 26*7 #days (maximum time point) +tau = 21 #days between doses +compartment = 2 #compartment drug is dosed into + +#set sampling and dosing event (ev) table +ev = eventTable(amount.units="nmol", time.units="days") +sample.points = c(seq(-7, tmax, 0.1), 10^(-3:0)) # sample time, increment by 0.1 +sample.points = sort(sample.points) +sample.points = unique(sample.points) +ev$add.sampling(sample.points) +ev$add.dosing(dose = dose.nmol, + nbr.doses = floor(tmax/tau)+1, + dosing.interval = tau, + dosing.to = compartment) + +#set the initial condition +init = model$init(param.as.double) + +#run the model +out = model$rxode$solve(param.as.double, ev, init) +out = model$rxout(out) +out = out %>% + mutate(Sfree.pct = S3/init["S3"], + Mfree.pct = M3/init["M3"], + dose.nmol = dose.nmol) + + +#modify output dataframe for plotting +out.plot = out %>% + select(time,D1,D3,M3,DM3,Mfree.pct) %>% + gather(cmt,value,-c(time)) + +#plot the results + g = ggplot(out.plot,aes(x=time,y=value,group=cmt,color=cmt)) + g = g + geom_line() + g = g + scale.y.log10() + g = xscale("w27",increment = 3) + print(g) +``` + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/ModelF/Task13_Example_Individual_Run_Plot.html b/ModelF/Task13_Example_Individual_Run_Plot.html new file mode 100644 index 0000000..ca2ea4d --- /dev/null +++ b/ModelF/Task13_Example_Individual_Run_Plot.html @@ -0,0 +1,352 @@ + + + + + + + + + + + + + + +Task13 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ + + + + + + + + + + + + + + + +
+
+
+
+
+ +
+ + + + + + + +
+

Preliminary stuff

+
# Sameed need this to get rxode working. I don't know why I keep having to run this.
+Sys.setenv(PATH = paste(Sys.getenv("PATH"), "C:/RBuildTools/3.4/bin/",
+            "C:/RBuildTools/3.4/mingw_64/bin", sep = ";"))
+Sys.setenv(BINPREF = "C:/RBuildTools/3.4/mingw_64/bin/")
+
+# To be called at the top of every Rmd file. Initialization code and some useful constants.
+source("ams_initialize_script.R")
+
## 
+## Attaching package: 'tidyr'
+
## The following object is masked from 'package:reshape2':
+## 
+##     smiths
+
## 
+## Attaching package: 'dplyr'
+
## The following object is masked from 'package:gridExtra':
+## 
+##     combine
+
## The following objects are masked from 'package:stats':
+## 
+##     filter, lag
+
## The following objects are masked from 'package:base':
+## 
+##     intersect, setdiff, setequal, union
+
+
+

Individual simulation examlpe with rxode

+
# Load model and parameters.
+model           = ivsc_4cmtct_shedct()
+drugs           = c("Pembrolizumab")
+filename        = paste0("../data/ModelF_", drugs,"_Params.xlsx")
+param.as.double = read.param.file(filename)
+df_param        = as.data.frame(t(param.as.double))
+
+#Set dosing parameters
+dose.nmol       = scale.mpk2nmol #1 mg/kg (mpk)
+tmax            = 26*7 #days (maximum time point)
+tau             = 21   #days between doses
+compartment     = 2    #compartment drug is dosed into
+
+#set sampling and dosing event (ev) table
+ev = eventTable(amount.units="nmol", time.units="days")
+sample.points = c(seq(-7, tmax, 0.1), 10^(-3:0)) # sample time, increment by 0.1
+sample.points = sort(sample.points)
+sample.points = unique(sample.points)
+ev$add.sampling(sample.points)
+ev$add.dosing(dose            = dose.nmol, 
+              nbr.doses       = floor(tmax/tau)+1, 
+              dosing.interval = tau,
+              dosing.to       = compartment)
+
+#set the initial condition
+init = model$init(param.as.double)
+    
+#run the model
+out  = model$rxode$solve(param.as.double, ev, init)
+out  = model$rxout(out)
+out  = out %>%
+        mutate(Sfree.pct = S3/init["S3"],
+               Mfree.pct = M3/init["M3"],
+               dose.nmol = dose.nmol)
+
+
+#modify output dataframe for plotting
+out.plot = out %>%
+    select(time,D1,D3,M3,DM3,Mfree.pct) %>%
+    gather(cmt,value,-c(time))
+  
+#plot the results
+  g = ggplot(out.plot,aes(x=time,y=value,group=cmt,color=cmt))
+  g = g + geom_line()
+  g = g + scale.y.log10()
+  g = xscale("w27",increment = 3)
+  print(g)
+
## Warning: Transformation introduced infinite values in continuous y-axis
+
## Warning: Removed 350 rows containing missing values (geom_path).
+

+
+ + + +
+
+ +
+ + + + + + + + diff --git a/ModelF/Users/steinanf/GitHub/TumorModeling_IMA/ModelF/ivsc_4cmtct_shedct.d/ivsc_4cmtct_shedct_.so b/ModelF/Users/steinanf/GitHub/TumorModeling_IMA/ModelF/ivsc_4cmtct_shedct.d/ivsc_4cmtct_shedct_.so new file mode 100755 index 0000000..b7711af Binary files /dev/null and b/ModelF/Users/steinanf/GitHub/TumorModeling_IMA/ModelF/ivsc_4cmtct_shedct.d/ivsc_4cmtct_shedct_.so differ diff --git a/ModelF/ivsc_4cmtct_shedct2.d/LHS_VARS.txt b/ModelF/ivsc_4cmtct_shedct2.d/LHS_VARS.txt new file mode 100644 index 0000000..c754739 --- /dev/null +++ b/ModelF/ivsc_4cmtct_shedct2.d/LHS_VARS.txt @@ -0,0 +1 @@ +D1 \ No newline at end of file diff --git a/ModelF/ivsc_4cmtct_shedct2.d/ODE_PARS.txt b/ModelF/ivsc_4cmtct_shedct2.d/ODE_PARS.txt new file mode 100644 index 0000000..6b21dca --- /dev/null +++ b/ModelF/ivsc_4cmtct_shedct2.d/ODE_PARS.txt @@ -0,0 +1 @@ +VD1 ka F k12D k21D VD2 k13D k31D VD3 kon1 koff1 keD1 kon3 koff3 keD3 ksynS1 kshedM1 k13S k31S VS3 VS1 keS1 ksynS3 kshedM3 keS3 ksynM1 k13M k31M VM3 VM1 keM1 ksynM3 keM3 kshedDM1 k13DS k31DS VDS3 VDS1 keDS1 kshedDM3 keDS3 k13DM k31DM VDM3 VDM1 keDM1 keDM3 \ No newline at end of file diff --git a/ModelF/ivsc_4cmtct_shedct2.d/STATE_VARS.txt b/ModelF/ivsc_4cmtct_shedct2.d/STATE_VARS.txt new file mode 100644 index 0000000..45fe637 --- /dev/null +++ b/ModelF/ivsc_4cmtct_shedct2.d/STATE_VARS.txt @@ -0,0 +1 @@ +AmtD0 AmtD1 D2 D3 S1 S3 M1 M3 DS1 DS3 DM1 DM3 \ No newline at end of file diff --git a/ModelF/ivsc_4cmtct_shedct2.d/error.txt b/ModelF/ivsc_4cmtct_shedct2.d/error.txt new file mode 100644 index 0000000..e69de29 diff --git a/ModelF/ivsc_4cmtct_shedct2.d/ivsc_4cmtct_shedct2.so b/ModelF/ivsc_4cmtct_shedct2.d/ivsc_4cmtct_shedct2.so new file mode 100755 index 0000000..de348c8 Binary files /dev/null and b/ModelF/ivsc_4cmtct_shedct2.d/ivsc_4cmtct_shedct2.so differ diff --git a/ModelF/ivsc_4cmtct_shedct2.d/model.txt b/ModelF/ivsc_4cmtct_shedct2.d/model.txt new file mode 100644 index 0000000..4c4da82 --- /dev/null +++ b/ModelF/ivsc_4cmtct_shedct2.d/model.txt @@ -0,0 +1,15 @@ + + D1 = AmtD1/VD1; + d/dt(AmtD0) = -ka *AmtD0; + d/dt(AmtD1) =(F*ka *AmtD0/VD1 - k12D*D1 + k21D*VD2/VD1*D2 - k13D *D1 + k31D *VD3 /VD1* D3 - kon1*D1*(S1+M1) + koff1*(DS1+DM1) - keD1 *D1)*VD1; + d/dt(D2) = - k21D*D2 + k12D*VD1/VD2*D1 ; + d/dt(D3) = - k31D *D3 + k13D *VD1 /VD3* D1 - kon3*D3*(S3+M3) + koff3*(DS3+DM3) - keD3 *D3; + d/dt(S1) = ksynS1 + kshedM1*M1 - k13S *S1 + k31S *VS3 /VS1* S3 - kon1*D1*S1 + koff1*DS1 - keS1 *S1 ; + d/dt(S3) = ksynS3 + kshedM3*M3 - k31S *S3 + k13S *VS1 /VS3* S1 - kon3*D3*S3 + koff3*DS3 - keS3 *S3; + d/dt(M1) = ksynM1 - kshedM1*M1 - k13M *M1 + k31M *VM3 /VM1* M3 - kon1*D1*M1 + koff1*DM1 - keM1 *M1 ; + d/dt(M3) = ksynM3 - kshedM3*M3 - k31M *M3 + k13M *VM1 /VM3* M1 - kon3*D3*M3 + koff3*DM3 - keM3 *M3; + d/dt(DS1) = kshedDM1*DM1 - k13DS*DS1 + k31DS*VDS3/VDS1*DS3 + kon1*D1*S1 - koff1*DS1 - keDS1*DS1 ; + d/dt(DS3) = kshedDM3*DM3 - k31DS*DS3 + k13DS*VDS1/VDS3*DS1 + kon3*D3*S3 - koff3*DS3 - keDS3*DS3 ; + d/dt(DM1) = - kshedDM1*DM1 - k13DM*DM1 + k31DM*VDM3/VDM1*DM3 + kon1*D1*M1 - koff1*DM1 - keDM1*DM1; + d/dt(DM3) = - kshedDM3*DM3 - k31DM*DM3 + k13DM*VDM1/VDM3*DM1 + kon3*D3*M3 - koff3*DM3 - keDM3*DM3; + diff --git a/latex/figures/Graphical_Abstract.png b/latex/figures/Graphical_Abstract.png new file mode 100644 index 0000000..c93c8c3 Binary files /dev/null and b/latex/figures/Graphical_Abstract.png differ diff --git a/latex/figures/Graphical_Abstract.pptx b/latex/figures/Graphical_Abstract.pptx new file mode 100644 index 0000000..0103f0f Binary files /dev/null and b/latex/figures/Graphical_Abstract.pptx differ diff --git a/poster/2018-05-30_Stein_etal_PAGE_AFTIR_Poster.pptx b/poster/2018-05-30_Stein_etal_PAGE_AFTIR_Poster.pptx new file mode 100644 index 0000000..cb4571f Binary files /dev/null and b/poster/2018-05-30_Stein_etal_PAGE_AFTIR_Poster.pptx differ diff --git a/results/Task04_Taccum.pdf b/results/Task04_Taccum.pdf index 2cce8db..eab69ac 100644 Binary files a/results/Task04_Taccum.pdf and b/results/Task04_Taccum.pdf differ diff --git a/results/Task05_AFTIR.Kd.pdf b/results/Task05_AFTIR.Kd.pdf index c4d93e0..c1b50c9 100644 Binary files a/results/Task05_AFTIR.Kd.pdf and b/results/Task05_AFTIR.Kd.pdf differ diff --git a/results/Task05_AFTIR.Kss.pdf b/results/Task05_AFTIR.Kss.pdf index 1bb245b..dc11bc9 100644 Binary files a/results/Task05_AFTIR.Kss.pdf and b/results/Task05_AFTIR.Kss.pdf differ diff --git a/results/Task05_AFTIR.Kssd.pdf b/results/Task05_AFTIR.Kssd.pdf index 63906ea..5b55eb6 100644 Binary files a/results/Task05_AFTIR.Kssd.pdf and b/results/Task05_AFTIR.Kssd.pdf differ diff --git a/results/Task05_TFTIR.Kd.pdf b/results/Task05_TFTIR.Kd.pdf index 253c7a0..4148255 100644 Binary files a/results/Task05_TFTIR.Kd.pdf and b/results/Task05_TFTIR.Kd.pdf differ diff --git a/results/Task05_TFTIR.Kss.pdf b/results/Task05_TFTIR.Kss.pdf index c3b8200..fa2826b 100644 Binary files a/results/Task05_TFTIR.Kss.pdf and b/results/Task05_TFTIR.Kss.pdf differ diff --git a/results/Task05_TFTIR.Kssd.pdf b/results/Task05_TFTIR.Kssd.pdf index 99a02ac..fa44136 100644 Binary files a/results/Task05_TFTIR.Kssd.pdf and b/results/Task05_TFTIR.Kssd.pdf differ diff --git a/results/Task09b_Compare_AFTIR_AFTIRstar.pdf b/results/Task09b_Compare_AFTIR_AFTIRstar.pdf index 3259084..3c4c4f8 100644 Binary files a/results/Task09b_Compare_AFTIR_AFTIRstar.pdf and b/results/Task09b_Compare_AFTIR_AFTIRstar.pdf differ diff --git a/results/Task09b_Compare_TFTIR_TFTIRstar.pdf b/results/Task09b_Compare_TFTIR_TFTIRstar.pdf index db8d12d..a8f974e 100644 Binary files a/results/Task09b_Compare_TFTIR_TFTIRstar.pdf and b/results/Task09b_Compare_TFTIR_TFTIRstar.pdf differ diff --git a/results/Task09c_PAGE_Figure_AFTIRstar.pdf b/results/Task09c_PAGE_Figure_AFTIRstar.pdf index e4e5399..a6ccbb1 100644 Binary files a/results/Task09c_PAGE_Figure_AFTIRstar.pdf and b/results/Task09c_PAGE_Figure_AFTIRstar.pdf differ diff --git a/results/Task09c_PAGE_Figure_TFTIRstar.pdf b/results/Task09c_PAGE_Figure_TFTIRstar.pdf index ae19ec2..dffda6a 100644 Binary files a/results/Task09c_PAGE_Figure_TFTIRstar.pdf and b/results/Task09c_PAGE_Figure_TFTIRstar.pdf differ