Why my plot is wrong if I decrease dt

I have a code, while I running the code is nothing error in the result. But if dt decrease (ex: 0.001) the plot is wrong.

I know the plot is wrong because I do that in excel. In excel is nothing wrong. I don't know what's wrong with my plot code

clc;clear;

%parameter
  delta=50;
  gamma=75;
  K1=10^-4;
  K2=5*10^-4;
  K3=10^-3;
  K4=5*10^-3;
  K5=10^-2;
  K6=5*10^-2;
  Ko=0.1;
  n=6;
  Oa=10;
  Pa=100;
  mu_1=10^-3;
  mu_2=10^-3;
  mu_3=10^-3;
  mu_4=10^-3;
  mu_5=10^-3;
  mu_6=10^-3;
  mu_o=10^-4;
  mu_p= 10^-5;
  
  %input for time
  t(1)=0;
  dt=0.001; %time interval
  t=0:dt:180; %time span

%component for microglia
Mg = 0.047;
lambdaMf = 5*10^-2; 
Fo = 3.36*10^-11;
KFo = 2.58*10^-11;
lambdaMa = 2.23*10^-2;
Ao = 0.14;
KAo = 10^-7;
beta = 10;
epsilon1 = 0.3333;
epsilon2 = 0.8;
lambdaM1Tb = 6*10^-3; 
Tb=10^-6;
KTb = 2.5*10^-7;
makro1 = 0.02;
makro2 = 0.02;
dM1= 0.015;
dM2= 0.015;
tmakro1 = Mg*(lambdaMf*(Fo/(Fo+KFo))+lambdaMa*(Ao/(Ao+KAo)))*(epsilon1*beta/beta*(epsilon1+epsilon2))-(lambdaM1Tb*(Tb/(Tb+KTb)*makro1))-(dM1*makro1);
tmakro2 = Mg*(lambdaMf*(Fo/(Fo+KFo))+lambdaMa*(Ao/(Ao+KAo)))*(epsilon2/beta*(epsilon1+epsilon2))+(lambdaM1Tb*(Tb/(Tb+KTb))*makro1)-(dM2*makro2);

%component drugs
Abi = 10^-6;
lambdaN = 8*10^-4; 
lambdaA = 8*10^-4; 
daboM = 2*10^-3; 
dabom = 10^-2; %clearance rate by macrophages 
makrofag1 = 0;
makrofag2 = 0;
mikro1 = 0.02;
mikro2 = 0.02;
teta = 0.9; 
h = 10;
Kaob = 7*10^-3;
AoB = 10^-8;
dtdab = Abi*(dt*(0.14*1.9*10^-4*exp(-1.9*10^-4)))+ (lambdaN*1)+(lambdaA*1)-(dabom*(makrofag1+(teta*makrofag2))+daboM*(mikro1+(teta*mikro2))*(1+h))*(AoB/(AoB+Kaob));

%initial condition with microglia
M1(1)=10;
M2(1)=0;
M3(1)=0;
M4(1)=0;
M5(1)=0;
M6(1)=0;
O(1)=0;
P(1)=0;

%initial condition without microglia
M12(1)=10;
M22(1)=0;
M32(1)=0;
M42(1)=0;
M52(1)=0;
M62(1)=0;
O2(1)=0;
P2(1)=0;

%initial condition with drugs depend time
M14(1)=10;
M24(1)=0;
M34(1)=0;
M44(1)=0;
M54(1)=0;
M64(1)=0;
O4(1)=0;
P4(1)=0;

%empty array with microglia
T=zeros(length(t)+1,1); %empty array for t
M1=zeros(length(t)+1,1); %empty array for M1
M2=zeros(length(t)+1,1); %empty array for M2
M3=zeros(length(t)+1,1); %empty array for M3
M4=zeros(length(t)+1,1); %empty array for M4
M5=zeros(length(t)+1,1); %empty array for M5
M6=zeros(length(t)+1,1); %empty array for M6
O=zeros(length(t)+1,1);  %empty array for O
P=zeros(length(t)+1,1);  %empty array for P

sumter=K2*M2+K3*M3+K4*M4+K5*M5;

%empty array without microglia
M12=zeros(length(t)+1,1); %empty array for M1
M22=zeros(length(t)+1,1); %empty array for M2
M32=zeros(length(t)+1,1); %empty array for M3
M42=zeros(length(t)+1,1); %empty array for M4
M52=zeros(length(t)+1,1); %empty array for M5
M62=zeros(length(t)+1,1); %empty array for M6
O2=zeros(length(t)+1,1); %empty array for O
P2=zeros(length(t)+1,1); %empty array for P

sumter2=K2*M22+K3*M32+K4*M42+K5*M52;

%empty array with drugs depend time
M14=zeros(length(t)+1,1); %empty array for M1
M24=zeros(length(t)+1,1); %empty array for M2
M34=zeros(length(t)+1,1); %empty array for M3
M44=zeros(length(t)+1,1); %empty array for M4
M54=zeros(length(t)+1,1); %empty array for M5
M64=zeros(length(t)+1,1); %empty array for M6
O4=zeros(length(t)+1,1); %empty array for O
P4=zeros(length(t)+1,1); %empty array for P

sumter2=K2*M24+K4*M34+K4*M44+K5*M54;

for j = 1:length(t)
    %with microglia
    T(j+1)=T(j)+dt;
    M1(j+1)=M1(j)+1./(1+exp(-T(j)));
    
    M1(j+1) = M1(j)+(dt*(delta*M1(j+1)*(1-(M1(j+1)/gamma))-2*K1*M1(j+1)*M1(j+1)-M1(j+1)*sumter(j+1))-((Oa-n)*K6*M1(j+1)*M6(j+1))-((Pa-Oa)*Ko*M1(j+1)*O(j+1))-(mu_1*M1(j+1)));
    M2(j+1) = M2(j)+(dt*(K1*M1(j)*M1(j)-K2*M1(j)*M2(j))-(mu_2*M2(j+1))-tmakro1-tmakro2);
    M3(j+1) = M3(j)+(dt*(K2*M1(j)*M2(j)-K3*M1(j)*M3(j))-mu_3*M3(j));
    M4(j+1) = M4(j)+(dt*(K3*M1(j)*M3(j)-K4*M1(j)*M4(j))-mu_4*M4(j));
    M5(j+1) = M5(j)+(dt*(K4*M1(j)*M4(j)-K5*M1(j)*M5(j))-mu_5*M5(j));
    M6(j+1) = M6(j)+(dt*(K5*M1(j)*M5(j)-K6*M1(j)*M6(j))-mu_6*M6(j));
    O(j+1)  = O(j)+(dt*(K6*M1(j)*M6(j)-Ko*M1(j)*O(j)-mu_o*O(j)));
    P(j+1)  = P(j)+(dt*(Ko*M1(j)*O(j)-mu_p*P(j)));
    
   
    %without microglia
    T(j+1)=T(j)+dt;
    M12(j+1) = M12(j)+1./(1+exp(-T(j)));
    
    M12(j+1) = M12(j)+(dt*(delta*M12(j+1)*(1-(M12(j+1)/gamma))-2*K1*M12(j+1)*M12(j+1)-M12(j+1)*sumter2(j+1))-((Oa-n)*K6*M12(j+1)*M62(j+1))-((Pa-Oa)*Ko*M12(j+1)*O2(j+1))-(mu_1*M12(j+1)));
    M22(j+1) = M22(j)+(dt*(K1*M12(j)*M12(j)-K2*M12(j)*M22(j))-(mu_2*M22(j+1)));
    M32(j+1) = M32(j)+(dt*(K2*M12(j)*M22(j)-K3*M12(j)*M32(j))-mu_3*M32(j));
    M42(j+1) = M42(j)+(dt*(K3*M12(j)*M32(j)-K4*M12(j)*M42(j))-mu_4*M42(j));
    M52(j+1) = M52(j)+(dt*(K4*M12(j)*M42(j)-K5*M12(j)*M52(j))-mu_5*M52(j));
    M62(j+1) = M62(j)+(dt*(K5*M12(j)*M52(j)-K6*M12(j)*M62(j))-mu_6*M62(j));
    O2(j+1)  = O2(j)+(dt*(K6*M12(j)*M62(j)-Ko*M12(j)*O2(j)-mu_o*O2(j)));
    P2(j+1)  = P2(j)+(dt*(Ko*M12(j)*O2(j)-mu_p*P2(j)));

%with drugs depend by time
    T(j+1)=T(j)+dt;
    M14(j+1) = M14(j)+1./(1+exp(-T(j)));
    
    M14(j+1) = M14(j)+(dt*(delta*M14(j+1)*(1-(M14(j+1)/gamma))-2*K1*M14(j+1)*M14(j+1)-M14(j+1)*sumter2(j+1))-((Oa-n)*K6*M14(j+1)*M64(j+1))-((Pa-Oa)*Ko*M14(j+1)*O4(j+1))-(mu_1*M14(j+1)));
    M24(j+1) = M24(j)+(dt*(K1*M14(j)*M14(j)-K2*M14(j)*M24(j))-(mu_2*M24(j+1))-tmakro1-tmakro2-dtdab);
    M34(j+1) = M34(j)+(dt*(K2*M14(j)*M24(j)-K3*M14(j)*M34(j))-mu_3*M34(j));
    M44(j+1) = M44(j)+(dt*(K3*M14(j)*M34(j)-K4*M14(j)*M44(j))-mu_4*M44(j));
    M54(j+1) = M54(j)+(dt*(K4*M14(j)*M44(j)-K5*M14(j)*M54(j))-mu_5*M54(j));
    M64(j+1) = M64(j)+(dt*(K5*M14(j)*M54(j)-K6*M14(j)*M64(j))-mu_6*M64(j));
    O4(j+1)  = O4(j)+(dt*(K6*M14(j)*M64(j)-Ko*M14(j)*O4(j)-mu_o*O4(j)));
    P4(j+1)  = P4(j)+(dt*(Ko*M14(j)*O4(j)-mu_p*P4(j)));

end

subplot (2,2,1)
plot(T,M1,'k', T,M12,'r',T,M14,'b','Linewidth',2)
legend ('M1 dengan with immune', 'M1 without immune', 'M1 with drugs');
%xticks ([60 120 180 240 300  360 420 480 540 600 660 720 780 840 900 960 1020 1080])
xticks ([30 60 90 120 150 180  ])
xlabel('time (day)','Fontsize',14,'FontWeight','bold')
ylabel('concentration (gr/mL)','Fontsize',14,'FontWeight','bold')
title ('M1')

subplot (2,2,2)
plot(T,M2,'k',T,M22,'r',T,M24,'b','Linewidth',2)
legend ('M2 with immune', 'M2 without immune','M2 with drug');
xticks ([30 60 90 120 150 180  ])
yticks ([ 0 2 4 6 8 10 12 14 16])
xlabel('time (day)','Fontsize',14,'FontWeight','bold')
ylabel('concentration (gr/mL)','Fontsize',14,'FontWeight','bold')
title ('M2')

subplot (2,2,3)
plot(T,M3,'k',T,M32,'r',T,M34,'b','Linewidth',2)
legend ('M3 with immune', 'M3 without immune','M3 with drug');
xticks ([30 60 90 120 150 180  ])
yticks ([0 0.5 1 1.5 2 2.5 3 3.5])
xlabel('time (day)','Fontsize',14,'FontWeight','bold')
ylabel('concentration (gr/mL)','Fontsize',14,'FontWeight','bold')
title ('M3')

subplot (2,2,4)
plot(T,M4,'k',T,M42,'r',T,M44,'b','Linewidth',2)
legend ('M4 with immune', 'M4 without immune','M4 with drug');
xticks ([30 60 90 120 150 180  ])
yticks ([0 0.1 0.2 0.3 0.4 0.5])
xlabel('time (day)','Fontsize',14,'FontWeight','bold')
ylabel('concentration (gr/mL)','Fontsize',14,'FontWeight','bold')
title ('M4')

subplot (2,2,1)
plot(T,M5,'k',T,M52,'r',T,M54,'b','Linewidth',2)
legend ('M5 with immune', 'M5 without immune','M5 with drug');
xticks ([30 60 90 120 150 180  ])
yticks ([0 0.05 0.1 0.15 0.2])
xlabel('time (day)','Fontsize',14,'FontWeight','bold')
ylabel('concentration (gr/mL)','Fontsize',14,'FontWeight','bold')
title ("M5")

subplot (2,2,2)
plot(T,M6,'k',T,M62,'r',T,M64,'b','Linewidth',2)
legend ('M6 with immune', 'M6 without imune','M6 with drug');
xticks ([30 60 90 120 150 180  ])
yticks ([0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045])
xlabel('time (day)','Fontsize',14,'FontWeight','bold')
ylabel('concentration (gr/mL)','Fontsize',14,'FontWeight','bold')
title ("M6")

subplot (2,2,3)
plot(T,O,'k',T,O2,'r',T,O4,'b','Linewidth',2)
legend ('O with immune', 'O without immune','O with drug');
xticks ([30 60 90 120 150 180  ])
yticks ([0 0.005 0.01 0.015 0.02])
xlabel('time (day)','Fontsize',14,'FontWeight','bold')
ylabel('concentration (gr/mL)','Fontsize',14,'FontWeight','bold')
title ("O")

subplot(2,2,4)
plot(T,P,'k',T,P2,'r',T,P4,'b','Linewidth',2)
legend ('P with immune', 'P without immune','P with drug');
xticks ([30 60 90 120 150 180  ])
%yticks ([0 2 4 6 8 10 12])
xlabel('time (day)','Fontsize',14,'FontWeight','bold')
ylabel('concentration (gr/mL)','Fontsize',14,'FontWeight','bold')
title ("P")

[Sam Chak]

Hi @cindyawati cindyawati

It seems like you are using the Forward Euler method. However, this method is known to have the propagation of errors at every time step over an extended period. Could you please share the original ODEs? We can obtain more reliable results using MATLAB ode45 solver, which is suitable for non-stiff systems without discontinuities.

[Naveed Salman]

Why have you used M1(j+1) and M2(j+1) on the right hand side when they are initialized to zero?

They are not contributing to the equation.

And you want to include there effect then you need to solve the simultneous equation by taking these unknown to right hand side.

    M1(j+1) = M1(j)+(dt*(delta*M1(j+1)*(1-(M1(j+1)/gamma))-2*K1*M1(j+1)*M1(j+1)-M1(j+1)*sumter(j+1))-((Oa-n)*K6*M1(j+1)*M6(j+1))-((Pa-Oa)*Ko*M1(j+1)*O(j+1))-(mu_1*M1(j+1)));
    M2(j+1) = M2(j)+(dt*(K1*M1(j)*M1(j)-K2*M1(j)*M2(j))-(mu_2*M2(j+1))-tmakro1-tmakro2);

[Torsten]

You don't seem to learn from answers already given:

https://uk.mathworks.com/matlabcentral/answers/2045850-how-to-calculate-first-and-second-derivative-for-concentration?s_tid=prof_contriblnk

but I repeat my doubts:

I don't think that setting

T(j+1)=T(j)+dt;

three times in the loop is correct.

Further

M6(j+1) and sumter(j+1) on the right-hand side are always zero in

M1(j+1) = M1(j)+(dt*(delta*M1(j+1)*(1-(M1(j+1)/gamma))-2*K1*M1(j+1)*M1(j+1)-M1(j+1)*sumter(j+1))-((Oa-n)*K6*M1(j+1)*M6(j+1))-((Pa-Oa)*Ko*M1(j+1)*O(j+1))-(mu_1*M1(j+1)));

M2(j+1) on the right-hand side is always zero in

M2(j+1) = M2(j)+(dt*(K1*M1(j)*M1(j)-K2*M1(j)*M2(j))-(mu_2*M2(j+1))-tmakro1-tmakro2);

M62(j+1) and sumter2(j+1) on the right-hand side are always zero in

M12(j+1) = M12(j)+(dt*(delta*M12(j+1)*(1-(M12(j+1)/gamma))-2*K1*M12(j+1)*M12(j+1)-M12(j+1)*sumter2(j+1))-((Oa-n)*K6*M12(j+1)*M62(j+1))-((Pa-Oa)*Ko*M12(j+1)*O2(j+1))-(mu_1*M12(j+1)));

M22(j+1) on the right-hand side is always zero in

M22(j+1) = M22(j)+(dt*(K1*M12(j)*M12(j)-K2*M12(j)*M22(j))-(mu_2*M22(j+1)));

M64(j+1) and sumter2(j+1) on the right-hand side are always zero in

M14(j+1) = M14(j)+(dt*(delta*M14(j+1)*(1-(M14(j+1)/gamma))-2*K1*M14(j+1)*M14(j+1)-M14(j+1)*sumter2(j+1))-((Oa-n)*K6*M14(j+1)*M64(j+1))-((Pa-Oa)*Ko*M14(j+1)*O4(j+1))-(mu_1*M14(j+1)));

M24(j+1) on the right-hand side is always zero in

M24(j+1) = M24(j)+(dt*(K1*M14(j)*M14(j)-K2*M14(j)*M24(j))-(mu_2*M24(j+1))-tmakro1-tmakro2-dtdab);

[cindyawati cindyawati]

thanks for your reply @Torsten . So, i must make empty array for M24(j+1) etc?

and if my code error why can show the graph? Is my graph getting error too because my code is error?

[Sam Chak]

Hi @cindyawati cindyawati

Could you show the Ordinary Differential Equations (ODEs) that are used to model the drug concentrations within the field of pharmacokinetics?

[cindyawati cindyawati]

Oh sorry I miss the reply your response @Sam Chak. This is my code ONLY for concentration Amyloid-beta from monomer until become fibril

clc;clear;

%input parameter
  delta=50;
  gamma=75;
  K1=10^-4;
  K2=5*10^-4;
  K3=10^-3;
  K4=5*10^-3;
  K5=10^-2;
  K6=5*10^-2;
  Ko=0.1;
  n=6;
  Oa=10;
  Pa=100;
  mu_1=10^-3;
  mu_2=10^-3;
  mu_3=10^-3;
  mu_4=10^-3;
  mu_5=10^-3;
  mu_6=10^-3;
  mu_o=10^-4;
  mu_p= 10^-5;
  
  
%input initial condition
M1(1)=10;
M2(1)=0;
M3(1)=0;
M4(1)=0;
M5(1)=0;
M6(1)=0;
O(1)=0;
P(1)=0;

%input for time
t(1)=0;
dt=0.01; %time interval
t=0:dt:180; %time span

%input empty array
T=zeros(length(t)+1,1); %empty array for t
M1=zeros(length(t)+1,1); %empty array for M1
M2=zeros(length(t)+1,1); %empty array for M2
M3=zeros(length(t)+1,1); %empty array for M3
M4=zeros(length(t)+1,1); %empty array for M3
M5=zeros(length(t)+1,1); %empty array for M3
M6=zeros(length(t)+1,1); %empty array for M3
O=zeros(length(t)+1,1);
P=zeros(length(t)+1,1);

sumter=K2*M2+K3*M3+K4*M4+K5*M5;
 
for j = 1:length(t)
    T(j+1)=T(j)+dt;
    M1(j+1)=M1(j)+1./(1+exp(-T(j))); 
    
    M1(j+1) = M1(j)+(dt*(delta*M1(j+1)*(1-(M1(j+1)/gamma))-2*K1*M1(j+1)*M1(j+1)-M1(j+1)*sumter(j+1))-((Oa-n)*K6*M1(j+1)*M6(j+1))-((Pa-Oa)*Ko*M1(j+1)*O(j+1))-(mu_1*M1(j+1)));
    M2(j+1) = M2(j)+(dt*(K1*M1(j)*M1(j)-K2*M1(j)*M2(j))-(mu_2*M2(j+1)));
    M3(j+1) = M3(j)+(dt*(K2*M1(j)*M2(j)-K3*M1(j)*M3(j))-mu_3*M3(j));
    M4(j+1) = M4(j)+(dt*(K3*M1(j)*M3(j)-K4*M1(j)*M4(j))-mu_4*M4(j));
    M5(j+1) = M5(j)+(dt*(K4*M1(j)*M4(j)-K5*M1(j)*M5(j))-mu_5*M5(j));
    M6(j+1) = M6(j)+(dt*(K5*M1(j)*M5(j)-K6*M1(j)*M6(j))-mu_6*M6(j));
    O(j+1)  = O(j)+(dt*(K6*M1(j)*M6(j)-Ko*M1(j)*O(j)-mu_o*O(j)));
    P(j+1)  = P(j)+(dt*(Ko*M1(j)*O(j)-mu_p*P(j)));
end

subplot (2,2,1)
%figure
%hold on
plot(T,M1,'r','Linewidth',2)
xticks ([30 60 90 120 150 180  ])
xlabel('Time (days)')
ylabel('Concentration(gr/mL)')
title ('M1')

subplot (2,2,2)
%figure
%hold on
plot(T,M2,'g','Linewidth',2)
xticks ([30 60 90 120 150 180  ])
xlabel('Time (days)')
ylabel('Concentration(gr/mL)')
title ('M2')

subplot (2,2,3)
%figure
%hold on
plot(T,M3,'b','Linewidth',2)
xticks ([30 60 90 120 150 180  ])
xlabel('Time (days)')
ylabel('Concentration(gr/mL)')
title ('M3')

subplot (2,2,4)
%figure
%hold on
plot(T,M4,'m','Linewidth',2)
xticks ([30 60 90 120 150 180  ])
xlabel('Time (days)')
ylabel('Concentration(gr/mL)')
title ('M4')

subplot (2,2,1)
%figure
%hold on
plot(T,M5,'r','Linewidth',2)
xticks ([30 60 90 120 150 180  ])
xlabel('Time (days)')
ylabel('Concentration(gr/mL)')
title ('M5')

subplot (2,2,2)
%hold on
plot(T,M6,'g','Linewidth',2)
xticks ([30 60 90 120 150 180  ])
xlabel('Time (days)')
ylabel('Concentration(gr/mL)')
title ('M6')

subplot (2,2,3)
%figure
%hold on
plot(T,O,'b','Linewidth',2)
xticks ([30 60 90 120 150 180  ])
xlabel('Time (days)')
ylabel('Concentration(gr/mL)')
title ('O')

subplot (2,2,4)
%figure
%hold on
plot(T,P,'m','Linewidth',3)
xticks ([30 60 90 120 150 180  ])
xlabel('Time (days)')
ylabel('Concentration(gr/mL)')
title ('P')

And Then I would combine with the drug and immune system (because there are is couple equation)

This is my code Already combine for drug and system imune

clc;clear;

%parameter
  delta=50;
  gamma=75;
  K1=10^-4;
  K2=5*10^-4;
  K3=10^-3;
  K4=5*10^-3;
  K5=10^-2;
  K6=5*10^-2;
  Ko=0.1;
  n=6;
  Oa=10;
  Pa=100;
  mu_1=10^-3;
  mu_2=10^-3;
  mu_3=10^-3;
  mu_4=10^-3;
  mu_5=10^-3;
  mu_6=10^-3;
  mu_o=10^-4;
  mu_p= 10^-5;
  
  %input for time
  t(1)=0;
  dt=0.001; %time interval
  t=0:dt:180; %time span

%component for SYSTEM IMMUNE
Mg = 0.047;
lambdaMf = 5*10^-2; 
Fo = 3.36*10^-11;
KFo = 2.58*10^-11;
lambdaMa = 2.23*10^-2;
Ao = 0.14;
KAo = 10^-7;
beta = 10;
epsilon1 = 0.3333;
epsilon2 = 0.8;
lambdaM1Tb = 6*10^-3; 
Tb=10^-6;
KTb = 2.5*10^-7;
makro1 = 0.02;
makro2 = 0.02;
dM1= 0.015;
dM2= 0.015;
tmakro1 = Mg*(lambdaMf*(Fo/(Fo+KFo))+lambdaMa*(Ao/(Ao+KAo)))*(epsilon1*beta/beta*(epsilon1+epsilon2))-(lambdaM1Tb*(Tb/(Tb+KTb)*makro1))-(dM1*makro1);
tmakro2 = Mg*(lambdaMf*(Fo/(Fo+KFo))+lambdaMa*(Ao/(Ao+KAo)))*(epsilon2/beta*(epsilon1+epsilon2))+(lambdaM1Tb*(Tb/(Tb+KTb))*makro1)-(dM2*makro2);

%component DRUGS
Abi = 10^-6;
lambdaN = 8*10^-4; 
lambdaA = 8*10^-4; 
daboM = 2*10^-3; 
dabom = 10^-2; 
makrofag1 = 0;
makrofag2 = 0;
mikro1 = 0.02;
mikro2 = 0.02;
teta = 0.9; 
h = 10;
Kaob = 7*10^-3;
AoB = 10^-8;
dtdab = Abi*(dt*(0.14*1.9*10^-4*exp(-1.9*10^-4)))+ (lambdaN*1)+(lambdaA*1)-(dabom*(makrofag1+(teta*makrofag2))+daboM*(mikro1+(teta*mikro2))*(1+h))*(AoB/(AoB+Kaob));

%initial condition for SYSTEM IMMUNE
M1(1)=10;
M2(1)=0;
M3(1)=0;
M4(1)=0;
M5(1)=0;
M6(1)=0;
O(1)=0;
P(1)=0;

%initial condition FOR NORMAL
M12(1)=10;
M22(1)=0;
M32(1)=0;
M42(1)=0;
M52(1)=0;
M62(1)=0;
O2(1)=0;
P2(1)=0;

%initial condition with DRUGS
M14(1)=10;
M24(1)=0;
M34(1)=0;
M44(1)=0;
M54(1)=0;
M64(1)=0;
O4(1)=0;
P4(1)=0;

%empty array FOR SYSTEM IMMUNE
T=zeros(length(t)+1,1); %empty array for t
M1=zeros(length(t)+1,1); %empty array for M1
M2=zeros(length(t)+1,1); %empty array for M2
M3=zeros(length(t)+1,1); %empty array for M3
M4=zeros(length(t)+1,1); %empty array for M4
M5=zeros(length(t)+1,1); %empty array for M5
M6=zeros(length(t)+1,1); %empty array for M6
O=zeros(length(t)+1,1);  %empty array for O
P=zeros(length(t)+1,1);  %empty array for P

sumter=K2*M2+K3*M3+K4*M4+K5*M5;

%empty array FOR NORMAL
M12=zeros(length(t)+1,1); %empty array for M1
M22=zeros(length(t)+1,1); %empty array for M2
M32=zeros(length(t)+1,1); %empty array for M3
M42=zeros(length(t)+1,1); %empty array for M4
M52=zeros(length(t)+1,1); %empty array for M5
M62=zeros(length(t)+1,1); %empty array for M6
O2=zeros(length(t)+1,1); %empty array for O
P2=zeros(length(t)+1,1); %empty array for P

sumter2=K2*M22+K3*M32+K4*M42+K5*M52;

%empty array FOR DRUGS
M14=zeros(length(t)+1,1); %empty array for M1
M24=zeros(length(t)+1,1); %empty array for M2
M34=zeros(length(t)+1,1); %empty array for M3
M44=zeros(length(t)+1,1); %empty array for M4
M54=zeros(length(t)+1,1); %empty array for M5
M64=zeros(length(t)+1,1); %empty array for M6
O4=zeros(length(t)+1,1); %empty array for O
P4=zeros(length(t)+1,1); %empty array for P

sumter2=K2*M24+K4*M34+K4*M44+K5*M54;

for j = 1:length(t)
    %with SYSTEM IMMUNE
    T(j+1)=T(j)+dt;
    M1(j+1)=M1(j)+1./(1+exp(-T(j)));
    
    M1(j+1) = M1(j)+(dt*(delta*M1(j+1)*(1-(M1(j+1)/gamma))-2*K1*M1(j+1)*M1(j+1)-M1(j+1)*sumter(j+1))-((Oa-n)*K6*M1(j+1)*M6(j+1))-((Pa-Oa)*Ko*M1(j+1)*O(j+1))-(mu_1*M1(j+1)));
    M2(j+1) = M2(j)+(dt*(K1*M1(j)*M1(j)-K2*M1(j)*M2(j))-(mu_2*M2(j+1))-tmakro1-tmakro2);
    M3(j+1) = M3(j)+(dt*(K2*M1(j)*M2(j)-K3*M1(j)*M3(j))-mu_3*M3(j));
    M4(j+1) = M4(j)+(dt*(K3*M1(j)*M3(j)-K4*M1(j)*M4(j))-mu_4*M4(j));
    M5(j+1) = M5(j)+(dt*(K4*M1(j)*M4(j)-K5*M1(j)*M5(j))-mu_5*M5(j));
    M6(j+1) = M6(j)+(dt*(K5*M1(j)*M5(j)-K6*M1(j)*M6(j))-mu_6*M6(j));
    O(j+1)  = O(j)+(dt*(K6*M1(j)*M6(j)-Ko*M1(j)*O(j)-mu_o*O(j)));
    P(j+1)  = P(j)+(dt*(Ko*M1(j)*O(j)-mu_p*P(j)));
    
   
    %FOR NORMAL SYSTEM
    T(j+1)=T(j)+dt;
    M12(j+1) = M12(j)+1./(1+exp(-T(j)));
    
    M12(j+1) = M12(j)+(dt*(delta*M12(j+1)*(1-(M12(j+1)/gamma))-2*K1*M12(j+1)*M12(j+1)-M12(j+1)*sumter2(j+1))-((Oa-n)*K6*M12(j+1)*M62(j+1))-((Pa-Oa)*Ko*M12(j+1)*O2(j+1))-(mu_1*M12(j+1)));
    M22(j+1) = M22(j)+(dt*(K1*M12(j)*M12(j)-K2*M12(j)*M22(j))-(mu_2*M22(j+1)));
    M32(j+1) = M32(j)+(dt*(K2*M12(j)*M22(j)-K3*M12(j)*M32(j))-mu_3*M32(j));
    M42(j+1) = M42(j)+(dt*(K3*M12(j)*M32(j)-K4*M12(j)*M42(j))-mu_4*M42(j));
    M52(j+1) = M52(j)+(dt*(K4*M12(j)*M42(j)-K5*M12(j)*M52(j))-mu_5*M52(j));
    M62(j+1) = M62(j)+(dt*(K5*M12(j)*M52(j)-K6*M12(j)*M62(j))-mu_6*M62(j));
    O2(j+1)  = O2(j)+(dt*(K6*M12(j)*M62(j)-Ko*M12(j)*O2(j)-mu_o*O2(j)));
    P2(j+1)  = P2(j)+(dt*(Ko*M12(j)*O2(j)-mu_p*P2(j)));

%with DRUGS
    T(j+1)=T(j)+dt;
    M14(j+1) = M14(j)+1./(1+exp(-T(j)));
    
    M14(j+1) = M14(j)+(dt*(delta*M14(j+1)*(1-(M14(j+1)/gamma))-2*K1*M14(j+1)*M14(j+1)-M14(j+1)*sumter2(j+1))-((Oa-n)*K6*M14(j+1)*M64(j+1))-((Pa-Oa)*Ko*M14(j+1)*O4(j+1))-(mu_1*M14(j+1)));
    M24(j+1) = M24(j)+(dt*(K1*M14(j)*M14(j)-K2*M14(j)*M24(j))-(mu_2*M24(j+1))-tmakro1-tmakro2-dtdab);
    M34(j+1) = M34(j)+(dt*(K2*M14(j)*M24(j)-K3*M14(j)*M34(j))-mu_3*M34(j));
    M44(j+1) = M44(j)+(dt*(K3*M14(j)*M34(j)-K4*M14(j)*M44(j))-mu_4*M44(j));
    M54(j+1) = M54(j)+(dt*(K4*M14(j)*M44(j)-K5*M14(j)*M54(j))-mu_5*M54(j));
    M64(j+1) = M64(j)+(dt*(K5*M14(j)*M54(j)-K6*M14(j)*M64(j))-mu_6*M64(j));
    O4(j+1)  = O4(j)+(dt*(K6*M14(j)*M64(j)-Ko*M14(j)*O4(j)-mu_o*O4(j)));
    P4(j+1)  = P4(j)+(dt*(Ko*M14(j)*O4(j)-mu_p*P4(j)));

end

subplot (2,2,1)
plot(T,M1,'k', T,M12,'r',T,M14,'b','Linewidth',2)
legend ('M1 with immune', 'M1 without immune', 'M1 with drugs');
%xticks ([60 120 180 240 300  360 420 480 540 600 660 720 780 840 900 960 1020 1080])
xticks ([30 60 90 120 150 180  ])
xlabel('time (day)','Fontsize',14,'FontWeight','bold')
ylabel('concentration (gr/mL)','Fontsize',14,'FontWeight','bold')
title ('M1')

subplot (2,2,2)
plot(T,M2,'k',T,M22,'r',T,M24,'b','Linewidth',2)
legend ('M2 with immune', 'M2 without immune','M2 with drug');
xticks ([30 60 90 120 150 180  ])
yticks ([ 0 2 4 6 8 10 12 14 16])
xlabel('time (day)','Fontsize',14,'FontWeight','bold')
ylabel('concentration (gr/mL)','Fontsize',14,'FontWeight','bold')
title ('M2')

subplot (2,2,3)
plot(T,M3,'k',T,M32,'r',T,M34,'b','Linewidth',2)
legend ('M3 with immune', 'M3 without immune','M3 with drug');
xticks ([30 60 90 120 150 180  ])
yticks ([0 0.5 1 1.5 2 2.5 3 3.5])
xlabel('time (day)','Fontsize',14,'FontWeight','bold')
ylabel('concentration (gr/mL)','Fontsize',14,'FontWeight','bold')
title ('M3')

subplot (2,2,4)
plot(T,M4,'k',T,M42,'r',T,M44,'b','Linewidth',2)
legend ('M4 with immune', 'M4 without immune','M4 with drug');
xticks ([30 60 90 120 150 180  ])
yticks ([0 0.1 0.2 0.3 0.4 0.5])
xlabel('time (day)','Fontsize',14,'FontWeight','bold')
ylabel('concentration (gr/mL)','Fontsize',14,'FontWeight','bold')
title ('M4')

subplot (2,2,1)
plot(T,M5,'k',T,M52,'r',T,M54,'b','Linewidth',2)
legend ('M5 with immune', 'M5 without immune','M5 with drug');
xticks ([30 60 90 120 150 180  ])
yticks ([0 0.05 0.1 0.15 0.2])
xlabel('time (day)','Fontsize',14,'FontWeight','bold')
ylabel('concentration (gr/mL)','Fontsize',14,'FontWeight','bold')
title ("M5")

subplot (2,2,2)
plot(T,M6,'k',T,M62,'r',T,M64,'b','Linewidth',2)
legend ('M6 with immune', 'M6 without imune','M6 with drug');
xticks ([30 60 90 120 150 180  ])
yticks ([0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045])
xlabel('time (day)','Fontsize',14,'FontWeight','bold')
ylabel('concentration (gr/mL)','Fontsize',14,'FontWeight','bold')
title ("M6")

subplot (2,2,3)
plot(T,O,'k',T,O2,'r',T,O4,'b','Linewidth',2)
legend ('O with immune', 'O without immune','O with drug');
xticks ([30 60 90 120 150 180  ])
yticks ([0 0.005 0.01 0.015 0.02])
xlabel('time (day)','Fontsize',14,'FontWeight','bold')
ylabel('concentration (gr/mL)','Fontsize',14,'FontWeight','bold')
title ("O")

subplot(2,2,4)
plot(T,P,'k',T,P2,'r',T,P4,'b','Linewidth',2)
legend ('P with immune', 'P without immune','P with drug');
xticks ([30 60 90 120 150 180  ])
%yticks ([0 2 4 6 8 10 12])
xlabel('time (day)','Fontsize',14,'FontWeight','bold')
ylabel('concentration (gr/mL)','Fontsize',14,'FontWeight','bold')
title ("P")

[Sam Chak]

Hi @cindyawati cindyawati

I wish to see the amyloid β math equations like this (an excerpt from the Mathematical Model on Alzheimer’s disease). I'm not a coder by profession, so it is rather difficult for me to decipher the math equations from the code alone.

Furthermore, seeing the equations allows us to countercheck if the equations in your code are described correctly or not.

[cindyawati cindyawati]

This is my equation for Amyloid beta that I use. and for the drugs code you are right, that is my code that i use (but i only use ODE not ODP) @Sam Chak

Answer 1

[Sam Chak]

Hi @cindyawati cindyawati

Please check if the differential equations and the initial values are enterer correctly or not.

% Call ode45 solver
tspan  = [0 180];                   % integration time interval
x0     = [10 0 0 0 0 0 0 0];        % initial values
[t, x] = ode45(@drugODE, tspan, x0);

% Plot results
subplot(221)
plot(t, x(:,5)), grid on, title('M_{5}')
subplot(222)
plot(t, x(:,6)), grid on, title('M_{6}')
subplot(223)
plot(t, x(:,7)), grid on, title('O')
subplot(224)
plot(t, x(:,8)), grid on, title('P')

% Differential equations
function dxdt = drugODE(t, x)
    % definitions
    dxdt    = zeros(8, 1);
    M1      = x(1);
    M2      = x(2);
    M3      = x(3);
    M4      = x(4);
    M5      = x(5);
    M6      = x(6);
    O       = x(7);
    P       = x(8);
    
    % Parameters
    delta   = 50;
    gamma   = 75;
    K1      = 1e-4;
    K2      = 5e-4;
    K3      = 1e-3;
    K4      = 5e-3;
    K5      = 1e-2;
    K6      = 5e-2;
    Ko      = 0.1;
    n       = 6;
    Oa      = 10;
    Pa      = 100;
    mu1     = 1e-3;
    mu2     = 1e-3;
    mu3     = 1e-3;
    mu4     = 1e-3;
    mu5     = 1e-3;
    mu6     = 1e-3;
    muo     = 1e-4;
    mup     = 1e-5;
    
    % ODEs
    dxdt(1) = delta*M1*(1 - M1/gamma) - 2*K1*M1^2 - M1*(K2*M2 + K3*M3 + K4*M4 + K5*M5) - (Oa - n)*K6*M1*M6 - (Pa - Oa)*Ko*M1*O - mu1*M1;
    dxdt(2) = K1*M1*M1 - K2*M1*M2 - mu2*M2;
    dxdt(3) = K2*M1*M2 - K3*M1*M3 - mu3*M3;
    dxdt(4) = K3*M1*M3 - K4*M1*M4 - mu4*M4;
    dxdt(5) = K4*M1*M4 - K5*M1*M5 - mu5*M5;
    dxdt(6) = K5*M1*M5 - K6*M1*M6 - mu6*M6;
    dxdt(7) = K6*M1*M6 - Ko*M1*O  - muo*O;
    dxdt(8) = Ko*M1*O             - mup*P;
end

[cindyawati cindyawati]

Thank you @Sam Chak. If I don't use the ODE45 is possible? I mean using euler for the drug, I modify the amyloid beta equation and put drug equation using euler.

[Sam Chak]

Hi @cindyawati

It is possible to use the user-defined Euler (ode1) solver if the step size is made arbitrarily small. As mentioned earlier, the Euler method is known for propagating errors at every time step. With the built-in ode45 solver, my only concern is ensuring that I've entered the equations and parameters correctly. In contrast, with Euler, one needs to ensure the correctness of the algorithm, not to mention the absence of tools to verify the accuracy of the results.

[cindyawati cindyawati]

oke thank you for your suggestion @Sam Chak