Normalization per trial?

[davidgwyrick]

Hello @DishengTang again! So I was revisiting this code for another project I'm working on and noticed that you calculate the normalization for each neuron pair per trial. In calculate_ccg_pair_single_trial() you divide by the triangle function and the firing_rates for that trial, but in the original formulation by Xiaoxuan Jia and in your own paper methods, the CCG between a pair of neurons is defined using the average firing rate. I don't think these are equivalent given how you're adding the CCG values across trials.

You could calculate the normalization for each pair in calculate_mean_ccg_corrected() after all the trials have been processed like this:

	# Compute firing rates for normalization
	firing_rates = np.mean(np.sum(spikeTrain,axis=2),axis=1)/ (T * self.time_bin)

	#Normalize by firing rates and triangle function
	normalization = (T - np.arange(self.window + 1)) * self.time_bin * np.sqrt(np.outer(firing_rates, firing_rates))[:, :, None]
	ccgs = np.divide(ccgs,normalization, where=normalization!=0,out=np.full(ccgs.shape,np.nan)) 
	ccg_jittered = np.divide(ccg_jittered,normalization, where=normalization!=0,out=np.full(ccg_jittered.shape,np.nan))

So in calculate_ccg_pair_single_trial() just comment out the per trial normalization.

	px, py = padded_st1[ind_A, :], padded_st2[ind_B, :]
	shifted = as_strided(px[self.window:], shape=(self.window + 1, T + self.window),
									strides=(-px.strides[0], px.strides[0]))
	return (shifted @ py)  # / ((T - np.arange(self.window + 1)) * self.time_bin * np.sqrt(firing_rates[ind_A] * firing_rates[ind_B]))

Let me know what you think. Maybe @myhannahchoi & @jiaxx can weigh in too.

[davidgwyrick]

Pretty close, but still different!

[DishengTang]

Hi @davidgwyrick thanks for catching this and for the clear suggestion!

You’re right: both Xiaoxuan’s original formulation and our paper Methods define the normalization using each neuron’s mean firing rate across trials, while the repo implementation was normalizing per trial and then averaging across trials. These are not strictly equivalent because per-trial normalization effectively re-weights trials, which can matter particularly for low-spiking neurons.

I’ve updated the code to match the Methods: the CCG normalization now uses mean firing rates over trials, and I also added optional random seed support so the jittered surrogates (and resulting sparse networks) are reproducible across runs. I ran a quick check on one example session/stimulus for the two variants: individual edge sets differ somewhat (expected for thresholded sparse networks), but the downstream summaries I inspected (motif profile and module partition metrics) were qualitatively consistent between them.

condition	weighted purity	weighted coverage	rand index
per_trial	0.493590	0.712384	0.156179
mean_over_trials	0.506711	0.698071	0.169860

Practical tip: when running this pipeline, it’s worth doing a small stability/power test for num_jitter—peak/interval detection can be noisy when spike counts are low, and edge sets can be sensitive to the number of jitter samples.

[davidgwyrick]

ya, it didn't seem like a big issue. Thanks for looking into it though!

how many jitter samples do you run? For larger datasets (I'm trying to run the method on >500 neurons), this is important because of computational time. I tried to convert things to jax/gpu but having some trouble with memory. It would be great to have a meeting with you and @myhannahchoi perhaps to discuss the results I have thus far to get your feedback. It involves mice on psilocybin.

[DishengTang]

We used 100 jitter samples, but later found that 10 does not significantly change the qualitative results if you have multiple sessions. I think a pilot run on maybe 100 neurons could provide some evidence and it won't take too much time especially if you can use multiple machines. Yeah I'd be happy to meet, my gmail is dishengtang3@gmail.com and maybe we can arrange an online meeting for this.

[davidgwyrick]

What do you mean by multiple sessions?
100 samples! woah ok. Evidence would be that the edge set is more stable? i.e. you find the same connected neurons between different jitter corrected instantiations? I have noticed that sometimes 10 jitter samples produces slightly different adjacency matrices...

[DishengTang]

Yes 10 samples can lead to some differences in the adjacency matrix, usually those are the connections with a significance level around the threshold (4-fold for example), the most significant connections would be usually more stable (such as Z-score of 6,7) even with a few jitter samples. This can also depend on number of trials/spikes, so in general it's advised to do a quick power test on your specific data and significance level.
I mean the downstream network analysis may not depend on every connection, so some fundamental properties could remain similar even with slight difference across jitter samples, especially if your findings are from multiple networks from different sessions (mice or something), then the average results across networks/sessions would be more reliable and consistent.