distributedActivity/runinitial.jl at main · SpikingNetwork/distributedActivity · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
function runinitial(p,w0Index,w0Weights,nc0)

    # copy param
    train_time = 20000.0
    # train_time = copy(p.train_time)
    dt = copy(p.dt)
    Nsteps = Int(train_time / dt) # network param
    # Nsteps = copy(p.Nsteps) # network param
    Ncells = copy(p.Ncells)
    Ne = copy(p.Ne)
    Ni = copy(p.Ni)
    taue = copy(p.taue) # neuron param
    taui = copy(p.taui)
    threshe = copy(p.threshe)
    threshi = copy(p.threshi)
    refrac = copy(p.refrac)
    vre = copy(p.vre)
    muemin = copy(p.muemin) # external input
    muemax = copy(p.muemax)
    muimin = copy(p.muimin)
    muimax = copy(p.muimax)
    tauedecay = copy(p.tauedecay) # synaptic time
    tauidecay = copy(p.tauidecay)
    maxrate = copy(p.maxrate)

    # set up variables
    mu = zeros(Ncells)
    mu[1:Ne] = (muemax-muemin)*rand(Ne) .+ muemin
    mu[(Ne+1):Ncells] = (muimax-muimin)*rand(Ni) .+ muimin

    thresh = zeros(Ncells)
    thresh[1:Ne] .= threshe
    thresh[(1+Ne):Ncells] .= threshi

    tau = zeros(Ncells)
    tau[1:Ne] .= taue
    tau[(1+Ne):Ncells] .= taui

    maxTimes = round(Int,maxrate*train_time/1000)
    times = zeros(Ncells,maxTimes)
    ns = zeros(Int,Ncells)

    forwardInputsE = zeros(Ncells) #summed weight of incoming E spikes
    forwardInputsI = zeros(Ncells)
    forwardInputsEPrev = zeros(Ncells) #as above, for previous timestep
    forwardInputsIPrev = zeros(Ncells)

    xedecay = zeros(Ncells)
    xidecay = zeros(Ncells)

    v = threshe*rand(Ncells) #membrane voltage

    lastSpike = -100.0*ones(Ncells) #time of last spike

    Nexam = 1000
    uavg = zeros(2*Nexam)
    utmp = zeros(Nsteps - Int(1000/p.dt),2*Nexam)
    t = 0.0
    r = zeros(Ncells)
    bias = zeros(Ncells)

    for ti=1:Nsteps
        if mod(ti,Nsteps/100) == 1  #print percent complete
            print("\r",round(Int,100*ti/Nsteps))
        end
        t = dt*ti;
        forwardInputsE .= 0.0;
        forwardInputsI .= 0.0;
        for ci = 1:Ncells
            xedecay[ci] += -dt*xedecay[ci]/tauedecay + forwardInputsEPrev[ci]/tauedecay
            xidecay[ci] += -dt*xidecay[ci]/tauidecay + forwardInputsIPrev[ci]/tauidecay
            synInput = xedecay[ci] + xidecay[ci]

            # excitatory (uavg)
            if ti > Int(1000/p.dt) && ci <= Nexam # 1000 ms
                uavg[ci] += (synInput + mu[ci]) / (Nsteps - Int(1000/p.dt)) # save synInput
            end

            # inhibitory (uavg)
            if ti > Int(1000/p.dt) && ci > Ne && ci <= Ne + Nexam # 1000 ms
                uavg[Nexam+ci-Ne] += (synInput + mu[ci]) / (Nsteps - Int(1000/p.dt)) # save synInput
            end

            # excitatory (ustd)
            if ti > Int(1000/p.dt) && ci <=Nexam
                utmp[ti - Int(1000/p.dt), ci] = synInput
            end

            # inhibitory (ustd)
            if ti > Int(1000/p.dt) && ci > Ne && ci <= Ne + Nexam # 1000 ms
                utmp[ti - Int(1000/p.dt), Nexam+ci-Ne] = synInput
            end

            bias[ci] = mu[ci]

            #not in refractory period
            if t > (lastSpike[ci] + refrac)
                v[ci] += dt*((1/tau[ci])*(bias[ci]-v[ci] + synInput))
                if v[ci] > thresh[ci]  #spike occurred
                    v[ci] = vre
                    lastSpike[ci] = t
                    ns[ci] = ns[ci]+1
                    if ns[ci] <= maxTimes
                        times[ci,ns[ci]] = t
                    end
                    for j = 1:nc0[ci]
                        if w0Weights[j,ci] > 0  #E synapse
                            forwardInputsE[w0Index[j,ci]] += w0Weights[j,ci]
                        elseif w0Weights[j,ci] < 0  #I synapse
                            forwardInputsI[w0Index[j,ci]] += w0Weights[j,ci]
                        end
                    end #end loop over synaptic projections
                end #end if(spike occurred)
            end #end not in refractory period
        end #end loop over neurons

        forwardInputsEPrev = copy(forwardInputsE)
        forwardInputsIPrev = copy(forwardInputsI)

    end #end loop over time
    print("\r")
    println("mean excitatory firing rate: ",mean(1000*ns[1:Ne]/train_time)," Hz")
    println("mean inhibitory firing rate: ",mean(1000*ns[(Ne+1):Ncells]/train_time)," Hz")

    ustd = std(utmp, dims=1)[:]
    ustd_mean = mean(ustd)

    return times, ns, ustd_mean

end