Solve 10 system of ODEs with separate functions and time dependent factors

Hi,

I am trying to solve this system of ODEs and I am using ode45. However, I do not know how to write the input Up(t) and Uf(t) which are essentially functions that are noise genereate + depend on previous values of Zp(t). See image for input. Also, I am not sure if my code is correct when handling the non-differential output (i.e. Zp(t), Vp(t)). Am I correct in creating a seperate variable X matrix? Also, if I wanted the ode solver to return Vp(t), how would I do so?

Thanks!

I have tried to solve this system using the following code but cant seem to generate the input:

clear all, close all, clc;

% filename = 'fitting.xlsx';
% fit = xlsread(filename);

fit = [-0.0500000000000000,-0.0350000000000000,-0.00650000000000000;54.5000000000000,57,31;81,97.5000000000000,136.500000000000;0.100000000000000,39,21;4.70000000000000,10.5000000000000,11.5000000000000;16,16,18;NaN,0,16.5000000000000;94.5000000000000,NaN,0;0.570000000000000,25.5000000000000,NaN;NaN,81,24.5000000000000;75.5000000000000,NaN,11.5000000000000;0,0,NaN;NaN,1,8.30000000000000;1,NaN,16.6000000000000;8.30000000000000,16.6000000000000,NaN];

Tspan = linspace(0,1.5,100);
[t, y] = ode45(@odefun,Tspan, zeros(1,10))

figure(1)
plot(t,y)


function dydt = odefun(t,y)
fit = [-0.0500000000000000,-0.0350000000000000,-0.00650000000000000;54.5000000000000,57,31;81,97.5000000000000,136.500000000000;0.100000000000000,39,21;4.70000000000000,10.5000000000000,11.5000000000000;16,16,18;NaN,0,16.5000000000000;94.5000000000000,NaN,0;0.570000000000000,25.5000000000000,NaN;NaN,81,24.5000000000000;75.5000000000000,NaN,11.5000000000000;0,0,NaN;NaN,1,8.30000000000000;1,NaN,16.6000000000000;8.30000000000000,16.6000000000000,NaN];
% Empty solution set
dydt = zeros(10,1);
% Empty functions
x = zeros(8,1);

% Constants from table 4
w_e = 75;
w_s = 30;
w_f = 75;
G_e = 5.17;
G_s = 4.45;
G_f = 57.1;
C_ep = 5;
C_pe = 25;
C_sp = 60;
e_0 = 2.5;
r = 0.56;

% Constants from fitting BA 19
C_ps_BA19 = fit(2,1);
C_pf_BA19 = fit(5,1);
C_fp_BA19 = fit(3,1);
C_fs_BA19 = fit(4,1);
C_ff_BA19 = fit(6,1);



% A DEF Pyramidal Neurons C
y_pDt = dydt(1);
x_pDt = dydt(2);

y_p = y(1);
x_p = y(2);

z_p = x(1);
v_p = x(2);

% B DEF Excitatory Neurons C
y_eDt = dydt(3);
x_eDt = dydt(4);

y_e = y(3);
x_e = y(4);

z_e = x(3);
v_e = x(4);

% C DEF Slow Inhibatory Interneurons C
y_sDt = dydt(5);
x_sDt = dydt(6);

y_s = y(5);
x_s = y(6);

z_s = x(5);
v_s = x(6);

% D DEF Fast Inhibatory Interneurons C
y_fDt = dydt(7);
x_fDt = dydt(8);
y_lDt = dydt(9);
x_lDt = dydt(10);

y_f = y(7);
x_f = y(8);
y_l = y(9);
x_l = y(10);

z_f = x(7);
v_f = x(8);

% Generate white noise
% variance = 5;
% AAA = randn(1, 0, variance);
AAA = sqrt(5)*randn(1,1);


% ODE Pyramidal Neurons
y_pDt = x_p;
x_pDt = G_e*w_e*z_p - 2*w_e*x_p - w_e^2*y_p;
z_p = 2*e_0./(1+exp(-r*v_p)) - e_0;
v_p = C_pe*y_e - C_ps_BA19*y_s - C_pf_BA19*y_f;

% Input noise


% ODE Excitatory Interneurons
y_eDt = x_e;
x_eDt = G_e*w_e*(z_e + AAA/C_pe) - 2*w_e*x_e - w_e^2*y_e;
z_p = 2*e_0./(1+exp(-r*v_e)) - e_0;
v_e = C_ep*y_p;

% ODE Slow Inhibatory Interneurons
y_sDt = x_s;
x_sDt = G_s*w_s*z_s - 2*w_s*x_s - w_s^2*y_s;
z_p = 2*e_0./(1+exp(-r*v_s)) - e_0;
v_e = C_sp*y_p;

% ODE Slow Inhibatory Interneurons
y_fDt = x_f;
x_fDt = G_f*w_f*z_f - 2*w_f*x_f - w_f^2*y_f;
y_lDt = x_l;
x_lDt = G_e*w_e*AAA - 2*w_e*x_l - w_e^2*y_l;
z_f = 2*e_0./(1+exp(-r*v_f)) - e_0;
v_f = C_fp_BA19*y_p - C_fs_BA19*y_s - C_ff_BA19*y_f + y_l;



% put back the solutions
dydt(1) = y_pDt;
dydt(2) = x_pDt;
dydt(3) = y_eDt;
dydt(4) = x_eDt;
dydt(5) = y_sDt;
dydt(6) = x_sDt;
dydt(7) = y_fDt;
dydt(8) = x_fDt;
dydt(9) = y_lDt;
dydt(10) = x_lDt;



end

Answer 1

[James Tursa]

This isn't going to work. You can't use random values in the derivative function for ode45. The derivative calls that ode45 makes to your function must be consistent ... i.e. the same exact derivative should result when called repeatedly with the same inputs. This won't happen if you have randomness in your function, e.g. this line

AAA = sqrt(5)*randn(1,1);

The way it is now, ode45 will be thoroughly confused as to what the derivative is at any given point since repeated calls near the same point won't give consistent results. You will either get an error or a wrong result.

[Ali Najmaldin]

Thanks for your help. Could you point me in the right direction. What ODE solver should I be looking at? How should I go about implementing the randomness element (the input Is noise + the function Zp(t))

And if I want to implement a function that has no differential (doesnt have dydt in this case), would I just create a new vector for functions?

Thanks

[Torsten]

https://en.wikipedia.org/wiki/Stochastic_differential_equation