A biologically plausible decision-making model based on interacting neural populations.

1. # %% [markdown]
2. # # Notebook of a biophysically plausible decision-making model based on interacting cortical columns
3. # by Emre Baspinar, CNRS, NeuroPSI, Laboratory of Computational Neuroscience, Paris-Saclay
4. #
5. # This notebook contains the implementation of our biophysically realistic AdEx mean-field model proposed for a reward-driven consequential decision-making task. The implementation is in Python 3.8.5. The code writing and the preparation of the Jupyter notebook was done by Emre Baspinar. The model design was contributed by Emre Baspinar, Gloria Cecchini (U. of Barcelona), Rubén Moreno Bote (Pompeu Fabra University, Barcelona), Ignasi Cos (U. of Barcelona) and Alain Destexhe (CNRS, NeuroPSI, Paris-Saclay). The related manuscript [1] is a joint work of Emre Baspinar, Gloria Cecchini, Michael DePass (U. of Barcelona), Marta Andujar (U. of Rome), Pierpaolo Pani (U. of Rome), Stefano Ferraina (U. of Rome), Rubén Moreno Bote, Ignasi Cos and Alain Destexhe. See [1] for more information.
6. #
7. # This notebook is licensed by Creative Commons Attribution 4.0 International Public License (CC BY 4.0). Please cite [1] and [7] if you use this notebook in your work.
8. #
9. # Please feel free to contact for questions, comments and feedback.
10. #
11. # Contact: [email hidden]
12. #
13. #
14. #
15. # This notebook provides the implementation of the decision-making AdEx mean-field model which was presented in [1]. It provides the implementation, an example simulation of the behavioral experiments related to the decision-making paradigm considered in [1, 2], and three case studies. The results of the example simulation are saved automatically in the notebook folder in case of that the user wants to analyse the results. The default values for the simmulations are provided and they can be easily changed by the user. For more description regarding the model and the technical details, please see *ModelSummary.pdf* and *README.pdf*.
17. # %% [markdown]
18. # # Glossary
19. #
20. # <u>MainNotebook.ipnyb</u>: This is the main notebook to run the simulations and the case studies. It contains explanations related to the model and to the simulations. This file uses the .py files given below.
21. #
22. # <u>cell_library.py</u>: It contains the parameters of the biophysical cell properties of the neurons. Do not change unless you add new cell types.
23. #
24. # <u>DiffOperator.py</u>: It contains the stochastic differential equations of the AdEx mean-field equations corresponding to both cortical columns.
25. #
26. # <u>NeuronConnectivity.py</u>: It contains the functions which we use to load the transfer functions of Regular Spiking (RS) and Fast Spikin (FS) cells by using the method explained in~\cite{zerlaut2018modeling}. The transfer functions and their parameters are based on a fitting to experimental data, therefore the parameters should be kept fixed. Do not change this file.
27. #
28. # <u>SDEIntegrator.py</u>: This is the Euler-Maruyama integrator. It integrates the SDEs found in DiffOperator.py in time.
29. #
30. # <u>syn_and_connec_library.py</u>: It contains the connectivity and synaptic properties of the neurons.
31. #
32. # <u>theoretical\_tools.py</u>: It contains implementation of some analytical functions appearing in the AdEx mean-field equations.
33. #
34. # In addition to these .py files, there are two .npy files in *data* folder: FS-cell_CONFIG1_fit.npy and RS-cell_CONFIG1_fit.npy. They contain the fitted parameters of RS and FS cell transfer functions to the experimentally obtained RS and FS cell response vlaues. Do not change the location of these files. Finally, *showcaseData* folder contains the data which MainNotebook.ipnyb uses for the case studies.
35. #
36. # This notebook uses Matplotlib, Nump, Random, Os and Scipy.io libraries.
38. # %% [markdown]
39. # ## Initialization
41. # %%
42. import matplotlib.pyplot as plt
43. import numpy as np
44. import random
45. import os
46. import scipy.io
47. # import operator
48. from SDEIntegrator import RegulatoryPsi, TimeStepping
49. from NeuronConnectivity import ReformatSynParameters, LoadTransferFunctions
51. # %% [markdown]
52. # ## Introduction
53. # Our model contains two pools where each pool is composed of an excitatory and an inhibitory population. Each population is modeled as a mean-field adaptive exponential integrate-and-fire (AdEx) model, which was derived in [3] from the AdEx network model [4] by following the *master equation formalism* [5, 6], by including adaptive dynamics as well. This model was proposed for only one excitatory and one inhibitory population, that is, for only one pool. We extend this model to a double pool architecture, where now also cross connections and correlations between the cross-populations must be considered. The challenge, besides the increased dimension of the mean field model, is to implement the cross-connections between the pools, since it requires the derivatives of cross-population covariances. Yet, it allows to take into account biophysical mechanisms which involve in the decision making process and which were not taken into account previously.
55. # %% [markdown]
56. # ## Experimental task
57. # Experimental task was applied on healthy human subjects. Each session consists of a number of episodes. Each episode is composed of a fixed number of trials. The number of trials are determined by the horizon number. If it is 0, then we have 1 trial per episode, and if it is 1, we have 2 trials per episode. In each trial, the subject is asked to choose between two stimuli, of which each one is a partially filled vertical bar. The amount of the filling is different for each bar. The total height of the bars are identical. Depending on the choice the subject makes and the implemented strategy for every episode, the subject has a reward or loss for the next trial. The reward is a fixed amount of gain added to the filled parts of both stimuli, and the loss is the same amount of gain subtracted from the filled parts of both stimuli. Each bar is located either on the right hand side or on the left hand side of the screen, and their positions might be randomly exchanged between two trials. The sum of the chosen stimuli with the reward or the loss obtained at the end of each trial of each episode is accepted as the total result of the episode. More precisely, if we denote the chosen stimulus plus the reward or loss at the end of the $i$th trial by $A_i$, then the total result of an episode with $N\in \mathbb{Z}$ trial is
58. # \begin{equation}
59. # \text{Total reward} = \sum\limits_{i=0}^N = A_i.
60. # \end{equation}
61. # The goal of the subject is to maximize the total result in each episode, that is, to be rewarded as many times as possible in each episode by making decisions in accordance with the implemented strategy. The behavioral performance of the subject is scored at the end of each episode as a scaled measure of the total reward between 0 and 1.
63. # %% [markdown]
64. # |  |
65. # |:--:|
66. # | <b>Fig. 1: Example episode with 2 trials. Implemented strategy is to choose small and large stimuli in the first and second trials, respectively. Subject makes the bad choice in the first trial and receives a loss. In the second she/he makes the good choice, thus receives a reward. The total result is the sum of the chosen stimuli (after their reward or loss), which are highlighted by the orange dashed rectangles. Corresponding performance results is the total result normalized by the sum of the stimuli corresponding to the maximum total result, i.e., by the double rewarded total result. The total result is not explicitly shown to the participant. </b> |
68. # %% [markdown]
69. # ## Double-column AdEx mean field model
70. # Departure point of this study is the AdEx mean-field model which was proposed in [3]. We extend to a bicolumnar setting the model equations and the simulation framework provided in [3], which was for only a single cortical column.
71. #
72. # The extended model equations are as follows:
73. # \begin{equation}
74. # \begin{split}
75. # T \partial_t v_\alpha & = (F_\alpha -v_\alpha) + \frac{1}{2} C_{\xi\eta}\partial_{\xi\eta}F_\alpha +\sigma \omega_\alpha\\
76. # T\partial_t C_{\alpha\beta} & =\delta_{\alpha\beta}A_{\alpha\alpha} + (F_\alpha-v_\alpha)(F_\beta-v_\beta) + C_{\beta\xi}\partial_\xi F_\alpha + C_\alpha\xi\partial_\xi F_\eta-2 C_{\alpha\beta}\\
77. # \partial_t W_\alpha & = -\frac{W_\alpha}{\tau_w} + (\delta_{\alpha e_A}+\delta_{\alpha e_B})(b v_\alpha + a(\mu_V(v_{e_A}, v_{e_B}, W_\alpha) - E_L)),\hspace{12cm}(1)
78. # \end{split}
79. # \end{equation}
80. # where
81. # \begin{equation}
82. # \begin{split}
83. # F_{e_A}=F_{e_A}(\tilde{v}_{e_A}, \tilde{v}_{i_A}, W_{e_A}), & \quad\quad
84. # F_{i_A}=F_{e_A}(\tilde{v}_{e_A}, \tilde{v}_{i_A}, W_{i_A}),\\
85. # F_{e_B}=F_{e_B}(\tilde{v}_{e_B}, \tilde{v}_{i_B}, W_{e_B}), &\quad\quad
86. # F_{i_B}=F_{e_A}(\tilde{v}_{e_B}, \tilde{v}_{i_B}, W_{i_B}),
87. # \end{split}
88. # \end{equation}
89. # are the population transfer functions with the corresponding population indicated as subindex and where $\omega_\alpha=w_\alpha(t)$ is a white Gaussian noise generated independently for each time instant and for each population. More precisely,
90. # $$
91. # \mathbb{E}[\omega_\alpha(t)]=0,\;\text{for all}\; t\geq 0,\quad\text{and}\quad \mathbb{E}[\omega_\alpha(t) \omega_\beta(t^\prime)]=\delta_{\alpha\beta}\delta_{tt^\prime}.
92. # $$
93. # Here $\sigma>0$ shows the noise intensity level. Function $A_{\alpha\beta}$ is defined as follows:
94. # $$
95. # A_{\alpha\beta}=\delta_{\alpha\beta}\frac{N_\alpha}{F_\alpha(1/T-F_\beta)},
96. # $$
97. # with $N_\alpha$ denoting the number of neurons in population $\alpha$. Here $T$ is the time scale parameter both for the firing rate and for the covariance variables appearing in (1). We refer to [Section 2.3.1, 3] for a detailed explanation of $\mu_V$ appearing in the last line of (1). The time dependency of the terms are tacitly understood from (3). We use the notation given by $\partial_\alpha = \frac{\partial}{\partial v_\alpha}$ and $\partial_{\alpha\beta} = \frac{\partial^2}{\partial {v_\alpha}\partial_{v_\beta}}$ for the partial derivatives. Here $\delta$ is the Dirac delta function, $E_L$ is a constant representing the reversal leakage potential, and finally $\tau_w$ is a time scale-like parameter for $W_\alpha$. The latter determines the evolution time scale of $W_\alpha$ approximately. We denote the total excitatory input firing rate for each transfer function by using the tilde notation. We express the probability of cross-connections with $p_c$. Cross-connections are always from excitatory populations, therefore, $p_c$ can be thought of as the ratio of exciatory neurons to the total number of neurons in the whole double pool network. Then we write the total excitatory input firing rates for the transfer functions as follows:
98. # \begin{equation}
99. # \begin{split}
100. # \tilde{v}_{e_A}(t) & = v_{e_A}(t)+ v_{\text{AI}}+ \lambda^A_{w_r}(v_A(t), v_B(t))+w^e_c \Big ( v_{e_B}(t) + v_{\text{AI}}+\lambda^B_{w_r}(v_A(t), v_B(t), t) \Big ) p_c N_{e_B},\\
101. # \tilde{v}_{i_A}(t) & = v_{i_A}(t)+ \lambda^A_{w_r}(v_A(t), v_B(t)) + w^i_c \Big ( v_{e_B}(t) + v_{\text{AI}}+\lambda^B_{w_r}(v_A(t), v_B(t), t) \Big ) p_c N_{e_B},\\
102. # \tilde{v}_{e_B}(t) & = v_{e_B}(t) + v_{\text{AI}}+\lambda^B_{w_r}(v_A(t), v_B(t))+ w^e_c \Big ( v_{e_A}(t) + v_{\text{AI}}+\lambda^A_{w_r}(v_A(t), v_B(t)) \Big ) p_c N_{e_A},\\
103. # \tilde{v}_{i_B}(t) & = v_{i_B}(t)+ \lambda^B_{w_r}(v_A(t), v_B(t)) + w^i_c \Big ( v_{e_A}t)+v_{\text{AI}}+\lambda^A_{w_r}(v_A(t), v_B(t), t) \Big ) p_c N_{e_B},\hspace{10cm} (2)
104. # \end{split}
105. # \end{equation}
106. # where $\lambda$ represents the firing rate of the regulatory pool, which is explained in the next section. In (3), all terms for which the time dependency is not explicit are constants.
108. # %% [markdown]
109. # | |
110. # |:--:|
111. # | <b>Fig. 2: Mean-field model. Each cortical column is represented as one pool of excitatory ($e_A$, $e_B$) and inhibitory ($i_A$, $i_B$) neuronal populations. Excitatory (blue) and inhibitory (red) in-pool connection weights are denoted by $w_+=w_-=1$, respectively. Cross-pool connections model intercolumnar long range connections in the prefrontal cortex, thus they are only excitatory. The cross-pool connections towards excitatory and inhibitory populations are weighted by $0<w_c^e=w_c^i<1$. Here $V_{\text{ext}}=v_{\text{AI}}$ represents the based-drive potential keeping the model in the asynchronous irregular (awake) state. Finally, $\lambda$ represents the regulatory mechanism introducing the bias to the input stimuli as explained in the next section.</b> |
113. # %% [markdown]
114. # ## Regulatory mechanism
115. #
116. # We interpret $v_A$ and $v_B$ as the stimuli with large quantity sample and the small quantity sample, respectively. In the experiment task conducted on human subjects, one of them appears on the left hand side of the screen and the other one appears on the right hand side. Those positions of the stimuli might be shuffled at each trial. The decision is made for choosing either $v_A$ or $v_B$, and it depends mostly on the regulatory pool. The regulatory pool weights both stimuli in accordance with the decision to be made. In this way, it promotes one of the stimulus, making more likely that the pool sensitive to the promoted stimulus has a higher firing rate compared to the firing rate of the other pool. This results in that the pool with the higher firing rate dominates the other pool and the decision is made for the promoted stimulus. We find the time instant where the decision is made as the moment where the difference between the firing rates of the excitatory populations of two pools exceeds a prefixed value.
117. #
118. # We model the regulatory pool evolving during the $i$th trial of the $E$th episode in terms of the following stochastic differential system:
119. # \begin{equation}
120. # \begin{cases}
121. # \tau_\psi \frac{d\psi^E_i(t)}{dt} = -4 \psi^E_i(t) \Big ( \psi^E_i(t)-1 \Big )\Big ( \psi^E_i(t)-1/2 \Big )+\frac{\sigma_r}{(c_0 t)^2}\zeta_i\\
122. # \psi_i^E(0)= \phi^{E-1}_i,\hspace{15cm}(2)
123. # \end{cases}
124. # \end{equation}
125. # where $\tau_\psi$ is the time scale parameter, $\psi$ is the firing rate of the regulatory pool evolving in time independently for each trial, $\zeta_i = \zeta_i(t)$ is white Gaussian noise whose intensity level is determined by a temporally scaled version of $\sigma>0$, and finally $c_0>0$ is a constant. Here $E$ and $i$ denote the episode number and the trial number within the corresponding episode, respectively. The time dependent denominator appearing in front of this noise term introduces a strong stochastic behavior to the regulatory pool at the beginning of its corresponding trial, and this allows us to take into account that the system learns the strategy by exploring it initially over the episodes. Moreover, $c_0$ determines how much the model has the flexibility to deviate from the captured strategy once it learns the strategy. This is important to model the deviations occurring due to either perceptual difficulty to distinguish two stimuli or the human factor which makes the subject sometimes not follow the strategy to explore what happens. The evolution process (2) is restarted at the beginning of each trial by setting the initial condition to the value determined by the reward $\phi^E_i$ (to be explained later) corresponding to the same trial but of the precedent episode, i.e., of the $(E-1)$th episode. In this way, the reward mechanism provides for each trial a feedback to the regulatory pool at the end of the $(E-1)$th episode. This feedback determines towards which decision the regulatory pool will be biased in each trial of the $(E-1)$th episode. The final time of the evolution given by (2) is the final time of each trial and it is the same for every trial.
126. #
127. #
128. # The regulatory pool activation is transmitted to Pool $A$ and Pool $B$ populations via $\lambda$ function appearing in (2). Its explicit form is given by
129. # \begin{equation}\label{eq:lambdaExpression}
130. # \begin{split}
131. # \lambda^A_{w_r}(v_A(t), v_B(t), t) = & \psi\,v_A(t) + w_r(1-\psi(t))\,v_B(t),\\
132. # \lambda^B_{w_r}(v_A(t), v_B(t), t) = & \psi\,v_B(t) + w_r(1-\psi(t))\,v_A(t),
133. # \end{split}
134. # \end{equation}
135. # where $w_r$ is a constant (which we choose to be 1 for now).
137. # %% [markdown]
138. # ## Parameters and initial conditions
139. # Please refer to [1, 3] for more details about the parameters. Here we set the parameters in a way that the model is always in asynchronous irregular (AI) state, which is the awake state. We focus on rather AI state here, but it is possible to set the model to up-and-down state (sleep/anesthesia state) by changing the parameters; see [1] for proper parameters corresponding to up-and-down state.
140. #
141. # <b>Do not forget to choose if SI units are considered or not! Note that the default parameters given below are without SI units!</b>
142. #
143. # <b>The default AdEx parameters are for the simulations in this notebook, which are related to the decision-making tasks considered in [1]. Please do not change them unless you would like to use the notebook for other simulations than the ones included here.</b>
145. # %%
146. # AdEx system parameters
148. siUnits=int(input("Biophysical parameters without --> 0 or with --> 1 SI units?"))
150. if siUnits==0:
152. ## without SI units (in the poster-paper, we use without SI units)
153. aRS, bRS, aFS, bFS, tauwRS, tauwFS = 4, 40, 0, 0, 5000, 1e9 # a, b and tauw parameters of adaptation mechanism for RS and FS cells
154. Ntot = 20000 # total number of neurons in the corresponding network
155. pc = 0.80 # ratio: Ne/Ntot, where Ne is the number of excitatory cells (RS cells)
156. Ne = Ntot * pc/2 # number of excitatory cells in one pool (RS cells)
157. Ni = Ntot * (1-pc)/2 # number of inhibitory cells in one pool (FS cells)
158. vAI = 5 # base potential keeping the population pairs in asynchronous irreguler (AI) state (awake state)
159. wce = 2.5e-4 # cross-connection weight from excitatory to excitatory population
160. wci = 2.5e-4 # cross-connection weight from excitatory to inhibitory population
161. sigma = 0.01 # intrinsic noise level (for the noise appearing in firing rate equations in (1))
162. El = -65 # reversal leakage potential
163. Qe, Qi = 1.5, 5 # excitatory and inhibitory quantal conductances
164. Te, Ti = 5, 5 # decay (refractory) time scale of excitatory and inhibitory synapses
165. Gl = 10 # leak conductance
166. Ee, Ei = 0, -80 # excitatory and inhibitory reversal potentials
167. tF, dt = 15, 0.05 # final time and time step of the time integration via Euler-Maruyama
168. T = 5 # time scale for the fast variables
169. tauPsi = T # time scale of the regulatory mechanism
170. sigma_r = 0.01 # Extrinsic noise level (noise in the regulatory mechanism)
171. c0 = 1.0 # c0 constant determining the decay rate of the extrinsic noise
173. elif siUnits==1:
175. ## with SI units
176. aRS, bRS, aFS, bFS, tauwRS, tauwFS = 4e-9, 40e-12, 0, 0, 5, 1e6 # a, b and tauw parameters of adaptation mechanism for RS and FS cells
177. Ntot = 20000 # total number of neurons in the corresponding network
178. pc = 0.80 # ratio: Ne/Ntot, where Ne is the number of excitatory cells (RS cells)
179. Ne = Ntot * pc/2 # number of excitatory cells in one pool (RS cells)
180. Ni = Ntot * (1-pc)/2 # number of inhibitory cells in one pool (FS cells)
181. vAI = 5 # base potential keeping the population pairs in asynchronous irreguler (AI) state (awake state)
182. wce = 2.5e-4 # cross-connection weight from excitatory to excitatory population
183. wci = 2.5e-4 # cross-connection weight from excitatory to inhibitory population
184. sigma = 0.01 # intrinsic noise level (for the noise appearing in firing rate equations in (1))
185. El = -65e-3 # reversal leakage potential
186. Qe, Qi = 1.5e-9, 5e-9 # excitatory and inhibitory quantal conductances
187. Te, Ti = 5e-3, 5e-3 # decay (refractory) time scale of excitatory and inhibitory synapses
188. Gl = 10e-9 # leak conductance
189. Ee, Ei = 0, -80e-3 # excitatory and inhibitory reversal potentials
190. tF, dt = 6, 0.05 # final time and time step of the time integration via Euler-Maruyama
191. T = 5e-3 # time scale for the fast variables
192. tauPsi = 10*T # time scale of the regulatory mechanism
193. sigma_r = 0.01 # Extrinsic noise level (noise in the regulatory mechanism)
194. c0 = 45 # c0 constant determining the decay rate of the extrinsic noise
197. # We gather all parameters together in the vector "params"
198. params = np.zeros(28)
199. params[0] = aRS
200. params[1] = bRS
201. params[2] = aFS
202. params[3] = bFS
203. params[4] = tauwRS
204. params[5] = tauwFS
205. params[6] = Ntot
206. params[7] = pc
207. params[8] = Ne
208. params[9] = Ni
209. params[10] = vAI
210. params[11] = wce
211. params[12] = wci
212. params[13] = sigma
213. params[14] = El
214. params[15] = Qe
215. params[16] = Qi
216. params[17] = Te
217. params[18] = Ti
218. params[19] = Gl
219. params[20] = Ee
220. params[21] = Ei
221. params[22] = tF
222. params[23] = dt
223. params[24] = T
224. params[25] = tauPsi
225. params[26] = sigma_r
226. params[27] = c0
228. # %% [markdown]
229. # We set below the initial conditions ad hoc to the decision-making tasks which we consider in [1].
230. #
231. # <b> These are the default initial conditions for the simulations in this notebook, which are related to the decision-making tasks considered in [1]. Do not change them unless you would like to use the notebook for other simulations than the ones included here.
233. # %%
234. # Initial conditions for state variables
235. v0 = 1. # v_eA : V[0]
236. v1 = 30. # v_iA : V[1]
237. v2 = 0.5 # C_eAeA : V[2]
238. v3 = 0.5 # C_eAiA = C_iAeA : V[3]
239. v4 = 0.5 # C_iAiA : V[4]
240. v5 = 1.e-10 # W_eA : V[5]
241. v6 = 0. # W_iA : V[6]
242. v7 = 1. # v_eB : V[7]
243. v8 = 30. # v_iB : V[8]
244. v9 = 0.5 # C_eBeB : V[9]
245. v10 = 0.5 # C_eBiB : V[10]
246. v11 = 0.5 # C_iBiB : V[11]
247. v12 = 1.e-10 # W_eB : V[12]
248. v13 = 0. # W_iB : V[13]
249. v14 = 0.05 # C_eAeB : V[14]
250. v15 = 0.05 # C_eAiB : V[15]
251. v16 = 0.05 # C_iAeB : V[16]
252. v17 = 0 # C_iAiB : V[17]
254. psi0 = 0.5 # initial condition for the regulatory mechanism at the beginning of the first trial in each episode
256. V0 = [v0, v1, v2, v3, v4, v5, v6, v7, v8,\
257. v9, v10, v11, v12, v13, v14, v15, v16, v17]
259. # %% [markdown]
260. # ## Single-trial simulation setup
261. # We first initialize the stimuli A and B, in a way that one of them is weaker in comparison to the other.
263. # %%
264. # Define the stimuli for one trial
265. from scipy.special import comb
267. def smoothstep(x, x_min=0, x_max=1, N=1):
268. x = np.clip((x - x_min) / (x_max - x_min), 0, 1)
270. result = 0
271. for n in range(0, N + 1):
272. result += comb(N + n, n) * comb(2 * N + 1, N - n) * (-x) ** n
274. result *= x ** (N + 1)
276. return result
278. t = np.linspace(0, tF, int(tF/dt)+1) # define the whole time interval
281. ampA = 1 # amplitude of stimulus A
282. ampB = 10 # amplitude of stimulus B
283. t0 = 2 # initial time instant of the stimuli
286. stimulusA = ampA*(smoothstep(t-t0)-smoothstep(t-t0-tF))
287. stimulusB = ampB*(smoothstep(t-t0)-smoothstep(t-t0-tF))
289. # %% [markdown]
290. # Here we load the transfer functions which were derived for both regularly spiking (RS, excitatory) and fast spiking (FS, inhibitory) neural populations via the semi-analytical approach given in [6]. Once they are loaded, we can run a one-trial simulation.
292. # %%
293. ## Build the simulation setup: load the transfer functions for both RS & FS populations
294. TF1, TF2 = LoadTransferFunctions('RS-cell', 'FS-cell', 'CONFIG1')
295. lambdaA, lambdaB, psi = RegulatoryPsi(psi0, stimulusA, stimulusB, params) # generate the regulated stimuli for the trial
297. # Time integration for one single trial
298. state = TimeStepping(V0, lambdaA, lambdaB, TF1, TF2, params)
300. # %% [markdown]
301. # Now we can plot the evolution of firing rates of excitatory populations in both A and B pools. Those plots are example of one trial simulations, in which the reward is not included yet. The population which wins the competition makes the decision in favor of the stimulus which it is selective to. We say that the population having the higher firing rate wins the competition when the difference betweenthe firing rates exceed a certain threshold. The threshold is arbitrary and fixed to a constant at the beginnning of the whole experiment session. In the next step we will consider multiple episodes, where each epsiode contains one trial. In those simulations, we will introduce also the reward mechanism.
303. # %%
304. # %matplotlib widget
305. # Plot Pool A (red) and B (blue) excitatory population firing rates
306. plt.figure(figsize=(6, 4))
307. plt.plot(t, state[:, 0], 'r') # plotting v_{eA} (red)
308. plt.plot(t, state[:, 7], 'b') # plotting v_{eB} (blue)
309. # plt.margins(x=0.001, y=-0.001)
311. # %% [markdown]
312. # Now it is time to introduce the reward mechanism. This will allow us to extend our simulation framework to multiple episodes, where each episode contains one (Horizon 0) or two (Horizon 1) trials, depending on which horizon we would like to study.
314. # %% [markdown]
315. # ## Reward mechanism
316. # Reward mechanism is the key module introducing an online learning to the model. It allows the model to capture the strategy maximizing the gain through each episode. Once the strategy is learned, the system makes the decisions in the way which maximizes the final gain, which is the total reward of one episode. This is the key mechanism endowing the model with working memory.
317. #
318. # Reward mechanism is updated through a discrete evolution where the temporal variable is the trial number. In other words, the reward function value remains constant during a trial and it is updated at the end of the trial. This updated value at the end of the trial is fed to the same trial corresponding to the episode coming after as the initial condition of the regulatory function $\psi$; see Eq. (2). We use the notation $M^E_i$ to denote the mean value of the stimuli corresponding to $i$the trial of the $E$th episode. We write $N$ with $N-1$ expressing the number of the last trial in one episode to refer to the total reward obtained at the end of the episode, although $N$th trial does not exist. Then we write the evolution for the reward function $\phi$ as follows:
319. # \begin{equation}\label{eq:rewardMechanismEqn}
320. # \begin{cases}
321. # \phi^{E+1}_i = \phi^E_i + k ( M^E_{i+1} - M^E_i ) (2\psi_i^E(T)-1)(\phi^E_i -1 )^2(\phi_i^E)^2,\\
322. # \phi^0_i = C,\quad \text{for all}\quad i\in \{0, 1,\dots, N-1 \},\quad\quad\quad\quad\quad\quad\quad\quad\quad (3)
323. # \end{cases}
324. # \end{equation}
325. # with $C$ denoting a constant, which is fixed to $0.5$ for the first episode in our framework. This system initiates the reward mechanism for each trial series over the episodes separately, yet the trials are not independent due to the coupling effect arising from $( M^E_{i+1} - M^E_{i})$ factor.
326. #
327. # Reward mechanism is implemented as given below, in an integrated way to the simulation setup.
329. # %% [markdown]
330. # ## Model in action: simulations
331. #
332. # Each simulation session is composed of a certain number (nOfEpisodes) of episodes. Here we fix nOfEpisodes to 5 for simplicity although in our simulations we use higher values (e.g., 100). In Horizon 0, each episode has one single trial. In Horizon 1, each episode has two trials. See [1] for more details about the experiments.
334. # %%
335. # Parameters for episode initializations
336. nOfEpisodes = 5 # number of episodes in the whole session. In [1], it was selected as 100.
337. dSet = [0.01,0.05,0.1,0.15,0.2] # difficulty set of visual discriminization. Please do not change these values.
338. # difficulty is proportional to the difference between the amplitudes of the two stimuli
339. dMax = max(dSet) # max. value of the difficulty
340. nH =int(input("Enter the horizon number (0 or 1): "))# horizon number: n. of trials per episode = nH + 1
342. gain = 0.3 # reward gain. Please do not change this value.
344. if nH==0:
345. minStim = gain # min. value of the stimuli (for Horizon 0)
346. maxStim = minStim + 6 # max. value of the stimuli (for Horizon 0)
347. elif nH==1:
348. minStim = 2 * gain # min. value of the stimuli (for Horizon 1)
349. maxStim = minStim + 6 # max. value of the stimuli (for Horizon 1)
350. else:
351. print("Invalid horizon number! Choose either 0 or 1.")
353. reward = np.zeros((nH+1, nOfEpisodes)) # initialize reward list
354. decisionThreshold = 5 # the decision is made once the absolute value of the difference between excitatory population firing rates exceeds this value. Please do not change this value.
356. # Initialize the vectors which store the recorded results
357. performance = np.zeros(nOfEpisodes) # performance vector saving performance results of each episode
358. decisionTimeList = np.zeros(nOfEpisodes * (nH+1)) # vector saving the decision times of each trial
359. difficultyList = np.zeros(nOfEpisodes * (nH+1)) # vector saving the difficulty values used in each episodeNow it is time to introduce the reward mechanism. This will allow us to extend our simulation framework to multiple episodes with each episode containing one (Horizon 0) or two (Horizon 1) trials, depending on which horizon we would like to study.
361. # %% [markdown]
362. # We perform each simulation session for several times (nOfIterations) to be able to obtain statistics allowing us to quantify the behavioral results produced by the model framework and to compare them to the experimental data. We fit the model to the experimental data by tuning two parameters: learning speed ($k$ as given in Eq. (3)) and noise decay rate ($c_0$ as given in Eq. (2)). The user can chose one fixed value for either one of those two parameters and vary the other by assigning a vector of several value. <b>Please do not put comma between varied values.</b> Example--> Enter the k values: 0 0.025 0.050 0.075 1
363. #
364. # <b>DEFAULT VALUES:</b>
365. #
366. # I. For Horizon 0 with fixed $k$ and varied $c_0$: $k=0.3$ and $c_0$ is varied between 1.0 and 1.9.
367. #
368. # II. For Horizon 0 with fixed $c_0$ and varied $k$: $c_0=1.1$ and $k$ is varied between 0.05 and 0.5.
369. #
370. # III. For Horizon 1 with fixed $k$ and varied $c_0$: $k=0.05$ and $c_0$ is varied between 1.0 and 2.0.
371. #
372. # IV. For Horizon 1 with fixed $c_0$ and varied $k$: $c_0=2.0$ and $k$ is varied between 0 and 0.10.
373. #
374. # These simulations were used to produce the results found in Figure 8 and 14 in [1].
376. # %%
377. # Iteration parameters
378. nOfIterations = int(input("Enter the number of iterations per each experiment session: ")) # number of how many times we perform each simulation session. In [1], it was selected 10.
379. fixedPar = int(input("Enter 0 to fix k (learning speed), enter 1 to fix c0 (decay rate): ")) # determine which parameter will be fixed and the other will be varied
381. if fixedPar==0:
382. learningSpeed=float(input("Enter fixed k (learning speed) value: "))
383. # nc0Val=int(input("Number of c0 values: ")) # enter for how many values of c0 will be the experiment performed
384. c0List=list(map(float, input("Enter the c0 values: ").strip().split())) # Enter the values as, e.g., 1 2 5 7, without any comma in between
385. lengthParameterList = len(c0List)
386. nc0Val = lengthParameterList
387. elif fixedPar==1:
388. c0=float(input("Enter fixed c0 (decay rate) value: "))
389. # nkVal=int(input("Number of k values: ")) # enter for how many values of k will be the experiment performed
390. learningSpeedList=list(map(float, input("Enter the k values: ").strip().split())) # Enter the values as, e.g., 1 2 5 7, without any comma in between
391. lengthParameterList = len(learningSpeedList)
392. nkVal = lengthParameterList
395. # %% [markdown]
396. # Now we can run the simulation for the values entered in the above code box. <b>Results will be saved automatically to the notebook folder with the parameters noted in the title.</b> These results can be used for the statistical analysis of the simulation results. <b>You should comment out the three np.save lines at the end of the code block to turn this option off.</b>
398. # %%
399. # Run the simulations by varying the parameter chosen in the previous command box
400. for parameterNo in range(lengthParameterList):
401. iteration = 0
403. if fixedPar==0:
404. c0 = c0List[parameterNo] # to let the regulator have time to evolve, a bit ad hoc!
405. params[27] = c0
406. elif fixedPar==1:
407. learningSpeed = learningSpeedList[parameterNo]
408. params[27] = c0
410. # Initiate the iteration
411. for iteration in range(nOfIterations):
413. # Keep the trace of which parameters are evaluated
414. print("Processing for...")
415. print("Parameter learning speed: ", learningSpeed)
416. print("Parameter c0: ", c0)
417. print("Iteration: ", iteration)
419. # Set the initial reward for each horizon in the first episode
420. for i in range(0, nH+1):
421. reward[i][0] = 0.5 # set the initial reward to 0.5 (neutral value)
423. # Simulate each episode
424. for episodeNo in range(0, nOfEpisodes):
426. d = dSet[random.randint(0,len(dSet)-1)] # assign randomly the difficulty for the current episode
427. difficultyList[episodeNo] = d # save the difficulty of the current episode
429. # Verify that the stimuli and rewards are compatible (to avoid negative stimuli etc.)
430. lowest = minStim + d/2 + (nH+1)*gain # generated stimuli cannot be lower than this value
431. highest = maxStim - d/2 - (nH+1)*gain # generated stimuli cannot be higher than this value
432. mean0Val = np.random.uniform(lowest, highest) # initial mean of the stimuli
433. minAchievable = mean0Val - d/2 - (nH+1)*gain # to check if this value is lower than "lowest" (Cond. A)
434. maxAchievable = mean0Val + d/2 + (nH+1)*gain # to chech if this value is higher than "higher" (Cond. B)
436. while (maxAchievable > maxStim) or (minAchievable < minStim): # if Cond. A or Cond. B, regenerate your stimuli
437. print('Mean value is out of bounds! Regenerating the mean value...')
438. mean0Val = np.random.uniform(lowest, highest)
439. minAchievable = mean0Val - d/2 - (nH+1)*gain
440. maxAchievable = mean0Val + d/2 + (nH+1)*gain
443. saveDecision = np.zeros(nH+1) # initialize the decision vector, it will save which stimuli are chosen throughout the trials
444. pointer = 0
446. mean0ValInitial = mean0Val # store mean value for nH=1 simulations
448. for trial in range(0,nH+1): # initiate the first trial of the current episode
449. shuffle = np.random.randint(0,2) # to randomly shuffle the positions of two stimuli
451. ampA = mean0Val - d/2 # weak stimulus amplitude
452. ampB = mean0Val + d/2 # strong stimulus amplitude
454. # Initialize the stimuli, and by shuffling them if "shuffle" is on
455. if shuffle==1:
456. ampA_0 = ampA
457. ampA = ampB
458. ampB = ampA_0
460. stimulusA = ampA*(smoothstep(t-t0)-smoothstep(t-t0-tF)) # weak stimulus (appears at t=2)
461. stimulusB = ampB*(smoothstep(t-t0)-smoothstep(t-t0-tF)) # strong stimulus (appears at t=2)
463. # Activity evolution (SDE system is integrated in time for each trial separately)
464. psi0=reward[trial][episodeNo] # set initial condition for the regulatory mechanism
465. lambdaA, lambdaB, psi = RegulatoryPsi(psi0, stimulusA, stimulusB, params) # generate the regulated stimuli for the trial
466. state = TimeStepping(V0, lambdaA, lambdaB, TF1, TF2, params)
467. veA, veB = state[:, 0], state[:, 7]
469. # Find the decision time
470. diff = np.abs(veA - veB) # compute the difference between the excitatory populations
471. diff2 = diff[int(2/dt):-1] # ignore the first instants of the trial since the activity is degenerate the beginning
472. decisionIndex = np.where(diff2>decisionThreshold) # identify the time instant at which the decision is made (bicolumnar competition is assumed to end at this instant)
475. # Save the results and update the reward: learning the strategy through episodes
476. if len(decisionIndex[0])>0:
477. decisionTime = (decisionIndex[0][0]) * dt # pick the first instant where the difference exceeds dec. threshold
478. decisionTimeList[episodeNo] = decisionTime
479. if episodeNo<nOfEpisodes-1:
480. if nH == 1:
481. if trial==0:
482. if psi[-1]<0.5:
483. ampA = ampA + gain
484. ampB = ampB + gain
485. saveDecision[trial] = min(ampA, ampB)
486. mean1Val = (ampA + ampB)/2 # mean after the gain
487. reward[trial][episodeNo+1] = reward[trial][episodeNo] +\
488. learningSpeed*(mean1Val - mean0Val)*\
489. (2*psi[-1]-1)*(reward[trial][episodeNo]-1)**2\
490. *(reward[trial][episodeNo])**2
491. mean0Val = mean1Val # update mean value for the next trial
492. else:
493. ampA = ampA - gain
494. ampB = ampB - gain
495. saveDecision[trial] = max(ampA, ampB)
496. mean1Val = (ampA + ampB)/2 # mean after the gain
497. reward[trial][episodeNo+1] = reward[trial][episodeNo] +\
498. learningSpeed*(mean1Val - mean0Val)*(2*psi[-1]-1)*(reward[trial][episodeNo]-1)**2\
499. *(reward[trial][episodeNo])**2
500. mean0Val = mean1Val # update mean value for the next trial
501. elif trial==1:
502. if psi[-1]>0.5:
503. ampA = ampA + gain
504. ampB = ampB + gain
505. saveDecision[trial] = max(ampA, ampB)
506. mean1Val = (ampA + ampB)/2 # mean after the gain
507. reward[trial][episodeNo+1] = reward[trial][episodeNo] +\
508. learningSpeed*(mean1Val - mean0Val)*(2*psi[-1]-1)*(reward[trial][episodeNo]-1)**2\
509. *(reward[trial][episodeNo])**2
510. mean0Val = mean1Val # update mean value for the next trial
511. else:
512. ampA = ampA - gain
513. ampB = ampB - gain
514. saveDecision[trial] = min(ampA, ampB)
515. mean1Val = (ampA + ampB)/2 # mean after the gain
516. reward[trial][episodeNo+1] = reward[trial][episodeNo] +\
517. learningSpeed*(mean1Val - mean0Val)*(2*psi[-1]-1)*(reward[trial][episodeNo]-1)**2\
518. *(reward[trial][episodeNo])**2
519. mean0Val = mean1Val # update mean value for the next trial
520. elif nH == 0:
521. if psi[-1]<0.5:
522. ampA = ampA + gain
523. ampB = ampB + gain
524. saveDecision[trial] = min(ampA, ampB)
525. performance[episodeNo] = 1
526. mean1Val = (ampA + ampB)/2 # mean after the gain
527. reward[trial][episodeNo+1] = reward[trial][episodeNo] +\
528. learningSpeed*(mean1Val - mean0Val)*(2*psi[-1]-1)*(reward[trial][episodeNo]-1)**2\
529. *(reward[trial][episodeNo])**2
530. else:
531. ampA = ampA - gain
532. ampB = ampB - gain
533. saveDecision[trial] = max(ampA, ampB)
534. performance[episodeNo] = 0
535. mean1Val = (ampA + ampB)/2 # mean after the gain
536. reward[trial][episodeNo+1] = reward[trial][episodeNo] +\
537. learningSpeed*(mean1Val - mean0Val)*(2*psi[-1]-1)*(reward[trial][episodeNo]-1)**2\
538. *(reward[trial][episodeNo])**2
539. else:
540. pointer = 1
541. decisionTimeList[episodeNo] = -1
543. # Save the final performance of the episode
544. if pointer == 0:
545. if nH == 1:
546. mxBound = 2*mean0ValInitial + 3*gain
547. mnBound = 2*mean0ValInitial - 3*gain
548. performance[episodeNo] = (sum(saveDecision)-mnBound)/(mxBound - mnBound)
549. if performance[episodeNo]>0.99:
550. performance[episodeNo]=1
551. elif performance[episodeNo]<0.01:
552. performance[episodeNo]=0
554. elif nH == 0:
555. '''DO NOTHING!'''
557. elif episodeNo!=nOfEpisodes-1:
558. print("Episode %1.0i is not valid!" % episodeNo)
559. performance[episodeNo] = -1
561. learningSpeedSave2 = learningSpeed*1000 # define again to comply with indentation of Python
562. c0Save = c0*1000
563. performanceArray = np.array(performance)
565. # Comment out these three lines to turn of automatic saving of the simulation results
566. np.save("performanceResults_Iter_%d_k_%d_c0_%d.npy" % (iteration, learningSpeedSave2, c0Save), performanceArray)
567. np.save("decisionTimeList_Iter_%d_k_%d_c0_%d.npy" % (iteration, learningSpeedSave2, c0Save), decisionTimeList)
568. np.save("difficultyList_Iter_%d_k_%d_c0_%d.npy" % (iteration, learningSpeedSave2, c0Save), difficultyList)
570. print('Simulation is over.')
572. # %% [markdown]
573. # If an episode is not valid (Episode X is not valid!), it means that there was no decision made in that episode. <b>This can happen especially when the values of $k$ or $c_0$ are chosen out of the suggested ranges in the default values above.</b> We ignore such episode in the performance results. Such cases model the experimental episodes in which the subject makes no decision within in the given time duration or makes decision before the authorized period of time and so on. This is why we usually simulate more episodes than we need, in this way we can eliminate the invalid episodes and use the ones ending up with a decision for our analysis.
575. # %% [markdown]
576. # ## Case study 1: changing noise decay rate $c_0$
577. # Noise decay rate $c_0$ appears in time evolution Eq. (2) of the regulatory mechanism, and it determines how much the model has the flexibility to deviate from the strategy which it captures. This flexibility models both the perceptual difficulties that the human subject undergoes in the experimental task and the exploratory nature of the subject which sometimes makes the subject to deviate from the strategy which the subject alread captured on purpose in order to explore, see what happens in such cases of deviation.
578. #
579. # Here we will have a look at how global measures of simulation performance results change with varying $c_0$. We fix $k=0.2$. Performance results are scaled between 0 and 1, with 0 denoting the worst case performance and 1 denoting the best case performance. In the best case, the decisions are made in complete agreement with the implemented strategy. In the worst case, there is no decision made in agreement with the implemented strategy. In the all cases in between, there is a partial overlap with the made decisions and the implemented strategy.
580. #
581. # We define the following global measures to quantify the performance results:
582. #
583. # **Performance deviation:** A performance deviation is a point with a performance lower than 1 in the set of performance samples.
584. #
585. # **Deviation cluster:** A deviation cluster is a set of <u>at least</u> 3 successive performance deviations with respective to the episode numbers.
586. #
587. # **Performance cluster:** Performance cluster is the set of performance samples in which there is no deviation cluster and whose last performance sample is also the last sample of the whole experiment.
588. #
589. # **In-cluster deviation:** A performance deviation is called in-cluster if it occurs in a performance cluster.
590. # Two measures are particularly important to quantify the behavioral performance in terms of strategy learning and flexibility of the model. Those are performance cluster, since it gives at which episode the model/subject captures the implemented strategy; and in-cluster deviation, since it provides information about how much flexible the model/subject is after capturing the strategy.
591. #
592. # We will consider only Horizon 0 simulations for the sake of simplicity, i.e., we have only one trial for each episode. Let us first load the performance data which was produced by using the same routine given in the previous section.
594. # %%
595. notebook_path = os.path.abspath("MainNotebook.ipynb") # find the actual path of the Python notebook
597. # Define "deviations" finding the first and the last episode numbers of all performance deviation ranges
598. def deviations(nums):
599. nums = sorted(set(nums))
600. gaps = [[s, e] for s, e in zip(nums, nums[1:]) if s+1 < e]
601. edges = iter(nums[:1] + sum(gaps, []) + nums[-1:])
602. return list(zip(edges, edges))
604. # Set the total number of episodes in each simulation
605. nOfEpisodes = 100 # since in human experiments it is 100.
606. episodeAxis = np.linspace(1,nOfEpisodes,nOfEpisodes, dtype=int) # for the x axes of performance result plots
608. ## Import the experimental data
609. # experimentData = scipy.io.loadmat('workspace_scriptGeneralH210528.mat')
611. ## Load and arrange the experimental human performance data
612. # subject7 = experimentData['Performance'][0][6]
613. # subject7 = subject7[~np.isnan(subject7)]
614. # subject7 = subject7[(1 >= subject7) & (subject7 >= 0)]
615. # subject7 = subject7.flatten()
616. # subject7 = subject7[0:nOfEpisodes]
618. # Load the simulation performance data
619. c0Set = [1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9] # simulated c0 values
621. inClusterSet = [] # initialize an array to save number of in-cluster performance deviations for future plots
623. for c0Val in c0Set:
624. c0Iter = c0Val*1000 # to avoid rounding in the title of the data file (only for Python convention)
625. path = os.path.join(os.path.dirname(notebook_path), \
626. "showcaseData/showcase_1_H0/varying_c0/performanceResults_Iter_0_k_200_c0_%d.npy" %c0Iter) # set the data file path
628. simulationData = np.load(path) # load the simulation result of the corresponding c0
629. simulationData = simulationData[0:nOfEpisodes] # take only the first 100 episodes!
631. # Plot the performance results episode-by-episode
632. print('SIMULATION FOR c0: %.1f' %c0Val)
633. fig = plt.figure(figsize=(6,4))
634. plt.scatter(episodeAxis, simulationData)
635. plt.xlabel("episode no")
636. plt.ylabel("performance")
637. plt.show()
640. # Episode numbers of the performance deviations...
641. # print('Episode numbers of EXPERIMENT cluster gaps:', np.where(subject7 < 1)) #...of the human experiments
642. print('Episode numbers SIMULATION cluster gaps:', np.where(simulationData < 1)) #...of the simulations
643. print('Number of performance deviations:', len(np.where(simulationData < 1)[0])) # total number of performance deviations (simulations)
645. deviationsAll = np.where(simulationData < 1) # performance deviation episode numbers
646. devRanges = deviationsAll[0] # jsut for Python convention
647. devRanges = deviations(devRanges) # performance deviation ranges to find deviation blocks (thus initial episode number of performance blocks)
648. # print(devRanges) # the first and last episode numbers of deviation blocks
650. # Arrange the performance deviation ranges
651. firstEpisode = []
652. lastEpisode = []
654. for element in devRanges:
655. firstEpisode.append(element[0]) # keep the first
656. lastEpisode.append(element[1]) # keep the second
658. # Find the last deviation block, the last episode number of this block is the initial episode number of the performance block
659. devBlockIndices = list(np.where(np.subtract(lastEpisode, firstEpisode)>=2))
661. # Check if there is any deviation block. If not, the whole experiment is actually a performance block! Save also the in-cluster performance deviations.
662. if len(devBlockIndices[0]>0):
663. initialPerformanceBlock = devBlockIndices[-1][0]
664. devRanges = list(devRanges)
665. initialPerformanceBlock = devRanges[initialPerformanceBlock][1]
666. numberOfInClusterDeviations = len(np.where(deviationsAll>initialPerformanceBlock)[0]) # find the number of in-cluster performance deviations
667. inClusterSet.append(numberOfInClusterDeviations) # save the number of in-cluster performance deviations for future plots
668. # Print now also the initial episode number of the performance blocks
669. print('Initial episode number of performance block:', initialPerformanceBlock)
670. print('Number of in-cluster performance deviations:', numberOfInClusterDeviations)
671. else:
672. noDeviationBlock = len(np.where(simulationData < 1)[0][:])
673. inClusterSet.append(noDeviationBlock)
674. print('No deviation block! Performance block starts from the beginning.')
675. print('Number of in-cluster performance deviations:', noDeviationBlock)
678. print('')
679. print('')
680. print('')
683. # %% [markdown]
684. # As expected, increasing $c_0$ results in a decreasing trend in the number of in-cluster deviations as can be seen in the following plot:
686. # %%
687. fig = plt.figure(figsize=(6,4))
688. plt.plot(c0Set, inClusterSet)
689. plt.xlabel("$c_0$")
690. plt.ylabel("number of in-cluster devations")
691. plt.show()
693. # %% [markdown]
694. # ## Case study 2: changing learning speed $k$
695. # Learning speed parameter $k$ appears in episode-based discrete evolution Eq. (3) of the reward mechanism. It determines (approximately) the initial episode number of the performance cluster, which was defined in Case study 1. In other words, it determines how fast the model captures the implemented strategy. Here we will see the effect of varying $k$, where we fix $c_0=2$.
696. #
697. # An important remark is that the learning speed parameter is related to in-cluster deviations, although this relation is not explicit. In the cases where the initial episode number of performance cluster is too large, since the number of episodes for in-cluster deviations is not large, the effect of $c_0$ on the number of in-cluster episodes is outweighted by the effect of learning speed $k$. Therefore, extreme values for $k$ might produce misleading results in the simulations.
698. #
699. # Here we perform the data obtained via the same type of simulations as in Case study 1, with the only difference of that now we fix $c_0$ and vary $k$. We consider only Horizon 0 case as before for the sake of simplicity.
701. # %%
702. notebook_path = os.path.abspath("MainNotebook.ipynb") # find the actual path of the Python notebook
704. # Define "deviations" finding the first and the last episode numbers of all performance deviation ranges
705. def deviations(nums):
706. nums = sorted(set(nums))
707. gaps = [[s, e] for s, e in zip(nums, nums[1:]) if s+1 < e]
708. edges = iter(nums[:1] + sum(gaps, []) + nums[-1:])
709. return list(zip(edges, edges))
711. # Set the total number of episodes in each simulation
712. nOfEpisodes = 100 # since in human experiments it is 100.
713. episodeAxis = np.linspace(1,nOfEpisodes,nOfEpisodes, dtype=int) # for the x axes of performance result plots
715. ## Import the experimental data
716. # experimentData = scipy.io.loadmat('workspace_scriptGeneralH210528.mat')
718. ## Load and arrange the experimental human performance data
719. # subject7 = experimentData['Performance'][0][6]
720. # subject7 = subject7[~np.isnan(subject7)]
721. # subject7 = subject7[(1 >= subject7) & (subject7 >= 0)]
722. # subject7 = subject7.flatten()
723. # subject7 = subject7[0:nOfEpisodes]
725. # Load the simulation performance data
726. kSet = [0.05, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.45, 0.50] # simulated k values
728. inClusterSet = [] # initialize an array to save number of in-cluster performance deviations for future plots
730. for kVal in kSet:
731. kIter = kVal*1000 # to avoid rounding in the title of the data file (only for Python convention)
732. path = os.path.join(os.path.dirname(notebook_path), \
733. "showcaseData/showcase_2_H0/varying_k/performanceResults_Iter_0_k_%d.npy" %kIter) # set the data file path
734. simulationData = np.load(path) # load the simulation result of the corresponding k
735. simulationData = simulationData[0:nOfEpisodes] # take only the first 100 episodes!
737. # Plot the performance results episode-by-episode
738. print('SIMULATION FOR k: %.2f' %kVal)
739. fig = plt.figure(figsize=(6,4))
740. plt.scatter(episodeAxis, simulationData)
741. plt.xlabel("episode no")
742. plt.ylabel("performance")
743. plt.show()
745. # Episode numbers of the performance deviations...
746. # print('Episode numbers of EXPERIMENT cluster gaps:', np.where(subject7 < 1)) #...of the human experiments
747. print('Episode numbers SIMULATION cluster gaps:', np.where(simulationData < 1)) #...of the simulations
748. print('Number of performance deviations:', len(np.where(simulationData < 1)[0])) # total number of performance deviations (simulations)
750. deviationsAll = np.where(simulationData < 1) # performance deviation episode numbers
751. devRanges = deviationsAll[0] # jsut for Python convention
752. devRanges = deviations(devRanges) # performance deviation ranges to find deviation blocks (thus initial episode number of performance blocks)
753. # print(devRanges) # the first and last episode numbers of deviation blocks
755. # Arrange the performance deviation ranges
756. firstEpisode = []
757. lastEpisode = []
759. for element in devRanges:
760. firstEpisode.append(element[0]) # keep the first
761. lastEpisode.append(element[1]) # keep the second
763. # Find the last deviation block, the last episode number of this block is the initial episode number of the performance block
764. devBlockIndices = list(np.where(np.subtract(lastEpisode, firstEpisode)>=2))
766. # Check if there is any deviation block. If not, the whole experiment is actually a performance block! Save also the in-cluster performance deviations.
767. if len(devBlockIndices[0]>0):
768. initialPerformanceBlock = devBlockIndices[-1][0]
769. devRanges = list(devRanges)
770. initialPerformanceBlock = devRanges[initialPerformanceBlock][1]
771. numberOfInClusterDeviations = len(np.where(deviationsAll>initialPerformanceBlock)[0]) # find the number of in-cluster performance deviations
772. inClusterSet.append(numberOfInClusterDeviations) # save the number of in-cluster performance deviations for future plots
773. # Print now also the initial episode number of the performance blocks
774. print('Initial episode number of performance block:', initialPerformanceBlock)
775. print('Number of in-cluster performance deviations:', numberOfInClusterDeviations)
776. else:
777. noDeviationBlock = len(np.where(simulationData < 1)[0][:])
778. inClusterSet.append(noDeviationBlock)
779. print('No deviation block! Performance block starts from the beginning.')
780. print('Number of in-cluster performance deviations:', noDeviationBlock)
783. print('')
784. print('')
785. print('')
788. # %% [markdown]
789. # As expected, increasing $k$ results in a decreasing trend in the initial number of performance block, that is, the model learns more and more quickly as $k$ increases; see the following plot:
791. # %%
792. plt.figure(figsize=(6, 4))
793. plt.plot(kSet, inClusterSet)
794. plt.xlabel("$k$")
795. plt.ylabel("number of in-cluster devations")
796. plt.show()
798. # %% [markdown]
799. # ## Case study 3: effect of decision threshold on reaction time
800. # Each cortical column is composed of a pair of excitatory and inhibitory neuronal populations, and in our model framework, each column makes the decision in favor of the stimulus which it is sensitive to. A column is represented as pool of pair of excitatory and inhibitory populations of mean-field AdEx equations in the model. Decision is made through bicolumnar competition between those two pool. It is made by the winning pool at the instant where the difference between the excitatory population firing rates exceeds a certain threshold, which we call decision threshold and which we fix to a convenient value. The decision threshold is directly related to the reaction time of the model, which is the time duration between the moment where the stimuli are shown on the screen, and the decision is made by the model.
801. #
802. # Here we will see this relation between the decision threshold and the reaction time by focusing on single trial. We simulate a single trial via:
804. # %%
805. # Define the stimuli for one trial
806. from scipy.special import comb
808. def smoothstep(x, x_min=0, x_max=1, N=1):
809. x = np.clip((x - x_min) / (x_max - x_min), 0, 1)
811. result = 0
812. for n in range(0, N + 1):
813. result += comb(N + n, n) * comb(2 * N + 1, N - n) * (-x) ** n
815. result *= x ** (N + 1)
817. return result
819. t = np.linspace(0, tF, int(tF/dt)+1) # define the whole time interval
822. ampA = 1 # amplitude of stimulus A
823. ampB = 10 # amplitude of stimulus B
824. t0 = 2 # initial time instant of the stimuli
827. stimulusA = ampA*(smoothstep(t-t0)-smoothstep(t-t0-tF))
828. stimulusB = ampB*(smoothstep(t-t0)-smoothstep(t-t0-tF))
830. ## Build the simulation setup: load the transfer functions for both RS & FS populations
831. TF1, TF2 = LoadTransferFunctions('RS-cell', 'FS-cell', 'CONFIG1')
832. lambdaA, lambdaB, psi = RegulatoryPsi(psi0, stimulusA, stimulusB, params) # generate the regulated stimuli for the trial
834. # Time integration for one single trial
835. state = TimeStepping(V0, lambdaA, lambdaB, TF1, TF2, params)
837. # Plot Pool A (red) and B (blue) excitatory population firing rates
838. fig = plt.figure(figsize=(6,4))
839. plt.plot(t, state[:, 0], 'r') # plotting v_{eA} (red)
840. plt.plot(t, state[:, 7], 'b') # plotting v_{eB} (blue)
841. plt.xlabel("time (s)")
842. plt.ylabel("exc. population firing rates (red-eA, blue-eB)")
843. # plt.margins(x=0.001, y=-0.001)
844. plt.show()
846. # %%
847. veA, veB = state[:, 0], state[:, 7]
848. decisionThreshold = 3 # the decision is made once the absolute value of the difference between excitatory population firing rates exceeds this value.
851. # Find the reaction time
852. diff = np.abs(veA - veB) # compute the difference between the excitatory populations
853. diff2 = diff[int(2/dt):-1] # ignore the first instants of the trial since the activity is degenerate the beginning
854. decisionIndex = np.where(diff2>decisionThreshold) # identify the time instant at which the decision is made (bicolumnar competition is assumed to end at this instant)
855. decisionTime = (decisionIndex[0][0]) * dt # pick the first instant where the difference exceeds dec. threshold
856. # decisionTimeList[episodeNo] = decisionTime
857. print('Decision is made at: %.2f s' % decisionTime)
859. # %% [markdown]
860. # Let us see now how the trend of decision time looks like with respect to varying decision threshold.
862. # %%
863. inThreshold = 2 # initial decision threshold value
864. finThreshold = 10 # final decision threshold value
865. nOfSamples = 20 # number of samples of threshold values between the initial and final values
866. decisionThresholdList = np.linspace(inThreshold, finThreshold, nOfSamples) # create the list of threshold samples
868. decisionTimeList = [] # initialize reaction time list
870. # Find the reaction time for each threshold value
871. diff = np.abs(veA - veB) # compute the difference between the excitatory populations
872. diff2 = diff[int(2/dt):-1] # ignore the first instants of the trial since the activity is degenerate the beginning
874. for decisionThreshold in decisionThresholdList:
875. decisionIndex = np.where(diff2>decisionThreshold) # identify the time instant at which the decision is made (bicolumnar competition is assumed to end at this instant)
876. decisionTime = (decisionIndex[0][0]) * dt # pick the first instant where the difference exceeds dec. threshold
877. decisionTimeList.append(decisionTime)
879. # Plot the reaction time with respect to increasing decision threshold values
880. fig = plt.figure(figsize=(6,4))
881. plt.plot(decisionTimeList, decisionTimeList)
882. plt.xlabel("decision threshold")
883. plt.ylabel("reaction time (s)")
884. # plt.margins(x=0.001, y=-0.001)
885. plt.show()
888. # %% [markdown]
889. # As expected, reaction time increases as the decision threshold increases. Moreover, the relation has a rather linear character.
891. # %% [markdown]
892. # ## Bibliography
893. # [1]: E. Baspinar, G. Cecchini, M. DePass, M. Andujar, P. Pani, S. Ferraina, R. Moreno-Bote, I. Cos, A. Destexhe, "A biologically plausible decision-making model based on interacting cortical columns", bioRxiv, 2023.
894. #
895. # [2]: G. Cecchini, M. DePass, E. Baspinar, M. Andujar, S. Ramawat, P. Pani, A. Destexhe, R. Moreno-Bote, I. Cos, "A theoretical formalization of consequence-based decision-making", bioRxiv, 2023.
896. #
897. # [3]: M. di Volo, A. Romagnoni, C. Capone, A. Destexhe, "Biologically realistic mean-field models of conductance-based networks of spiking neurons with adaptation", Neural Computation, vol. 31, no. 4, pp. 653-680, 2019.
898. #
899. # [4]: R. Brette, W. Gerstner, "Adaptive exponential integrate-and-fire model as an effective description of neuronal activity", Journal of Neurophysiology, vol. 94, no. 5, pp. 3637-3642, 2005.
900. #
901. # [5]: S. El Boustani, A. Destexhe, "A master equation formalism for macroscopic modeling of asynchronous irregular activity states", Neural Computation, vol. 21, no. 1, pp.46-100, 2009.
902. #
903. # [6]: Y. Zerlaut, S. Chemla, F. Chavane, A. Destexhe, "Modeling mesoscopic cortical dynamics using a mean-field model of conductance-based networks of adaptive exponential integrate-and-fire neurons", Journal of Computational Neuroscience, vol. 44, no. 1, pp. 45-61, 2018.
904. #
905. # [7]: E. Baspinar, G. Cecchini, R. Moreno-Bote, I. Cos, A. Destexhe, "Jupyter notebook of a biophysically plausible decision-making model based on interacting cortical columns (1.0.0)", Zenodo, 2023, https://doi.org/10.5281/zenodo.7682309

MainNotebook.ipynb at commit 55fd876, under CC-BY-4.0 · at the source