diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml new file mode 100644 index 0000000..1203161 --- /dev/null +++ b/.github/workflows/ci.yml @@ -0,0 +1,66 @@ +name: CI + +on: + push: + branches: + - main + pull_request: + +permissions: + actions: write + contents: read + +jobs: + test: + runs-on: ubuntu-latest + + steps: + # Step 1: Check out the repository + - name: Checkout repository + uses: actions/checkout@v3 + + # Step 2: Set up Julia + - name: Set up Julia + uses: julia-actions/setup-julia@latest + with: + version: '1.9' # Specify the Julia version you need + + # Step 3: Set up cache for dependencies + - name: Setup cache + uses: julia-actions/cache@v2 + + # Step 4: Install dependencies (registered and unregistered) + - name: Install Dependencies + run: | + julia -e ' + import Pkg; + # Add unregistered dependency from GitHub + Pkg.develop(Pkg.PackageSpec(url="https://github.com/JuliaSNN/SpikingNeuralNetworks.jl")); + # Install all dependencies listed in Project.toml + Pkg.update(); + Pkg.resolve(); + Pkg.instantiate(); + ' + + # Step 5: Run tests + # Step 5: Run tests with force_latest_compatible_version=true + - name: Run tests + run: | + julia -e ' + import Pkg; + # Add unregistered dependency from GitHub + Pkg.develop(Pkg.PackageSpec(url="https://github.com/JuliaSNN/SpikingNeuralNetworks.jl")); + # Install all dependencies listed in Project.toml + Pkg.update(); + Pkg.resolve(); + Pkg.instantiate(); + Pkg.activate("."); + Pkg.test() + ' + # Step 6: Upload test results on failure + - name: Upload test results + if: failure() + uses: actions/upload-artifact@v3 + with: + name: test-results + path: test/ diff --git a/Project.toml b/Project.toml index 81648c0..799714b 100644 --- a/Project.toml +++ b/Project.toml @@ -1 +1,23 @@ +name = "SNNExamples" +uuid = "276b8b60-e50f-4651-bb99-b251f196bd69" + [deps] +Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" +DrWatson = "634d3b9d-ee7a-5ddf-bec9-22491ea816e1" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" + +[extras] +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +MLJ = "add582a8-e3ab-11e8-2d5e-e98b27df1bc7" +Reexport = "189a3867-3050-52da-a836-e630ba90ab69" + +[targets] +test = ["Test"] + +[compat] +SpikingNeuralNetworks = "0.1.0" +MLJ = "0.20.7" +Reexport = "1.2.2" diff --git a/README.md b/README.md index 1796d21..0c7e888 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,3 @@ # SNNExamples +![CI](https://github.com/russelljjarvis/SNNExamples/actions/workflows/ci.yml/badge.svg) Collection of examples from the JuliaSNN organization diff --git a/experiments/Mongillo2008_WorkingMemory/STP_WM_Mongillo2008.jl b/experiments/Mongillo2008_WorkingMemory/STP_WM_Mongillo2008.jl index 9c8e82b..fe8815b 100644 --- a/experiments/Mongillo2008_WorkingMemory/STP_WM_Mongillo2008.jl +++ b/experiments/Mongillo2008_WorkingMemory/STP_WM_Mongillo2008.jl @@ -9,8 +9,8 @@ using StatsBase using Distributions using LaTeXStrings -## -include("../../parameters/Mongillo_WM2008.jl") + +include("../../parameters/mongillo_WM2008.jl") model, assemblies = Mongillo2008(n_assemblies=2) peak_rate = 2kHz diff --git a/experiments/PotjansDiesmann2014_CorticalCircuit/multilayer_potjans-diesmann2015.jl b/experiments/PotjansDiesmann2014_CorticalCircuit/multilayer_potjans-diesmann2015.jl index afb9ed4..b037dd5 100644 --- a/experiments/PotjansDiesmann2014_CorticalCircuit/multilayer_potjans-diesmann2015.jl +++ b/experiments/PotjansDiesmann2014_CorticalCircuit/multilayer_potjans-diesmann2015.jl @@ -10,6 +10,85 @@ using ProgressMeter using Plots using SpikingNeuralNetworks using SNNUtils +using JLD2 +using Distributions + + +""" +Model parameters pasted here from original 2015 paper. + +https://pmc.ncbi.nlm.nih.gov/articles/PMC3920768/ + +Layer-specific external inputs +External inputs L2/3 L4 L5 L6 +Total 1606 2111 1997 2915 + + +Parameter specification +Populations and inputs +Name L2/3e L2/3i L4e L4i L5e L5i L6e L6i Th +Population size, N 20 683 5834 21 915 5479 4850 1065 14 395 2948 902 +External inputs, kext (reference) 1600 1500 2100 1900 2000 1900 2900 2100 n/a +External inputs, kext (layer independent) 2000 1850 2000 1850 2000 1850 2000 1850 n/a +Background rate, νbg 8 Hz +Connectivity + + from (columns are presynaptic, rows are post synaptic) + L2/3e L2/3i L4e L4i L5e L5i L6e L6i Th +to L2/3e 0.101 0.169 0.044 0.082 0.032 0.0 0.008 0.0 0.0 + L2/3i 0.135 0.137 0.032 0.052 0.075 0.0 0.004 0.0 0.0 + L4e 0.008 0.006 0.050 0.135 0.007 0.0003 0.045 0.0 0.0983 + L4i 0.069 0.003 0.079 0.160 0.003 0.0 0.106 0.0 0.0619 + L5e 0.100 0.062 0.051 0.006 0.083 0.373 0.020 0.0 0.0 + L5i 0.055 0.027 0.026 0.002 0.060 0.316 0.009 0.0 0.0 + L6e 0.016 0.007 0.021 0.017 0.057 0.020 0.040 0.225 0.0512 + L6i 0.036 0.001 0.003 0.001 0.028 0.008 0.066 0.144 0.0196 +Name Value Description +w ± δw 87.8 ± 8.8 pA Excitatory synaptic strengths +g –4 Relative inhibitory synaptic strength +de ± δde 1.5 ± 0.75 ms Excitatory synaptic transmission delays +di ± δdi 0.8 ± 0.4 ms Inhibitory synaptic transmission delays +Neuron model +Name Value Description +τm 10 ms Membrane time constant +τref 2 ms Absolute refractory period +τsyn 0.5 ms Postsynaptic current time constant +Cm 250 pF Membrane capacity +Vreset −65 mV Reset potential +θ −50 mV Fixed firing threshold +νth 15 Hz Thalamic firing rate during input perio + + +Populations Nine; 8 cortical populations and 1 thalamic population +Topology — +Connectivity Random connections +Neuron model Cortex: Leaky integrate-and-fire, fixed voltage threshold, fixed absolute refractory period (voltage clamp), thalamus: Fixed-rate Poisson +Synapse model Exponential-shaped postsynaptic currents +Plasticity — +Input Cortex: Independent fixed-rate Poisson spike trains +Measurements Spike activity, membrane potentials +Populations +Type Elements +Cortical network iaf neurons, 8 populations (2 per layer), type specific size N +Th Poisson, 1 population, size Nth +Connectivity +Type Random connections with independently chosen pre- and postsynaptic neurons; see Table 5 for probabilities +Weights Fixed, drawn from Gaussian distribution +Delays Fixed, drawn from Gaussian distribution multiples of computation stepsize +Neuron and synapse model +Name iaf neuron +Type Leaky integrate-and-fire, exponential-shaped synaptic current inputs +Subthreshold dynamics Inline graphic + V(t) = Vreset else + Inline graphic +Spiking If Inline graphic + 1. set t* = t, 2. emit spike with time stamp t* +Input +Type Description +Background Independent Poisson spikes to iaf neurons (Table 5) +Measurements +Spiking activity and membrane potentials from a subset of neurons in every population +""" """ @@ -30,7 +109,8 @@ function potjans_neurons(scale=1.0) neurons = Dict{Symbol, SNN.AbstractPopulation}() for (k, v) in ccu if occursin("E", String(k)) - neurons[k] = AdEx(N = v, param=LKD2014SingleExp.AdEx, name=string(k)) + # PD model only consists of LIF neurons. + neurons[k] = IF(N = v, param=LKD2014SingleExp.PV, name=string(k)) else neurons[k] = IF(N = v, param=LKD2014SingleExp.PV, name=string(k)) end @@ -42,23 +122,43 @@ end Define Potjans parameters for neuron populations and connection probabilities """ function potjans_conn(Ne) + # adding static synaptic weight parameters + # based on PSPs that are more accurate to PD based on the PyNN + # model + # Descrepency between PyNN model and original paper + # original paper + # w ± δw 87.8 ± 8.8 pA Excitatory synaptic strengths + # PyNN model: syn_weight_dist = Normal(0.15, 0.1) - function j_from_name(pre, post, g) + # g –4 Relative inhibitory synaptic strength + # + # + function j_from_name(pre, post) if occursin("E", String(pre)) && occursin("E", String(post)) - return g.jee + syn_weight_dist = Normal(87.8, 8.8) # Unit is pA + + return rand(syn_weight_dist) elseif occursin("I", String(pre)) && occursin("E", String(post)) - return g.jei + syn_weight_dist = Normal(0.15, 0.1) + return -4.0#*rand(syn_weight_dist) + elseif occursin("E", String(pre)) && occursin("I", String(post)) - return g.jie + syn_weight_dist = Normal(87.8, 8.8) # Unit is pA + + return rand(syn_weight_dist) elseif occursin("I", String(pre)) && occursin("I", String(post)) - return g.jii + syn_weight_dist = Normal(0.15, 0.1) + return -4.0#*rand(syn_weight_dist) else throw(ArgumentError("Invalid pre-post combination: $pre-$post")) end end + + layer_names = [:E23, :I23, :E4, :I4, :E5, :I5, :E6, :I6] + # Replace static matrix with a regular matrix for `conn_probs` # ! the convention is j_post_pre. This is how the matrices `w` are built. Are you using that when defining the parameters? conn_probs = Float32[ @@ -72,33 +172,64 @@ function potjans_conn(Ne) 0.0364 0.001 0.0034 0.0005 0.0277 0.008 0.0658 0.1443 ] - # !Assign dimensions to these parameters, and some reference - pree = 0.1 - g = 1.0 - tau_meme = 10.0 # (ms) - # Synaptic strengths for each connection type - K = round(Int, Ne * pree) - sqrtK = sqrt(K) - - # !Same, the convention is j_post_pre - je = 2.0 / sqrtK * tau_meme * g - ji = 2.0 / sqrtK * tau_meme * g - jee = 0.15 * je - jie = je - jei = -0.75 * ji - jii = -ji - g_strengths = dict2ntuple(SNN.@symdict jee jie jei jii) - + conn_j = zeros(Float32, size(conn_probs)) for pre in eachindex(layer_names) for post in eachindex(layer_names) - conn_j[post, pre ] = j_from_name(layer_names[pre], layer_names[post], g_strengths) + + conn_j[post, pre ] = j_from_name(layer_names[pre], layer_names[post]) + end end return layer_names, conn_probs, conn_j end +function unused_volumetric_part_of_model() + + total_cortical_thickness = 1500.0 + N_full = [20683, 5834, 21915, 5479, 4850, 1065, 14395, 2948] + N_E_total = N_full[1]+N_full[3]+N_full[5]+N_full[7] + dimensions_3D = Dict( + "x_dimension"=> 1000, + "z_dimension"=> 1000, + "total_cortical_thickness"=> total_cortical_thickness, + + # Have the thicknesses proportional to the numbers of E cells in each layer + "layer_thicknesses"=> Dict( + "L23"=> total_cortical_thickness*N_full[1]/N_E_total, + "L4" => total_cortical_thickness*N_full[3]/N_E_total, + "L5" => total_cortical_thickness*N_full[5]/N_E_total, + "L6" => total_cortical_thickness*N_full[7]/N_E_total, + "thalamus" => 100 + ) + ) + + net_dict = Dict{String, Any}( + "PSP_e"=> 0.15, + # Relative standard deviation of the postsynaptic potential. + "PSP_sd"=> 0.1, + # Relative inhibitory synaptic strength (in relative units). + "g"=> -4, + # Rate of the Poissonian spike generator (in Hz). + "bg_rate"=> 8., + # Turn Poisson input on or off (True or False). + "poisson_input"=> true, + # Delay of the Poisson generator (in ms). + "poisson_delay"=> 1.5, + # Mean delay of excitatory connections (in ms). + "mean_delay_exc"=> 1.5, + # Mean delay of inhibitory connections (in ms). + "mean_delay_inh"=> 0.75, + # Relative standard deviation of the delay of excitatory and + # inhibitory connections (in relative units). + "rel_std_delay"=> 0.5 + ) + dimensions_3D,net_dict + + +end + """ Main function to setup Potjans layer with memory-optimized connectivity @@ -111,8 +242,15 @@ function potjans_layer(scale) inh_pop = filter(x -> occursin("I", String(x)), keys(neurons)) Ne = trunc(Int32, sum([neurons[k].N for k in exc_pop])) Ni = trunc(Int32, sum([neurons[k].N for k in inh_pop])) - @show exc_pop, Ne, Ni layer_names, conn_probs, conn_j = potjans_conn(Ne) + ## + # Added synaptic delays from distribution informed by PyNN/PD model + ## + # de ± δde 1.5 ± 0.75 ms Excitatory synaptic transmission delays + # di ± δdi 0.8 ± 0.4 ms Inhibitory synaptic transmission delays + + delay_dist_exc = Normal(1.5, 0.75) # Unit is ms + delay_dist_inh = Normal(0.8, 0.4) # Unit is ms ## Create the synaptic connections based on the connection probabilities and synaptic weights assigned to each pre-post pair connections = Dict() @@ -124,24 +262,32 @@ function potjans_layer(scale) J = conn_j[j, i] sym = J>=0 ? :ge : :gi μ = abs(J) - s = SNN.SpikingSynapse(neurons[pre], neurons[post], sym; μ = μ, p=p, σ=0) + if J>=0 + s = SNN.SpikingSynapse(neurons[pre], neurons[post], sym; μ = μ, p=p, σ=0, delay_dist=delay_dist_exc) + else + s = SNN.SpikingSynapse(neurons[pre], neurons[post], sym; μ = -μ, p=p, σ=0, delay_dist=delay_dist_inh) + + end connections[Symbol(string(pre,"_", post))] = s end end ## Create the Poisson stimulus for each population + ## The order of these rates needs to be carefully checked as I have changed the order of populations from the + # PyNN model to conveniently group excitatory and inhibitory connections. + full_mean_rates = [0.971, 2.868, 4.746, 5.396, 8.142, 9.078, 0.991, 7.523] stimuli = Dict() - for pop in exc_pop - νe = 2.5kHz + for (ind,pop) in enumerate(exc_pop) + νe = full_mean_rates[ind]kHz post = neurons[pop] s = SNN.PoissonStimulus(post, :ge; param = νe, cells=:ALL, μ=1.f0, name="PoissonE_$(post.name)") stimuli[Symbol(string("PoissonE_", pop))] = s end - return merge_models(neurons, connections, stimuli) + return merge_models(neurons,connections, stimuli),neurons,connections,stimuli end -model = potjans_layer(0.01) +model,neurons,connections,stimuli = potjans_layer(0.01) duration = 15000ms SNN.monitor([model.pop...], [:fire]) SNN.monitor([model.pop...], [:v], sr=200Hz) @@ -150,12 +296,13 @@ SNN.raster(model.pop, [10s, 15s]) Trange = 5s:10:15s frE, interval, names_pop = SNN.firing_rate(model.pop, interval = Trange) -plot(interval, mean.(frE), label=hcat(names_pop...), xlabel="Time [ms]", ylabel="Firing rate [Hz]", legend=:topleft) +display(plot(interval, mean.(frE), label=hcat(names_pop...), xlabel="Time [ms]", ylabel="Firing rate [Hz]", legend=:topleft)) ## -vecplot(model.pop.E23, :v, neurons =1, r=0s:15s,label="soma") -layer_names, conn_probs, conn_j = potjans_conn(4000) -pj = heatmap(conn_j, xticks=(1:8,layer_names), yticks=(1:8,layer_names), aspect_ratio=1, color=:bluesreds, title="Synaptic weights", xlabel="Presynaptic", ylabel="Postsynaptic", size=(500,500), clims=(-maximum(abs.(conn_j)), maximum(abs.(conn_j)))) -pprob=heatmap(conn_probs, xticks=(1:8,layer_names), yticks=(1:8,layer_names), aspect_ratio=1, color=:viridis, title="Connection probability", xlabel="Presynaptic", ylabel="Postsynaptic", size=(500,500)) -plot(pprob, pj, layout=(1,2), size=(1000,500), margin=5Plots.mm) +display(vecplot(model.pop.E23, :v, neurons =1, r=0s:15s,label="soma")) + +#layer_names, conn_probs, conn_j = potjans_conn(4000) +#pj = heatmap(conn_j, xticks=(1:8,layer_names), yticks=(1:8,layer_names), aspect_ratio=1, color=:bluesreds, title="Synaptic weights", xlabel="Presynaptic", ylabel="Postsynaptic", size=(500,500), clims=(-maximum(abs.(conn_j)), maximum(abs.(conn_j)))) +#pprob=heatmap(conn_probs, xticks=(1:8,layer_names), yticks=(1:8,layer_names), aspect_ratio=1, color=:viridis, title="Connection probability", xlabel="Presynaptic", ylabel="Postsynaptic", size=(500,500)) +#plot(pprob, pj, layout=(1,2), size=(1000,500), margin=5Plots.mm) ## diff --git a/experiments/Quaresima2023_SequenceRecognition/tripod_quaresima2023.jl b/experiments/Quaresima2023_SequenceRecognition/tripod_quaresima2023_broken_or_depricated.jl similarity index 100% rename from experiments/Quaresima2023_SequenceRecognition/tripod_quaresima2023.jl rename to experiments/Quaresima2023_SequenceRecognition/tripod_quaresima2023_broken_or_depricated.jl diff --git a/src/SNNExamples.jl b/src/SNNExamples.jl new file mode 100644 index 0000000..e69de29 diff --git a/test/Tonic_NMNIST_Stimulus.jl b/test/Tonic_NMNIST_Stimulus.jl new file mode 100644 index 0000000..b339ece --- /dev/null +++ b/test/Tonic_NMNIST_Stimulus.jl @@ -0,0 +1,98 @@ +module Tonic_NMNIST_Stimulus + using JLD2 + using PyCall + using DataFrames + + using Random + using Plots + """ + NMNIST data set is a 3D labelled data set played out onto a 400 by 400 matrix of pixels over time. For any time instant there can be positive and negative spikes. + Here the 2D pixel matrix is collapsed into 1D so that it the matrix can be more easily projected onto the synapses of a biological network population. + """ + + """ + Call Pythons tonic module to get NMNIST a spiking data set of numerals. + """ + + function getNMNIST(ind)#, cnt, input_shape) + pushfirst!(PyVector(pyimport("sys")."path"), "") + nmnist_module = pyimport("batch_nmnist_motions") + dataset::PyObject = nmnist_module.NMNIST("./") + A = zeros((35,35)) + I = LinearIndices(A) + pop_stimulation= Vector{Int32}([])#Vector{UInt32}([]) + l = -1 + events,l = dataset._dataset[ind]#._dataset(ind)#;limit=15) + l = convert(Int32,l) + x = Vector{Int32}([e[1] for e in events]) + y = Vector{Int32}([e[2] for e in events]) + ts = Vector{Float32}([e[3]/100000.0 for e in events]) + p = Vector{Int8}([e[4] for e in events]) + for (x_,y_) in zip(x,y) + push!(pop_stimulation,Int32(I[CartesianIndex(convert(Int32,x_+1),convert(Int32,y_+1))])) + end + (ts,p,l,pop_stimulation) + end + + function spike_packets_by_labels(cached_spikes) + ts = [cached_spikes[i][1] for (i,l) in enumerate(cached_spikes)] + p = [cached_spikes[i][2] for (i,l) in enumerate(cached_spikes)] + labels = [cached_spikes[i][3] for (i,l) in enumerate(cached_spikes)] + pop_stimulation = [cached_spikes[i][4] for (i,l) in enumerate(cached_spikes)] + + # Create DataFrame + df = DataFrame( + ts = ts, + p = p, + labels = labels, + pop_stimulation = pop_stimulation + ) + + # Sort the DataFrame by `labels` + sorted_df = sort(df, :labels) + + # Group by `labels` to create indexed/enumerated groups + grouped_df = groupby(sorted_df, :labels) + unique_labels = unique(labels) + labeled_packets = Dict() + for filter_label in unique_labels + time_and_offset=[] + population_code=[] + next_offset = 0 + group = grouped_df[filter_label+1] + for row in eachrow(group) + append!(population_code,row["pop_stimulation"]) + append!(time_and_offset,next_offset.+row["ts"]) + windowb = minimum(row["ts"]) + windowa = maximum(row["ts"]) + next_offset = windowa-windowb + + + labeled_packets[filter_label] = (population_code,time_and_offset) + + end + end + labeled_packets + end + + """ + For a given numeral characters label plot the corresponding spike packets. + """ + function spike_character_plot(filter_label,population_code,time_and_offset) + p1 = Plots.scatter() + scatter!(p1, + population_code, + time_and_offset, + ms = 1, # Marker size + ylabel = "Time (ms)", + xlabel = "Neuron Index", + title = "Spiking Activity with Distinct Characters", + legend=false + ) + plot(p1) + savefig("label$filter_label.png") + end +end + + + diff --git a/test/bistable_network_with_stimulus.jl b/test/bistable_network_with_stimulus.jl new file mode 100644 index 0000000..cd26aed --- /dev/null +++ b/test/bistable_network_with_stimulus.jl @@ -0,0 +1,161 @@ +using Test + +using DrWatson +using Revise +using SpikingNeuralNetworks +SpikingNeuralNetworks.@load_units; +using SNNUtils +using Plots +using Statistics + +# Define each of the network recurrent assemblies +function define_network(N = 800) + # Number of neurons in the network + N = N + # Create dendrites for each neuron + E = SNN.AdEx(N = N, param = SNN.AdExParameter(Vr = -55mV, At = 0mV, b=0, a=0)) + # Define interneurons + I = SNN.IF(; N = N ÷ 2, param = SNN.IFParameter(τm = 20ms, El = -50mV)) + # Define synaptic interactions between neurons and interneurons + E_to_I = SNN.SpikingSynapse(E, I, :he, p = 0.2, μ = 2.0) + E_to_E = SNN.SpikingSynapse(E, E, :he, p = 0.2, μ = 3.0)#, param = SNN.vSTDPParameter()) + I_to_I = SNN.SpikingSynapse(I, I, :hi, p = 0.2, μ = 1.0) + I_to_E = SNN.SpikingSynapse( + I, + E, + :hi, + p = 0.2, + μ = 1, + param = SNN.iSTDPParameterRate(r = 4Hz), + ) + norm = SNN.SynapseNormalization(E, [E_to_E], param = SNN.AdditiveNorm(τ = 10ms)) + + # Store neurons and synapses into a dictionary + pop = dict2ntuple(@strdict E I) + syn = dict2ntuple(@strdict I_to_E E_to_I E_to_E norm I_to_I) + # Return the network as a tuple + SNN.monitor([E, I], [:fire]) + (pop = pop, syn = syn) +end + +n_assemblies = 2 +## Instantiate the network assemblies and local inhibitory populations +subnets = Dict(Symbol("sub_$n") => define_network(400) for n = 1:n_assemblies) +# Add noise to each assembly +noise = Dict(Symbol("noise_$(i)") => SNN.PoissonStimulus(subnets[i].pop.E, :he, param=2.5kHz, cells=:ALL) for i in eachindex(subnets)) +# Create synaptic connections between the assemblies and the lateral inhibitory populations +syns = Dict{Symbol,Any}() +for i in eachindex(subnets) + for j in eachindex(subnets) + i == j && continue + push!( + syns, + Symbol("lateral_$(i)E_to_$(j)I") => SNN.SpikingSynapse( + subnets[i].pop.E, + subnets[j].pop.I, + :he, + p = 0.2, + μ = 2.25, + ), + ) + end +end + +# select only excitatory populations +input = function (t, param::PSParam) + id::Int = param.variables[:id] + n_assemblies::Int = param.variables[:n_assemblies] + if (t ÷ 1000)%n_assemblies+1 == id + return 5kHz * exp(-(t%1000)/200) + else + return 0 + end +end +# select only excitatory populations + + +stimuli = Dict{Symbol,Any}() +for n in 1:n_assemblies + push!(stimuli, Symbol("stim_$n")=> + SNN.PoissonStimulus( + subnets[Symbol("sub_$n")].pop.E, + :ge, + cells=:ALL, + param=PSParam( + rate=input, + variables = Dict( + :id => n, + :n_assemblies => n_assemblies)))) +end +network = SNN.merge_models(noise, subnets, syns, stimuli) + + +## Merge the models and run the simulation, the merge_models function will return a model object (syn=..., pop=...); the function has strong type checking, see the documentation. +train!(model = network, duration = 5000ms, pbar = true, dt = 0.125) + +## Create a model object with only the populations to run the analysis +SNN.monitor([network.pop...], [:fire]) + +# Define a time object to keep track of the simulation time, the time object will be passed to the train! function, otherwise the simulation will not create one on the fly. +time_keeper = SNN.Time() +train!(model = network, duration = 15000ms, time = time_keeper, pbar = true, dt = 0.125) + + +# Plot the raster plot of the network +SNN.raster(network.pop, [4s, 15s]) +## + +@test mean(network.syn.sub_2_I_to_E.W)!=0.0 +# define the time interval for the analysis + +exc_populations = SNN.filter_populations(populations, :E) + +# get the spiketimes of the excitatory populations and the indices of each population +spiketimes = SNN.spiketimes(exc_populations) +indices = SNN.population_indices(populations, :E) + +# calculate the firing rate of each excitatory population +interval = 0:5:SNN.get_time(time_keeper) +rates = map(eachindex(indices)) do i + rates, intervals = + SNN.firing_rate(spiketimes, interval = interval, pop = indices[i], τ = 10) + mean_rate = mean(rates) +end + +## Plot the firing rate of each assembly and the correlation matrix +p1 = plot() +for i in eachindex(rates) + plot!( + interval, + rates[i], + label = "Assembly $i", + xlabel = "Time (ms)", + ylabel = "Firing rate (Hz)", + xlims = (4_000, 15_000), + legend = :topleft, + ) +end +plot!() + +cor_mat = zeros(length(rates), length(rates)) +for i in eachindex(rates) + for j in eachindex(rates) + cor_mat[i, j] = cor(rates[i], rates[j]) + end +end +p2 = heatmap( + cor_mat, + c = :bluesreds, + clims = (-1, 1), + xlabel = "Assembly", + ylabel = "Assembly", + title = "Correlation matrix", + xticks = 1:3, + yticks = 1:3, +) +plot(p1, p2, layout = (2, 1), size = (600, 800), margin = 5Plots.mm) + +using Statistics, StatsBase +pcor = plot() +n = 800 +plot(0:5:5*n,autocor(rates[1], 0:n)) diff --git a/test/long_short_memory.jl b/test/long_short_memory.jl new file mode 100644 index 0000000..2176b25 --- /dev/null +++ b/test/long_short_memory.jl @@ -0,0 +1,62 @@ +using Plots +using DrWatson +using Revise +using SpikingNeuralNetworks +SNN.@load_units; +using SNNUtils + +## +# Define the network +network = let + # Number of neurons in the network + N = 1000 + # Create dendrites for each neuron + E = SNN.AdEx(N = N, param = SNN.AdExParameter(Vr = -60mV)) + # Define interneurons + I = SNN.IF(; N = N ÷ 4, param = SNN.IFParameter(τm = 20ms, El = -50mV)) + # Define synaptic interactions between neurons and interneurons + E_to_I = SNN.SpikingSynapse(E, I, :ge, p = 0.2, μ = 3.0) + LSSPparam = SNN.LSSPParameter(;long=SNN.vSTDPParameter(), short=SNN.STPParameter()) + E_to_E = SNN.SpikingSynapse(E, E, :ge, p = 0.2, μ = 0.5, param = LSSPparam) + I_to_I = SNN.SpikingSynapse(I, I, :gi, p = 0.2, μ = 4.0) + I_to_E = SNN.SpikingSynapse( + I, + E, + :gi, + p = 0.2, + μ = 1, + param = SNN.iSTDPParameterRate(r = 4Hz), + ) + norm = SNN.SynapseNormalization(E, [E_to_E], param = SNN.AdditiveNorm(τ = 30ms)) + # Store neurons and synapses into a dictionary + pop = dict2ntuple(@strdict E I) + syn = dict2ntuple(@strdict I_to_E E_to_I E_to_E norm I_to_I) + # Return the network as a tuple + (pop = pop, syn = syn) +end + +noise = SNN.PoissonStimulus(network.pop.E, :ge, param=2.8kHz, cells=:ALL) +model = SNN.merge_models(network=network, noise=noise) +SNN.monitor([model.pop...], [:fire, :v]) +# SNN.monitor([network.syn.E_to_E], [:x]) +# SNN.monitor([network.syn.E_to_E], [:v]) +# SNN.monitor([network.syn.E_to_E], [:ρ]) + +simtime = SNN.Time() +train!(model=model, duration = 5000ms, time = simtime, dt = 0.1f0, pbar = true) + + +SNN.vecplot(model.pop.network_E, [:v], neurons = 1:10, r = 800ms:4999ms) + +spiketimes = SNN.spiketimes(model.pop.network_E) +SNN.raster([model.pop...], [0s, 5s]) +rates, intervals = SNN.firing_rate(network.pop.E, interval=0:10:5s, τ=100ms) +plot(intervals,rates[1]) +plot(mean(rates, dims = 1)[1, :], legend = false) +## + +network.pop.E.records + +## The simulation achieves > 1.0 iteration per second on my M1 Pro machine. + + diff --git a/test/multilayer_potjans-diesmann2015.jl b/test/multilayer_potjans-diesmann2015.jl new file mode 100644 index 0000000..8e8a115 --- /dev/null +++ b/test/multilayer_potjans-diesmann2015.jl @@ -0,0 +1,306 @@ +using Revise +using DrWatson +using SpikingNeuralNetworks +SNN.@load_units +using SNNPlots +#gr() +using Statistics, Random +using StatsBase + +using Plots +#using JLD +using ProgressMeter + +using SpikingNeuralNetworks +using SNNUtils +using JLD2 +using Distributions +using Test + +using DataFrames +using Plots + +using Random + +function potjans_neurons(scale=1.0) + ccu = Dict( + :E23 => trunc(Int32, 20683 * scale), + :E4 => trunc(Int32, 21915 * scale), + :E5 => trunc(Int32, 4850 * scale), + :E6 => trunc(Int32, 14395 * scale), + :I6 => trunc(Int32, 2948 * scale), + :I23 => trunc(Int32, 5834 * scale), + :I5 => trunc(Int32, 1065 * scale), + :I4 => trunc(Int32, 5479 * scale), + :Th => trunc(Int32, 902 * scale) + ) + #potjans2015_param = SNN.IFCurrentParameter(; El = -49mV,τm=10ms,τabs=2ms,Vt=−50mV,Vr=-65mV) + + neurons = Dict{Symbol, SNN.AbstractPopulation}() + for (k, neuron_count) in pairs(ccu) + + neurons[k] = SNN.IF(; N = neuron_count, param = SNN.IFParameter(τm = 20ms, El = -50mV), name=string(k)) + #IFCurrent(N = v, param=potjans2015_param, name=string(k)) + + end + return neurons +end + +""" +Define Potjans parameters for neuron populations and connection probabilities +# adding static synaptic weight parameters +# based on PSPs that are more accurate to PD based on the PyNN +# model +# Descrepency between PyNN model and original paper +# original paper +# w ± δw 87.8 ± 8.8 pA Excitatory synaptic strengths +# PyNN model: syn_weight_dist = Normal(0.15, 0.1) + +# g –4 Relative inhibitory synaptic strength +""" +function potjans_conn(Ne) + + function j_from_name(pre, post) + if occursin("E", String(pre)) && occursin("E", String(post)) + syn_weight_dist = Normal(87.8, 8.8) # Unit is pA + return rand(syn_weight_dist)*10^-6#pA + + elseif occursin("I", String(pre)) && occursin("E", String(post)) + #syn_weight_dist = Normal(0.15, 0.1) + return -4.0#mV + elseif occursin("E", String(pre)) && occursin("I", String(post)) + syn_weight_dist = Normal(87.8, 8.8)*10^-6# # Unit is pA + return rand(syn_weight_dist)mV + + elseif occursin("I", String(pre)) && occursin("I", String(post)) + #syn_weight_dist = Normal(0.15, 0.1) + return -4.0#mV + elseif occursin("Th", String(pre)) && occursin("I", String(post)) + syn_weight_dist = Normal(87.8, 8.8)*10^-6# # Unit is pA + return rand(syn_weight_dist)mV + + elseif occursin("Th", String(pre)) && occursin("E", String(post)) + syn_weight_dist = Normal(87.8, 8.8)*10^-6# # Unit is pA + return rand(syn_weight_dist)mV + + + else + throw(ArgumentError("Invalid pre-post combination: $pre-$post")) + end + end + pre_layer_names = [:E23, :I23, :E4, :I4, :E5, :I5, :E6, :I6, :Th] + post_layer_names = [:E23, :I23, :E4, :I4, :E5, :I5, :E6, :I6] + + conn_probs = Float32[ + 0.1009 0.1689 0.0437 0.0818 0.0323 0.0 0.0076 0.0 0.0 + 0.1346 0.1371 0.0316 0.0515 0.0755 0.0 0.0042 0.0 0.0 + 0.0077 0.0059 0.0497 0.135 0.0067 0.0003 0.0453 0.0 0.0983 + 0.0691 0.0029 0.0794 0.1597 0.0033 0.0 0.1057 0.0 0.0619 + 0.1004 0.0622 0.0505 0.0057 0.0831 0.3726 0.0204 0.0 0.0 + 0.0548 0.0269 0.0257 0.0022 0.06 0.3158 0.0086 0.0 0.0 + 0.0156 0.0066 0.0211 0.0166 0.0572 0.0197 0.0396 0.225 0.0512 + 0.0364 0.001 0.0034 0.0005 0.0277 0.008 0.0658 0.144 0.0196 + ] + + + conn_j = zeros(Float32, size(conn_probs)) + for pre in eachindex(pre_layer_names) + for post in eachindex(post_layer_names) + + conn_j[post, pre] = j_from_name(pre_layer_names[pre], post_layer_names[post]) + + end + end + + return pre_layer_names,post_layer_names, conn_probs, conn_j +end + + +""" +Main function to setup Potjans layer with memory-optimized connectivity +""" +function potjans_layer(scale) + + ## Create the neuron populations + neurons = potjans_neurons(scale) + exc_pop = filter(x -> occursin("E", String(x)), keys(neurons)) + inh_pop = filter(x -> occursin("I", String(x)), keys(neurons)) + Ne = trunc(Int32, sum([neurons[k].N for k in exc_pop])) + Ni = trunc(Int32, sum([neurons[k].N for k in inh_pop])) + pre_layer_names,post_layer_names, conn_probs, conn_j = potjans_conn(Ne) + ## + # Added synaptic delays from distribution informed by PyNN/PD model + ## + # de ± δde 1.5 ± 0.75 ms Excitatory synaptic transmission delays + # di ± δdi 0.8 ± 0.4 ms Inhibitory synaptic transmission delays + delay_dist_exc = Normal(1.5, 0.75) # Unit is ms + delay_dist_inh = Normal(0.8, 0.4) # Unit is ms + connections = Dict() + conn_map_pre_density = Dict() + conn_map_post_density = Dict() + cnt=0 + for i in eachindex(pre_layer_names) + for j in eachindex(post_layer_names) + pre = pre_layer_names[i] + post = post_layer_names[j] + p = conn_probs[j, i] + @show(p) + J = conn_j[j, i] + sym = J>=0 ? :ge : :gi + μ = J + #@show(μ) + if J>=0 + s = SNN.SpikingSynapse(neurons[pre], neurons[post], sym; μ = μ, p=p, σ=0, delay_dist=delay_dist_exc) + + else + s = SNN.SpikingSynapse(neurons[pre], neurons[post], sym; μ = μ, p=p, σ=0, delay_dist=delay_dist_inh) + end + keyed = Symbol(string(pre,post)) + if haskey(conn_map_pre_density, keyed) + conn_map_pre_density[keyed] += length(s.I) + conn_map_post_density[keyed] += length(s.J) + + else + conn_map_pre_density[keyed] = length(s.I) # or set to some default value + conn_map_post_density[keyed] = length(s.J) + end + conn_map_pre_density[Symbol(string(pre,post))] += length(s.I) + connections[Symbol(string(i,"_",pre,"_",j,"_",post))] = s + end + end + + # Name L2/3e L2/3i L4e L4i L5e L5i L6e L6i Th + # Population size, N 20 683 5834 21 915 5479 4850 1065 14 395 2948 902 + # External inputs, kext (reference) 1600 1500 2100 1900 2000 1900 2900 2100 n/a + # External inputs, kext (layer independent) 2000 1850 2000 1850 2000 1850 2000 1850 n/a + # Background rate, νbg 8 Hz + νe = 8Hz # background stimulation + stimuli = Dict() + μ=0.0525 + peak_rate = 2kHz + stim_parameters = Dict(:decay=>1ms, :peak=>peak_rate, :start=>peak_rate) + + param = PSParam(rate=attack_decay, + variables=variables) + SNN.PoissonStimulus(model.pop.E, :ge, μ=1pF, cells=assembly.cells, param=param, name=string(assembly.name)) + + SE23 = SNN.PoissonStimulus(neurons[:E23], :ge, nothing; cells=[i for i in range(1,neurons[:E23].N)],μ=μ, param = νe,N=1600) + stimuli[Symbol(string("Poisson_SE23", :E23))] = SE23 + SI23 = SNN.PoissonStimulus(neurons[:I23], :ge, nothing; cells=[i for i in range(1,neurons[:I23].N)],μ=μ, param = νe,N=1500) + stimuli[Symbol(string("Poisson_SI23", :I23))] = SI23 + SE4 = SNN.PoissonStimulus(neurons[:E4], :ge, nothing; cells=[i for i in range(1,neurons[:E4].N)],μ=μ, param = νe,N=2100) + stimuli[Symbol(string("Poisson_SE4", :E4))] = SE4 + SI4 = SNN.PoissonStimulus(neurons[:I4], :ge, nothing; cells=[i for i in range(1,neurons[:I4].N)],μ=μ, param = νe,N=1900) + stimuli[Symbol(string("Poisson_SI4", :I4))] = SI4 + SE5 = SNN.PoissonStimulus(neurons[:E5], :ge, nothing; cells=[i for i in range(1,neurons[:E5].N)],μ=μ, param = νe,N=2000) + stimuli[Symbol(string("Poisson_SE5", :E5))] = SE5 + SI5 = SNN.PoissonStimulus(neurons[:I5], :ge, nothing; cells=[i for i in range(1,neurons[:I5].N)],μ=μ, param = νe,N=1900) + stimuli[Symbol(string("Poisson_SI5", :I5))] = SI5 + SE6 = SNN.PoissonStimulus(neurons[:E6], :ge, nothing; cells=[i for i in range(1,neurons[:E6].N)],μ=μ, param = νe,N=2900,name=string("SE6")) + stimuli[Symbol(string("Poisson_SE6", SE6))] = SE6 + SI6 = SNN.PoissonStimulus(neurons[:I6], :ge, nothing; cells=[i for i in range(1,neurons[:I6].N)],μ=μ,param = νe,N=2100,name=string("I6")) + stimuli[Symbol(string("Poisson_SI6", :I6))] = SI6 + for (ind,pop) in enumerate(exc_pop) + #νe = 4Hz # background stimulation + post = neurons[pop] + #@show(μ=1/100) + s = SNN.PoissonStimulus(post, :ge; param = νe, cells=:ALL, μ=μ, name="PoissonE_$(post.name)") + stimuli[Symbol(string("PoissonE_", pop))] = s + end + for (ind,pop) in enumerate(inh_pop) + #νe = 4Hz # background stimulation + post = neurons[pop] + s = SNN.PoissonStimulus(post, :ge; param = νe, cells=:ALL, μ=μ, name="PoissonI_$(post.name)") + stimuli[Symbol(string("PoissonI_", pop))] = s + end + + return merge_models(neurons,connections, stimuli),neurons,connections,stimuli,conn_map_pre_density,conn_map_post_density,pre_layer_names,post_layer_names + +end + +if !isfile("data_connectome.jld") + model,neurons,connections,stimuli,conn_map_pre_density,conn_map_post_density,pre_layer_names,post_layer_names = potjans_layer(0.125) + pre_synapses = [] + post_synapses = [] + layer_names = [] + + for (k,s) in pairs(connections) + push!(pre_synapses,s.I) + push!(post_synapses,s.J) + push!(layer_names,k) + + end + @save "data_connectome.jld" pre_synapses post_synapses layer_names conn_map_pre_density conn_map_post_density pre_layer_names post_layer_names neurons + #@save "just_model.jld" model + +else + @load "data_connectome.jld" pre_synapses post_synapses layer_names conn_map_pre_density conn_map_post_density pre_layer_names post_layer_names neurons + #@load "just_model.jld" model + +end + +model,neurons,connections,stimuli,conn_map_pre_density,conn_map_post_density,pre_layer_names,post_layer_names = potjans_layer(0.125) +@save "data_connectome.jld" pre_synapses post_synapses layer_names conn_map_pre_density conn_map_post_density pre_layer_names post_layer_names neurons +#@save "just_model.jld" model + +function pre_process_for_sankey(pre_layer_names,post_layer_names,conn_map_pre_density) + labels = [] + connections = [] + for i in eachindex(pre_layer_names) + for j in eachindex(post_layer_names) + pre = pre_layer_names[i] + post = post_layer_names[j] + push!(labels,string(pre)) + push!(connections,(i,j,conn_map_pre_density[Symbol(string(pre,post))])) + end + end + @save "sankey_data.jld" pre_layer_names post_layer_names connections labels +end + +#include("sankey_only.jl") +#sankey_applied(true) +sankeyed=false +if sankeyed + sources = pre_synapses[1] # Indices of source nodes + targets = post_synapses[1] # Indices of target nodes + values_ = [1 for i in range(1,length(sources))] # Values or densities for ribbons + layer_names = [layer_names[ind] for (ind,post) in enumerate(post_synapses) if length(post) != 0] + pre_synapses = [pre for pre in pre_synapses if length(pre) != 0] + post_synapses = [post for post in post_synapses if length(post) != 0] + + pre_process_for_sankey(pre_layer_names,post_layer_names,conn_map_pre_density) + + p = SNNPlots.doparallelCoords_optimized(pre_synapses,post_synapses;figure_name="parallelCoordinatesPlot.png") +end + #p = SNNPlots.doparallelCoords(pre_synapses,post_synapses) +#Plots.plot(p) +#Plots.savefig("parallelCoordinatesPlot.png") + + +#@snn_kw struct PoissonStimulusInterval{R=Float32, } <: PoissonStimulusParameter +# rate::Vector{R} = fill(0.0, N) +# intervals::Vector{Vector{R}} +#end + +param = PoissonStimulusInterval(rate,[[]]) +TotalNNeurons = sum([v.N for v in values(neurons)]) +duration = 1500ms +SNN.monitor([model.pop...], [:fire]) +SNN.monitor([model.pop...], [:v], sr=400Hz) +SNN.sim!(model=model; duration = duration, pbar = true, dt = 0.125) +gui(SNN.raster(model.pop, [10s, 15s])) +gui(SNN.raster(model.pop, [14.5s, 15s])) +gui(SNN.raster(model.pop.E23, [14.5s, 15s])) + +#gui(plot(interval, mean.(frE), label=hcat(names_pop...), xlabel="Time [ms]", ylabel="Firing rate [Hz]", legend=:topleft)) +gui(vecplot(model.pop.E23, :v, neurons =1:200, r=0s:0.0125s,label="soma")) +gui(vecplot(model.pop.Th, :v, neurons =1:50, r=0s:0.125s,label="soma")) +gui(vecplot(model.pop.I4, :v, neurons =1:200, r=0s:0.125s,label="soma")) +gui(vecplot(model.pop.E4, :v, neurons =1:200, r=0s:0.125s,label="soma")) +#gui(vecplot(model.pop.E23, :v, neurons =1:200, r=0s:0.0125s,label="soma")) +layer_names, conn_probs, conn_j = potjans_conn(4000) +pj = heatmap(conn_j, xticks=(1:8,layer_names), yticks=(1:8,layer_names), aspect_ratio=1, color=:bluesreds, title="Synaptic weights", xlabel="Presynaptic", ylabel="Postsynaptic", size=(500,500), clims=(-maximum(abs.(conn_j)), maximum(abs.(conn_j)))) +gui(pj) +Trange = 0s:0.25:1.5s +frE, interval, names_pop = SNN.firing_rate(model.pop, interval = Trange) +gui(scatter([i for i in range(1,9)],mean.(frE), labels=names_pop, xlabel="Time [ms]", ylabel="Firing rate [Hz]", legend=:topright)) diff --git a/test/multiple_competing_populations.jl b/test/multiple_competing_populations.jl new file mode 100644 index 0000000..09a14c0 --- /dev/null +++ b/test/multiple_competing_populations.jl @@ -0,0 +1,126 @@ +using Revise +using SpikingNeuralNetworks +using DrWatson +SNN.@load_units; +using SNNUtils +using Plots +using Statistics + +# Define each of the network recurrent assemblies +function define_network(N, name) + # Number of neurons in the network + N = N + # Create dendrites for each neuron + E = SNN.AdEx(N = N, param = SNN.AdExParameter(Vr = -55mV, At = 1mV, b=0, a=0), name="Exc_$name") + # Define interneurons + I = SNN.IF(; N = N ÷ 4, param = SNN.IFParameter(τm = 20ms, El = -50mV), name="Inh_$name") + # Define synaptic interactions between neurons and interneurons + E_to_I = SNN.SpikingSynapse(E, I, :he, p = 0.2, μ = 1.0) + E_to_E = SNN.SpikingSynapse(E, E, :he, p = 0.2, μ = 0.5)#, param = SNN.vSTDPParameter()) + I_to_I = SNN.SpikingSynapse(I, I, :hi, p = 0.2, μ = 1.0) + I_to_E = SNN.SpikingSynapse( + I, + E, + :hi, + p = 0.2, + μ = 1, + param = SNN.iSTDPParameterRate(r = 4Hz, η=0.02), + ) + norm = SNN.SynapseNormalization(E, [E_to_E], param = SNN.AdditiveNorm(τ = 10ms)) + + # Store neurons and synapses into a dictionary + pop = SNN.@symdict E I + syn = SNN.@symdict I_to_E E_to_I E_to_E norm I_to_I + noise = SNN.PoissonStimulus(E, :he, param=4.5kHz, cells=:ALL) + # Return the network as a tuple + SNN.monitor([E, I], [:fire]) + SNN.merge_models(pop, syn, noise=noise, silent=true) +end + +n_assemblies = 4 +## Instantiate the network assemblies and local inhibitory populations +subnets = Dict(Symbol("sub_$n") => define_network(400, n) for n = 1:n_assemblies) +# Add noise to each assembly +# Create synaptic connections between the assemblies and the lateral inhibitory populations +syns = Dict{Symbol,Any}() +for i in eachindex(subnets) + for j in eachindex(subnets) + i == j && continue + push!( + syns, + Symbol("$(i)E_to_$(j)I_lateral") => SNN.SpikingSynapse( + subnets[i].pop.E, + subnets[j].pop.I, + :he, + p = 0.2, + μ = 1.25, + ), + ) + end +end + +## Merge the models and run the simulation, the merge_models function will return a model object (syn=..., pop=...); the function has strong type checking, see the documentation. +network = SNN.merge_models(subnets, syns) + +# Define a time object to keep track of the simulation time, the time object will be passed to the train! function, otherwise the simulation will not create one on the fly. +time_keeper = SNN.Time() +train!(model = network, duration = 5000ms, time = time_keeper, pbar = true, dt = 0.125) + +time_keeper = SNN.Time() +SNN.clear_records([network.pop...]) +train!(model = network, duration = 15000ms, time = time_keeper, pbar = true, dt = 0.125) + +# Plot the raster plot of the network +SNN.raster([network.pop...], [1s, 15s]) +## + +# define the time interval for the analysis +interval = 0:20:SNN.get_time(time_keeper) + +# select only excitatory populations +exc_populations = SNN.filter_populations(populations, :E) +exc_populations.sub_1_E + +# get the spiketimes of the excitatory populations and the indices of each population +spiketimes = SNN.spiketimes(exc_populations) +indices = SNN.population_indices(populations, :E) + +# calculate the firing rate of each excitatory population +rates = map(eachindex(indices)) do i + rates, intervals = + SNN.firing_rate(exc_populations, interval = interval, pop = indices[i], τ = 50) + mean_rate = mean(rates) +end + +## Plot the firing rate of each assembly and the correlation matrix +p1 = plot() +for i in eachindex(rates) + plot!( + interval, + rates[i], + label = "Assembly $i", + xlabel = "Time (ms)", + ylabel = "Firing rate (Hz)", + xlims = (10_000, 15_000), + legend = :topleft, + ) +end +plot!() + +cor_mat = zeros(length(rates), length(rates)) +for i in eachindex(rates) + for j in eachindex(rates) + cor_mat[i, j] = cor(rates[i], rates[j]) + end +end +p2 = heatmap( + cor_mat, + c = :bluesreds, + clims = (-1, 1), + xlabel = "Assembly", + ylabel = "Assembly", + title = "Correlation matrix", + xticks = 1:3, + yticks = 1:3, +) +plot(p1, p2, layout = (2, 1), size = (600, 800), margin = 5Plots.mm) diff --git a/test/network_soma.jl b/test/network_soma.jl new file mode 100644 index 0000000..e165db1 --- /dev/null +++ b/test/network_soma.jl @@ -0,0 +1,72 @@ +using Plots +using DrWatson +using Revise +using SpikingNeuralNetworks +SNN.@load_units; +using SNNUtils +using Statistics + +# Define the network +network = let + # Number of neurons in the network + N = 2000 + # Create dendrites for each neuron + E = SNN.AdEx(N = N, param = SNN.AdExParameter(Vr = -60mV)) + # Define interneurons + I = SNN.IF(; N = N ÷ 4, param = SNN.IFParameter(τm = 20ms, El = -50mV)) + # Define synaptic interactions between neurons and interneurons + E_to_I = SNN.SpikingSynapse(E, I, :ge, p = 0.2, μ = 1.0) + E_to_E = SNN.SpikingSynapse(E, E, :ge, p = 0.2, μ = 0.5, param = SNN.vSTDPParameter()) + I_to_I = SNN.SpikingSynapse(I, I, :gi, p = 0.2, μ = 4.0) + I_to_E = SNN.SpikingSynapse( + I, + E, + :gi, + p = 0.2, + μ = 1, + param = SNN.iSTDPParameterRate(r = 4Hz), + ) + norm = SNN.SynapseNormalization(E, [E_to_E], param = SNN.AdditiveNorm(τ = 30ms)) + + # Store neurons and synapses into a dictionary + pop = dict2ntuple(@strdict E I) + syn = dict2ntuple(@strdict I_to_E E_to_I E_to_E norm I_to_I) + # Return the network as a tuple + # Return the network as a tuple + (pop = pop, syn = syn) +end + +# Create background for the network simulation +noise = SNN.PoissonStimulus(network.pop.E, :ge, param=2.8kHz, cells=:ALL) + +# Combine all +model_soma = SNN.merge_models(network, noise=noise) + +# +@info "Initializing network" +simtime = SNN.Time() +SNN.monitor([network.pop...], [:fire]) +SNN.monitor([network.pop...], [:v], sr=200Hz) + +train!(model=model_soma, duration = 2500ms, time = simtime, dt = 0.125f0, pbar = true) +## +plots = map(1:10) do i + st = spiketimes(model.pop.E)[i] + τ = 400ms + sr = 10ms + ac = SNN.autocorrelogram(st, τ=τ) + histogram(ac,bins=-τ:sr:τ) +end +plot(plots..., layout = (2,5), size = (800, 400), legend = false) + +## +SNN.raster(network.pop, [11s, 15s]) +SNN.vecplot(network.pop.E, :ge, neurons = 1, r = 800ms:24s) +SNN.vecplot(network.pop.E, :v, neurons = 1:1, r=10s:0.125:15s) +rates, intervals = SNN.firing_rate(network.pop.E, interval=0:10:15s, τ=20ms) +plot(intervals, mean(rates,dims=1)[1,:], xlabel="Time [s]", ylabel="Firing rate [Hz]", legend=false) +## + +v, r = SNN.interpolated_record(network.pop.E, :v) + +r \ No newline at end of file diff --git a/test/plasticity_test.jl b/test/plasticity_test.jl new file mode 100644 index 0000000..627ef2f --- /dev/null +++ b/test/plasticity_test.jl @@ -0,0 +1,75 @@ +using Plots +using DrWatson +using Revise +using SpikingNeuralNetworks +SNN.@load_units; +using SNNUtils + +# Define the network +network = let + # Number of neurons in the network + N = 1000 + # Create dendrites for each neuron + E = SNN.AdEx(N = N, param = SNN.AdExParameter(Vr = -60mV)) + # Define interneurons + I = SNN.IF(; N = N ÷ 4, param = SNN.IFParameter(τm = 20ms, El = -50mV)) + # Define synaptic interactions between neurons and interneurons + E_to_I = SNN.SpikingSynapse(E, I, :ge, p = 0.2, μ = 3.0) + E_to_E = SNN.SpikingSynapse(E, E, :ge, p = 0.2, μ = 0.5)#, param = SNN.vSTDPParameter()) + I_to_I = SNN.SpikingSynapse(I, I, :gi, p = 0.2, μ = 4.0) + I_to_E = SNN.SpikingSynapse( + I, + E, + :gi, + p = 0.2, + μ = 1, + param = SNN.iSTDPParameterRate(r = 4Hz), + ) + norm = SNN.SynapseNormalization(E, [E_to_E], param = SNN.AdditiveNorm(τ = 30ms)) + + # Store neurons and synapses into a dictionary + pop = dict2ntuple(@strdict E I) + syn = dict2ntuple(@strdict I_to_E E_to_I E_to_E norm I_to_I) + # Return the network as a tuple + # Return the network as a tuple + (pop = pop, syn = syn) +end + +# Create background for the network simulation +noise = SNN.PoissonStimulus(network.pop.E, :ge, param=1.8kHz, cells=:ALL) +cellA = 23 +cellB = 58 +W = zeros(network.pop.E.N, network.pop.E.N) +W[cellB, cellA] = 5 + +measured_syn = SNN.SpikingSynapse( + network.pop.E, + network.pop.E, + :ge, + w = W, + param = SNN.vSTDPParameter(), +) + +## Combine all +model = SNN.merge_models(network, noise=noise, measured_syn=measured_syn) + + +@info "Initializing network" +simtime = SNN.Time() +SNN.monitor([network.pop...], [:fire]) + +train!(model=model, duration = 5000ms, time = simtime, dt = 0.125f0, pbar = true) + + + +spiketimes = SNN.spiketimes(network.pop.E) +SNN.raster([network.pop...], [1s, 2s]) +# SNN.vecplot(network.pop.E, [:ge], neurons = 1:10, r = 800ms:4999ms) +rates = SNN.firing_rate(network.pop.E, 100ms) +# plot(rates[1:100,1:end]', legend=false) +plot(mean(rates, dims = 1)[1, :], legend = false) +## + +network.pop.E.records + +## The simulation achieves > 1.0 iteration per second on my M1 Pro machine. diff --git a/test/potjans_nmnist.jl b/test/potjans_nmnist.jl new file mode 100644 index 0000000..82f63a6 --- /dev/null +++ b/test/potjans_nmnist.jl @@ -0,0 +1,385 @@ +using Revise +using DrWatson +using SpikingNeuralNetworks +SNN.@load_units +import SpikingNeuralNetworks: AdExParameter, IFParameter +using Statistics, Random +using Plots +#using SparseArrays +#using ProgressMeter +using Plots +using SpikingNeuralNetworks +#using SNNUtils +using JLD2 +using Distributions +include("Tonic_NMNIST_Stimulus.jl") +using .Tonic_NMNIST_Stimulus + +function potjans_neurons(scale=1.0) + ccu = Dict( + :E23 => trunc(Int32, 20683 * scale), + :E4 => trunc(Int32, 21915 * scale), + :E5 => trunc(Int32, 4850 * scale), + :E6 => trunc(Int32, 14395 * scale), + :I6 => trunc(Int32, 2948 * scale), + :I23 => trunc(Int32, 5834 * scale), + :I5 => trunc(Int32, 1065 * scale), + :I4 => trunc(Int32, 5479 * scale), + :Th => trunc(Int32, 902 * scale) + ) + potjans2015_param = SNN.IFCurrentParameter(; El = -49mV,τm=10ms,τabs=2ms,Vt=−50mV,Vr=-65mV) + + neurons = Dict{Symbol, SNN.AbstractPopulation}() + for (k, v) in ccu + + neurons[k] = IFCurrent(N = v, param=potjans2015_param, name=string(k)) + + end + return neurons +end + +""" +Define Potjans parameters for neuron populations and connection probabilities +""" +function potjans_conn(Ne) + # adding static synaptic weight parameters + # based on PSPs that are more accurate to PD based on the PyNN + # model + # Descrepency between PyNN model and original paper + # original paper + # w ± δw 87.8 ± 8.8 pA Excitatory synaptic strengths + # PyNN model: syn_weight_dist = Normal(0.15, 0.1) + + # g –4 Relative inhibitory synaptic strength + # + # + function j_from_name(pre, post) + if occursin("E", String(pre)) && occursin("E", String(post)) + syn_weight_dist = Normal(87.8, 8.8) # Unit is pA + return rand(syn_weight_dist)pA + + elseif occursin("I", String(pre)) && occursin("E", String(post)) + #syn_weight_dist = Normal(0.15, 0.1) + return -4.0pA + elseif occursin("E", String(pre)) && occursin("I", String(post)) + syn_weight_dist = Normal(87.8, 8.8) # Unit is pA + return rand(syn_weight_dist)pA + + elseif occursin("I", String(pre)) && occursin("I", String(post)) + #syn_weight_dist = Normal(0.15, 0.1) + return -4.0pA + elseif occursin("Th", String(pre)) && occursin("I", String(post)) + syn_weight_dist = Normal(87.8, 8.8) # Unit is pA + return rand(syn_weight_dist)pA + + elseif occursin("Th", String(pre)) && occursin("E", String(post)) + syn_weight_dist = Normal(87.8, 8.8) # Unit is pA + return rand(syn_weight_dist)pA + + + else + throw(ArgumentError("Invalid pre-post combination: $pre-$post")) + end + end + + + + pre_layer_names = [:E23, :I23, :E4, :I4, :E5, :I5, :E6, :I6, :Th] + post_layer_names = [:E23, :I23, :E4, :I4, :E5, :I5, :E6, :I6] + + + # Replace static matrix with a regular matrix for `conn_probs` + # ! the convention is j_post_pre. This is how the matrices `w` are built. Are you using that when defining the parameters? + # L2/3e L2/3i L4e L4i L5e L5i L6e L6i Th + + """ + from (columns are presynaptic, rows are post synaptic) + L2/3e L2/3i L4e L4i L5e L5i L6e L6i Th +to L2/3e 0.101 0.169 0.044 0.082 0.032 0.0 0.008 0.0 0.0 + L2/3i 0.135 0.137 0.032 0.052 0.075 0.0 0.004 0.0 0.0 + L4e 0.008 0.006 0.050 0.135 0.007 0.0003 0.045 0.0 0.0983 + L4i 0.069 0.003 0.079 0.160 0.003 0.0 0.106 0.0 0.0619 + L5e 0.100 0.062 0.051 0.006 0.083 0.373 0.020 0.0 0.0 + L5i 0.055 0.027 0.026 0.002 0.060 0.316 0.009 0.0 0.0 + L6e 0.016 0.007 0.021 0.017 0.057 0.020 0.040 0.225 0.0512 + L6i 0.036 0.001 0.003 0.001 0.028 0.008 0.066 0.144 0.0196 + + """ + + + + + conn_probs = Float32[ + 0.1009 0.1689 0.0437 0.0818 0.0323 0.0 0.0076 0.0 0.0 + 0.1346 0.1371 0.0316 0.0515 0.0755 0.0 0.0042 0.0 0.0 + 0.0077 0.0059 0.0497 0.135 0.0067 0.0003 0.0453 0.0 0.0983 + 0.0691 0.0029 0.0794 0.1597 0.0033 0.0 0.1057 0.0 0.0619 + 0.1004 0.0622 0.0505 0.0057 0.0831 0.3726 0.0204 0.0 0.0 + 0.0548 0.0269 0.0257 0.0022 0.06 0.3158 0.0086 0.0 0.0 + 0.0156 0.0066 0.0211 0.0166 0.0572 0.0197 0.0396 0.225 0.0512 + 0.0364 0.001 0.0034 0.0005 0.0277 0.008 0.0658 0.144 0.0196 + ] + + + conn_j = zeros(Float32, size(conn_probs)) + for pre in eachindex(pre_layer_names) + for post in eachindex(post_layer_names) + + conn_j[post, pre] = j_from_name(pre_layer_names[pre], post_layer_names[post]) + + end + end + + return pre_layer_names,post_layer_names, conn_probs, conn_j +end + + +""" +Main function to setup Potjans layer with memory-optimized connectivity +""" +function potjans_layer(scale) + + ## Create the neuron populations + neurons = potjans_neurons(scale) + exc_pop = filter(x -> occursin("E", String(x)), keys(neurons)) + inh_pop = filter(x -> occursin("I", String(x)), keys(neurons)) + Ne = trunc(Int32, sum([neurons[k].N for k in exc_pop])) + Ni = trunc(Int32, sum([neurons[k].N for k in inh_pop])) + pre_layer_names,post_layer_names, conn_probs, conn_j = potjans_conn(Ne) + ## + # Added synaptic delays from distribution informed by PyNN/PD model + ## + # de ± δde 1.5 ± 0.75 ms Excitatory synaptic transmission delays + # di ± δdi 0.8 ± 0.4 ms Inhibitory synaptic transmission delays + + delay_dist_exc = Normal(1.5, 0.75) # Unit is ms + delay_dist_inh = Normal(0.8, 0.4) # Unit is ms + + ## Create the synaptic connections based on the connection probabilities and synaptic weights assigned to each pre-post pair + connections = Dict() + stdp_param = STDPParameter(A_pre =5e-1, + A_post=-5e-1, + τpre =20ms, + τpost =15ms) + + for i in eachindex(pre_layer_names) + for j in eachindex(post_layer_names) + pre = pre_layer_names[i] + post = post_layer_names[j] + p = conn_probs[j, i] + J = conn_j[j, i] + sym = J>=0 ? :ge : :gi + μ = abs(J) + if J>=0 + s = SNN.SpikingSynapse(neurons[pre], neurons[post], sym; μ = μ, p=p, σ=0,param=stdp_param, delay_dist=delay_dist_exc) + else + s = SNN.SpikingSynapse(neurons[pre], neurons[post], sym; μ = -μ, p=p, σ=0, delay_dist=delay_dist_inh) + + end + connections[Symbol(string(pre,"_", post))] = s + end + end + + ## Create the Poisson stimulus for each population + ## The order of these rates needs to be carefully checked as I have changed the order of populations from the + # PyNN model to conveniently group excitatory and inhibitory connections. + #full_mean_rates = [0.971, 2.868, 4.746, 5.396, 8.142, 9.078, 0.991, 7.523] + stimuli = Dict() + for (ind,pop) in enumerate(exc_pop) + νe = 8Hz # background stimulation + #full_mean_rates[ind]kHz + post = neurons[pop] + ## + # TODO replace all, with something that better matches the number of poisson spike generator sources and targets from the Potjan's model. + ## + s = SNN.PoissonStimulus(post, :ge; param = νe, cells=:ALL, μ=1.f0, name="PoissonE_$(post.name)") + stimuli[Symbol(string("PoissonE_", pop))] = s + end + for (ind,pop) in enumerate(inh_pop) + νe = 8Hz # background stimulation + #full_mean_rates[ind]kHz + post = neurons[pop] + ## + # TODO replace all, with something that better matches the number of poisson spike generator sources and targets from the Potjan's model. + ## + + s = SNN.PoissonStimulus(post, :ge; param = νe, cells=:ALL, μ=1.f0, name="PoissonE_$(post.name)") + stimuli[Symbol(string("PoissonE_", pop))] = s + end + + return merge_models(neurons,connections, stimuli),neurons,connections,stimuli +end + +#if !isfile("cached_potjans_model.jld2") + model,neurons_,connections_,stimuli_,net_dict = potjans_layer(0.085) + @save "cached_potjans_model.jld2" model neurons_ connections_ stimuli_ net_dict +#else + @load "cached_potjans_model.jld2" model neurons_ connections_ stimuli_ net_dict +#end +#@show(connections_.vals) + +before_learnning_weights = model.syn[1].W +@show(before_learnning_weights) +ΔTs = -100:1:100ms +ΔWs = zeros(Float32, length(ΔTs)) +Threads.@threads for i in eachindex(ΔTs) + ΔT = ΔTs[i] + #spiketime = [2000ms, 2000ms+ΔT] + #neurons = [[1], [2]] + #inputs = SpikeTime(spiketime, neurons) + w = zeros(Float32, 2,2) + w[1, 2] = 1f0 + st = Identity(N=max_neurons(inputs)) + stim = SpikeTimeStimulusIdentity(st, :g, param=inputs) + syn = SpikingSynapse( st, st, nothing, w = w, param = stdp_param) + + model = merge_models(pop=st, stim=stim, syn=syn, silent=true) + SNN.monitor(model.pop..., [:fire]) + SNN.monitor(model.syn..., [:tpre, :tpost]) + train!(model=model, duration=3000ms, dt=0.1ms) + ΔWs[i] = model.syn[1].W[1] - 1 +end + +duration = 15000ms +SNN.monitor([model.pop...], [:fire]) +#SNN.monitor([model.pop...], [:v], sr=200Hz) +SNN.sim!(model=model, duration=3000ms, dt=0.125, pbar = true)#, dt = 0.125) + +#= +# Example data: Spike times and trial start times +s#pike_times = [0.1, 0.15, 0.2, 0.4, 0.6, 0.65, 0.8, 1.0, 1.1, 1.3, 1.5] +t#rial_starts = [0.0, 1.0, 2.0] # Start times of each trial +num_trials = length(trial_starts) + +# Parameters +bin_width = 0.1 # Bin width in seconds +time_window = (0.0, 1.0) # Time window relative to each trial start + +# Create bins +edges = collect(time_window[1]:bin_width:time_window[2]) +num_bins = length(edges) - 1 + +# Initialize a matrix to hold spike counts +spike_counts = zeros(num_trials, num_bins) + +# Bin spikes for each trial +for (i, trial_start) in enumerate(trial_starts) + aligned_spikes = spike_times .- trial_start # Align spikes to trial start + filtered_spikes = aligned_spikes[aligned_spikes .>= time_window[1] .& aligned_spikes .< time_window[2]] + spike_counts[i, :] = histcounts(filtered_spikes, edges) +end + +# Compute average spike activity (optional) +average_spike_activity = mean(spike_counts, dims=1) + +# Plot heatmap +heatmap( + edges[1:end-1] .+ bin_width / 2, # Bin centers + 1:num_trials, # Trial indices + spike_counts, + xlabel="Time (s)", + ylabel="Trial", + color=:viridis, + title="Spike Activity Heatmap" +) +=# +SNN.train!(model=model, duration=3000ms, dt=0.125, pbar = true)#, dt = 0.125) + +#SNN.sim!(model=model; duration = duration, pbar = true, dt = 0.125) +display(SNN.raster(model.pop, [1.0s, 2s])) +savefig("without_stimulus.png") + + + + +after_learnning_weights0 = model.syn[1].W + +training_order = shuffle(0:60000-1) +cached_spikes=[] +for i in 1:40 + push!(cached_spikes,Tonic_NMNIST_Stimulus.getNMNIST(training_order[i])) +end + +""" +Filter for spike packets that belong to provided labels that pertain to single numeral characters. +""" + + +if !isfile("labeled_packets.jld2") + labeled_packets = Tonic_NMNIST_Stimulus.spike_packets_by_labels(cached_spikes) + filter_label = 0 + (population_code,time_and_offset) = labeled_packets[filter_label] + for filter_label in range(0,9) + (population_code_,time_and_offset_) = labeled_packets[filter_label] + append!(population_code,population_code_) + append!(time_and_offset,time_and_offset_) + end + @save "labeled_packets.jld2" population_code time_and_offset +else + @load "labeled_packets.jld2" population_code time_and_offset +end + + + +p1 = Plots.scatter() +scatter!(p1, + time_and_offset, + population_code, + ms = 0.5, # Marker size + ylabel = "Neuron Index" , + xlabel ="Time (ms)", + title = "Spiking Activity with Distinct Characters", + legend=false +) +display(plot(p1)) +savefig("stimulus.png") + +neurons_as_nested_array = [ Vector{Int64}([n]) for n in population_code] +inputs = SpikeTime(time_and_offset,neurons_as_nested_array) + +st = neurons_[:E4] #Identity(N=max_neurons(inputs)) +w = ones(Float32,neurons_[:E4].N,max_neurons(inputs))*15 + + +st = Identity(N=max_neurons(inputs)) +stim = SpikeTimeStimulusIdentity(st, :g, param=inputs) + + +syn = SpikingSynapse( st, neurons_[:E4], nothing, w = w)#, param = stdp_param) +model2 = merge_models(pop=[st,model], stim=[stim,stimuli_], syn=[syn,connections_], silent=false) + +duration = 15000ms +SNN.monitor([model2.pop...], [:fire]) +SNN.monitor([model2.pop...], [:v], sr=200Hz) +SNN.train!(model=model2; duration = duration, pbar = true, dt = 0.125) +#display(SNN.raster(model2.pop, [0s, 15s])) + +after_learnning_weights1 = model.syn[1].W + +@show(mean(before_learnning_weights)) +@show(mean(after_learnning_weights0)) +@show(mean(after_learnning_weights1)) + +#mean(model2.syn[1].W) + +SNN.spiketimes(model.pop[1]) + +#x, y, y0 = SNN._raster(model2.pop.pop_2_E5,[1.95s, 2s])) + +display(SNN.raster(model2.pop, [1.75s, 2s])) + +savefig("with_stimulus.png") + +Trange = 0:10:15s +frE, interval, names_pop = SNN.firing_rate(model2.pop, interval = Trange) +plot(mean.(frE), label=hcat(names_pop...), xlabel="Time [ms]", ylabel="Firing rate [Hz]", legend=:topleft) +savefig("firing_rate.png") + +## + +vecplot(model2.pop.E4, :v, neurons =1, r=0s:15s,label="soma") +savefig("vector_vm_plot.png") +layer_names, conn_probs, conn_j = potjans_conn(4000) +pj = heatmap(conn_j, xticks=(1:8,layer_names), yticks=(1:8,layer_names), aspect_ratio=1, color=:bluesreds, title="Synaptic weights", xlabel="Presynaptic", ylabel="Postsynaptic", size=(500,500), clims=(-maximum(abs.(conn_j)), maximum(abs.(conn_j)))) +pprob=heatmap(conn_probs, xticks=(1:8,layer_names), yticks=(1:8,layer_names), aspect_ratio=1, color=:viridis, title="Connection probability", xlabel="Presynaptic", ylabel="Postsynaptic", size=(500,500)) +plot(pprob, pj, layout=(1,2), size=(1000,500), margin=5Plots.mm) diff --git a/test/runtests.jl b/test/runtests.jl new file mode 100644 index 0000000..f64958b --- /dev/null +++ b/test/runtests.jl @@ -0,0 +1,6 @@ +using SpikingNeuralNetworks +using Test +#SNN.@load_units + +include("multilayer_potjans-diesmann2015.jl") +include("bistable_network_with_stimulus.jl") \ No newline at end of file diff --git a/test/sankey_only.jl b/test/sankey_only.jl new file mode 100644 index 0000000..c00da13 --- /dev/null +++ b/test/sankey_only.jl @@ -0,0 +1,49 @@ +using PlotlyJS +using ElectronDisplay +using JLD2 +""" +Assuming you have already created a connectome and you have formated the date into three seperate dense vectors, which are encoded as the 1,2 and 3rd elements of the List of Lists represented by the loaded variable +"connections" Make a Sankey diagram from the given information about pre synaptic and post synaptic connection densities. +Network layer sources (pre-synaptic densities), network layer targets (post synaptic densities), and ribbon thickness (values). +""" +function sankey_applied(from_jld=true) + if from_jld + @load "sankey_data.jld" connections labels + else + throw("implement connectome creating code here.") + end + if isa(connections, Vector{<:Vector}) + connections = hcat(connections...) + @save "sankey_data.jld" connections + end + sankey_trace = sankey( + arrangement = "snap", + node = attr( + label = labels, + pad = 15, + thickness = 20, + line = attr(color = "black", width = 0.5) + ), + #if isa(connections, Vector{<:Any}) + + link = attr( + source = [i[1] for i in connections], + target = [i[2] for i in connections], + value = [i[3] for i in connections] + ) + + ) + plt = PlotlyJS.plot(sankey_trace) + ElectronDisplay.display(plt) # opens a new Electron window +end +#= +else + #@show(typeof(connections)) + #@show([i[1] for i in connections]) + link = attr( + source = connections[:, 1], # First column for source + target = connections[:, 2], # Second column for target + value = connections[:, 3] # Third column for value + ) +end +=# \ No newline at end of file diff --git a/test/two_competing_populations.jl b/test/two_competing_populations.jl new file mode 100644 index 0000000..5402adf --- /dev/null +++ b/test/two_competing_populations.jl @@ -0,0 +1,71 @@ +using DrWatson +using Revise +using Plots +using SpikingNeuralNetworks +SNN.@load_units; +using SNNUtils + +# Define the network +function define_network(N = 800) + # Number of neurons in the network + N = N + # Create dendrites for each neuron + E = SNN.AdEx(N = N, param = SNN.AdExParameter(Vr = -60mV)) + # Define interneurons + I = SNN.IF(; N = N ÷ 4, param = SNN.IFParameter(τm = 20ms, El = -50mV)) + # Define synaptic interactions between neurons and interneurons + E_to_I = SNN.SpikingSynapse(E, I, :he, p = 0.2, μ = 3.0) + E_to_E = SNN.SpikingSynapse(E, E, :he, p = 0.2, μ = 0.5)#, param = SNN.vSTDPParameter()) + I_to_I = SNN.SpikingSynapse(I, I, :hi, p = 0.2, μ = 1.0) + I_to_E = SNN.SpikingSynapse( + I, + E, + :hi, + p = 0.2, + μ = 10., + param = SNN.iSTDPParameterRate(r = 4Hz), + ) + norm = SNN.SynapseNormalization(E, [E_to_E], param = SNN.AdditiveNorm(τ = 10ms)) + + # Store neurons and synapses into a dictionary + pop = dict2ntuple(@strdict E I) + syn = dict2ntuple(@strdict I_to_E E_to_I E_to_E norm I_to_I) + # Return the network as a tuple + (pop = pop, syn = syn) +end + + +network1 = define_network(800) +network2 = define_network(800) +inter = ( + E1_to_I2 = SNN.SpikingSynapse(network1.pop.E, network2.pop.I, :he, p = 0.2, μ = 4.25), + E2_to_I1 = SNN.SpikingSynapse(network2.pop.E, network1.pop.I, :he, p = 0.2, μ = 4.25) +) + +noise = ( + A = SNN.PoissonStimulus(network1.pop.E, :he, param=4.5kHz, cells=:ALL, μ=1.7f0,), + B = SNN.PoissonStimulus(network2.pop.E, :he, param=4.5kHz, cells=:ALL, μ=1.7f0,), +) + +model = merge_models( noise, inter; n1=network1, n2=network2) + +## @info "Initializing network" +SNN.monitor([model.pop...], [:fire, :v, :he, :ge]) +train!(model = model, duration = 15000ms, pbar = true, dt = 0.125ms) +SNN.raster([model.pop...], [14s, 15s]) + +train!(model = model, duration = 1000ms, pbar = true, dt = 0.125ms) +SNN.vecplot(model.pop.n1_E, [:v, :he, :ge], r=0.01s:0.001:1s, dt=0.125ms) + + +## +using Statistics +rate1, intervals = SNN.firing_rate(no_noise.pop.network1_E, τ = 10ms) +rate2, intervals = SNN.firing_rate(no_noise.pop.network2_E, τ = 10ms) +r1 = mean(rate1) +r2 = mean(rate2) +cor(r1, r2) +plot(r1, label = "Network 1", xlabel = "Time (s)", ylabel = "Firing rate (Hz)") +plot!(r2, label = "Network 2", title = "Correlation: $(cor(r1, r2))") +plot!(xlims = (100, 500)) +## diff --git a/test/working_NMNIST_reading.j b/test/working_NMNIST_reading.j new file mode 100644 index 0000000..ae80b67 --- /dev/null +++ b/test/working_NMNIST_reading.j @@ -0,0 +1,70 @@ +using Plots +using JLD2 +@load "labeled_packets.jld2" population_code time_and_offset + +gui(Plots.scatter(time_and_offset,population_code)) + + +p1 = Plots.scatter() +scatter!(p1, + time_and_offset, + population_code, + ms = 0.5, # Marker size + ylabel = "Neuron Index" , + xlabel ="Time (ms)", + title = "Spiking Activity with Distinct Characters", + legend=false +) +display(plot(p1)) +savefig("stimulus.png") + +neurons_as_nested_array = [ Vector{Int64}([n]) for n in population_code] +inputs = SpikeTime(time_and_offset,neurons_as_nested_array) + +st = neurons_[:E4] #Identity(N=max_neurons(inputs)) +w = ones(Float32,neurons_[:E4].N,max_neurons(inputs))*15 + + +st = Identity(N=max_neurons(inputs)) +stim = SpikeTimeStimulusIdentity(st, :g, param=inputs) + + +syn = SpikingSynapse( st, neurons_[:E4], nothing, w = w)#, param = stdp_param) +model2 = merge_models(pop=[st,model], stim=[stim,stimuli_], syn=[syn,connections_], silent=false) + +duration = 15000ms +SNN.monitor([model2.pop...], [:fire]) +SNN.monitor([model2.pop...], [:v], sr=200Hz) +SNN.train!(model=model2; duration = duration, pbar = true, dt = 0.125) +#display(SNN.raster(model2.pop, [0s, 15s])) + +after_learnning_weights1 = model.syn[1].W + +@show(mean(before_learnning_weights)) +@show(mean(after_learnning_weights0)) +@show(mean(after_learnning_weights1)) + +#mean(model2.syn[1].W) + +SNN.spiketimes(model.pop[1]) + +#x, y, y0 = SNN._raster(model2.pop.pop_2_E5,[1.95s, 2s])) + +display(SNN.raster(model2.pop, [1.75s, 2s])) + +savefig("with_stimulus.png") + +Trange = 0:10:15s +frE, interval, names_pop = SNN.firing_rate(model2.pop, interval = Trange) +plot(mean.(frE), label=hcat(names_pop...), xlabel="Time [ms]", ylabel="Firing rate [Hz]", legend=:topleft) +savefig("firing_rate.png") + +## + +vecplot(model2.pop.E4, :v, neurons =1, r=0s:15s,label="soma") +savefig("vector_vm_plot.png") +layer_names, conn_probs, conn_j = potjans_conn(4000) +pj = heatmap(conn_j, xticks=(1:8,layer_names), yticks=(1:8,layer_names), aspect_ratio=1, color=:bluesreds, title="Synaptic weights", xlabel="Presynaptic", ylabel="Postsynaptic", size=(500,500), clims=(-maximum(abs.(conn_j)), maximum(abs.(conn_j)))) +pprob=heatmap(conn_probs, xticks=(1:8,layer_names), yticks=(1:8,layer_names), aspect_ratio=1, color=:viridis, title="Connection probability", xlabel="Presynaptic", ylabel="Postsynaptic", size=(500,500)) +plot(pprob, pj, layout=(1,2), size=(1000,500), margin=5Plots.mm) + diff --git a/test/working_NMNIST_reading.jl b/test/working_NMNIST_reading.jl new file mode 100644 index 0000000..eb909c9 --- /dev/null +++ b/test/working_NMNIST_reading.jl @@ -0,0 +1,117 @@ +using Plots +using JLD2 +@load "labeled_packets.jld2" population_code time_and_offset + +gui(Plots.scatter(time_and_offset,population_code)) + + +p1 = Plots.scatter() +scatter!(p1, + time_and_offset, + population_code, + ms = 0.5, # Marker size + ylabel = "Neuron Index" , + xlabel ="Time (ms)", + title = "Spiking Activity with Distinct Characters", + legend=false +) +display(plot(p1)) +savefig("stimulus.png") + +neurons_as_nested_array = [ Vector{Int64}([n]) for n in population_code] +#= +inputs = SpikeTime(time_and_offset,neurons_as_nested_array) + +st = neurons_[:E4] #Identity(N=max_neurons(inputs)) +w = ones(Float32,neurons_[:E4].N,max_neurons(inputs))*15 + + +st = Identity(N=max_neurons(inputs)) +stim = SpikeTimeStimulusIdentity(st, :g, param=inputs) + + +syn = SpikingSynapse( st, neurons_[:E4], nothing, w = w)#, param = stdp_param) +model2 = merge_models(pop=[st,model], stim=[stim,stimuli_], syn=[syn,connections_], silent=false) + +duration = 15000ms +SNN.monitor([model2.pop...], [:fire]) +SNN.monitor([model2.pop...], [:v], sr=200Hz) +SNN.train!(model=model2; duration = duration, pbar = true, dt = 0.125) +#display(SNN.raster(model2.pop, [0s, 15s])) + +after_learnning_weights1 = model.syn[1].W + +@show(mean(before_learnning_weights)) +@show(mean(after_learnning_weights0)) +@show(mean(after_learnning_weights1)) + +#mean(model2.syn[1].W) + +SNN.spiketimes(model.pop[1]) + +#x, y, y0 = SNN._raster(model2.pop.pop_2_E5,[1.95s, 2s])) + +display(SNN.raster(model2.pop, [1.75s, 2s])) + +savefig("with_stimulus.png") + +Trange = 0:10:15s +frE, interval, names_pop = SNN.firing_rate(model2.pop, interval = Trange) +plot(mean.(frE), label=hcat(names_pop...), xlabel="Time [ms]", ylabel="Firing rate [Hz]", legend=:topleft) +savefig("firing_rate.png") + +## + +vecplot(model2.pop.E4, :v, neurons =1, r=0s:15s,label="soma") +savefig("vector_vm_plot.png") +layer_names, conn_probs, conn_j = potjans_conn(4000) +pj = heatmap(conn_j, xticks=(1:8,layer_names), yticks=(1:8,layer_names), aspect_ratio=1, color=:bluesreds, title="Synaptic weights", xlabel="Presynaptic", ylabel="Postsynaptic", size=(500,500), clims=(-maximum(abs.(conn_j)), maximum(abs.(conn_j)))) +pprob=heatmap(conn_probs, xticks=(1:8,layer_names), yticks=(1:8,layer_names), aspect_ratio=1, color=:viridis, title="Connection probability", xlabel="Presynaptic", ylabel="Postsynaptic", size=(500,500)) +plot(pprob, pj, layout=(1,2), size=(1000,500), margin=5Plots.mm) + +using Plots +using SpikingNeuralNetworks +SNN.@load_units +import SpikingNeuralNetworks: AdExParameter +using Statistics, Random + + + + +## + +# Parameters +T = 10.0s # Total time +# Example spike trains as floating-point times +t_pre = sort(rand(10000) .* T) # Post-synaptic spikes (50 random times in [0, T]) +# t_pre = sort(rand(30) .* T) # Post-synaptic spikes (50 random times in [0, T]) +t_post = t_pre .+rand(length(t_pre))*100 .-50 #.-200 # Pre-synaptic spikes (75 random times in [0, T]) + + + +r, cv = compute_covariance_density(Float32.(t_pre), Float32.(t_post), T, sr=40Hz, τ=300ms) +bar(r, cv, legend=false, fill=true, xlabel="ΔT(ms)", ylabel="C(τ)", size=(500, 300), alpha=0.5, color=:black, margin=5Plots.mm) +plot!(frame=:origin, yticks=:none, ) +## + +# +z = [] +tpre +for x in reverse(2:length(tpre)) + if tpre[x] - tpre[x-1] < 50ms + # if rand() < 0.5 + push!(z, x) + # end + end +end +z +for x in z + popat!(tpre, x) +end + + +taus = autocorrelogram(tpre, τ= 400ms) +histogram(taus, bins=100, legend=false, fill=true, xlabel="ΔT(mAs)", ylabel="Autocorrelogram", size=(500, 300), alpha=0.5, color=:black, margin=5Plots.mm, frame=:origin, yticks=:none) +=# + + diff --git a/tests/bistable_network_with_stimulus.jl b/tests/bistable_network_with_stimulus.jl new file mode 100644 index 0000000..5816b4c --- /dev/null +++ b/tests/bistable_network_with_stimulus.jl @@ -0,0 +1,159 @@ +using DrWatson +using Revise +using SpikingNeuralNetworks +SNN.@load_units; +using SNNUtils +using Plots +using Statistics + +# Define each of the network recurrent assemblies +function define_network(N = 800) + # Number of neurons in the network + N = N + # Create dendrites for each neuron + E = SNN.AdEx(N = N, param = SNN.AdExParameter(Vr = -55mV, At = 0mV, b=0, a=0)) + # Define interneurons + I = SNN.IF(; N = N ÷ 2, param = SNN.IFParameter(τm = 20ms, El = -50mV)) + # Define synaptic interactions between neurons and interneurons + E_to_I = SNN.SpikingSynapse(E, I, :he, p = 0.2, μ = 2.0) + E_to_E = SNN.SpikingSynapse(E, E, :he, p = 0.2, μ = 3.0)#, param = SNN.vSTDPParameter()) + I_to_I = SNN.SpikingSynapse(I, I, :hi, p = 0.2, μ = 1.0) + I_to_E = SNN.SpikingSynapse( + I, + E, + :hi, + p = 0.2, + μ = 1, + param = SNN.iSTDPParameterRate(r = 4Hz), + ) + norm = SNN.SynapseNormalization(E, [E_to_E], param = SNN.AdditiveNorm(τ = 10ms)) + + # Store neurons and synapses into a dictionary + pop = dict2ntuple(@strdict E I) + syn = dict2ntuple(@strdict I_to_E E_to_I E_to_E norm I_to_I) + # Return the network as a tuple + SNN.monitor([E, I], [:fire]) + (pop = pop, syn = syn) +end + +n_assemblies = 2 +## Instantiate the network assemblies and local inhibitory populations +subnets = Dict(Symbol("sub_$n") => define_network(400) for n = 1:n_assemblies) +# Add noise to each assembly +noise = Dict(Symbol("noise_$(i)") => SNN.PoissonStimulus(subnets[i].pop.E, :he, param=2.5kHz, cells=:ALL) for i in eachindex(subnets)) +# Create synaptic connections between the assemblies and the lateral inhibitory populations +syns = Dict{Symbol,Any}() +for i in eachindex(subnets) + for j in eachindex(subnets) + i == j && continue + push!( + syns, + Symbol("lateral_$(i)E_to_$(j)I") => SNN.SpikingSynapse( + subnets[i].pop.E, + subnets[j].pop.I, + :he, + p = 0.2, + μ = 2.25, + ), + ) + end +end + +# select only excitatory populations +input = function (t, param::PSParam) + id::Int = param.variables[:id] + n_assemblies::Int = param.variables[:n_assemblies] + if (t ÷ 1000)%n_assemblies+1 == id + return 5kHz * exp(-(t%1000)/200) + else + return 0 + end +end +# select only excitatory populations + + +stimuli = Dict{Symbol,Any}() +for n in 1:n_assemblies + push!(stimuli, Symbol("stim_$n")=> + SNN.PoissonStimulus( + subnets[Symbol("sub_$n")].pop.E, + :ge, + cells=:ALL, + param=PSParam( + rate=input, + variables = Dict( + :id => n, + :n_assemblies => n_assemblies)))) +end +network = SNN.merge_models(noise, subnets, syns, stimuli) + + +## Merge the models and run the simulation, the merge_models function will return a model object (syn=..., pop=...); the function has strong type checking, see the documentation. +train!(model = network, duration = 5000ms, pbar = true, dt = 0.125) + +## Create a model object with only the populations to run the analysis +SNN.monitor([network.pop...], [:fire]) + +# Define a time object to keep track of the simulation time, the time object will be passed to the train! function, otherwise the simulation will not create one on the fly. +time_keeper = SNN.Time() +train!(model = network, duration = 15000ms, time = time_keeper, pbar = true, dt = 0.125) + + +# Plot the raster plot of the network +SNN.raster([network.pop...], [4s, 15s]) +## + +network.syn.sub_2_I_to_E.W +# define the time interval for the analysis + +exc_populations = SNN.filter_populations(populations, :E) + +# get the spiketimes of the excitatory populations and the indices of each population +spiketimes = SNN.spiketimes(exc_populations) +indices = SNN.population_indices(populations, :E) + +# calculate the firing rate of each excitatory population +interval = 0:5:SNN.get_time(time_keeper) +rates = map(eachindex(indices)) do i + rates, intervals = + SNN.firing_rate(spiketimes, interval = interval, pop = indices[i], τ = 10) + mean_rate = mean(rates) +end + +## Plot the firing rate of each assembly and the correlation matrix +p1 = plot() +for i in eachindex(rates) + plot!( + interval, + rates[i], + label = "Assembly $i", + xlabel = "Time (ms)", + ylabel = "Firing rate (Hz)", + xlims = (4_000, 15_000), + legend = :topleft, + ) +end +plot!() + +cor_mat = zeros(length(rates), length(rates)) +for i in eachindex(rates) + for j in eachindex(rates) + cor_mat[i, j] = cor(rates[i], rates[j]) + end +end +p2 = heatmap( + cor_mat, + c = :bluesreds, + clims = (-1, 1), + xlabel = "Assembly", + ylabel = "Assembly", + title = "Correlation matrix", + xticks = 1:3, + yticks = 1:3, +) +plot(p1, p2, layout = (2, 1), size = (600, 800), margin = 5Plots.mm) + +using Statistics, StatsBase +pcor = plot() +n = 800 +plot(0:5:5*n,autocor(rates[1], 0:n)) diff --git a/tests/long_short_memory.jl b/tests/long_short_memory.jl new file mode 100644 index 0000000..2176b25 --- /dev/null +++ b/tests/long_short_memory.jl @@ -0,0 +1,62 @@ +using Plots +using DrWatson +using Revise +using SpikingNeuralNetworks +SNN.@load_units; +using SNNUtils + +## +# Define the network +network = let + # Number of neurons in the network + N = 1000 + # Create dendrites for each neuron + E = SNN.AdEx(N = N, param = SNN.AdExParameter(Vr = -60mV)) + # Define interneurons + I = SNN.IF(; N = N ÷ 4, param = SNN.IFParameter(τm = 20ms, El = -50mV)) + # Define synaptic interactions between neurons and interneurons + E_to_I = SNN.SpikingSynapse(E, I, :ge, p = 0.2, μ = 3.0) + LSSPparam = SNN.LSSPParameter(;long=SNN.vSTDPParameter(), short=SNN.STPParameter()) + E_to_E = SNN.SpikingSynapse(E, E, :ge, p = 0.2, μ = 0.5, param = LSSPparam) + I_to_I = SNN.SpikingSynapse(I, I, :gi, p = 0.2, μ = 4.0) + I_to_E = SNN.SpikingSynapse( + I, + E, + :gi, + p = 0.2, + μ = 1, + param = SNN.iSTDPParameterRate(r = 4Hz), + ) + norm = SNN.SynapseNormalization(E, [E_to_E], param = SNN.AdditiveNorm(τ = 30ms)) + # Store neurons and synapses into a dictionary + pop = dict2ntuple(@strdict E I) + syn = dict2ntuple(@strdict I_to_E E_to_I E_to_E norm I_to_I) + # Return the network as a tuple + (pop = pop, syn = syn) +end + +noise = SNN.PoissonStimulus(network.pop.E, :ge, param=2.8kHz, cells=:ALL) +model = SNN.merge_models(network=network, noise=noise) +SNN.monitor([model.pop...], [:fire, :v]) +# SNN.monitor([network.syn.E_to_E], [:x]) +# SNN.monitor([network.syn.E_to_E], [:v]) +# SNN.monitor([network.syn.E_to_E], [:ρ]) + +simtime = SNN.Time() +train!(model=model, duration = 5000ms, time = simtime, dt = 0.1f0, pbar = true) + + +SNN.vecplot(model.pop.network_E, [:v], neurons = 1:10, r = 800ms:4999ms) + +spiketimes = SNN.spiketimes(model.pop.network_E) +SNN.raster([model.pop...], [0s, 5s]) +rates, intervals = SNN.firing_rate(network.pop.E, interval=0:10:5s, τ=100ms) +plot(intervals,rates[1]) +plot(mean(rates, dims = 1)[1, :], legend = false) +## + +network.pop.E.records + +## The simulation achieves > 1.0 iteration per second on my M1 Pro machine. + + diff --git a/tests/multiple_competing_populations.jl b/tests/multiple_competing_populations.jl new file mode 100644 index 0000000..09a14c0 --- /dev/null +++ b/tests/multiple_competing_populations.jl @@ -0,0 +1,126 @@ +using Revise +using SpikingNeuralNetworks +using DrWatson +SNN.@load_units; +using SNNUtils +using Plots +using Statistics + +# Define each of the network recurrent assemblies +function define_network(N, name) + # Number of neurons in the network + N = N + # Create dendrites for each neuron + E = SNN.AdEx(N = N, param = SNN.AdExParameter(Vr = -55mV, At = 1mV, b=0, a=0), name="Exc_$name") + # Define interneurons + I = SNN.IF(; N = N ÷ 4, param = SNN.IFParameter(τm = 20ms, El = -50mV), name="Inh_$name") + # Define synaptic interactions between neurons and interneurons + E_to_I = SNN.SpikingSynapse(E, I, :he, p = 0.2, μ = 1.0) + E_to_E = SNN.SpikingSynapse(E, E, :he, p = 0.2, μ = 0.5)#, param = SNN.vSTDPParameter()) + I_to_I = SNN.SpikingSynapse(I, I, :hi, p = 0.2, μ = 1.0) + I_to_E = SNN.SpikingSynapse( + I, + E, + :hi, + p = 0.2, + μ = 1, + param = SNN.iSTDPParameterRate(r = 4Hz, η=0.02), + ) + norm = SNN.SynapseNormalization(E, [E_to_E], param = SNN.AdditiveNorm(τ = 10ms)) + + # Store neurons and synapses into a dictionary + pop = SNN.@symdict E I + syn = SNN.@symdict I_to_E E_to_I E_to_E norm I_to_I + noise = SNN.PoissonStimulus(E, :he, param=4.5kHz, cells=:ALL) + # Return the network as a tuple + SNN.monitor([E, I], [:fire]) + SNN.merge_models(pop, syn, noise=noise, silent=true) +end + +n_assemblies = 4 +## Instantiate the network assemblies and local inhibitory populations +subnets = Dict(Symbol("sub_$n") => define_network(400, n) for n = 1:n_assemblies) +# Add noise to each assembly +# Create synaptic connections between the assemblies and the lateral inhibitory populations +syns = Dict{Symbol,Any}() +for i in eachindex(subnets) + for j in eachindex(subnets) + i == j && continue + push!( + syns, + Symbol("$(i)E_to_$(j)I_lateral") => SNN.SpikingSynapse( + subnets[i].pop.E, + subnets[j].pop.I, + :he, + p = 0.2, + μ = 1.25, + ), + ) + end +end + +## Merge the models and run the simulation, the merge_models function will return a model object (syn=..., pop=...); the function has strong type checking, see the documentation. +network = SNN.merge_models(subnets, syns) + +# Define a time object to keep track of the simulation time, the time object will be passed to the train! function, otherwise the simulation will not create one on the fly. +time_keeper = SNN.Time() +train!(model = network, duration = 5000ms, time = time_keeper, pbar = true, dt = 0.125) + +time_keeper = SNN.Time() +SNN.clear_records([network.pop...]) +train!(model = network, duration = 15000ms, time = time_keeper, pbar = true, dt = 0.125) + +# Plot the raster plot of the network +SNN.raster([network.pop...], [1s, 15s]) +## + +# define the time interval for the analysis +interval = 0:20:SNN.get_time(time_keeper) + +# select only excitatory populations +exc_populations = SNN.filter_populations(populations, :E) +exc_populations.sub_1_E + +# get the spiketimes of the excitatory populations and the indices of each population +spiketimes = SNN.spiketimes(exc_populations) +indices = SNN.population_indices(populations, :E) + +# calculate the firing rate of each excitatory population +rates = map(eachindex(indices)) do i + rates, intervals = + SNN.firing_rate(exc_populations, interval = interval, pop = indices[i], τ = 50) + mean_rate = mean(rates) +end + +## Plot the firing rate of each assembly and the correlation matrix +p1 = plot() +for i in eachindex(rates) + plot!( + interval, + rates[i], + label = "Assembly $i", + xlabel = "Time (ms)", + ylabel = "Firing rate (Hz)", + xlims = (10_000, 15_000), + legend = :topleft, + ) +end +plot!() + +cor_mat = zeros(length(rates), length(rates)) +for i in eachindex(rates) + for j in eachindex(rates) + cor_mat[i, j] = cor(rates[i], rates[j]) + end +end +p2 = heatmap( + cor_mat, + c = :bluesreds, + clims = (-1, 1), + xlabel = "Assembly", + ylabel = "Assembly", + title = "Correlation matrix", + xticks = 1:3, + yticks = 1:3, +) +plot(p1, p2, layout = (2, 1), size = (600, 800), margin = 5Plots.mm) diff --git a/tests/network_soma.jl b/tests/network_soma.jl new file mode 100644 index 0000000..e165db1 --- /dev/null +++ b/tests/network_soma.jl @@ -0,0 +1,72 @@ +using Plots +using DrWatson +using Revise +using SpikingNeuralNetworks +SNN.@load_units; +using SNNUtils +using Statistics + +# Define the network +network = let + # Number of neurons in the network + N = 2000 + # Create dendrites for each neuron + E = SNN.AdEx(N = N, param = SNN.AdExParameter(Vr = -60mV)) + # Define interneurons + I = SNN.IF(; N = N ÷ 4, param = SNN.IFParameter(τm = 20ms, El = -50mV)) + # Define synaptic interactions between neurons and interneurons + E_to_I = SNN.SpikingSynapse(E, I, :ge, p = 0.2, μ = 1.0) + E_to_E = SNN.SpikingSynapse(E, E, :ge, p = 0.2, μ = 0.5, param = SNN.vSTDPParameter()) + I_to_I = SNN.SpikingSynapse(I, I, :gi, p = 0.2, μ = 4.0) + I_to_E = SNN.SpikingSynapse( + I, + E, + :gi, + p = 0.2, + μ = 1, + param = SNN.iSTDPParameterRate(r = 4Hz), + ) + norm = SNN.SynapseNormalization(E, [E_to_E], param = SNN.AdditiveNorm(τ = 30ms)) + + # Store neurons and synapses into a dictionary + pop = dict2ntuple(@strdict E I) + syn = dict2ntuple(@strdict I_to_E E_to_I E_to_E norm I_to_I) + # Return the network as a tuple + # Return the network as a tuple + (pop = pop, syn = syn) +end + +# Create background for the network simulation +noise = SNN.PoissonStimulus(network.pop.E, :ge, param=2.8kHz, cells=:ALL) + +# Combine all +model_soma = SNN.merge_models(network, noise=noise) + +# +@info "Initializing network" +simtime = SNN.Time() +SNN.monitor([network.pop...], [:fire]) +SNN.monitor([network.pop...], [:v], sr=200Hz) + +train!(model=model_soma, duration = 2500ms, time = simtime, dt = 0.125f0, pbar = true) +## +plots = map(1:10) do i + st = spiketimes(model.pop.E)[i] + τ = 400ms + sr = 10ms + ac = SNN.autocorrelogram(st, τ=τ) + histogram(ac,bins=-τ:sr:τ) +end +plot(plots..., layout = (2,5), size = (800, 400), legend = false) + +## +SNN.raster(network.pop, [11s, 15s]) +SNN.vecplot(network.pop.E, :ge, neurons = 1, r = 800ms:24s) +SNN.vecplot(network.pop.E, :v, neurons = 1:1, r=10s:0.125:15s) +rates, intervals = SNN.firing_rate(network.pop.E, interval=0:10:15s, τ=20ms) +plot(intervals, mean(rates,dims=1)[1,:], xlabel="Time [s]", ylabel="Firing rate [Hz]", legend=false) +## + +v, r = SNN.interpolated_record(network.pop.E, :v) + +r \ No newline at end of file diff --git a/tests/plasticity_test.jl b/tests/plasticity_test.jl new file mode 100644 index 0000000..627ef2f --- /dev/null +++ b/tests/plasticity_test.jl @@ -0,0 +1,75 @@ +using Plots +using DrWatson +using Revise +using SpikingNeuralNetworks +SNN.@load_units; +using SNNUtils + +# Define the network +network = let + # Number of neurons in the network + N = 1000 + # Create dendrites for each neuron + E = SNN.AdEx(N = N, param = SNN.AdExParameter(Vr = -60mV)) + # Define interneurons + I = SNN.IF(; N = N ÷ 4, param = SNN.IFParameter(τm = 20ms, El = -50mV)) + # Define synaptic interactions between neurons and interneurons + E_to_I = SNN.SpikingSynapse(E, I, :ge, p = 0.2, μ = 3.0) + E_to_E = SNN.SpikingSynapse(E, E, :ge, p = 0.2, μ = 0.5)#, param = SNN.vSTDPParameter()) + I_to_I = SNN.SpikingSynapse(I, I, :gi, p = 0.2, μ = 4.0) + I_to_E = SNN.SpikingSynapse( + I, + E, + :gi, + p = 0.2, + μ = 1, + param = SNN.iSTDPParameterRate(r = 4Hz), + ) + norm = SNN.SynapseNormalization(E, [E_to_E], param = SNN.AdditiveNorm(τ = 30ms)) + + # Store neurons and synapses into a dictionary + pop = dict2ntuple(@strdict E I) + syn = dict2ntuple(@strdict I_to_E E_to_I E_to_E norm I_to_I) + # Return the network as a tuple + # Return the network as a tuple + (pop = pop, syn = syn) +end + +# Create background for the network simulation +noise = SNN.PoissonStimulus(network.pop.E, :ge, param=1.8kHz, cells=:ALL) +cellA = 23 +cellB = 58 +W = zeros(network.pop.E.N, network.pop.E.N) +W[cellB, cellA] = 5 + +measured_syn = SNN.SpikingSynapse( + network.pop.E, + network.pop.E, + :ge, + w = W, + param = SNN.vSTDPParameter(), +) + +## Combine all +model = SNN.merge_models(network, noise=noise, measured_syn=measured_syn) + + +@info "Initializing network" +simtime = SNN.Time() +SNN.monitor([network.pop...], [:fire]) + +train!(model=model, duration = 5000ms, time = simtime, dt = 0.125f0, pbar = true) + + + +spiketimes = SNN.spiketimes(network.pop.E) +SNN.raster([network.pop...], [1s, 2s]) +# SNN.vecplot(network.pop.E, [:ge], neurons = 1:10, r = 800ms:4999ms) +rates = SNN.firing_rate(network.pop.E, 100ms) +# plot(rates[1:100,1:end]', legend=false) +plot(mean(rates, dims = 1)[1, :], legend = false) +## + +network.pop.E.records + +## The simulation achieves > 1.0 iteration per second on my M1 Pro machine. diff --git a/tests/two_competing_populations.jl b/tests/two_competing_populations.jl new file mode 100644 index 0000000..5402adf --- /dev/null +++ b/tests/two_competing_populations.jl @@ -0,0 +1,71 @@ +using DrWatson +using Revise +using Plots +using SpikingNeuralNetworks +SNN.@load_units; +using SNNUtils + +# Define the network +function define_network(N = 800) + # Number of neurons in the network + N = N + # Create dendrites for each neuron + E = SNN.AdEx(N = N, param = SNN.AdExParameter(Vr = -60mV)) + # Define interneurons + I = SNN.IF(; N = N ÷ 4, param = SNN.IFParameter(τm = 20ms, El = -50mV)) + # Define synaptic interactions between neurons and interneurons + E_to_I = SNN.SpikingSynapse(E, I, :he, p = 0.2, μ = 3.0) + E_to_E = SNN.SpikingSynapse(E, E, :he, p = 0.2, μ = 0.5)#, param = SNN.vSTDPParameter()) + I_to_I = SNN.SpikingSynapse(I, I, :hi, p = 0.2, μ = 1.0) + I_to_E = SNN.SpikingSynapse( + I, + E, + :hi, + p = 0.2, + μ = 10., + param = SNN.iSTDPParameterRate(r = 4Hz), + ) + norm = SNN.SynapseNormalization(E, [E_to_E], param = SNN.AdditiveNorm(τ = 10ms)) + + # Store neurons and synapses into a dictionary + pop = dict2ntuple(@strdict E I) + syn = dict2ntuple(@strdict I_to_E E_to_I E_to_E norm I_to_I) + # Return the network as a tuple + (pop = pop, syn = syn) +end + + +network1 = define_network(800) +network2 = define_network(800) +inter = ( + E1_to_I2 = SNN.SpikingSynapse(network1.pop.E, network2.pop.I, :he, p = 0.2, μ = 4.25), + E2_to_I1 = SNN.SpikingSynapse(network2.pop.E, network1.pop.I, :he, p = 0.2, μ = 4.25) +) + +noise = ( + A = SNN.PoissonStimulus(network1.pop.E, :he, param=4.5kHz, cells=:ALL, μ=1.7f0,), + B = SNN.PoissonStimulus(network2.pop.E, :he, param=4.5kHz, cells=:ALL, μ=1.7f0,), +) + +model = merge_models( noise, inter; n1=network1, n2=network2) + +## @info "Initializing network" +SNN.monitor([model.pop...], [:fire, :v, :he, :ge]) +train!(model = model, duration = 15000ms, pbar = true, dt = 0.125ms) +SNN.raster([model.pop...], [14s, 15s]) + +train!(model = model, duration = 1000ms, pbar = true, dt = 0.125ms) +SNN.vecplot(model.pop.n1_E, [:v, :he, :ge], r=0.01s:0.001:1s, dt=0.125ms) + + +## +using Statistics +rate1, intervals = SNN.firing_rate(no_noise.pop.network1_E, τ = 10ms) +rate2, intervals = SNN.firing_rate(no_noise.pop.network2_E, τ = 10ms) +r1 = mean(rate1) +r2 = mean(rate2) +cor(r1, r2) +plot(r1, label = "Network 1", xlabel = "Time (s)", ylabel = "Firing rate (Hz)") +plot!(r2, label = "Network 2", title = "Correlation: $(cor(r1, r2))") +plot!(xlims = (100, 500)) +##