jl_code_for_sigmaInvestigation.jl 4.77 KB
Newer Older
Narges Chinichian's avatar
Narges Chinichian committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

# For explanations check the notebook
#checking
cd(@__DIR__)
using Pkg
Pkg.activate(@__DIR__)

using DelimitedFiles
using LightGraphs
using NetworkDynamics
using OrdinaryDiffEq
using Plots
using BenchmarkTools


G = readdlm(joinpath(@__DIR__, "weights.txt"), '\t', Float64, '\n')
g = SimpleDiGraph(G)
18
N = nv(g)
Narges Chinichian's avatar
Narges Chinichian committed
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
synaptic_weights = G
synaptic_weights ./= max(synaptic_weights...)

#N = 30
#H = G[1:N,1:N]
#H ./= maximum(H)
#g = SimpleDiGraph(H)
#synaptic_weights = H

# FHN parameters
const ϵ = 0.1        # time separation
const a = 0.45        # bifurcation parameter
const b = 0.9         # bifurcation parameter
#const σ = 0.2       # coupling strength

# STP params
const τᴰ = 200.    # depression time
const τᶠ = 1500.    # facilitation time
const u₀ = -2.0        # baseline activity
const y₀ = 0.2        # baseline utilization

# Speed optimization parameters
const α = 1. / 3.
const  = a * ϵ
const  = b * ϵ
const τᴰinv = 1. / τᴰ
const τᶠinv = 1. / τᶠ
const y₀τᶠinv = y₀ * τᶠinv

#we only need random x0 for the first step
49
x0 = Vector(vec([rand(N) .* 4 .- 2 rand(N) .* 4 .- 2 rand(N) .* 0.1 rand(N) .* 0.1 .+ 1]'))
Narges Chinichian's avatar
Narges Chinichian committed
50
51
#println(x0)
for sig in [10,5,2,1.0,0.9,0.8,0.7,0.6,0.5,0.4,0.3,0.2,0.1,0.05,0]
Narges Chinichian's avatar
Narges Chinichian committed
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
    global x0
    σ = sig
    # Weight parameters
    edgep = Array{Float64}(undef, 2, ne(g))
    for (i, edge) in zip(1:ne(g), collect(edges(g)))
        edgep[1, i] = synaptic_weights[edge.src, edge.dst]
        edgep[2, i] = synaptic_weights[edge.dst, edge.src]
    end
    edgep .*= σ # multiply by coupling strength beforehand to save computations later

    @inline function fhn_stp!(, ξ, e_s, e_d, p, t)
        # u, v, x, y = ξ[1:4]
        # FHN model
        [1] = ξ[1] - ξ[1]^3 * α - ξ[2] + .8
        [2] = ξ[1] * ϵ +  -  * ξ[2]
        # STP model
        [3] = τᴰinv - τᴰinv * ξ[3] + ξ[3] * ξ[4] * (u₀ - ξ[1])
        [4] = y₀τᶠinv - ξ[4] * τᶠinv - y₀ * (ξ[4] - 1.) * (ξ[1] - u₀)
70
        # coupling
Narges Chinichian's avatar
Narges Chinichian committed
71
72
73
        @inbounds for e in e_s # edges for which vertex is source
            [1] += e[1]
        end
lindnemi's avatar
lindnemi committed
74
75
76
        #@inbounds for e in e_d
        #    dξ[1] += e[2] # edges for which vertex is destination
        #end
Narges Chinichian's avatar
Narges Chinichian committed
77
78
79
80
81
82
83
        nothing
    end

    @inline Base.@propagate_inbounds function electrical_edge!(e, v_s, v_d, p, t)
        # The coupling is not symmetric wrt. change of source and destination
        # Hence we compute the flow in both directions
        e[1] =  v_d[3] * v_d[4] * (v_d[1] - v_s[1]) * p[1]
lindnemi's avatar
lindnemi committed
84
        #e[2] =  v_s[3] * v_s[4] * (v_s[1] - v_d[1]) * p[2]
Narges Chinichian's avatar
Narges Chinichian committed
85
86
87
        nothing
    end

lindnemi's avatar
lindnemi committed
88
    electricaledge = StaticEdge(f! = electrical_edge!, dim = 1)
Narges Chinichian's avatar
Narges Chinichian committed
89
90
91
92
93
    odevertex = ODEVertex(f! = fhn_stp!, dim = 4, sym = [:u, :v, :x, :y])

    fhn_network! = network_dynamics(odevertex, electricaledge, g);

    p = (nothing, edgep) # network_dynamics uses tuple syntax for parameters
Narges Chinichian's avatar
Narges Chinichian committed
94
    tspan = (0., 100000.)
Narges Chinichian's avatar
Narges Chinichian committed
95

Narges Chinichian's avatar
Narges Chinichian committed
96
    prob  = ODEProblem(fhn_network!, x0, tspan, p)
Narges Chinichian's avatar
Narges Chinichian committed
97

98
    @time sol = solve(prob, Tsit5(),reltol=1e-6, saveat=90000:.2:100000);
Narges Chinichian's avatar
Narges Chinichian committed
99
100
101
102
103
104
105
106
107
108
109
110
    #we only need the values for u from the sol
    #u values are saved in indices 1,5,9,13,...
    #sol is of the size 1508(377*4)*5001

    # u_noinput = plot(sol, tspan = tspan, vars = idx_containing(fhn_network!, :u), legend=false,color = "darkred")
    # savefig(u_noinput,"u_noinput_$σ.png")
    # v_noinput = plot(sol, tspan = tspan, vars = idx_containing(fhn_network!, :v), legend=false, color = "lightblue")
    # savefig(v_noinput,"v_noinput_$σ.png")

    # plot(sol,  vars = idx_containing(fhn_network!, :x), legend=false, color = "darkred")
    # xy = plot!(sol,  vars = idx_containing(fhn_network!, :y), legend=false, color = "lightblue")
    # savefig(xy,"xy_$σ.png")
Narges Chinichian's avatar
Narges Chinichian committed
111
112
    writedlm( "/Volumes/004915257822556/sigma-output-compactInfo/reltol_rev_100t_u_v_x_y_$σ.csv",  sol[:,end-1000:end], ',')
    x0 = sol[:,end]
Narges Chinichian's avatar
Narges Chinichian committed
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
    # tstart = 1
    # plot(sol[tstart:end], vars=(1,2))
    # plot!(sol[tstart:end], vars=(5,6))
    # limit_cicle = plot!(sol[tstart:end], vars=(33,34))
    # savefig(limit_cicle,"limit_cicle_$σ.png")

    # hm = heatmap(Array(sol)[1:4:end,1:end])
    # savefig(hm,"heatmap_$σ.png")

    # anim = @animate for t in 250:4:2000
    #     display_time = round(.2t, digits=3)
    #     scatter(sol[1:4:end, t], ylims = [-2, 2], title="t=$display_time", legend=false)
    # end
    # gif(anim,"oscil_$σ.gif", fps=12)

    # using FFTW

    # spectrum= (abs.(fft(sol[1,250:1000] .- mean(sol[1,250:1000]))).^2)
    # findall(x -> x == maximum(spectrum), spectrum)

    # vline([-0.813, 0.813]./2, label = "T = 2.46", linewidth = 2.5, color=:red)
    # plspec = plot!( collect(-75:0.2:75) ./ 150, abs.(fft(sol[1,250:1000] .- mean(sol[1,250:1000]))).^2, color = :darkblue, xlabel="Hz")
    # xticks!([-1, -0.813, -0.5,0.,0.5,0.813,1.] ./2)
    # savefig(plspec,"spect_$σ.png")
Narges Chinichian's avatar
Narges Chinichian committed
137
    #println(x0)
Narges Chinichian's avatar
Narges Chinichian committed
138

139
140

end