Adjacency matrix for ODE45

Hi

I am trying to model a systema as below:

where kij is :

I tried to model two sets of states manually first, without Aij (the adjacency matrix) as below:

clear all;clc;

%%

global b r I x0 alpha beta gamma

b=3.0;

r=0.006;

I=2.8;

x0=-1.6;

alpha=1;

beta=10;

gamma=0.5;

initial_condition=[-1,0.2,0.4,0.8,0,0.2,0.1,0.2];

tspan=[0 3000];

[t,y]=ode45(@EquationSys,tspan,initial_condition);

plot(t,y(:,1))

hold on

%%

function dy=EquationSys(t,y)

global b r I x0 alpha beta gamma

x1=y(1); %-1

y1=y(2);%0.2

z1=y(3);%0.4

x2=y(4);%0.8

y2=y(5);%0

z2=y(6);%0.2

k12=y(7);%0.1

k21=y(8);%0.2

dy=[y1-x1^3+b*x1^2-z1+I+k12*(x2-x1);

1-5*x1^2-y1;

r*(4*(x1-x0)-z1);

y2-x2^3+b*x2^2-z2+I+k21*(x1-x2);

1-5*x2^2-y2;

r*(4*(x2-x0)-z2);

k12*(alpha*exp(-beta*(x1 - x2)^2) - gamma*(k12+1));

k21*(alpha*exp(-beta*(x2 - x1)^2) - gamma*(k21+1));

];

end

%

%%%%%%%%%%%%%%%%%

For some reason, I want to add Aij to this model and vectorization the code and I programmed the code as below:

clear all;clc;

% HR Parameters

global b r x0 I N alpha beta gamma Adj_Matrix

b=3.0;

r=0.006;

I=2.8;

x0=-1.6;

alpha=1;

beta=10;

gamma=0.5;

%%

% Simulation parameters

T = 3000; % Total time

time_span = [0 T]; % Time vector

%%

% Number of neurons

% n=2; %Number of neurons on each edges

N = 2; % Number of all neurons

%% Adjacency Matrix

% Initialize the matrix with zeros

Adj_Matrix = zeros(N, N);

Adj_Matrix(1,2)=1;

Adj_Matrix(2,1)=1;

%%

% Flatten initial conditions and initial coupling strengths into a single vector

initial_conditions = [-1,0.8,0.2,0,0.4,0.2,0,0.1,0.2,0];

%%

% Solve the ODE for the coupled system with STDP

[time, sol] = ode45(@ hr_equations, time_span, initial_conditions);

%x, y, z, and k from the solution

x = sol(:, 1:N);

y = sol(:, N+1:2*N);

z = sol(:, 2*N+1:3*N);

k = sol(:, 3*N+1:end);

%%

% % Plotting results

%

% % If you want to visualize all the neurons' behavior together

plot(time,x(:,1))

hold on

%%

function dstate_dt = hr_equations(time, state)

global b r x0 I N alpha beta gamma Adj_Matrix

% Reshape the state vector

x = state(1:N);

y = state(N+1:2*N);

z = state(2*N+1:3*N);

k = reshape(state(3*N+1:end), N, N);

% HR Model Equations

dx = y-x.^3+b*x.^2-z+I+sum((Adj_Matrix.*k).*(x'-x),2);

dy = 1-5*x.^2-y;

dz = r*(4*(x-x0)-z);

dk = zeros(N, N);

for i = 1:N

for j = 1:N

if i ~= j

dk(i,j) = k(i,j) * (alpha * exp(-beta * (x(i) - x(j))^2) - gamma * (k(i,j) + 1));

end

end

end

% Flatten the derivatives to return as a single column vector

dstate_dt = [dx; dy; dz; reshape(dk, [], 1)];

end

However, after running this code and comparing two plots I found that there is slight differences between two results. I cannot find the issue. Can you help me in this regard?

Thanks

Answer 1

[Torsten]

k = [0,0.1,0.2,0];
reshape(k,2,2)

ans = 2x2

0 0.2000 0.1000 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

I think you want the reshaped matrix to be

reshape(k,2,2).'

ans = 2x2

0 0.1000 0.2000 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[N/A]

Thanks! it works.