dPHZ* PYO titration Dap estimate from Potential Step
09_16_18
library(knitr)
library(kableExtra)
knitr::opts_knit$set(root.dir = '/Users/scottsaunders/OneDrive - California Institute of Technology/Documents/Summer 2018/IDA_biofilms/09_11_18_DapEstimation')
knitr::opts_chunk$set(tidy.opts=list(width.cutoff=60),tidy=TRUE)
knitr::opts_chunk$set(echo = TRUE, message=F,warning=F,fig.align="center")
suppressMessages(library('tidyverse'))
suppressMessages(library('cowplot'))
suppressMessages(library('broom'))
suppressMessages(library(modelr))
suppressMessages(library(viridis))
theme_set(theme_bw())Experiment
The following analysis is based on data acquired during the 09/04/18 IDA biofilm experiment. During this experiment I grew two dPHZ* biofilms and exposed them to PYO. I let the biofilms equilibrate during a transfer step after each soak, as usual. For this experiment I also tried to collect usable potential step scans. The scans were 30 seconds with three steps. First the baseline step was at -170mV, then after 10 seconds we jumped to -300mV and after another 10 seconds it jumped back to -170mV. Everything that follows is only from the first reductive step (from -170 to -300) which should be reducing PYO.
Also note that IDA A and B echem was done using separate potentiostats. IDA A was collected on the CHI 760b and IDA B was collected on the CHI 760e. Apparently the newer 760e can collect data points slightly faster (~0.2ms vs. 0.8ms). I totally messed up IDA B by accidentally unclipping the reference during a run, but I wanted to compare the data to see if the collection rate had a big influence.
Results
files1 <- paste("idaB", dir(path = "/Users/scottsaunders/OneDrive - California Institute of Technology/Documents/Summer 2018/IDA_biofilms/09_11_18_DapEstimation/idaB",
pattern = "*.txt"), sep = "/")
# reads each file (from above) into a dataframe that is
# stored in another data frame with the filename
df1 <- data_frame(filename = files1) %>% mutate(file_contents = map(filename,
~read_csv(., col_names = c("a", "b", "c", "d"))))
potStepIdaB <- unnest(df1) %>% separate(filename, c("subDir",
"filename"), sep = "/") %>% separate(filename, c("strain",
"IDA", "reactor", "PHZadded", "molecule", "condition", "echem",
"rep"), sep = "_") %>% mutate(PHZaddedInt = as.integer(str_extract(PHZadded,
"^[0-9]+"))) %>% mutate(reactor = "M", t = a, current = b) %>%
select(-subDir, -a, -b, -c, -d) %>% filter(echem == "potStep")
files2 <- paste("raw090418potStepIDAa", dir(path = "/Users/scottsaunders/OneDrive - California Institute of Technology/Documents/Summer 2018/IDA_biofilms/09_11_18_DapEstimation/raw090418potStepIDAa",
pattern = "*.txt"), sep = "/")
# reads each file (from above) into a dataframe that is
# stored in another data frame with the filename
df2 <- data_frame(filename = files2) %>% mutate(file_contents = map(filename,
~read_csv(., col_names = c("t", "current"))))
potStepCsigFigs <- unnest(df2) %>% separate(filename, c("subDir",
"filename"), sep = "/") %>% separate(filename, c("strain",
"IDA", "reactor", "PHZadded", "molecule", "condition", "echem",
"rep"), sep = "_") %>% mutate(PHZaddedInt = as.integer(str_extract(PHZadded,
"^[0-9]+"))) %>% mutate(reactor = "C") %>% select(-subDir)
# must export with enough sig figs!!!!???
## make subsets for the first 300ms of the reduction steps
idaAsubsetSoak <- potStepCsigFigs %>% filter(condition == "soak" &
rep == "s5.txt")
idaAsubsetTran <- potStepCsigFigs %>% filter(condition == "tran" &
rep == "s7.txt")
idaAsubset <- bind_rows(idaAsubsetSoak, idaAsubsetTran) %>% filter(t >
10 & t <= 10.5) %>% mutate(t = t - 10) %>% mutate(t0.5 = t^-0.5)
idaBsubsetSoak <- potStepIdaB %>% filter(condition == "soak" &
rep == "s5.txt")
idaBsubsetTran <- potStepIdaB %>% filter(condition == "tran" &
rep == "s7.txt")
idaBsubset <- bind_rows(idaBsubsetSoak, idaBsubsetTran) %>% filter(t >
10 & t <= 10.5) %>% mutate(t = t - 10) %>% mutate(t0.5 = t^-0.5)Potential Step Current Decays (first 300ms)
IDA A
ggplot(idaAsubset %>% filter(t >= 0.002), aes(x = t, y = current,
color = PHZaddedInt)) + geom_point() + scale_color_viridis() +
facet_wrap(~condition, scales = "free") + labs(title = "IDA A",
y = "Current (A)", x = "time (sec)")
Examining the raw data shows that the current decay does in fact seem sensitive to the concentration of PYO that the electrode was soaked in. You can see the challenge of quantifying the transfer step since the noise is quite large compared to the signal.
IDA B
ggplot(idaBsubset %>% filter(t >= 0.002), aes(x = t, y = current,
color = PHZaddedInt)) + geom_point() + scale_color_viridis() +
facet_wrap(~condition, scales = "free") + labs(title = "IDA B",
y = "Current (A)", x = "time (sec)")
The IDA B data is kinda messed up (not sure why 200uM is negative - maybe I changed the params?), but the soak steps look reasonably similar to IDA A.
Integration to get Concentration
Integration should yield the total amount of charge passed through the electrode over the time frame. In theory, converting charge to mol PYO should give the concentration of PYO at the electrode (if we know the relevant volume).
Integration is performed using trapezoid method.
idaAint <- idaAsubset %>% filter(t >= 0.002 & t < 0.3) %>% arrange(t) %>%
group_by(PHZadded, PHZaddedInt, rep, condition) %>% mutate(int = (lead(t) -
t) * (lead(abs(current)) + abs(current))/2) %>% mutate(cumInt = cumsum(int))
idaBint <- idaBsubset %>% filter(t > 0.002 & t < 0.3) %>% arrange(t) %>%
group_by(PHZadded, PHZaddedInt, rep, condition) %>% mutate(int = (lead(t) -
t) * (lead(abs(current)) + abs(current))/2) %>% mutate(cumInt = cumsum(int))Chroncoulometry Plots
IDA A
ggplot(idaAint, aes(x = t, y = cumInt, color = PHZaddedInt)) +
geom_point() + scale_color_viridis() + facet_wrap(~condition,
scales = "free") + labs(x = "time (sec)", y = "Charge (Coulombs)",
title = "Accumulated Charge over time")
Accumulated Charge again seems to scale pretty well with added PYO.
IDA B
ggplot(idaBint, aes(x = t, y = cumInt, color = PHZaddedInt)) +
geom_point() + scale_color_viridis() + facet_wrap(~condition,
scales = "free") + labs(x = "time (sec)", y = "Charge (Coulombs)",
title = "Accumulated Charge over time")
Calculate concentration from charge, compare to PYO added
F=96485.3329
A=0.039 #cm^2
z=0.001 #cm (10um)
idaAintSumm <- idaAint %>%
group_by(PHZadded,PHZaddedInt,rep,condition) %>%
summarise(sumInt=sum(int,na.rm=T)) %>%
mutate(PYOmol=sumInt/(2*F)) %>% #convert integration (coulombs) to moles PYO
mutate(PYOconcMperCm3=PYOmol/(A*z)) %>%
mutate(PYOconcMperL=PYOconcMperCm3*1000)
idaBintSumm <- idaBint %>%
group_by(PHZadded,PHZaddedInt,rep,condition) %>%
summarise(sumInt=sum(int,na.rm=T)) %>%
mutate(PYOmol=sumInt/(2*F)) %>% #convert integration (coulombs) to moles PYO
mutate(PYOconcMperCm3=PYOmol/(A*z)) %>%
mutate(PYOconcMperL=PYOconcMperCm3*1000)\[C=\frac{Q}{nFAz} \]
\(F= 9.6485333\times 10^{4}\) Coulombs/mol (Faraday’s constant)
\(A= 0.039\) cm^2 (electrode area)
\(z= 10\) um (biofilm height/working volume)
So this calculated C value will be sensitive to amount of time we integrate (affects Q) and our value of z. For the soak step, it seems we may have actually chosen good values, because the calculated C is similar to the known amount of added PYO. Estimating C is useful in general, but will be very important for the Dap estimates. See further discussion below of what estimated C might mean…
IDA A
ggplot(idaAintSumm, aes(x = PHZaddedInt, y = PYOconcMperL, color = condition)) +
geom_point() + geom_line() + facet_wrap(~condition, scales = "free") +
labs(title = "IDA A calculated [PYO] vs. added PYO", y = "Calculated [PYO] (M)",
x = "added [PYO] to soak (uM)")
This plot shows that accumulated charge (over this time scale) is a shockingly good proxy for concentration. The estimated concentration in Molar is very close to the amount of PYO (uM) added to the soak. Also the transfer step looks similar to the other electrochemical signals we measure…concentration seems to build up and then saturate. The question is, do we believe that the working volume of the electrode is the same in both cases? In the soak step the volume likely includes the biofilm and the solution, while in the transfer step only the biofilm has PYO in it…so this calculated concentration still isn’t quite how much is in the biofilm, but it’s how much is in the closest 10um to the electrode?
idaAechemPYO <- idaAintSumm %>% filter(condition == "tran") %>%
pull(PYOconcMperL)
tibble(transferStep = c(10, 25, 50, 75, 100, 150, 200), Measured_PYO_Solution = c(0.116,
0.13, 0.134, 0.188, 0.196, 0.342, 0.494), Calculated_PYO_Biofilm = idaAechemPYO *
10^6) %>% kable(digits = 12, format = "html", caption = "IDA A: Comparing PYO measured at electrode (echem) vs. PYO measured in Solution (LC-Ms) in uM") %>%
kable_styling(bootstrap_options = "striped", full_width = FALSE)| transferStep | Measured_PYO_Solution | Calculated_PYO_Biofilm |
|---|---|---|
| 10 | 0.116 | 4.514011 |
| 25 | 0.130 | 2.579644 |
| 50 | 0.134 | 5.954776 |
| 75 | 0.188 | 5.715927 |
| 100 | 0.196 | 2.731453 |
| 150 | 0.342 | 3.012939 |
| 200 | 0.494 | 3.904269 |
That said, if we believe the echem estimates for PYO concentration, they are significantly higher than what we measure in solution by LC-MS, so that’s encouraging. It suggests that the biofilms are still retaining PYO in the uM range at equilibrium.
IDA B
ggplot(idaBintSumm, aes(x = PHZaddedInt, y = PYOconcMperL, color = condition)) +
geom_point() + geom_line() + facet_wrap(~condition, scales = "free") +
labs(title = "IDA B calculated [PYO] vs. added PYO", y = "Calculated [PYO] (M)",
x = "added [PYO] to soak (uM)")
idaBechemPYO <- idaBintSumm %>% filter(condition == "tran") %>%
pull(PYOconcMperL)
tibble(transferStep = c(5, 25, 200), Measured_PYO_Solution = c(0.114,
0.138, 0.488), Calculated_PYO_Biofilm = idaBechemPYO * 10^6) %>%
kable(digits = 12, format = "html", caption = "IDA B: Comparing PYO measured at electrode (echem) vs. PYO measured in Solution (LC-Ms) in uM") %>%
kable_styling(bootstrap_options = "striped", full_width = FALSE)| transferStep | Measured_PYO_Solution | Calculated_PYO_Biofilm |
|---|---|---|
| 5 | 0.114 | 4.665094 |
| 25 | 0.138 | 4.847746 |
| 200 | 0.488 | 2.406896 |
IDA A IDA B5 0.16 0.114 10 0.116 0.12 25 0.13 0.138 50 0.134 0.136 75 0.188 0.124 100 0.196
150 0.342
200 0.494 0.488
Linear fits of Current to get Dap
From the Cottrell Equation: \[I(t)=\frac{nFAC \sqrt{D_{ap}}}{\sqrt{\pi t}}\]
\[I=\frac{nFAC \sqrt{D_{ap}}}{\sqrt{\pi}}t^{-1/2}\]
for I vs. t^-(1/2) \[slope=m=\frac{nFAC \sqrt{D_{ap}}}{\sqrt{\pi}}\]
\[D_{ap}= \pi (\frac{m}{nFAC})^2\]
or for I/A vs. T^-(1/2)
\[D_{ap}= \pi (\frac{m}{nFC})^2\]
This expression shows that the estimated Dap is extremely sensitive to the estimated C (scales with the square). Consequently, if we underestimate concentration by ~30% we would calculate a Dap 2x larger, and overestimating by 30% would decrease Dap by ~2x.
Potential Step Current Decays (300ms) vs. t^-0.5, fitting for slopes.
IDA A
ggplot(idaAsubset %>% filter(condition == "soak" & rep == "s5.txt") %>%
filter(t >= 0.001), aes(x = t^-(0.5), y = current/A, color = PHZaddedInt)) +
geom_point() + scale_color_viridis() + geom_smooth(data = idaAsubset %>%
filter(condition == "soak" & rep == "s5.txt") %>% filter(t >=
0.01 & t <= 0.2), aes(x = t^-0.5, y = current/A, group = PHZaddedInt),
color = "black", method = "lm") + labs(title = "IDA A")
IDA B
ggplot(idaBsubset %>% filter(condition == "soak" & rep == "s5.txt") %>%
filter(t >= 0.001), aes(x = t^-(0.5), y = current/A, color = PHZaddedInt)) +
geom_point() + scale_color_viridis() + geom_smooth(data = idaBsubset %>%
filter(condition == "soak" & rep == "s5.txt") %>% filter(t >=
0.01 & t <= 0.2), aes(x = t^-0.5, y = current/A, group = PHZaddedInt),
color = "black", method = "lm") + labs(title = "IDA A")
Dap Estimates
lmAsubset <- idaAsubset %>% filter(t >= 0.01 & t <= 0.2) %>%
group_by(PHZadded, PHZaddedInt, condition, rep) %>% do(tidy(lm(current/A ~
t0.5, .)))
lmBsubset <- idaBsubset %>% filter(t >= 0.01 & t <= 0.2) %>%
group_by(PHZadded, PHZaddedInt, condition, rep) %>% do(tidy(lm(current/A ~
t0.5, .)))
dapIdaA <- left_join(lmAsubset, idaAintSumm, by = c("PHZadded",
"PHZaddedInt", "condition", "rep")) %>% ungroup() %>% mutate(Dap = pi *
(estimate/(2 * F * PYOconcMperCm3))^2) %>% mutate(DapStdError = pi *
(std.error/(2 * F * PYOconcMperCm3))^2)
dapIdaB <- left_join(lmBsubset, idaBintSumm, by = c("PHZadded",
"PHZaddedInt", "condition", "rep")) %>% ungroup() %>% mutate(Dap = pi *
(estimate/(2 * F * PYOconcMperCm3))^2) %>% mutate(DapStdError = pi *
(std.error/(2 * F * PYOconcMperCm3))^2)IDA A
ggplot(dapIdaA %>% filter(term == "t0.5"), aes(x = factor(PHZaddedInt),
y = Dap, fill = condition)) + geom_col(position = "dodge",
color = "black") + geom_errorbar(aes(ymin = Dap - DapStdError,
ymax = Dap + DapStdError), position = "dodge") + labs(title = "Dap estimated from current decay (IDA A)",
x = "[PYO] soak (uM)")
Here the estimated Dap for the soak step is actually lower than the transfer step, which is not consistent with other measurements of Fernanda’s measurements. So the relationship between the soak and transfer steps did not match, but the estimates for Dap are somewhat close (within about an order of magnitude). See the table at the end.
IDA B
ggplot(dapIdaB %>% filter(term == "t0.5"), aes(x = factor(PHZaddedInt),
y = Dap, fill = condition)) + geom_col(position = "dodge",
color = "black") + geom_errorbar(aes(ymin = Dap - DapStdError,
ymax = Dap + DapStdError), position = "dodge") + labs(title = "Dap estimated from current decay (IDA A)",
x = "[PYO] soak (uM)")
Linear fits of Charge to get Dap
From the integrated Cottrell Equation: \[Q(t)=\frac{2nFAC \sqrt{D_{ap}t}}{\sqrt{\pi}}\]
\[Q(t)=\frac{2nFAC \sqrt{D_{ap}}}{\sqrt{\pi}}t^{-1/2}\]
for Q vs. t^(1/2) \[slope=m=\frac{2nFAC \sqrt{D_{ap}}}{\sqrt{\pi}}\]
\[D_{ap}= \pi (\frac{m}{2nFAC})^2\] The scaling of C is the same as chronoamperometry. Therefore this Dap is also sensitive to our C estimate.
Chronocoulometry plots vs. t^0.5
Try to just fit the later linear portion of the plot (50ms to 250ms).
IDA A
ggplot(idaAint %>% filter(t <= 0.25), aes(x = t^0.5, y = cumInt,
color = PHZaddedInt)) + geom_point() + scale_color_viridis() +
geom_smooth(data = idaAint %>% filter(t >= 0.05 & t <= 0.25),
method = "lm", color = "black", aes(x = t^0.5, y = cumInt,
group = PHZaddedInt)) + facet_wrap(~condition, scales = "free")
IDA B
ggplot(idaBint %>% filter(t <= 0.25), aes(x = t^0.5, y = cumInt,
color = PHZaddedInt)) + geom_point() + scale_color_viridis() +
geom_smooth(data = idaBint %>% filter(t >= 0.05 & t <= 0.25),
method = "lm", color = "black", aes(x = t^0.5, y = cumInt,
group = PHZaddedInt)) + facet_wrap(~condition, scales = "free")
Dap estimates
lmAsubsetQ <- idaAint %>% filter(t >= 0.05 & t <= 0.25) %>% group_by(PHZadded,
PHZaddedInt, condition, rep) %>% do(tidy(lm(cumInt ~ t0.5,
.)))
lmBsubsetQ <- idaBint %>% filter(t >= 0.05 & t <= 0.25) %>% group_by(PHZadded,
PHZaddedInt, condition, rep) %>% do(tidy(lm(cumInt ~ t0.5,
.)))
dapQidaA <- left_join(lmAsubsetQ, idaAintSumm, by = c("PHZadded",
"PHZaddedInt", "condition", "rep")) %>% ungroup() %>% mutate(Dap = pi *
(estimate/(4 * F * A * PYOconcMperCm3))^2) %>% mutate(DapStdError = pi *
(std.error/(4 * F * A * PYOconcMperCm3))^2)
dapQidaB <- left_join(lmBsubsetQ, idaBintSumm, by = c("PHZadded",
"PHZaddedInt", "condition", "rep")) %>% ungroup() %>% mutate(Dap = pi *
(estimate/(2 * F * PYOconcMperCm3))^2) %>% mutate(DapStdError = pi *
(std.error/(4 * F * A * PYOconcMperCm3))^2)IDA A
ggplot(dapQidaA %>% filter(term == "t0.5"), aes(x = factor(PHZaddedInt),
y = Dap, fill = condition)) + geom_col(position = "dodge",
color = "black") + geom_errorbar(aes(ymin = Dap - DapStdError,
ymax = Dap + DapStdError), position = "dodge") + labs(title = "",
x = "[PYO] soak (uM)")
Now the transfer step shows a lower Dap than the soak step, which makes more sense. However, the estimates are at least an order of magnitude lower than the current estimate above.
IDA B
ggplot(dapQidaB %>% filter(term == "t0.5"), aes(x = factor(PHZaddedInt),
y = Dap, fill = condition)) + geom_col(position = "dodge",
color = "black") + geom_errorbar(aes(ymin = Dap - DapStdError,
ymax = Dap + DapStdError), position = "dodge") + labs(title = "",
x = "[PYO] soak (uM)")
Conclusions
Let’s look at all of our different estimates of Dap compared to Dphys (cm^2/s).
tibble(method = c("Scott Potential Step Current", "Fernanda Potential Step",
"Scott GC vs SWV", "Scott Potential Step Charge", "Scott Dphys estimate (z=10um)"),
soak = c(5e-07, 4.7e-06, 8.5e-06, 5.6e-08, NA), transfer = c(2e-06,
7e-07, 1.4e-06, 4.5e-08, 2.5e-09)) %>% kable(digits = 12,
format = "html") %>% kable_styling(bootstrap_options = "striped",
full_width = FALSE)| method | soak | transfer |
|---|---|---|
| Scott Potential Step Current | 5.0e-07 | 2.0e-06 |
| Fernanda Potential Step | 4.7e-06 | 7.0e-07 |
| Scott GC vs SWV | 8.5e-06 | 1.4e-06 |
| Scott Potential Step Charge | 5.6e-08 | 4.5e-08 |
| Scott Dphys estimate (z=10um) | NA | 2.5e-09 |
In some regards the estimates of Dap are pretty consistent. This could be encouraging, because none of the Dap measurements are as low as the Dphys measurement. It is important to remember that the parameters effect these estimates a lot, so to believe these differences we need to think about how we can be completely confident in our choice of parameters.
As mentioned in the last notebook, I think the best way to do this is to repeat some of these measurements on an blank IDA (in Vitro). We should be able to get reasonable values for PYO diffusing in solution (Dap=Dphys).