-
Notifications
You must be signed in to change notification settings - Fork 4
/
run_brain_control.m
120 lines (82 loc) · 2.56 KB
/
run_brain_control.m
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
function run_brain_control
% Local control design and dynamical simulation of brain network
clear;clc;
addpath(['.' filesep 'Brainnet']);
%% Data Preparation
data=load('Coactivation_matrix');
alpha_p=1935;
alpha_h=2852;
beta=150;
gamma=100;
Y=data.Coactivation_matrix;
%=================================
% take part of the network for illustration purpose
% comment out this block to use the whole network
n = 50;
Y=Y(1:n,1:n);
%================================
n=size(Y,1);
actnode=1:n;
Q=[100*ones(n,1);1e-8*ones(n,1);ones(n,1)];
R=ones(length(actnode),1)*1e-6;
dataMat=beta*Y+alpha_p*eye(n);
Wup=InfoDistanceUpperbdd(dataMat);
%% Build Linear System and SIN
radius=10;
l=cell(n,1);
for i=1:n
[nlist, ~] = ucs_geodesic_k(Wup,i,radius);
l{i}=sparse(1,n); l{i}(nlist)=1;
end
L=vertcat(l{:});
x0_h=[0.1*ones(n,1); zeros(n,1);];
x0_p=[0.2*ones(n,1); zeros(n,1); zeros(n,1)];
dt = 1.e-3; T = 1; N = ceil(T/dt);
xs_h = zeros(N+1,2*n);
xs_p = zeros(N+1,3*n);
xs_pc = zeros(N+1,3*n);
f_h = @(t,x)brain_eqns(t,x,Y,alpha_h,gamma,beta);
f_p = @(t,x,x_ref,u)controlled_brain_eqns(t,x,actnode,Y,alpha_p,gamma,beta,x_ref,u);
u0=zeros(length(actnode),1);
x_h=x0_h;
x_p=x0_p;
x_pc=x0_p;
for nt = 1:N
xs_h(nt,:) = x_h';
xs_p(nt,:) = x_p';
xs_pc(nt,:) = x_pc';
t = (nt-1)*dt
sys=BrainLinearSys(x_p,Y,alpha_p,gamma,beta);
% [Klocal1,S,e] = lqr(sys.A,sys.B(:,actnode),diag(Q),diag(R));
% K1=full(Klocal1); % Global control strategy
Klocal=optcontrol_brain_local(sys,Q,R,L,actnode); % Local control strategy
u = -Klocal*[x_pc(1:n)-x_h(1:n);x_pc(n+1:2*n);x_pc(2*n+1:3*n)];
x_p = rk4_step_ctrl(f_p,t,x_p,x_h,u0,dt);
x_pc = rk4_step_ctrl(f_p,t,x_pc,x_h,u,dt);
x_h = rk4_step_ref(f_h,t,x_h,dt);
end
plot_num=7;
figure;
ts = (0:N).'*dt;
plot(ts,xs_h(:,plot_num)); hold on;
plot(ts,xs_p(:,plot_num)); hold on;
plot(ts,xs_pc(:,plot_num)); hold on;
ylabel('x (mV)');
xlabel('time (s)');
legend('healthy','pathological','treated');
end
%% ODE solvers
function x = rk4_step_ref(f,t,x0,dt)
k1 = dt*feval(f,t,x0);
k2 = dt*feval(f,t+0.5*dt,x0 + k1/2);
k3 = dt*feval(f,t+0.5*dt,x0 + k2/2);
k4 = dt*feval(f,t+dt,x0 + k3);
x = x0 + (k1 + 2*k2 + 2*k3 + k4)/6;
end
function x = rk4_step_ctrl(f,t,x0,x_ref,K,dt)
k1 = dt*feval(f,t,x0,x_ref,K);
k2 = dt*feval(f,t+0.5*dt,x0 + k1/2,x_ref,K);
k3 = dt*feval(f,t+0.5*dt,x0 + k2/2,x_ref,K);
k4 = dt*feval(f,t+dt,x0 + k3,x_ref,K);
x = x0 + (k1 + 2*k2 + 2*k3 + k4)/6;
end