Non linear ODE system matlab fourier analysis problem

[calum.g]

Hello all.

I am studying a system and want to investigate how the frequency of y(3) varies under different conditions. However, my the fft I perform on it tells me the frequency is zero, which must be incorrect. I have tried a stack of things but can't see what the problem is. I'm relatively new to the MATLAB environment.

The main code:

close all; clear all;
%--------------------------------------%
Time control and initial conditions
timeinmin=600; % Time in minutes
tfin=timeinmin*60; % Time in seconds
%--------------------------------------%
Initial Levels and system Saturation
totIKK=0.08; % Total IKKn at start
totNFkB=0.08; % Total NFkB cytoplasmic at start
T=1; % Saturation of the system
yint=[totNFkB,0,0,0,0,0,totIKK,0,0,0,0,totIKK,totNFkB,T]; % Initial conditions
%--------------------------------------%
ODEs
[t,y]=ode45(@nfkb,[0 tfin],yint); % calls Model
%--------------------------------------%
Extract Data from the ODE
fs = 1; % The sampling freqency in Hz
time = 0:1/fs:tfin; % Define time of interesrt
Datapoints = tfin*fs; % # of Data points being used
yi = interp1(t,y(:,3),time);
%plot(time/60,yi,'o',t/60,y(:,3)) % test to ensure it is working
%dbstop
%--------------------------------------%
Fourier Analysis from
NFFT = 2^nextpow2(tfin); % Next power of 2 from length of y
Y = fft(yi,NFFT)/Datapoints;
f = fs/2*linspace(0,1,NFFT/2+1);
% Plot single-sided amplitude spectrum.
plot(f,2*abs(Y(1:NFFT/2+1)))
title('Single-Sided Amplitude Spectrum of y(t)')
xlabel('Frequency (Hz)')
ylabel('Amplitude - |Y(f)|')
% Get value of frequency
[B,IX] = sort(2*abs(Y(1:NFFT/2+1)));
BFloor=0.1; %BFloor is the minimum amplitude value (ignore small values)
Amplitudes=B(B>=BFloor) %find all amplitudes above the BFloor
Frequencies=f(IX(1+end-numel(Amplitudes):end)) %frequency of the peaks

The ODEs:

function dydt=nfkb(t,y)
parameters;
dydt=[kd1*(y(13)-(y(3)+y(5))/kv-y(11)-y(1))-ka1*y(2)*y(1)-ki1*y(1)+...
kdegc*(y(13)-(y(3)+y(5))/kv-y(11)-y(1))+ke1f*y(3)+kdegpin*y(11);
kd1*(y(13)-(y(3)+y(5))/kv-y(11)-y(1))-ka1*y(2)*y(1)-ki2*y(2)+...
kdegf*y(2)+ktria*y(6)-kc2*y(8)*y(2);
kd1*y(5)-ka1*y(4)*y(3)+kv*ki1*y(1)-kv*ke1f*y(3);
kd1*y(5)-ka1*y(4)*y(3)+kv*ki2*y(2)-kv*ke2*y(4)-kdegf*y(4);
ka1*y(4)*y(3)-kd1*y(5)-kv*ke1c*y(5);
kitria*(y(3).^h)/((y(3).^h)+(k.^h))-kdegtia*y(6);
kp*(y(12)-y(8)-y(7))*kbA20/(kbA20+y(10)*y(14))-y(14)*ka*y(7);
y(14)*ka*y(7)-ki*y(8);
kitra*(y(3).^h)/((y(3).^h)+(k.^h))-kdegta*y(9);
ktra*y(9)-kda*y(10);
kc3*y(8)*(y(13)-(y(3)+y(5))/kv-y(11)-y(1))-kdegpin*y(11);
0;
0;
0];

The parameters:

%rates
kv=3.3; % C:N ratio
kp=0.0006; % Conversion rate of IKKi
ka=0.004; %activation rate of IKK
ki=0.003; % Usage rate of IKK
kda=0.0045; % Deg rate of A20 protein
kbA20=0.0018; % A20 concentration
kc2=0.074; % phosopho of IkB
kc3=0.37; % phospho of IkB-NFkB
kdegpin=0.1; % splitting of Ikb-NFkB due to phospho
ka1=0.5; % nuclear IkB-NFkB formation
kd1=0.0005; % dissisociation of nuclear IkB-NFkB
kdegf=0.0005; % Deg rate of nuclear IkB
kdegc=0.000022; % dissisociation of cyto IkB-NFkB
ki1=0.0026; % nuclear entry of NF-kB
ke1f=52*10^(-6); % nuclear exit of NF-kB
ke1c=0.01; % Transport of IkB-NFkB out of nuc.
ki2=0.00067; % nuclear entry of IkB
ke2=3.35*10^(-4); % Nuclear exit of IkB
h=2; % hill coefficient of IkB transcription
k=0.065; % dissasociation constant for IkB tanscription with NFkB
kitria=1.4*10^(-7); % trans rate of IkB
ktria=0.5; % translation of IkB
kdegtia=0.0003; % Deg rate of IkB RNA
kitra=1.4*10^(-7); % Transcription rate of A20 RNA
ktra=0.5; % transcription rate of A20
kdegta=0.00048; % deg rate of A20 RNA

Thank you in advance for any insights.

 

[Valkarie]

It's possible that the frequency of y(3) is indeed zero, and that you're not seeing any variation because it's not changing. If you can provide more information about the system you're studying and the conditions you're testing, it may help us to troubleshoot the issue.