Primary cortical cells from rat embryo were grown separate plates and transduced with a Per2::Luc AAV virus at DIV5. Then the cells were left maturing till I could detect regular network bursting (from DIV7-8 on). The circadian experiments were performed over the weekend from DIV9-12. Bioluminescence was performed in a lumicycler, sampling every 12min. HD-MEA network activity was measured every 2 hrs for 5-10min.
I synchronized both the 24-well plate for the luminescence and the 6-well HD-MEA plate in parallel using a 30min Forskolin incubation.
Per2::Luciferase Reporter in 24-well
The bioluminescence is scaled across wells. The first and second row showed lower signal which is probably due to differences in PMT sensitivity. Nevertheless, oscillations can be measured also in the control conditions, showing that a simple feeding is able to synchronize the clock.
library(readxl)
library(dplyr)
library(ggplot2)
luc <- readRDS("~/dev/r-projects/general/circadian_hdmea/Bioluc_n1.rds")
luc$group <- factor(luc$group, levels=c("Control", "1uM Fsk", "5uM Fsk", "10uM Fsk"))
ggplot(luc, aes(V1, scaledVal, color=file_name)) + geom_line() + facet_grid(group~Col, scales="free")+ geom_vline(xintercept = c(24, 48), linetype="dashed")+ xlab("hours after synchronization")+
theme_minimal()
We can also detrend the signal to take into account loss in signal over time
fit <- lm(scaledVal ~ V1 + as.factor(well_id), data=luc)
luc$resid_signal <- resid(fit)
ggplot(luc, aes(V1, resid_signal, color=file_name)) + geom_line() + facet_grid(group~Col, scales="free")+ geom_vline(xintercept = c(24, 48), linetype="dashed")+ xlab("hours after synchronization")+
theme_minimal()
Here, we plot the development of burst frequency (how many times the cultures burst within the 5min recording) in Hz over the course of 60hrs.
hdmea <- readRDS("~/dev/r-projects/general/circadian_hdmea/hdmea_n1.rds")
ggplot(hdmea, aes(x=as.numeric(V1), y=burst_frequency_hz, color=as.factor(well_number), fill=as.factor(well_number))) +
stat_summary(geom='line', fun=mean, size=1, alpha=0.7) +
theme_minimal() + facet_wrap(condition~well_number, scales="free") + expand_limits(y=0) + geom_vline(xintercept = c(24, 48), linetype="dashed") + xlab("hours after synchronization")
hdmea$group <- hdmea$condition
hdmea$group <- ifelse(hdmea$group == "10uM_Fsk", "10uM Fsk", "Control")
hdmea$scaledVal <- ave(hdmea$burst_frequency_hz, as.factor(hdmea$well_number), FUN = scale)
#hdmea$scaledVal <- ave(hdmea$burst_peak_firing_rate_90th_percentile_hz, as.factor(hdmea$well_number), FUN = scale)
fit <- lm(scaledVal ~ V1 + as.factor(well_number), data=hdmea)
summary(fit)
Call:
lm(formula = scaledVal ~ V1 + as.factor(well_number), data = hdmea)
Residuals:
Min 1Q Median 3Q Max
-2.4802 -0.5627 -0.2150 0.6689 2.2837
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.739e-01 2.081e-01 2.278 0.023854 *
V1 -1.478e-02 3.823e-03 -3.865 0.000152 ***
as.factor(well_number)2 -6.798e-18 2.377e-01 0.000 1.000000
as.factor(well_number)3 -2.903e-16 2.377e-01 0.000 1.000000
as.factor(well_number)4 2.516e-17 2.377e-01 0.000 1.000000
as.factor(well_number)5 -9.812e-19 2.377e-01 0.000 1.000000
as.factor(well_number)6 -1.913e-16 2.377e-01 0.000 1.000000
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.9656 on 191 degrees of freedom
Multiple R-squared: 0.07253, Adjusted R-squared: 0.04339
F-statistic: 2.489 on 6 and 191 DF, p-value: 0.02425hdmea$resid_signal <- resid(fit)
luc_long <- luc %>%
mutate(hour = floor(as.numeric(V1))) %>%
group_by(group, well_id, hour) %>%
summarise(signal = mean(resid_signal, na.rm = TRUE), .groups = "drop") %>%
mutate(dataset = "luc", well = well_id)
# reshape hdmea
hdmea_long <- hdmea %>%
mutate(hour = round(as.numeric(V1))) %>%
group_by(group, well_number, hour) %>%
summarise(signal = mean(resid_signal, na.rm = TRUE), .groups = "drop") %>%
mutate(dataset = "hdmea", well = as.character(well_number))
# stack (not join)
joined_long <- bind_rows(luc_long, hdmea_long)
summary_long <- joined_long %>%
group_by(dataset, group, hour) %>%
summarise(
mean_signal = mean(signal, na.rm = TRUE),
se_signal = sd(signal, na.rm = TRUE) / sqrt(n()),
.groups = "drop"
)
library(ggplot2)
# distinct colors for datasets
dataset_colors <- c("luc" = "#E41A1C", # red
"hdmea" = "#377EB8") # blue
summary_long$group <- factor(summary_long$group, levels=c("Control", "1uM Fsk", "5uM Fsk", "10uM Fsk"))
ggplot() +
# mean ± SE ribbon (per dataset)
geom_ribbon(
data = summary_long,
aes(x = hour,
ymin = mean_signal - se_signal,
ymax = mean_signal + se_signal,
fill = dataset),
alpha = 0.2,
inherit.aes = FALSE
) +
geom_line(
data = summary_long,
aes(hour, mean_signal, color = dataset, group = dataset),
size = 1.2,
inherit.aes = FALSE
) +
# manual colors for datasets only
scale_color_manual(values = dataset_colors, aesthetics = "color", guide = "legend") +
scale_fill_manual(values = dataset_colors) +
facet_wrap(.~ group, scales = "free_y") +
geom_vline(xintercept = c(24, 48, 72), linetype = "dashed") +
ylab("Residual signal (Z-score)") +
theme_minimal() + ggtitle("Burst Frequency and Per2::Luc-Activity")