.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/binary_rate_simulation/binary_sim_vs_thy.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note Click :ref:`here ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_binary_rate_simulation_binary_sim_vs_thy.py: Binary rates: simulation vs mean-field ====================================== .. image:: ../../../../examples/binary_rate_simulation/binary.png :width: 1000 :alt: Plot of simulated and estimated rates Here we simulate an E-I network of binary neurons, calculate the mean-field estimate of the firing rates, and plot them together. .. GENERATED FROM PYTHON SOURCE LINES 12-18 .. code-block:: default import nnmt import numpy as np import matplotlib.pyplot as plt .. GENERATED FROM PYTHON SOURCE LINES 19-23 Helper functions ---------------- First, we define the functions used to perform the simulation. .. GENERATED FROM PYTHON SOURCE LINES 23-109 .. code-block:: default def _sort_queue(q): """ Sorts a 2d array of id and update time point by update time points. Sorts the queue of update time points id_0, t_next_0 id_1, t_next_1 ... by the update time point t_next. Parameters ---------- q : array 2d array of neuron ids and update time points. Returns ------- array Sorted array of ids and update time points. """ return q[np.argsort(q[:, 1], kind='stable')] def update_poisson(J, S, T, tau, thetas): """ Simulates a network of binary neurons. Evolves the initial network state `S` in time by drawing exponentially distributed update times (Poisson process) with time constant `tau` until time `T` is reached, using connectivity matrix `J` and thresholds `thetas`. Parameters ---------- J : array Connectivity matrix. S : array Initial state of each neuron. T : float Simulation time. tau : float Time constant of all binary neurons. thetas : [array|list] Thresholds of neurons. Returns ------- array Update times. array State of each neuron at each update time point. array Mean population activity at each update time point. """ # check dimensions of J and S assert J.shape[0] == J.shape[1] == S.shape[0] # get total number of neurons N = S.shape[0] # create update queue update_queue = np.empty((N, 2), dtype=np.float32) update_queue[:, 0] = np.arange(N) # draw first update time point for every neuron # from an exponential distribution with mean tau update_queue[:, 1] = np.random.exponential(tau, N) # sort queue ascendingly according to update time update_queue = _sort_queue(update_queue) # initialize storage lists ts, m, Ss = [], [], [] t = 0 while t < T: # select neuron with next update time point i, t = update_queue[0, :] i = int(i) # input to neuron i h_i = J[i, :].dot(S) # update state of neuron i S[i] = np.heaviside(h_i-thetas[i], 0) # draw new update time update_queue[0, 1] += np.random.exponential(tau) update_queue = _sort_queue(update_queue) m.append(S.mean()) ts.append(t) Ss.append(S.copy()) return np.array(ts), np.array(Ss).T, np.array(m) .. GENERATED FROM PYTHON SOURCE LINES 110-112 Here we define the functions that construct the network properties in a format needed to perform the simulation. .. GENERATED FROM PYTHON SOURCE LINES 112-174 .. code-block:: default def fixed_indegree_connectivity(N, J, K): """ Constructs a fixed indegree matrix for the given network parameters. Parameters ---------- N : array of ints Number of neurons in each population. J : array of floats Weight matrix. K : array of ints Indegree matrix. Returns ------- array Connectivity matrix. """ W = np.zeros((N.sum(), N.sum()), dtype=np.float32) # list of neurons (each one gets a unique number) neurons = np.arange(N.sum()) population_ix = N.cumsum() for pre_ix, pre_pop in enumerate( np.array_split(neurons, population_ix[:-1])): for post_ix, post_pop in enumerate( np.array_split(neurons, population_ix[:-1])): for post_neuron in post_pop: pre_neurons = np.random.choice(pre_pop[pre_pop!=post_neuron], size=K[post_ix][pre_ix], replace=False) W[post_neuron, pre_neurons] = J[post_ix][pre_ix] return W def neuron_thresholds(N, theta): """ Creates a list of thresholds for each neuron. Parameters ---------- N : array Numbers of neurons in each population. theta : array Threshold for each population Returns ------- array Threshold for each individual neuron numbered from 0 to N-1. """ thresholds = np.zeros(N.sum()) population_ix = np.append([0], N.cumsum()) for i in range(len(population_ix[:-1])): thresholds[population_ix[i]: population_ix[i+1]] = theta[i] return thresholds .. GENERATED FROM PYTHON SOURCE LINES 175-179 Parameter definition -------------------- Here we define the network parameters .. GENERATED FROM PYTHON SOURCE LINES 179-194 .. code-block:: default # numbers of neurons in each population and their respective name N = np.array([1000, 1000]) label = ['E', 'I'] # indegree matrix K = np.array([[150, 200], [350, 200]]) # weight matrix J = np.array([[0.1, -0.2], [0.1, -0.2]]) network_params = {'J': J, 'K': K, 'N': N} .. GENERATED FROM PYTHON SOURCE LINES 195-201 Calculation of mean-field estimate ---------------------------------- Lets define the network model. We decided to use the ``Plain`` model, as we just want to load the parameters into the models dicts and do not want to calculate any dependent parameters from them. .. GENERATED FROM PYTHON SOURCE LINES 201-205 .. code-block:: default # here we have to copy the dict due to a bug in ``nnmt.models.Network`` network = nnmt.models.Plain(dict(network_params)) .. GENERATED FROM PYTHON SOURCE LINES 206-209 We could set the firing threshold of the neurons directly, but here we decided to calculate the threshold using a balanced condition for some expected rates .. GENERATED FROM PYTHON SOURCE LINES 209-213 .. code-block:: default expected_rates = [0.7, 0.4] theta = nnmt.binary.balanced_threshold(network, expected_rates) network.network_params['theta'] = theta .. GENERATED FROM PYTHON SOURCE LINES 214-215 We calculate the mean-field estimate of the rates .. GENERATED FROM PYTHON SOURCE LINES 215-217 .. code-block:: default rates_thy = nnmt.binary.mean_activity(network) .. GENERATED FROM PYTHON SOURCE LINES 218-223 Simulation ---------- For simulating the network we require a concrete realization of the network connectivity .. GENERATED FROM PYTHON SOURCE LINES 223-225 .. code-block:: default W = fixed_indegree_connectivity(**network_params) .. GENERATED FROM PYTHON SOURCE LINES 226-228 And we are going to use a list that defines the thresholds of each neuron separately .. GENERATED FROM PYTHON SOURCE LINES 228-231 .. code-block:: default thresholds = neuron_thresholds(network.network_params['N'], network.network_params['theta']) .. GENERATED FROM PYTHON SOURCE LINES 232-234 The simulation additionally requires to define the time constant of the neurons .. GENERATED FROM PYTHON SOURCE LINES 234-236 .. code-block:: default sim_params = {'tau': 1., 'thetas': thresholds} .. GENERATED FROM PYTHON SOURCE LINES 237-238 Then we can run the simulation .. GENERATED FROM PYTHON SOURCE LINES 238-246 .. code-block:: default # simulation time T = 10 # initial state S = np.zeros(W.shape[0], dtype=np.float32) # simulate t, S, m = update_poisson(W, S, T=T, **sim_params) .. GENERATED FROM PYTHON SOURCE LINES 247-248 Now we calculate the mean rates for each population .. GENERATED FROM PYTHON SOURCE LINES 248-251 .. code-block:: default rates_sim = np.array([pop.mean(axis=0) for pop in np.array_split(S, N.cumsum()[:-1])]) .. GENERATED FROM PYTHON SOURCE LINES 252-254 Plotting -------- .. GENERATED FROM PYTHON SOURCE LINES 254-270 .. code-block:: default fig, axs = plt.subplots(1, len(N), figsize=(6, 2.5)) for i in range(len(N)): axs[i].plot(t, rates_sim[i], '-', color='k', label='sim') axs[i].plot(t, rates_thy[i]*np.ones_like(t), '--', color='gray', label='thy') axs[i].set_xlim(0, T) axs[i].set_ylim(-0.1, 1.1) axs[i].set_xlabel('time') axs[i].set_ylabel('mean activity') axs[i].legend(loc=0) axs[i].set_title(f'{label[i]} population') plt.tight_layout() plt.savefig('binary.png', dpi=600) plt.show() .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 0.000 seconds) .. _sphx_glr_download_auto_examples_binary_rate_simulation_binary_sim_vs_thy.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: binary_sim_vs_thy.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: binary_sim_vs_thy.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_