Abstract

In this project, we analyze the effect of visual stimuli on the neural activity of neurons in the visual cortex of mice. In particular, determine the nature of the relationship between the stimuli and neural activity, and we build a predictor model to predict the outcome of trials based on visual stimuli and neural activity. The data used in this project comes from the novel study by Steinmetz et al. (2019). Throughout the analysis, this study addresses some of the unique features and challenges in working with this dataset. It was found that the mean firing rate of each trial has an interactive effect with the levels of visual stimuli presented to the left and/or right side. Furthermore, the predictive model was found to perform moderately well given the data at hand. Further research using a more robust set of factors are suggested in order to make concrete conclusions regarding the neural activity as a result of visual stimuli.

Introduction

In this project, we study the effect of visual stimuli on the activity of neurons in the visual cortex of mice. The dataset used in this project is a subset of the data collected by Steinmetz et al. (2019) in a study where 10 mice were randomly presented with some visual stimuli presented on the left and/or right side over 39 sessions each comprising of several hundred trials, and a reward was awarded to the mice based on the success of their resulting action. In particular, only use data from five sessions (Sessions 1 to 5) from two mice - Cori (Sessions 1 to 3) and Frossman (Sessions 4 to 5), and we specifically focus on the spike trains of neurons in the visual cortex in a 0.4 second window from onset of the stimuli.

The main goal of this project is to understand how the neural activity in the visual cortex is affected by the different levels of stimuli on the left and right, and how this information can be used to predict the outcome of the trial. As such, the two specific questions of interest are:

1. How do neurons in the visual cortex respond to the stimuli presented on the left and right?

2. How can the outcome of each trial be predicted using the neural activities and stimuli?

Background

The study conducted by Steinmetz et al. (2019) was novel as a result of the development and application of the Neuropixels probes which enabled the researchers to record neural activity at a much higher resolution and frequency than was possible before. In their study, experiments were conducted on 10 mice over the span of 39 sessions consisting of few hundred trials. The trials consisted of two screens placed on either side of the mouse, who was secured in place with the Neuropixels probes inserted into various regions of the left hemisphere of the brain. During each trial, visual stimuli were shown on either or both screens in varying levels of contrasts (0%, 25%, 50%, and 100%). The task of the mouse was to turn a wheel in the direction from the higher contrast to the lower contrast. For example: If the right screen displayed a higher contrast compared to the left screen, the trial was considered a success if the mouse turned the wheel to the left (counter-clockwise), moving the contrast point from the right to the center. Upon a successful trial, the mouse was rewarded with a sip of water via an attached delivery tube. On trials where the two contrasts were equal, a reward was randomly awarded irrespective of the action of the mouse. From start to end, each trial lasted 1.5 seconds. Throughout the entire process of the trial, the Neuropixels probes were recording the neural activity in the inserted brain regions.

Throughout the 39 session with 10 different mice, the Neuropixels probes were inserted into various regions of the left hemisphere, in different depths and angles. This allowed for a robust dataset to be built covering each of the 42 regions of the brain.

In this project, we only focus on the neural activity in the visual cortex of the brain. The dataset we use is a subset of the original dataset collected by Steinmetz et al. (2019), consisting of the spike trains spanning 0.4 seconds after stimulus onset, for trials conducted on the two mice Cori and Frossman over five sessions. In order to address the first question of interest, we conduct ANOVA analysis on our data to determine the nature of the relationship between neural activity and visual stimuli. In particular, we test whether the relation between left and right contrast and neural activity is additive or interactive. For the second question of interest, our goal is to build a prediction model that is optimized for the specificity and sensitivity of the prediction.

Descriptive analysis

In this project, the primary statistical unit is each individual trial, where a random visual stimuli was presented to a mouse, following which the resulting neural activity and outcome of trial were recorded. In our dataset, the five variables available for each trial are:

• feedback_type: type of the feedback, 1 for success and -1 for failure
• contrast_left: contrast of the left stimulus
• contrast_right: contrast of the right stimulus
• time: centers of the time bins for spks
• spks: numbers of spikes of neurons in the visual cortex in time bins defined in time

To get a better understanding of the data we are dealing with, we can visualize a sample spike train from an individual trial. In Fig 1 below, the x-axis represents the bins of 0.01 seconds each (spanning 0.4 seconds total), and the y-axis represents the 178 neurons that were observed in Session 1. The shaded cells correspond to the activity of each neuron during each 0.01s time interval, measured by the number of spikes. We can see for this particular trial, the number of spikes ranged between 0 to 4.

temp = reshape::melt(session[[1]]$spks[[1]])
names(temp) = c("neuron", "time", "spikes")

ggplot(temp, aes(time, neuron, fill = spikes)) +
  geom_tile() +
  scale_fill_gradient(low="white", high="blue4") +
  theme_linedraw() +
  labs(x = "Time Bins", y = "Neurons",title = "Sample Spike Train", subtitle = "Session 1, Trial 1")

Fig 1: Sample Spike Train

In addition, we also observe that a majority of the neurons are mostly inactive, and there is no obvious discernible pattern we can identify from the spike train. That leads us to one of the most interesting characteristics of this data. Since the Neuropixels probes were inserted in different parts of the brain in each session (Steinmetz et al., 2019), i.e. different parts of the visual cortex in our case, the number of neurons (see Tab 1) and the actual neurons recorded vary from session to session. This gives us a wider coverage of the visual cortex, however it makes handling and analyzing the data more challenging.

count.neurons = numeric(5)
names(count.neurons) = paste("Session", 1:5)

for (ID in 1:5) {
  count.neurons[ID] = dim(session[[ID]]$spks[[1]])[1]
}

kbl(t(count.neurons), caption = "Tab 1: Number of Neurons per session") %>%
  kable_styling(full_width = F)

In order to proceed with our analysis, we need to first define our outcome variable corresponding to the neural activity in the visual cortex. We can not use the activity of each individual neuron, since the number of neurons are different across sessions, and the neurons themselves are not uniquely identified. For this project, we using the Mean Firing Rate (MFR) as the outcome variable. The MFR for each trial is defined as:

\[ MFR = \frac{\text{sum of FR across all neurons}}{\text{total number of neurons}} \]

where FR is the Firing Rate of each neuron over the 0.4s window of each trial, and is defined as:

\[ FR = \frac{\text{total number of spikes}}{\text{0.4}} \]

The MFR should be a justifiable measure of the neural activity since it is essentially the average neural activity for each trial. Using this as our outcome variable condenses the dimension of our data, and averages out any noise or outliers in the spikes.

As such, we build the dataset for our project with the MFR for each trial. In addition, we encode the feedback_type variable such that 1 corresponds to success and 0 corresponds to failure, and we include a factor variable for trial_type, which identifies each trial by the combination of left and right contrasts.

kbl(head(d1), caption = "Tab 2: First 6 observations of dataset")

To understand the distribution of Mean Firing Rate, we can observe the plots in Fig 2.

d = ggplot(aes(mfr_all),data = d1) +
  geom_density(fill="steelblue", alpha=0.8) +
  labs(x = "Mean Firing Rate", y = "Density", title = "(a) Density Plot of MFR")

bp_mfr = ggplot(d1, aes(session, mfr_all)) +
  geom_boxplot(aes(fill = name)) +
  labs(x = "Session", y = "MFR", title = "(b) MFR by Session")

plot_mfr = ggplot(d1, aes(1:nrow(d1), mfr_all)) +
  geom_point(aes(color = session), shape = "x", alpha = 0.8) +
  labs(x = "Trial Index", y = "MFR", title = "(c) MFR by Trial")

grid.arrange(d, bp_mfr, plot_mfr, layout_matrix = matrix(c(1, 2,
                   1, 3), ncol = 2, byrow = TRUE))

Fig 2: Distribution of MFR

From Fig 2(a), we observe that MFR appears to be mostly normally distributed, with a slight positive skew. This is consistent with our initial observation from the sample spike train (Fig 1) that most neurons nearly inactive, with very few showing high activity. From Fig 2(b), we observe that the MFR appears to be decreasing as the sessions progress. In fact, as we observe from Fig 2(c), the MFR also presents a decreasing trend with respect to the trials. This may be due to the mice getting more tired as the trials proceed, and therefore presenting less activity in their visual cortex. However, that does not explain the difference between the MFR of Cori and Frossman, which may just be due to the inherent difference between the two mice. We will need to take these findings into account when building our models.

ggplot(d1, aes(contrast_right, contrast_left, fill = mfr_all)) + 
  geom_tile() +
  facet_wrap(~session) +
  scale_fill_gradient(low="blue4", high="orange") +
  labs(x = "Right Contrast", y = "Left Contrast", fill = "MFR", title = "MFR by Combination of Contrasts across Sessions")

Fig 3: Heatmap of MFR by Combination of Contrasts across Sessions

In Fig 3, we again see the decreasing trend of MFR across sessions, with the overall shade of the heatmaps getting darker. In-spite of that we can still assess if there appears to be any relation between MFR and combinations of contrasts. Intuitively, a visual stimuli with more contrast should result in more neural activity, and as such we would expect to see lighter shades towards the right and top edges of the heatmaps and darker shades towards the bottom left. However, we don’t see that pattern reflected in the heatmaps. We will look into this further after our primary analysis.

Inferential Analysis

For our analysis, we define a mixed effect model with two interacting fixed-effect factors for left contrast and right contrast, and a random intercept for each session. The left and right contrasts are taken as fixed effects since their levels were selected by design, and we are only concerned with analyzing the effect on neural activity at these levels of visual stimuli. The each session is taken as a random effect since the five sessions we are studying are a subset of the 39 total sessions in the experiment, and we want to be able to generalize our findings beyond just the five sessions in our data. The defined full model is as follows:

\[ Y_{ijkl} = \mu_{...} + \alpha_i + \beta_j + (\alpha\beta)_{ij} + \gamma_k + \epsilon_{ijkl}, \\ l = 1, \dots, n,\ k=1, \dots, 5,\ j=1,\dots, 4,\ i=1,\dots, 4 \]

where

• The population mean of the MFR is \(\mu_{...}\)

• The fixed-effect of left contrast with level \(i\) is \(\alpha_i\), with constraint \(\sum \alpha_i =0\)

• The fixed-effect of right contrast with level \(j\) is \(\beta_j\), with constraint \(\sum \beta_j = 0\)

• The interaction effect of left contrast \(i\) and right contrast \(j\) is \((\alpha \beta)_{ij}\), with constraints \(\sum_i(\alpha\beta)_{ij} = \sum_j(\alpha\beta)_{ij} = 0\)

• The random intercept corresponding to the \(k\)th session is \(\gamma_k\), with constraint \(\gamma_k \sim_{i.i.d.} N(0,\sigma_{\gamma}^2)\)

• The random error term \(\epsilon_{ijkl} \sim_{i.i.d.} N(0, \sigma^2)\)

• All random variables are mutually independent.

• \(Y_{ijkl}\) represents the mean firing corresponding to the \(l\)th trial in the \(k\)th session with right contrast \(j\) and left contrast \(i\).

options(contrasts = c("contr.treatment", "contr.poly"))

fit0.1 = lmer(mfr_all ~ contrast_left*contrast_right + (1 | session), d1)
# summary(fit0.1)

kbl(anova(fit0.1), caption = "Tab 3: ANOVA table for fixed effects")

Upon fitting the proposed model, we observe (Tab 3) that all three fixed-effect terms are significant at a level of 5%. We are not reporting the individual fixed-effect coefficients since the list would be long, and it is not of interest for us to see how MFR is affected by different contrast levels. Instead, we will now test whether the relation between MFR and the visual stimuli indeed has an interactive effect, oi is it simply additive. To that end, we fit a reduced model without the interaction term, and conduct a test for \(H_0: (\alpha \beta)_{ij} = 0\) at a 5% level of significance.

fit0.2 = lmer(mfr_all ~ contrast_left + contrast_right + (1 | session), d1)
# summary(fit0.2)
# anova(fit0.2)

kbl(anova(fit0.1, fit0.2), caption = "Tab 4: Output of Anova test. Full Model: fit0.1, Reduced Model: fit0.2")

With a p-value of 0.04, we reject \(H_0\) at a 5% level of significance. Therefore, we conclude at 5% level of significance that the visual stimuli and MFR has an interactive effect.

Sensitivity analysis

Model Diagnostics

From our initial descriptive analysis of the data, we had decided to include a random intercept for each session into our model. However, the validity of that choice needs to be verified. As such, we fit another reduced model without the random-effect term, and compare that with our full model to test \(H_0: \gamma_k = 0\) at 5% level of significance. From the output of the test (tab 5) we see a p-value much smaller than 0.05. Therefore, we reject \(H_0\), and conclude at 5% level of significance that the random intercept is necessary in our model.

fit0.3 = aov(mfr_all ~ contrast_left*contrast_right, d1)
# summary(fit0.3)

kbl(anova(fit0.1, fit0.3), caption = "Tab 5: Output of Anova test. Full Model: fit0.1, Reduced Model: fit0.3")

Tab 5: Output of Anova test. Full Model: fit0.1, Reduced Model: fit0.3

	npar	AIC	BIC	logLik	deviance	Chisq	Df	Pr(>Chisq)
fit0.3	17	3787.622	3874.097	-1876.811	3753.622	NA	NA	NA
fit0.1	18	2349.779	2441.341	-1156.890	2313.779	1439.843	1	0

In order to validate our model further, we can perform model diagnostics using some diagnostic plots. The Residual vs. Fitted Values plot (Fig 4) shows no obvious signs of violation of our assumptions of homoskedasticity. And the Normal QQ Plot (Fig 5) presents a mostly straight line, showing our assumption of Normality is justified.

plot(fit0.1, xlab = "Fitted Values", ylab = "Residuals", main = "Residuals vs. Fitted Values")

Fig 4: Residual vs. Fitted Values Plot

qqnorm(resid(fit0.1), main = "Residual QQ Plot")

Fig 5: Normal QQ Plot

As such, our current model appears to be adequate with no obvious pitfalls. However, we can explore some alternative models to address some of the concerns and observations we had during our initial descriptive analysis.

Alternate Models

Mean Scaled Firing Rate

for (ID in 1:5) {
  n.neurons = dim(session[[ID]]$spks[[1]])[1]
  
  s[[ID]]$sfr_neurons = lapply(s[[ID]]$fr_neurons, function(x) {(x - min(x))/max(x)})
  s[[ID]]$msfr_all = sapply(s[[ID]]$sfr_neurons, function(x) {sum(x)/n.neurons})
}

d1$msfr_all = c(s[[1]]$msfr_all, s[[2]]$msfr_all, s[[3]]$msfr_all, s[[4]]$msfr_all, s[[5]]$msfr_all)

As we had identified previously, the MFR appears to have a decreasing trend as trial progress across all sessions. Even though we are accounting for the random effects of sessions in our model, the effect of trials remain. One possible way to address that may be to scale the firing rates of all neurons for each individual trial, and then take the mean of the scaled firing rates (MSFR) as the target variable. The scaling can be done using various methods, and for our purposes we chose a min/max scaler. The distribution of the scaled variable can be seen in Fig 6. We observe that the overall distribution of MSFR (Fig 6(a)) is comparable to that of MFR, being a slightly positively skewed distribution. Whereas from Fig 6(b) and Fig 6(c) we observe that while there is still a difference in MSFR across sessions, the decreasing trend that was present across trials within each session is now much reduced.

d = ggplot(aes(msfr_all),data = d1) +
  geom_density(fill="steelblue", alpha=0.8) +
  labs(x = "Mean Scaled Firing Rate", y = "Density", title = "(a) Density Plot of MSFR")

bp_msfr = ggplot(d1, aes(session, msfr_all)) +
  geom_boxplot(aes(fill = name)) +
  labs(x = "Session", y = "MSFR", title = "(b) MSFR by Session")

plot_msfr = ggplot(d1, aes(1:nrow(d1), msfr_all)) +
  geom_point(aes(color = session), shape = "x", alpha = 0.8) +
  labs(x = "Trial Index", y = "MSFR", title = "(c) MSFR by Trial")

grid.arrange(d, bp_msfr, plot_msfr, layout_matrix = matrix(c(1, 2,
                   1, 3), ncol = 2, byrow = TRUE))

Fig 6: Distribution of Mean Scaled Firing Rate

Observing the heatmaps in Fig 7, it appears that there is now slightly more discernible difference in the neural activity for higher contrast stimuli (compared to Fig 3). However, the pattern is still not what we would intuitively expect. For example: there are clearly some high activity for lower contrast levels in sessions 1 and 4.

ggplot(d1, aes(contrast_right, contrast_left, fill = msfr_all)) + 
  geom_tile() +
  facet_wrap(~session) +
  scale_fill_gradient(low="blue4", high="orange") +
  labs(x = "Right Contrast", y = "Left Contrast", fill = "MSFR", title = "MSFR by Combination of Contrasts across Sessions")

Fig 7: Heatmap of MSFR by Combination of Contrasts across Sessions

We now fit the same model as before, with two interacting fixed-effect factors for the left and right contrasts, and a random effect term for the sessions, but this the the target variable is the mean scaled firing rate (MSFR). Interestingly, none of the fixed effect terms came up as significant in the ANOVA decomposition (Tab 6). In addition, we also tried a reduced model without the random intercept, and even that did not result in any significant fixed effect terms. It appears that our original model with MFR was clearly more significant compared to the model with MSFR.

fit1 = lmer(msfr_all ~ (1 | session) + contrast_left*contrast_right, d1)
# summary(fit1)

kbl(anova(fit1), caption = "Tab 6: ANOVA table for fixed effects")

Tab 6: ANOVA table for fixed effects

	Sum Sq	Mean Sq	NumDF	DenDF	F value	Pr(>F)
contrast_left	0.0031617	0.0010539	3	1176.239	1.822946	0.1411432
contrast_right	0.0029853	0.0009951	3	1176.144	1.721253	0.1608171
contrast_left:contrast_right	0.0072519	0.0008058	9	1176.370	1.393751	0.1858372

fit2 = lmer(msfr_all ~ (1 | session) + contrast_left + contrast_right, d1)
summary(fit2)
anova(fit2)
anova(fit1, fit2)

Clustering Neurons

In all our previous analysis, we condensed the information of the activity of all neurons for each trial into one mean metric. This is definitely useful in reducing the dimension of our model, however in the process we may be losing a lot of information. For example, let us look at the distribution of the scaled firing rates of all neurons for one particular trial. From Fig 8 we observe that in this particular trial, there are a large number of neurons with very low scaled firing rates, some with moderate scaled firing rates, and even fewer with very high rates. We note that the categorization of scaled firing rates described was done so arbitrarily in reference to Fig 8. However, there is such a distribution of scaled firing rates for all trials, and therefore simply taking the mean scaled firing rate per trial may not be an appropriate statistic for our research.

hist(s[[1]]$sfr_neurons[[44]], col = "steelblue", xlab = " Scaled Firing Rate", main = "Distribution of Scaled Firing Rate of all neurons: Session 1 Trial 44")

Fig 8: Sample Distribution of per Neuron Scaled Firing Rate

To address this issue, we will have to use more granular data as a measure for the neural activity of the visual cortex. However, we can not use the data of every individual neuron, since the number of neurons vary between sessions, and the neurons are not uniquely identified; not to mention it would significantly increase the dimension of our model. To strike a balance, what we propose is to aggregate the scaled firing rate of each neuron with respect to the trial type (based on combination of left and right contrast). This will give us 16 aggregated (mean) statistics of scaled firing rates for each neuron. Using this data, we can implement some clustering algorithm to identify each neuron of each session into some unique cluster. We can then use the mean scaled firing rates of each cluster of neurons as the target variables for our model.

To that end, after aggregating the scaled firing rates of each neuron based on trial type, we implemented a k-means clustering algorithm on the neurons. From the scree plot (Fig 9), the elbow was determined to be at k=3. As such, we applied the clustering model on our data and clustered the neurons into three groups.

d2.1 = d2 %>% select(3:18)

wss = 0

for (i in 1:15) {
  km.out = kmeans(d2.1, centers = i, nstart = 20, iter.max = 50)
  wss[i] = km.out$tot.withinss
}

plot(1:15, wss, type = "b", 
     xlab = "Number of Clusters", 
     ylab = "Within groups sum of squares",
     main = "Scree plot")

Fig 9: Scree plot of K-Means Clustering

k = 3

km.neurons = kmeans(d2.1, centers = k, nstart = 20, iter.max = 50)
d2$cluster = km.neurons$cluster

table("Session" = d2$session, "Cluster" = d2$cluster)

Next, we computed the mean scaled firing rates (MSFR) of each cluster of neurons, and based on the overall MSFR values, the clusters were labelled as ‘low’, ‘med’, and ‘high’. A sample of the updated dataset is as follows in Tab 7.

for (ID in 1:5) {
  n.neurons = dim(session[[ID]]$spks[[1]])[1]
  
  c = d2$neuron[which(d2$session == ID & d2$cluster == 1)]
  s[[ID]]$msfr_low = colMeans(N[[ID]]$M[c,])
  
  c = d2$neuron[which(d2$session == ID & d2$cluster == 3)]
  s[[ID]]$msfr_med = colMeans(N[[ID]]$M[c,])
  
  c = d2$neuron[which(d2$session == ID & d2$cluster == 2)]
  s[[ID]]$msfr_high = colMeans(N[[ID]]$M[c,])
}

d1$msfr_low = c(s[[1]]$msfr_low, s[[2]]$msfr_low, s[[3]]$msfr_low, s[[4]]$msfr_low, s[[5]]$msfr_low)
d1$msfr_med = c(s[[1]]$msfr_med, s[[2]]$msfr_med, s[[3]]$msfr_med, s[[4]]$msfr_med, s[[5]]$msfr_med)
d1$msfr_high = c(s[[1]]$msfr_high, s[[2]]$msfr_high, s[[3]]$msfr_high, s[[4]]$msfr_high, s[[5]]$msfr_high)

kbl(head(d1), caption = "Tab 7: Head of Updated Dataset")

Tab 7: Head of Updated Dataset

session	date	name	contrast_left	contrast_right	trial_type	feedback_type	mfr_all	msfr_all	msfr_low	msfr_med	msfr_high
1	2016-12-14	Cori	1	0	L1-R0	1	6.193820	0.1303962	0.2182018	0.5789474	0.0535161
1	2016-12-14	Cori	0	0.5	L0-R0.5	1	4.087079	0.0961666	0.1776961	0.5401070	0.0222442
1	2016-12-14	Cori	1	0.5	L1-R0.5	1	5.870787	0.0939326	0.1766667	0.4400000	0.0285714
1	2016-12-14	Cori	0	0	L0-R0	1	4.087079	0.1021770	0.1666667	0.5965909	0.0304622
1	2016-12-14	Cori	0.5	1	L0.5-R1	0	5.533708	0.1006129	0.1723485	0.5082645	0.0339954
1	2016-12-14	Cori	0	0	L0-R0	0	3.553371	0.0789638	0.1550926	0.4191919	0.0168067

We then fit three separate models for the MSFR of each cluster with two interacting fixed effect terms for the left and right contrasts and a random effect term for the sessions. However, none of the three models showed any significant fixed effect terms in their ANOVA decomposition. This could be because our original model with MFR is just a more appropriate model, or it could also be because once the time-dependent trend was removed, there is no longer a significant association the visual stimuli and neural activity. The latter seems counter-intuitive, and further research is necessary to confidently determine the cause of what we observed.

fit3.1 = lmer(msfr_low ~ (1 | session) + contrast_left*contrast_right, d1)
fit3.2 = lmer(msfr_low ~ (1 | session) + contrast_left + contrast_right, d1)

anova(fit3.1, fit3.2)
anova(fit3.1)

fit4.1 = lmer(msfr_med ~ contrast_left*contrast_right + (1 | session), d1)
fit4.2 = lmer(msfr_med ~ contrast_left + contrast_right + (1 | session), d1)

anova(fit4.1, fit4.2)
anova(fit4.1)

fit5.1 = lmer(msfr_high ~ contrast_left*contrast_right + (1 | session), d1)
fit5.2 = lmer(msfr_high ~ contrast_left + contrast_right + (1 | session), d1)

anova(fit5.1, fit5.2)
anova(fit5.1)

Predictive Modeling

In this part of the project, our goal is to build a well-performing predictive model that accurately predicts the outcome (feedback type) of each trial. To do so, we will first start by visualizing the distribution of outcomes in our dataset. From Fig 10, we observe that the mice have a greater chance of success when the difference between the left and right contrast is greater. In particular, we see a greater probability of success for more contrast on the right side.

d1_prob <- d1 %>%
  group_by(session, contrast_right, contrast_left) %>%
  summarize(prob = mean(feedback_type == 1))

ggplot(d1_prob, aes(x = contrast_right, y = contrast_left, fill = prob)) +
  geom_tile() +
  scale_fill_gradient(low = "blue4", high = "orange") +
  labs(x = "Right Contrast", y = "Left Contrast", fill = "Probability of Success", title = "Probability of Success by Combination of Contrasts")

Fig 10: Probability of Success by Combination of Contrasts

In Fig 11, we observe a roughly similar pattern with some variation. This raises the question of whether to include a term for sessions in our predictive models.

ggplot(d1_prob, aes(x = contrast_right, y = contrast_left, fill = prob)) +
  geom_tile() +
  scale_fill_gradient(low = "blue4", high = "orange") +
  facet_wrap(~session) +
  labs(x = "Right Contrast", y = "Left Contrast", fill = "Probability of Success", title = "Probability of Success by Combination of Contrasts across Sessions")

Fig 11: Probability of Success by Combination of Contrasts across Sessions

In order to perform the predictive analysis, we will use the logistic regression model since we are dealing with a binary classification. In particular, our full model is:

\[ {logit}(\mathbb{E}[y_i|X_i])= \beta_0 + \sum_{i=2}^{5} \beta_{1i}X_{1i} + \sum_{j=2}^{4} \beta_{2j}X_{2j} + \sum_{k=2}^{4} \beta_{3k}X_{3k} + \sum_{j=2}^{4} \sum_{k=2}^{4}\beta_{4jk}(X_{2j}X_{3k}) \]

where \(X_{1i}\) is the indicator variable for the \(i\)th session, \(X_{2j}\) is the indicator variable for the \(j\)th level of left contrast, \(X_{3k}\) is the indicator variable for the \(k\)th level of right contrast, and \((X_{2j}X_{3k})\) is the indicator variable for the interaction of left and right contrast. All indices start from 2 since the base case is considered to be the reference class.

d3.train = d1[-(1:100), ]
d3.test = d1[(1:100), ]

Before moving on to model fitting, we first set aside the data of the first 100 trials of Session 1 as our test data, and we will build our model using the remaining dataset. We fit two logistic regression models, one as the full model defined above, and the other as a reduced model without the session variable. We then conduct a likelihood ratio test for \(H_0: \beta_{1i} = 0\) at 5% level of significance. From Tab 8, we can see the p-value is 0. Therefore, we reject \(H_0\) at 5% level of significance, and conclude that the session variable can not be dropped from our model.

fit6.1 = glm(feedback_type ~ session + contrast_left*contrast_right + mfr_all, data = d3.train, family="binomial")

summary(fit6.1)

fit6.2 = glm(feedback_type ~ contrast_left*contrast_right + mfr_all, data = d3.train, family="binomial")

summary(fit6.2)

kbl(lrtest(fit6.1, fit6.2), caption = "Tab 8: Likelihood Ratio Test") %>%
  kable_styling(full_width = F)

To further validate our model, we plot the ROC curves for both the models, and then determine the respective threshold value for each model to use in our final prediction on the test data. From Fig 11 and Fig 12, we observe that the ROC curves of both models are comparable, with the full model appearing to have a slightly larger AUC. Furthermore, the threshold value for the full model is 0.624, and that for the reduced model is 0.673. If the predicted/fitted value of our model is greater than the threshold, the observation will be assigned as a success, and otherwise as a failure.

#Full Model
predicted_prob = predict(fit6.1, d3.train, type = "response")

result = prediction(predicted_prob, d3.train$feedback_type) %>%
  performance(measure = "tpr", x.measure = "fpr")

d3.train.plot = data.frame(x = result@x.values[[1]],
                       y = result@y.values[[1]], 
                       p = result@alpha.values[[1]])

roc = ggplot(data = d3.train.plot) +
  geom_path(aes(x = x, y = y)) + 
  xlab(result@x.name) +
  ylab(result@y.name) +
  theme_bw()

dist_vec = d3.train.plot$x^2 + (1 - d3.train.plot$y)^2
opt_pos = which.min(dist_vec)

roc1 = roc + 
  geom_point(data = d3.train.plot[opt_pos, ], 
             aes(x = x, y = y), col = "red") +
  annotate("text", 
           x = d3.train.plot[opt_pos, ]$x + 0.1,
           y = d3.train.plot[opt_pos, ]$y,
           label = paste("p =", round(d3.train.plot[opt_pos, ]$p, 3)))

roc1 + labs(title = "ROC Curve of Full model")

Fig 11: ROC Curve of Full Model

#Reduced Model
predicted_prob = predict(fit6.2, d3.train, type = "response")

result = prediction(predicted_prob, d3.train$feedback_type) %>%
  performance(measure = "tpr", x.measure = "fpr")

d3.train.plot = data.frame(x = result@x.values[[1]],
                       y = result@y.values[[1]], 
                       p = result@alpha.values[[1]])

roc = ggplot(data = d3.train.plot) +
  geom_path(aes(x = x, y = y)) + 
  xlab(result@x.name) +
  ylab(result@y.name) +
  theme_bw()

dist_vec = d3.train.plot$x^2 + (1 - d3.train.plot$y)^2
opt_pos = which.min(dist_vec)

roc2 = roc + 
  geom_point(data = d3.train.plot[opt_pos, ], 
             aes(x = x, y = y), col = "red") +
  annotate("text", 
           x = d3.train.plot[opt_pos, ]$x + 0.1,
           y = d3.train.plot[opt_pos, ]$y,
           label = paste("p =", round(d3.train.plot[opt_pos, ]$p, 3)))

roc2 + labs(title = "ROC Curve of Reduced model")

Fig 12: ROC Curve of Reduced Model

With the threshold values determined, we now fit our model to the test data to get some predictions. Comparing the predicted values with the actual values, we obtain the following confusion matrix for the two model.

predicted_prob = predict(fit6.1, d3.test, type = "response")
prediction = as.integer(predicted_prob > 0.624)
confusion_mat1 = addmargins(table(d3.test$feedback_type, prediction))
names(dimnames(confusion_mat1)) <- c("True Feedback", "Predicted Feedback")
colnames(confusion_mat1) <- c("Fail", "Success", "Total")
rownames(confusion_mat1) <- c("Fail", "Success", "Total")

kbl(confusion_mat1, caption = "Tab 9: Confusion Matrix of Full Model") %>%
    add_header_above(c(" ", "Predicted Feedback" = 2, " "))

predicted_prob = predict(fit6.2, d3.test, type = "response")
prediction = as.integer(predicted_prob > 0.73)
confusion_mat2 = addmargins(table(d3.test$feedback_type, prediction))
names(dimnames(confusion_mat2)) <- c("True Feedback", "Predicted Feedback")
colnames(confusion_mat2) <- c("Fail", "Success", "Total")
rownames(confusion_mat2) <- c("Fail", "Success", "Total")

kbl(confusion_mat2, caption = "Tab 10: Confusion Matrix of reduced Model") %>%
    add_header_above(c(" ", "Predicted Feedback" = 2, " "))

Performance of full model:

• Sensitivity: \(\frac{63}{74} = 0.8514\)

• Specificity: \(\frac{9}{26} = 0.3462\)

Performance of reduced model:

• Sensitivity: \(\frac{47}{74} = 0.6351\)

• Specificity: \(\frac{11}{26} = 0.4231\)

Comparing the predictive performance of the two models, we determine that the full model is the better performer, since the gain in specificity in the full model is small compared to the loss in sensitivity.

Discussion

In this project we explored how neural activity in the visual cortex of mice is affected by visual stimuli presented to them, and found that the mean firing rate has an interactive effect with the levels of stimuli on the left and right side of the mice. However, the same result was not obtained when the target variable was the mean scaled firing rate of all neurons or the scaled mean firing rate of different clusters of neurons. Further analysis is required to determine an alternate statistic which is more robust than the mean firing rate. Furthermore, we built a predictive model for predicting the outcome each trial, and came up with a model with around 85% sensitivity and 35% specificity. For both our inferential and predictive analyses, we used a dataset that is a small subset of the original dataset of the experiment by Steinmetz, et al. (2019), and we believe further analysis using more features like the identity of unique neurons can improve the quality of out findings. For future research, we suggest looking at a dataset with more sessions across more mice recorded for a longer time window starting before stimulus onset. We also suggest incorporating a reliable way to uniquely identify the neurons to better understand the modulation of neural activity as a result of visual stimuli.