Optimal Transport Theory to Extract Spiking Motifs

1. # %%
2. %load_ext autoreload
3. %autoreload 2
5. # %%
6. import torch
7. import matplotlib.pyplot as plt
8. import numpy as np
9. from dataset_generation import generate_dataset, sm_generative_model
10. from wassa_plots import plot_SM, plot_colored_raster
12. # %% [markdown]
13. # # Generative model for raster plots and spiking motifs
15. # %% [markdown]
16. # List of parameters for this generative model
18. # %%
19. dataset_parameters = {
20. 'seed' : 81,
22. 'N_pre' : 100, # number of neurons
23. 'N_timesteps' : 1000, # number of timesteps for the raster plot (in ms)
24. 'N_samples' : 60, # total number of samples in the dataset
26. 'N_delays' : 51, # number of timesteps in spiking motifs, must be a odd number for convolutions
27. 'N_SMs' : 3, # number of structured spiking motifs
28. 'avg_fr' : 20, # average firing rate of the neurons (in Hz)
29. 'std_fr' : .1, # standard deviation for the firing rates of the different neurons
30. 'overlapping_sms' : True, # possibility to have overlapping sequences
32. 'temporal_jitter' : .1, # temporal jitter for the spike generation in motifs
33. 'dropout_proba' : 0, # probabilistic participations of the different neurons to the spiking motif
34. 'additive_noise' : .1, # percentage of background noise/spontaneous activity
35. 'warping_coef' : 1, # coefficient for time warping
37. }
38. dataset_parameters.update({'N_involved':dataset_parameters['N_pre']*torch.ones(dataset_parameters['N_SMs'])}) # number of neurons involved in the spiking motif
39. dataset_parameters.update({'frs' : torch.normal(dataset_parameters['avg_fr'], dataset_parameters['std_fr'], size=(dataset_parameters['N_pre'],)).abs()})
40. dataset_parameters.update({'freq_sms' : 4*torch.ones(dataset_parameters['N_SMs'])}) # frequency of apparition of the different spiking motifs (in Hz)
41. dataset_parameters.update({'coefficient_variation' : 1*torch.ones(dataset_parameters['N_pre'])}) # coefficient of variation of the FRs of each neuron (within VS. without a motif)
43. # %%
44. sm = sm_generative_model(dataset_parameters)
46. # %% [markdown]
47. # The generative model will draw raster plots generated from both a spontaneous (random) activity and from spiking motifs. There are *parameters.N_SM* spiking motifs represented in the next plot. There are given as a spiking probability for each presynatpic neuron and each timestep (or delay).
49. # %%
50. plot_SM(sm.SMs);
52. # %% [markdown]
53. # The raster plot is produced with a sequence of apparition of the different motifs and random activity when no motif occurs.
55. # %%
56. number_of_trials = 100
57. input_rp, output_rp = sm.draw_input(nb_trials=number_of_trials)
59. # %%
60. random_ind = torch.randint(dataset_parameters['N_samples'],[1])
61. plot_colored_raster(input_rp[random_ind].squeeze(0), output_rp[random_ind].squeeze(0), dataset_parameters['N_delays'], dataset_parameters['warping_coef']);
63. # %% [markdown]
64. # The not structured spiking activity of the different neurons is given by their FRs (*parameters.frs*). Then, it is given by the occurence of a motif (colored parts of the raster plot) and the structure of the motif. The occurence of a motif is given by a random trial defining its timing and its type (or color).
65. #
66. # The raster plot is obtained with a Bernoulli trial on the probability distribution representing the non structured (flat distribution) or the structured (probability distributions of the different motifs). When motifs overlap the formula used to merge the motifs and obtain the expected raster plot is the following:
67. # $$Bernoulli(p) + Bernoulli(q) = Bernoulli(p+q-p*q)$$
68. #
69. # It leads to Poisson distribution with an expected number of spike during the trial and their relative temporal dependencies that can vary with the number of motifs, their frequency of apparition (*parameters.freq_sm*) or the number of neurons involved in the motif (*parameters.N_involved*)
71. # %%
72. def moving_average(x, w):
73. return np.convolve(x, np.ones(w), 'valid') / w
75. fig, ax = plt.subplots(figsize = (12, 1.61803))
76. ax.plot(moving_average(input_rp[0].sum(axis=(0)),1))
77. ax.set_xlim(0, dataset_parameters['N_timesteps'])
78. ax.set_xlabel('Time (a. u.)')
79. ax.set_ylabel('Spike number')
80. ax.set_title('Moving average of the number of spike per timebin');
82. # %% [markdown]
83. # The instantaneous firing rate sumed over the different trials shows that the distribution is Poisson. And we can check the effective number of spikes per neuron compared to the expected firing rate. Firing rate of the different neurons can be easily changed, default is chosen by following a normal distribution around an average firing rate.
85. # %%
86. fig, ax = plt.subplots()
87. ax.bar(np.arange(dataset_parameters['N_pre']),1e3*input_rp.mean(axis=(0,2)))
88. ax.plot(dataset_parameters['frs'], 'ro', label='expected firing rate')
89. ax.set_xlabel('@ Neuron')
90. ax.set_ylabel('Firing rate (in Hz)')
91. ax.legend();
92. ax.set_title('Average number of spikes per neuron');
94. # %% [markdown]
95. # We can also check the number of occurence of the motifs and compare it to the expected number of occurence.
97. # %%
98. fig, ax = plt.subplots()
99. ax.bar(np.arange(dataset_parameters['N_SMs']),output_rp.sum(axis=(0,2)))
100. ax.plot(dataset_parameters['freq_sms']*number_of_trials*dataset_parameters['N_timesteps']*1e-3, 'ro', label='expected number')
101. ax.set_xlabel('Spiking motif #')
102. ax.set_ylabel('Number of occurence')
103. ax.set_title('Occurence of the different motifs');
104. ax.legend();
106. # %% [markdown]
107. # ## No overlap of the motifs
108. #
109. # If you set the following parameter to False, then the overlap of the motifs is prevented.
111. # %%
112. sm.opt['overlapping_sms'] = False
113. input_rp, output_rp = sm.draw_input(nb_trials=100)
115. # %%
116. random_ind = torch.randint(dataset_parameters['N_samples'],[1])
117. plot_colored_raster(input_rp[random_ind].squeeze(0), output_rp[random_ind].squeeze(0), dataset_parameters['N_delays'], dataset_parameters['warping_coef']);
119. # %%
120. fig, ax = plt.subplots()
121. ax.bar(np.arange(dataset_parameters['N_SMs']),output_rp.sum(axis=(0,2)))
122. ax.plot(dataset_parameters['freq_sms']*dataset_parameters['N_timesteps']*number_of_trials*1e-3, 'ro', label='expected number')
123. ax.set_xlabel('Spiking motif #')
124. ax.set_ylabel('Number of occurence')
125. ax.set_title('Occurence of the different motifs');
126. ax.legend();
128. # %% [markdown]
129. # ## Save a dataset
130. #
131. # You can also easily save a dataset with a defined number of trials (*total_number_of_trials*). It will split the dataset into a training (80% of the trials) and a testing (20% of the trials) of the trials sets.
133. # %%
134. total_number_of_trials = 1000
135. sm, dataset_input, dataset_output = generate_dataset(dataset_parameters, record_path = './synthetic_data/')
137. # %%
138. random_ind = torch.randint(dataset_parameters['N_samples'],[1])
139. plot_colored_raster(input_rp[random_ind].squeeze(0), output_rp[random_ind].squeeze(0), dataset_parameters['N_delays'], dataset_parameters['warping_coef']);
141. # %% [markdown]
142. # ## Alteration of the motifs
143. #
144. # There are different types of noise that can be added to the motifs:
145. # - temporal jitter: the timing of a spiking can be shifted in time. It is modeled by convolving the precise spiking motifs with a gaussian kernel (its standard deviation (*parameters.temporal_jitter*) gives the amount of jitter).
146. # - stochatistic activation of the neurons: neurons can also not participate in the generation of a spike at a specific spatio temporal location. This participation is obtained with a random trial on which neurons will be activated with a probability of activation for each neuron (*parameters.dropout_proba*)
147. # - additive noise: percentage of background activity given by (*parameters.additive_noise*).
148. # - time warping: the kernel is streched in time with a specific coefficient (parameters.warping_coef)
150. # %% [markdown]
151. # Because we maintain a constant firing rates for the different trials and the way the generative model is defined, these types of alteration will occur only at the motif level and then, a background activity of 100% leads to the same type of unstructured activity as a dropout probability of 1 for each of the neurons or an infinite temporal jitter.
153. # %%

TUTO_generation_synthetic_datasets.ipynb at commit bb3d841, no license · at the source