Data overview and fits

First, let’s look at all of the SWV signals vs. timestamp (in min) for all of the different experiments. This is only the transfer scans from i1, since i2 is problematic.

ggplot(dap_dataset %>% filter(reactor=='transfer' & electrode_from_swv=='i1') , 
       aes(x = minutes_from_swv, signal_from_swv, color = factor(run))) + 
  geom_point() + 
  facet_wrap(treatment~exp_id, scales = 'free')

Overall it looks pretty good! Now let’s reset the start of each run to ~0 min (exactly 0.6 min because that’s the average time it takes between immersion and starting the scan). Let’s also go ahead and exclude the psoralen data right now and get rid of that one control dataset that was taken with SWVslow parameters (exp 3, run 4).

Once that’s done I can plot all of the datasets fit with the semi-infinite diffusion model:

dap_test <- dap_dataset %>%  
  filter(electrode_from_swv=='i1' & treatment=='control') %>% 
  filter(exp_id!='3' | run!=4)

dap_test_norm <- dap_test %>% 
  group_by(exp_id,reactor,run) %>% 
  mutate(min_time = min(minutes_from_swv)) %>% 
  mutate(norm_time = minutes_from_swv - min_time + 0.6)

ggplot(dap_test_norm %>% filter(reactor=='transfer'), 
       aes(x = norm_time, y = signal_from_swv, color = exp_id)) + 
  geom_point() +
  geom_smooth(method='nls',
              formula=y~b*(x)^-0.5+a,
              method.args=list(start=c(b=0.1,a=1e-7)),
              se=F)+
  facet_wrap(exp_id~run,scales='free')

Unlike for the other biofilm/blank datasets this original model actually fits reasonably well. Let’s see how the finite diffusion model fits:

e = exp(1)

ggplot(dap_test_norm %>% filter(reactor=='transfer'), 
       aes(x = norm_time, y = signal_from_swv, color = exp_id )) + 
  geom_point() +
  geom_smooth(method='nls',
              formula=y~(a*(x)^-0.5)*(2+e^(b*(x)^-0.5))+c,
              method.args=list(start=c(a=0.1,b=0.1,c=0)),
              se=F)+
    facet_wrap(exp_id~run,scales='free')

It looks like it does fit slightly better, although the difference is not that great. Let’s proceed to estimate \(D_m\) with both models for now, although ultimately I will just use the no flux / finite model because we already estimated a \(t_s\) for that one.

Also notice that experiment one looks a little different than the other two experiments. The points decay very quickly, and they never approach zero because of a weird high background. I always thought that dataset was a little weird, so we’ll throw that out later.

Estimate model coefficients

Now lets perform the nls fits and look at the parameter estimates for the two models:

nls_no_flux <- dap_test_norm %>% 
  filter(reactor=='transfer') %>% 
  group_by(exp_id,run) %>% 
  do(tidy(nls(.,formula = signal_from_swv~(a*(norm_time)^-0.5)*(2+e^(b/norm_time))+c,start=list(a=0.1,b=0.1,c=0)))) %>% 
  mutate(model='no_flux')

nls_og <- dap_test_norm %>% 
  filter(reactor=='transfer') %>% 
  group_by(exp_id,run) %>% 
  do(tidy(nls(.,formula = signal_from_swv~(a*(norm_time)^-0.5)+c,start=list(a=0.1,c=0)))) %>% 
  mutate(model='og')

nls_fits <- rbind(nls_no_flux,nls_og)

ggplot(nls_fits %>% filter(term!='c'),aes(x=model,y=estimate,color=factor(run),shape = exp_id))+geom_point()+
  facet_wrap(~term,scales='free')+ylim(0,NA)

Ok, this is very similar to what I saw before. Recall that if the OG and no_flux model had identical fits, they should estimate \(a\) with a 4x difference, which is about what we see because the fits are similar.

Calculate \(D_m\)

With those estimates let’s now get our functions to calculate \(D_m\) from \(a\). Recall that for the no flux / finite model \(D_m = \frac{I_0^2 D_{ap} t_s}{4 \pi a^2}\) and for the original semi infinite model \(D_m = \frac{I_0^2 D_{ap} t_s}{\pi a^2}\).

dm_from_a_noFlux <- function(estimate, i_0, dap, t_s=0.1){
  
  dm <- ( i_0^2 * dap * t_s ) / ( 4*pi* estimate^2)
  
  dm
}

dm_from_a_OG <- function(estimate, i_0, dap, t_s=0.1){
  
  dm <- ( i_0^2 * dap * t_s ) / ( pi* estimate^2)
  
  dm
}

soak_i0 <- dap_test_norm %>% 
  mutate(i0_soak = signal_from_swv) %>% 
  filter(rep==0 & reactor=='soak') %>% 
  select(exp_id,reactor, run, i0_soak) 

transfer_i0 <- dap_test_norm %>% 
  mutate(i0_transfer = signal_from_swv) %>% 
  filter(rep==0 & reactor=='transfer') %>% 
  select(exp_id,reactor, run, i0_transfer) 

dm_estimates_soak <- left_join(nls_fits %>% filter(term=='a'), soak_i0, by = c('exp_id','run'))

dm_estimates <- left_join(dm_estimates_soak, transfer_i0, by = c('exp_id','run'))

dm_estimates_wDap <- left_join(dm_estimates, 
                               dap_estimates %>% 
                                 filter(treatment=='control') %>% 
                                 select(exp_id, run, r.squared, dap, dap_high, dap_low), 
                               by = c('exp_id','run'))

dm_estimates_og <- dm_estimates_wDap %>% 
  filter(model=='og') %>% 
  mutate(dm_soak = dm_from_a_OG(estimate = estimate, i_0 = i0_soak, dap = dap, t_s = 0.008)) %>% 
  mutate(dm_transfer = dm_from_a_OG(estimate = estimate, i_0 = i0_transfer, dap = dap, t_s = 0.008))

dm_estimates_noFlux <- dm_estimates_wDap %>% 
  filter(model=='no_flux') %>% 
  mutate(dm_soak = dm_from_a_noFlux(estimate = estimate, i_0 = i0_soak, dap = dap, t_s = 0.008)) %>% 
  mutate(dm_transfer = dm_from_a_noFlux(estimate = estimate, i_0 = i0_transfer, dap = dap, t_s = 0.008)) 

dm_estimates_tidy <- rbind(dm_estimates_og,dm_estimates_noFlux) %>% 
  rowwise() %>% 
  mutate(mean = mean(c(dm_soak, dm_transfer)))

mean_biofilm_dap <- mean(dm_estimates_tidy$dap, na.rm = T)

ggplot(dm_estimates_tidy, aes(x = exp_id, shape = factor(run))) + 
  geom_point(aes(y = dm_soak),color = 'red') +
  geom_point(aes(y = dm_transfer), color = 'blue') +
  geom_hline(yintercept = mean_biofilm_dap, linetype = 2) +
  facet_wrap(~model) + 
  scale_y_log10(limits = c(1e-9,5e-6))

Comparing Biofilm to Blank IDA

With that type of visualization in hand let’s go back and get the \(D_m\) and \(D_{ap}\) data for the blank IDA. We can organize it in the same way and then show the blank and biofilm datasets on the same plot / scale.

d_blank <- left_join(df_2blank_estimates %>% filter(biofilm=='blank' & model=='no_flux'),
                     dap_estimates %>% filter(treatment=='blank'), 
                     by = c('biofilm'='exp_id')) %>% 
  mutate(dm = dm_soak) %>% 
  gather(key=d_cat,val=d_est,dap,dm) %>% 
  mutate(ymax = if_else(d_cat=='dap', dap_high, dm_soak)) %>% 
  mutate(ymin = if_else(d_cat=='dap', dap_low, dm_soak)) %>% 
  mutate(exp_id = treatment) %>% 
  distinct(d_est,.keep_all = T)

#dap_estimates

d_combined_wBlank <- bind_rows(d_combined, d_blank)

ggplot(d_combined_wBlank, 
       aes(x = d_cat, y = d_est, fill = exp_id)) + 
  geom_pointrange(aes(ymin = ymin, ymax = ymax), position = position_jitter(height = 0,width=0.2),shape=21)+
  facet_wrap(~biofilm, scales='free') + 
  scale_y_log10(limits = c(1e-9,1e-4))+
  scale_fill_viridis(discrete = T) +
  scale_x_discrete(breaks = c('dap','dm'), labels=c(expression(D[ap]),expression(D[m]))) +
  labs(x='',y=expression("Estimate ("~cm^2 / sec~")"))+
  guides(fill=F) + 
  theme(axis.text.x = element_text( size=14) )

#ggsave(filename = "../Figures/01_30_19_estimated_Dap_Dm_combined.pdf", width = 6, height = 4)

Here you can now get perspective on how well our pipeline estimates \(D_m\) and \(D_{ap}\) for the blank IDA, where we think \(D_{ap} \approx D_m\), while also seeing how the biofilm compares.

This figure summarizes a lot of our conclusions, and I plan for it to be a panel in a figure for the paper.

You can see that:

• the blank \(D_m \approx D_{ap}\)
• The biofilm \(D_{ap}\) is significantly less than the blank.
• The biofilm \(D_m\) is significantly lower than \(D_{ap}\)