- #1

calum.g

- 11

- 0

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.