First pass

So first let’s just guess some parameters that give a curve that’s close to our dataset. After playing around and re-reading what the parameters mean I started using A1 = 0.02; lrc1 = 16.5; A2 = 0.02; lrc2 = 15.7:

df_1 <- df %>% filter(exp_num == 2 & run == 2) %>% filter(time >= 
    8e-09) %>% filter(bg_sub > 0)

df_guess <- bind_cols(df_1, biexp = SSbiexp(df_1$time, 0.02, 
    16.5, 0.02, 15.7))

ggplot(df_guess, aes(x = time)) + geom_path(aes(y = biexp)) + 
    geom_point(aes(y = bg_sub), shape = 21) + geom_hline(yintercept = 0, 
    linetype = 2, color = "light gray")

I’m glossing over this, but I actually did a lot of guess and check to get to this point! Importantly the lrc parameters are roughly \(\log(\frac{1}{half-life})\) (FYI that’s a natural log base e…, and half life is in seconds). The A1 and A2 values add up to equal the max height of the data. So A1 = 0.02, A2 = 0.02 gives a height of 0.04.

The other two essential data processing things are that you need to exclude all the data < 0…but just after zero there’s also a rise in signal and we need to exclude that too so that we only observe the decay. That’s why I excluded data < 8 ns. Also, some samples have values that go below zero, which throws an error, so I filtered those out too.

Ok, now we’re ready to fit the model using nonlinear least squares (nls()).

# ggplot(df_1, aes(x = time, y = bg_sub)) + geom_point(shape
# = 21) + geom_smooth(method='nls', formula=y~SSbiexp(x, A1,
# lrc1, A2, lrc2), method.args=list(start=c(A1 = 0.02, lrc1 =
# 16.5, A2 = 0.02, lrc2 = 15.7)), se=F)

A1 = 0.02
lrc1 = 16.5
A2 = 0.02
lrc2 = 15.7

mod_1 <- nls(bg_sub ~ SSbiexp(time, A1, lrc1, A2, lrc2), data = df_1)

summary(mod_1)

## 
## Formula: bg_sub ~ SSbiexp(time, A1, lrc1, A2, lrc2)
## 
## Parameters:
##       Estimate Std. Error t value Pr(>|t|)    
## A1   4.003e-02  4.397e-04   91.03   <2e-16 ***
## lrc1 1.777e+01  2.110e-02  842.11   <2e-16 ***
## A2   1.077e-02  2.370e-04   45.47   <2e-16 ***
## lrc2 1.493e+01  3.608e-02  413.93   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0006689 on 233 degrees of freedom
## 
## Number of iterations to convergence: 0 
## Achieved convergence tolerance: 6.319e-07

Seems reasonable…idk. Let’s look at the actual fit graphically:

df_mod_1 <- df_1 %>% add_predictions(mod_1)

ggplot(df_mod_1, aes(x = time)) + geom_path(aes(y = pred)) + 
    geom_point(aes(y = bg_sub), shape = 21)

Looks beautiful! Let’s move forward to the whole dataset.

All datasets

Obviously I iterated back and forth between the previous section and this one before it actually worked for the whole thing…that’s how I found out about the negative value issue. Anyway, we’re going to select the subset of the data we’re interested in, which is experiment 2. Notice that again I filter out negative values, times before 8ns, and we are only looking at runs 1 - 9, because 10 is essentially noise in the background subtraction.

I’m using the same parameters shown above.

df_exp_2 <- df %>% filter(exp_num == 2 & run %in% 1:9) %>% filter(time >= 
    8e-09) %>% filter(bg_sub > 0) %>% group_by(run, quencher_eq) %>% 
    nest()

fit_biexp_bg_sub <- function(df) {
    
    A1 = 0.02
    lrc1 = 16.5
    A2 = 0.02
    lrc2 = 15.7
    
    mod <- nls(bg_sub ~ SSbiexp(time, A1, lrc1, A2, lrc2), data = df)
    
    mod
}

df_exp_2_models <- df_exp_2 %>% mutate(models = map(data, fit_biexp_bg_sub))

df_exp_2_preds <- df_exp_2_models %>% mutate(preds = map2(data, 
    models, add_predictions)) %>% unnest(preds)

df_exp_2_ests <- df_exp_2_models %>% mutate(ests = map(models, 
    tidy, conf.int = T)) %>% unnest(ests)


df_exp_2_preds %>% head()

## # A tibble: 6 x 17
##     run quencher_eq  date   rep compressed    time intensity material
##   <dbl>       <dbl> <dbl> <dbl> <chr>        <dbl>     <dbl> <chr>   
## 1     1          NA 82919    13 compresse… 8.01e-9    0.0862 biofilm 
## 2     1          NA 82919    13 compresse… 8.21e-9    0.0856 biofilm 
## 3     1          NA 82919    13 compresse… 8.41e-9    0.0849 biofilm 
## 4     1          NA 82919    13 compresse… 8.71e-9    0.0835 biofilm 
## 5     1          NA 82919    13 compresse… 9.01e-9    0.0818 biofilm 
## 6     1          NA 82919    13 compresse… 9.21e-9    0.0807 biofilm 
## # … with 9 more variables: material_id <dbl>, ruphen_conc <dbl>,
## #   quencher_conc <dbl>, quencher_id <chr>, comment <lgl>, exp_num <dbl>,
## #   min <dbl>, bg_sub <dbl>, pred <dbl>

Ok, we fit all the models individually to the different runs in the experiment and added the predictions from the model (e.g. best fit line), and we also have the parameter estimates:

# write_csv(df_exp_2_ests,
# '2019_09_27_spectroscopy_fits_bg_sub.csv')

df_exp_2_ests %>% kable() %>% kable_styling() %>% scroll_box(height = "400px")

Let’s see how the fits look:

ggplot(df_exp_2_preds, aes(x = time, group = run, color = quencher_eq)) + 
    geom_path(aes(y = pred), color = "black") + geom_point(aes(y = bg_sub), 
    alpha = 0.2, shape = 21) + xlim(0, 5e-07) + facet_wrap(~quencher_eq)

Nice! All of the fits seem totally reasonable. Note that the plot with quencher eq as NA is the biofilm only background. Overlaying them, we can see the same pattern as before:

ggplot(df_exp_2_preds, aes(x = time, group = run, color = quencher_eq)) + 
    geom_path(aes(y = pred)) + geom_point(aes(y = bg_sub), alpha = 0.2, 
    shape = 21) + xlim(0, 5e-07)

Notice the curves don’t overlap at all, except the gray one, which is the biofilm background.

Now, let’s look at the actual parameter estimates:

background_ests <- df_exp_2_ests %>% filter(is.na(quencher_eq))

# background_ests

ggplot(df_exp_2_ests, aes(x = quencher_eq, y = estimate)) + geom_hline(data = background_ests, 
    aes(yintercept = estimate), linetype = 2, color = "light gray") + 
    geom_pointrange(aes(ymin = conf.low, ymax = conf.high)) + 
    facet_wrap(~term, scales = "free")

Note that the dotted lines are the values for the biofilm only control.

Ok, so let’s break this down. To me this means that there are two components of different magnitudes - A1 is about 4 times larger than A2. The component described by A1 decays pretty fast with a half life of ~ 20 ns, and the the A2 signal decays much more slowly with a half life of ~ 150ns. You can see that the magnitudes of both A1 and A2 go down significantly as we add more quencher (note bars are 95% CI). Meanwhile, the rate constants do not really change significantly. There seem to be some statistically significant differences, but we’re talking about tiny effect sizes. Let’s see the stern volmer plot, but I think this supports what we originally thought - quenching seems to occur statically.

Just for fun let’s replot the fit data with those estimates of the lifetimes.

ggplot(df_exp_2_preds, aes(x = time, group = run, color = quencher_eq)) + 
    geom_vline(xintercept = c(2e-08, 1.5e-07), color = "light gray", 
        size = 2) + geom_path(aes(y = pred)) + xlim(0, 5e-07)

Seems reasonable!

First pass

let’s guess and check first:

df_raw <- df %>% filter(exp_num == 2 & quencher_eq == 0) %>% 
    filter(time >= 8e-09)

df_raw_guess <- bind_cols(df_raw, biexp = SSbiexp(df_raw$time, 
    0.02, 16.5, 0.02, 15.7))

ggplot(df_raw_guess, aes(x = time)) + geom_path(aes(y = biexp)) + 
    geom_point(aes(y = intensity), shape = 21)

This is actually the same guess we used above…it looks pretty bad, but right order of magnitude. Let’s see what happens when we fit and plot:

# ggplot(df_raw, aes(x = time, y = intensity)) +
# geom_point(shape = 21) + geom_smooth(method='nls',
# formula=y~SSbiexp(x, A1, lrc1, A2, lrc2),
# method.args=list(start=c(A1 = 0.05, lrc1 = 16.5, A2 = 0.05,
# lrc2 = 15.7)), se=F)

A1 = 0.02
lrc1 = 16.5
A2 = 0.02
lrc2 = 15.7

mod_raw_1 <- nls(intensity ~ SSbiexp(time, A1, lrc1, A2, lrc2), 
    data = df_raw)

df_mod_raw_1 <- df_raw %>% add_predictions(mod_raw_1)

ggplot(df_mod_raw_1, aes(x = time)) + geom_path(aes(y = pred)) + 
    geom_point(aes(y = intensity), shape = 21)

Looks perfect. This little exercise demonstrates how robust this fitting process is. Even with a bad guess, the model rapidly converges to what looks like a right answer. We really just need to give it starting parameters that are in the right order of magnitude.

All datasets

This time we’ll be able to fit all the datasets (runs 1 - 10), because they all have a decay dominated by the biofilm background.

df_exp_2_raw <- df %>% filter(exp_num == 2 & run %in% 1:10) %>% 
    filter(time >= 8e-09) %>% filter(intensity > 0) %>% group_by(run, 
    quencher_eq) %>% nest()

fit_biexp_raw <- function(df) {
    
    A1 = 0.02
    lrc1 = 16.5
    A2 = 0.02
    lrc2 = 15.7
    
    mod <- nls(intensity ~ SSbiexp(time, A1, lrc1, A2, lrc2), 
        data = df)
    
    mod
}

df_exp_2_raw_models <- df_exp_2_raw %>% mutate(models = map(data, 
    fit_biexp_raw))

df_exp_2_raw_preds <- df_exp_2_raw_models %>% mutate(preds = map2(data, 
    models, add_predictions)) %>% unnest(preds)

df_exp_2_raw_ests <- df_exp_2_raw_models %>% mutate(ests = map(models, 
    tidy, conf.int = T)) %>% unnest(ests)


df_exp_2_raw_preds %>% head()

## # A tibble: 6 x 17
##     run quencher_eq  date   rep compressed    time intensity material
##   <dbl>       <dbl> <dbl> <dbl> <chr>        <dbl>     <dbl> <chr>   
## 1     1          NA 82919    13 compresse… 8.01e-9    0.0862 biofilm 
## 2     1          NA 82919    13 compresse… 8.21e-9    0.0856 biofilm 
## 3     1          NA 82919    13 compresse… 8.41e-9    0.0849 biofilm 
## 4     1          NA 82919    13 compresse… 8.71e-9    0.0835 biofilm 
## 5     1          NA 82919    13 compresse… 9.01e-9    0.0818 biofilm 
## 6     1          NA 82919    13 compresse… 9.21e-9    0.0807 biofilm 
## # … with 9 more variables: material_id <dbl>, ruphen_conc <dbl>,
## #   quencher_conc <dbl>, quencher_id <chr>, comment <lgl>, exp_num <dbl>,
## #   min <dbl>, bg_sub <dbl>, pred <dbl>

Again we have predictions for our best fit line. And below we have our parameter estimates:

# write_csv(df_exp_2_raw_ests,
# '2019_09_27_spectroscopy_fits_raw.csv')

df_exp_2_raw_ests %>% kable() %>% kable_styling() %>% scroll_box(height = "400px")

Let’s take a look:

ggplot(df_exp_2_raw_preds, aes(x = time, group = run, color = quencher_eq)) + 
    geom_path(aes(y = pred), color = "black") + geom_point(aes(y = intensity), 
    alpha = 0.2, shape = 21) + xlim(0, 5e-07) + facet_wrap(~quencher_eq)

Looks perfect! Overlaying gives:

ggplot(df_exp_2_raw_preds, aes(x = time, group = run, color = quencher_eq)) + 
    geom_path(aes(y = pred)) + geom_point(aes(y = intensity), 
    alpha = 0.2, shape = 21) + xlim(0, 5e-07)

Ok, now let’s see what the actual parameter estimates are:

background_ests <- df_exp_2_raw_ests %>% filter(is.na(quencher_eq))

# background_ests

ggplot(df_exp_2_raw_ests, aes(x = quencher_eq, y = estimate)) + 
    geom_hline(data = background_ests, aes(yintercept = estimate), 
        linetype = 1, color = "light gray") + geom_hline(data = background_ests, 
    aes(yintercept = conf.low), linetype = 2, color = "light gray") + 
    geom_hline(data = background_ests, aes(yintercept = conf.high), 
        linetype = 2, color = "light gray") + geom_pointrange(aes(ymin = conf.low, 
    ymax = conf.high)) + facet_wrap(~term, scales = "free")

Let’s walk through it. Again, there are two components of different magnitude - A1 is about 10x larger than A2. A1 has a half life of about 15ns, and A2 has a half life of ~ 250ns, but there’s quite a bit of variation. A1 is probably the biofilm background signal, and interestingly it’s amplitude decreases significantly (~30%), while it’s lifetime barely changes. The A2 component may be the Ru signal - It’s amplitude decreases a lot (~80%), but there’s also a significant change in lifetime from about 250ns (lrc2 = 15) to about 100ns (lrc2 = 16). I’m not exactly sure what to make of this. I would say that the amplitude decrease is larger and more confident than the lifetime changes, but it may be consistent with a combination of static and dynammic quenching.

ggplot(df_exp_2_raw_preds, aes(x = time, group = run, color = quencher_eq)) + 
    geom_vline(xintercept = 1.5e-08, color = "light gray", size = 2) + 
    geom_rect(aes(xmin = 1e-07, xmax = 2.5e-07, ymin = 0, ymax = Inf), 
        fill = "light gray", color = NA) + geom_path(aes(y = pred)) + 
    xlim(0, 5e-07)

Yeah, when we visualize it like this it looks like there’s a pretty big range in that second lifetime…