Data import

First let’s import the datasets and exclude some of the problematic points that were obvious from the past analyses.

df_blank <- read_csv("../../11_28_18_blank_IDA/Processing/11_28_18_swv_gc_soak_processed.csv") %>% 
  select(reactor,electrode_from_swv,signal_from_swv,electrode_from_gc, signal_from_gc,echem_from_swv ) %>% 
  mutate(run = 1, rep = 0, treatment='blank', exp_id = 'blank')


df_3 <- read_csv("../../01_23_19_psoralen_nonequil_3/Processing/01_23_19_processed_swv_gc_all.csv") %>% 
  select(treatment,run,rep,reactor,electrode_from_swv,signal_from_swv,electrode_from_gc, signal_from_gc, echem_from_swv ) %>% 
  mutate(exp_id = '3') %>% 
  filter(run!=1 | treatment!='psoralen' | rep<14 ) %>% 
  filter(run!=4 | treatment!='control' | rep>1 )
  

df_2_control <- read_csv("../../01_17_19_psoralen_nonequil_2/Processing/01_17_19_swv_gc_control_dap_processed.csv")
df_2_psoralen <- read_csv("../../01_17_19_psoralen_nonequil_2/Processing/01_17_19_swv_gc_psoralen_dap_processed.csv")

df_2 <- rbind(df_2_control,df_2_psoralen) %>% 
  select(treatment,run, rep, reactor,electrode_from_swv,signal_from_swv,electrode_from_gc, signal_from_gc, echem_from_swv ) %>% 
  mutate(exp_id = '2')

df_1 <- read_csv("../../01_08_19_psoralen_nonequil/Processing/01_08_19_processed_swv_gc_signals.csv") %>% 
  select(treatment,run, rep, reactor,electrode_from_swv,signal_from_swv,electrode_from_gc, signal_from_gc , echem_from_swv) %>% 
  mutate(exp_id = "1")


df_all <- rbind(df_1,df_2,df_3,df_blank)

Examine Linear Fits

First, let’s look at the blank dataset.

ggplot(df_all %>% filter(electrode_from_swv=='i1' & electrode_from_gc=='i2' & treatment=='blank'), 
       aes(x = signal_from_swv, y = signal_from_gc)) + 
  geom_point() + 
  geom_smooth(method='lm')+
  facet_wrap(~exp_id)

You can see that the blank dataset is almost perfectly linear.

Now let’s look at all of the biofilm data. Below the data are broken down by experiment number and treatment. The colors represent each run within each experiment. The scales are free, so you can assess how well each dataset is fit by the linear model. Remember we are very interested to know whether the shape is linear or nonlinear.

ggplot(df_all %>% filter(electrode_from_swv=='i1' & electrode_from_gc=='i2' & treatment!='blank'), 
       aes(x = signal_from_swv, y = signal_from_gc, color = factor(run))) + 
  geom_point() + 
  geom_smooth(method='lm')+
  facet_wrap(treatment~exp_id,scales='free')+
  scale_color_viridis(discrete = T)

# ggplot(df_all %>% 
#          filter(electrode_from_swv=='i1' & electrode_from_gc=='i2') %>% 
#          filter(exp_id=='3' & run==4), 
#        aes(x = signal_from_swv, y = signal_from_gc, color = factor(run), shape = treatment)) + 
#   geom_point() + 
#   geom_smooth(method='lm') +
#   facet_wrap(reactor~exp_id,scales='free')

Overall, I think the data are fit pretty well by the linear models. There are definitely subsets that are not fit well, and we can break those into two categories: obviously bad data, and good data that doesn’t look linear. For experiment 3 psoralen treated biofilm, you can see that runs 2 and 3 do not look like nice monotonic trends. This is likely because the SWV scans were super weird for those runs - therefore we will ignore them later on. On the other hand, exp 1 run 1 for control and psoralen looks pretty clean, but they look quite nonlinear. There are several other runs in the dataset that look somewhat nonlinear, and it is interesting to think about what that might mean. For now though, let’s move forward with the assumption that most of the datasets can be reasonably fit with a linear model.

Given that we can fit with a linear model let’s now compare the slopes of the datasets, by looking at the same plot with fixed scales:

ggplot(df_all %>% filter(electrode_from_swv=='i1' & electrode_from_gc=='i2' & treatment!='blank'), 
       aes(x = signal_from_swv, y = signal_from_gc, color = factor(run))) + 
  geom_point() + 
  geom_smooth(method='lm')+
  facet_grid(treatment~exp_id)+
  scale_color_viridis(discrete = T)

We can see some clear differences in slopes across the datasets. However, keep in mind that exp 2 psoralen and experiment 3 run 4 both conditions were taken with SWVslow settings that will increase the slope due to the different scan rate. Still, it’s clear that there’s a variety of slopes, but also that none of the slopes seem like crazy outliers.

\(D_{ap}\) Estimation

Now, let’s grab the coefficients of those linear fits and calculate \(D_{ap}\) estimates using our known parameters.

Recall: \[ D_{ap} = \frac{(m A \psi)^2 }{S^2 \pi t_p }\]

Here’s our default \(D_{ap}\) calculator function to get from slope, \(m\), to our estimate.

dap_calc <- function(m, t_p=1/(2*300)){
  
  psi <-  0.91
  A <-  0.013 #cm^2
  S <-  18.4 #cm
  d_ap <- (m*A*psi)^2 / (S^2 * pi * t_p)
  
  d_ap
}

Now, we’ll fit all the subsets of data and grab metrics like \(R^2\) and confidence intervals for the slope coefficient. Let’s go ahead and plot all of the \(R^2\) values for the datasets to quantitatively assess the fits.

swv_gc_all <- df_all %>% 
  filter(electrode_from_swv=='i1' & electrode_from_gc=='i2' ) %>% 
  group_by(treatment,exp_id,run,echem_from_swv) %>% 
  do(fit = lm(signal_from_gc~signal_from_swv,data = .))

swv_gc_fit_tidy <- tidy(swv_gc_all,fit,conf.int=T)

swv_gc_fit_glance <- glance(swv_gc_all,fit) %>% 
  select(treatment,exp_id,run,r.squared,adj.r.squared)

swv_gc_fit <- left_join(swv_gc_fit_tidy %>% filter(term=='signal_from_swv'), swv_gc_fit_glance, by = c('treatment','exp_id','run','echem_from_swv'))

ggplot(swv_gc_fit, aes(x = exp_id, y = r.squared, shape=treatment,color = factor(run)))+
  geom_point()+
  scale_color_viridis(discrete = T)

You can see that all of the datasets have very high \(R^2\) values (1 is a perfect line), except for the problematic subsets we identified earlier (exp 3, psoralen, runs 2-3), and a very nonlinear subset (exp 1 control run 1)

Let’s move forward calculating \(D_{ap}\), but we’ll exclude all the models with \(R^2 < 0.9\). First calculate \(D_{ap}\) for all the data that was acquired with SWV fast settings, which has a pulse time of \(t_p=1/600\)

swv_gc_dap <- swv_gc_fit %>% 
  filter(r.squared>0.9) %>% 
  filter(echem_from_swv=='SWV' | echem_from_swv == 'SWV.txt') %>% 
  mutate(dap=dap_calc(m = estimate)) %>% 
  mutate(dap_high = dap_calc(m = conf.high)) %>% 
  mutate(dap_low = dap_calc(m = conf.low)) 
  

swv_gc_dap %>% 
  select(treatment, exp_id, run, echem_from_swv, estimate, conf.low, conf.high, r.squared, dap, dap_high, dap_low) %>% 
  kable() %>% 
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = F)

And now we’ll separately calculate \(D_{ap}\) for the SWV slow subsets, with \(t_p = 1/30\).

swvSlow_gc_dap <- swv_gc_fit %>% 
  filter(r.squared>0.9) %>% 
  filter(echem_from_swv=='SWVslow') %>% 
  mutate(dap=dap_calc(m = estimate,t_p= (1 / 30))) %>% 
  mutate(dap_high = dap_calc(m = conf.high,t_p= (1 / 30))) %>% 
  mutate(dap_low = dap_calc(m = conf.low,t_p= (1 / 30))) 

swvSlow_gc_dap %>% 
  select(treatment, exp_id, run, echem_from_swv, estimate, conf.low, conf.high, r.squared, dap, dap_high, dap_low) %>% 
  kable() %>% 
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),full_width = F)

Blank vs. Biofilm

Although there may not be a clear difference +/- psoralen, let’s examine the dataset in comparison to a dataset acquired with a blank (no biofilm) IDA. This dataset should report on the \(D_{ap}\) of free PYO in solution:

ggplot(swv_gc_dap_estimates, aes(x = treatment, y = dap, color = factor(run))) +
  geom_point() + 
  geom_pointrange(aes(ymin = dap_low, ymax = dap_high))+
  facet_grid(~exp_id,scales='free_x') +
  scale_color_viridis(discrete = T)+
  labs(y=expression(D[ap] (cm^2 / sec)),color='Run Number') 

# ggplot(swv_gc_dap_estimates %>% filter(treatment!='blank' & exp_id=='3'), aes(x = treatment, y = dap, color = factor(run),shape = factor(run))) +
#   geom_point() + 
#   geom_pointrange(aes(ymin = dap_low, ymax = dap_high))+
#   facet_wrap(~exp_id,scales = 'free')+
#   scale_color_viridis(discrete = T)

This plot definitely gives us a different view. It’s clear that the blank \(D_{ap}\) is much higher than the biofilm measurements, and that perspective also shows how relatively consistent our biofilm measurements are. Almost all of the biofilm measurements are between \(1\) and \(5 * 10^{-6} \text{cm}^2 / \text{sec}\), while the solution \(D_{ap}\) is nearly an order of magnitude higher.