How i use break in continues data in couple differential? it is possible?

Hello, this is my code for a diseases and I want to know the effect of drugs (blue line) if I stoped from days 30-60 and continue with microglia from days 60-180. It's use break code or anything else?

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.01; %time interval
t=0:dt:180; %time span
  
%component for microglia  ---> This is 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;
macro1 = 0.02;
macro2 = 0.02;
dM1= 0.015;
dM2= 0.015;
tmacro1 = (Mg*(lambdaMf*(Fo/(Fo+KFo))+lambdaMa*(Ao/(Ao+KAo)))*(epsilon1*beta/beta*(epsilon1+epsilon2)))-(lambdaM1Tb*(Tb/(Tb+KTb)*macro1))-(dM1*macro1);
tmacro2 = (Mg*(lambdaMf*(Fo/(Fo+KFo))+lambdaMa*(Ao/(Ao+KAo)))*(epsilon2/beta*(epsilon1+epsilon2)))+(lambdaM1Tb*(Tb/(Tb+KTb))*macro1)-(dM2*macro2);

%component drugs  --> This is component for drugs
Abi = 10^-6;
lambdaN = 8*10^-4; 
lambdaA = 8*10^-4; 
daboM = 2*10^-3; %clearance rate by microglia
dabom = 10^-2; %clearance rate by macrophages 
macrofag1 = 0;
macrofag2 = 0;
micro1 = 0.02;
micro2 = 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*(macrofag1+(teta*macrofag2))+daboM*(micro1+(teta*micro2))*(1+h))*(AoB/(AoB+Kaob));

%macro for macrofag and micro for microglia and each other is not same

%initial condition with microglia --> Initial condition for 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 drugs --> Initial condition for without drugs and 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 --> Intial 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 with microglia --> Array for 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 drugs --> Array for without Drugs and 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 --> 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 microglia (This component give the black line)
    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))-tmacro1-tmacro2);
    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 drugs (This component give the red line)
    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 (This component give the blue line)
    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))-tmacro1-tmacro2-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
 
figure
subplot (2,2,1)
plot( T,M12,'r',T,M1,'k',T,M14,'b','Linewidth',4) %M12 for without drugs, M1 for microglia, M14 for drugs
legend ('M1 without drugs','M1 with microglia', '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('days','Fontsize',14,'FontWeight','bold')
ylabel('(g/mL)','Fontsize',14,'FontWeight','bold')
title ('M1')

subplot (2,2,2)
plot(T,M22,'r',T,M2,'k',T,M24,'b','Linewidth',4) %M22 for without drugs, M2 for microglia, M24 for drugs
legend ('M2 without drugs','M2 with microglia','M2 with drugs');
xticks ([30 60 90 120 150 180  ])
yticks ([ 0 2 4 6 8 10 12 14 16])
xlabel('days','Fontsize',14,'FontWeight','bold')
ylabel('(g/mL)','Fontsize',14,'FontWeight','bold')
title ('M2')

subplot (2,2,3)
plot(T,M32,'r',T,M3,'k',T,M34,'b','Linewidth',4) %M32 for without drugs, M3 for microglia, M34 for drugs
legend ('M3 without drugs','M3 with microglia','M3 with drugs');
xticks ([30 60 90 120 150 180  ])
yticks ([0 0.5 1 1.5 2 2.5 3 3.5])
xlabel('days','Fontsize',14,'FontWeight','bold')
ylabel('(g/mL)','Fontsize',14,'FontWeight','bold')
title ('M3')

subplot (2,2,4)
plot(T,M42,'r',T,M4,'k',T,M44,'b','Linewidth',4) %M42 for without drugs, M4 for microglia, M44 for drugs
legend ('M4 without drugs','M4 with microglia', 'M4 with drugs');
xticks ([30 60 90 120 150 180  ])
yticks ([0 0.1 0.2 0.3 0.4 0.5])
xlabel('days','Fontsize',14,'FontWeight','bold')
ylabel('(g/mL)','Fontsize',14,'FontWeight','bold')
title ('M4')

figure
subplot (2,2,1)
plot(T,M52,'r',T,M5,'k',T,M54,'b','Linewidth',4) %M52 for without drugs, M1\5 for microglia, M54 for drugs
legend ('M5 without drugs','M5 with microglia', 'M5 with drugs');
xticks ([30 60 90 120 150 180  ])
yticks ([0 0.05 0.1 0.15 0.2])
xlabel('days','Fontsize',14,'FontWeight','bold')
ylabel('(g/mL)','Fontsize',14,'FontWeight','bold')
title ("M5")

subplot (2,2,2)
plot(T,M62,'r',T,M6,'k',T,M64,'b','Linewidth',4)%M62 for without drugs, M6 for microglia, M64 for drugs
legend ('M6 without drugs','M6 with microglia','M6 with drugs');
xticks ([30 60 90 120 150 180  ])
yticks ([0 0.01 0.02  0.03  0.04 ])
xlabel('days','Fontsize',14,'FontWeight','bold')
ylabel('(g/mL)','Fontsize',14,'FontWeight','bold')
title ("M6")

subplot (2,2,3)
plot(T,O2,'r',T,O,'k',T,O4,'b','Linewidth',4)%O2 for without drugs, O for microglia, O4 for drugs
legend ('O without drugs','O with microglia', 'O with drugs');
xticks ([30 60 90 120 150 180  ])
yticks ([0 0.005 0.01 0.015 0.02])
xlabel('days','Fontsize',14,'FontWeight','bold')
ylabel('(g/mL)','Fontsize',14,'FontWeight','bold')
title ("O")

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

[Sam Chak]

Hi @cindyawati cindyawati

Because your code is rather lengthy, I suggest that you annotate sections of the code so that users can easily identify which section of the drug delivery dynamics produces the blue curve.

Additionally, please specify which parameters describe the introduction of the "drug" and "microglia" in the drug delivery system.

By the way, is ode45 solver allowed to solve drug delivery question?

[cindyawati cindyawati]

Thank u for your advice @Sam Chak. I already edited and I hope that code can be understand. I suggest to using Euler for solve drug delivey but if using Euler is hard, I think use ode45 can be allow. But the main is to show the effect if the drug has been stopped from 0-60 days and from 60-180 continue only use microglia. Thx

[cindyawati cindyawati]

I'm already try using while loop but it's nothing happen. I dont know what is mistake. Can u help me @Sam Chak

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.01; %time interval
t=0:dt:180; %time span
t_end=10; %end time
  
%component for microglia  ---> This is 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;
macro1 = 0.02;
macro2 = 0.02;
dM1= 0.015;
dM2= 0.015;
tmacro1 = (Mg*(lambdaMf*(Fo/(Fo+KFo))+lambdaMa*(Ao/(Ao+KAo)))*(epsilon1*beta/beta*(epsilon1+epsilon2)))-(lambdaM1Tb*(Tb/(Tb+KTb)*macro1))-(dM1*macro1);
tmacro2 = (Mg*(lambdaMf*(Fo/(Fo+KFo))+lambdaMa*(Ao/(Ao+KAo)))*(epsilon2/beta*(epsilon1+epsilon2)))+(lambdaM1Tb*(Tb/(Tb+KTb))*macro1)-(dM2*macro2);

%component drugs  --> This is component for drugs
Abi = 10^-6;
lambdaN = 8*10^-4; 
lambdaA = 8*10^-4; 
daboM = 2*10^-3; %clearance rate by microglia
dabom = 10^-2; %clearance rate by macrophages 
macrofag1 = 0;
macrofag2 = 0;
micro1 = 0.02;
micro2 = 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*(macrofag1+(teta*macrofag2))+daboM*(micro1+(teta*micro2))*(1+h))*(AoB/(AoB+Kaob));

%macro for macrofag and micro for microglia and each other is not same

%initial condition with microglia --> Initial condition for 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 drugs --> Initial condition for without drugs and 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 --> Intial 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 with microglia --> Array for 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 drugs --> Array for without Drugs and 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 --> 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 microglia (This component give the black line)
    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))-tmacro1-tmacro2);
    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 drugs (This component give the red line)
    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 (This component give the blue line)
    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))-tmacro1-tmacro2-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)));

    while T(j+1)<t_end
        if M14(j+1)>M1(j+1)
            break;
        end
        if T(j+1)<t_end
            continue;
        end 
    end
end

 
figure
subplot (2,2,1)
plot( T,M12,'r',T,M1,'k',T,M14,'b','Linewidth',4) %M12 for without drugs, M1 for microglia, M14 for drugs
legend ('M1 without drugs','M1 with microglia', '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('days','Fontsize',14,'FontWeight','bold')
ylabel('(g/mL)','Fontsize',14,'FontWeight','bold')
title ('M1')

subplot (2,2,2)
plot(T,M22,'r',T,M2,'k',T,M24,'b','Linewidth',4) %M22 for without drugs, M2 for microglia, M24 for drugs
legend ('M2 without drugs','M2 with microglia','M2 with drugs');
xticks ([30 60 90 120 150 180  ])
yticks ([ 0 2 4 6 8 10 12 14 16])
xlabel('days','Fontsize',14,'FontWeight','bold')
ylabel('(g/mL)','Fontsize',14,'FontWeight','bold')
title ('M2')

subplot (2,2,3)
plot(T,M32,'r',T,M3,'k',T,M34,'b','Linewidth',4) %M32 for without drugs, M3 for microglia, M34 for drugs
legend ('M3 without drugs','M3 with microglia','M3 with drugs');
xticks ([30 60 90 120 150 180  ])
yticks ([0 0.5 1 1.5 2 2.5 3 3.5])
xlabel('days','Fontsize',14,'FontWeight','bold')
ylabel('(g/mL)','Fontsize',14,'FontWeight','bold')
title ('M3')

subplot (2,2,4)
plot(T,M42,'r',T,M4,'k',T,M44,'b','Linewidth',4) %M42 for without drugs, M4 for microglia, M44 for drugs
legend ('M4 without drugs','M4 with microglia', 'M4 with drugs');
xticks ([30 60 90 120 150 180  ])
yticks ([0 0.1 0.2 0.3 0.4 0.5])
xlabel('days','Fontsize',14,'FontWeight','bold')
ylabel('(g/mL)','Fontsize',14,'FontWeight','bold')
title ('M4')
figure
subplot (2,2,1)
plot(T,M52,'r',T,M5,'k',T,M54,'b','Linewidth',4) %M52 for without drugs, M1\5 for microglia, M54 for drugs
legend ('M5 without drugs','M5 with microglia', 'M5 with drugs');
xticks ([30 60 90 120 150 180  ])
yticks ([0 0.05 0.1 0.15 0.2])
xlabel('days','Fontsize',14,'FontWeight','bold')
ylabel('(g/mL)','Fontsize',14,'FontWeight','bold')
title ("M5")

subplot (2,2,2)
plot(T,M62,'r',T,M6,'k',T,M64,'b','Linewidth',4)%M62 for without drugs, M6 for microglia, M64 for drugs
legend ('M6 without drugs','M6 with microglia','M6 with drugs');
xticks ([30 60 90 120 150 180  ])
yticks ([0 0.01 0.02  0.03  0.04 ])
xlabel('days','Fontsize',14,'FontWeight','bold')
ylabel('(g/mL)','Fontsize',14,'FontWeight','bold')
title ("M6")

subplot (2,2,3)
plot(T,O2,'r',T,O,'k',T,O4,'b','Linewidth',4)%O2 for without drugs, O for microglia, O4 for drugs
legend ('O without drugs','O with microglia', 'O with drugs');
xticks ([30 60 90 120 150 180  ])
yticks ([0 0.005 0.01 0.015 0.02])
xlabel('days','Fontsize',14,'FontWeight','bold')
ylabel('(g/mL)','Fontsize',14,'FontWeight','bold')
title ("O")

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

[Sam Chak]

Hi Cindy,

Can you provide the differential equations of the drug delivery system and the input variables for the drug and microglia?

It would be easier for me to read the math than the code. An image would be enough, no need to type it out in LaTeX format. Once the drug and microglia variables are identified, I can likely introduce a time-dependent state to specify when to cut off the drug intake and when to activate the microglia.

[cindyawati cindyawati]

Thank u before @Sam Chak maybe this eq you can help me.

I have 4 eq. The first is for the normally disease

and then, I input the microglia to the eq (1)

and the last, I input the drug eq to the eq 1.

So, in the eq 1 have eq microglia and drug (combine drug and microglia). The code is work you can see the result on top but the problem is how to show the effect of the drug when i stoped.

[Sam Chak]

Hi @cindyawati cindyawati

I have simulated the drug delivery system using the ode45 solver. However, I still cannot identify which input variables represent the drug and the microglia. I believe the system is incomplete, as I can see there are still 3 differential equations yet to be included in the system below. You have variables like tmacro1, tmacro2, dtdab, and I believe they are related to these 3 differential equations.

Can you please include them in the drugDeliverySystem(t, x) function?

dxdt(9) = ...

dxdt(10)= ...

dxdt(11)= ...

tspan   = [0, 180];                                 % simulation time
M10     = 10;
M20     = 0;
M30     = 0;
M40     = 0;
M50     = 0;
M60     = 0;
O0      = 0;
P0      = 0;
x0      = [M10; M20; M30; M40; M50; M60; O0; P0];   % initial condition
[t, x]  = ode45(@drugDeliverySystem, tspan, x0);    % solve the system
plot(t, x), grid on, xlim([0, 180])                 % plot results
xlabel('Days')

function dxdt = drugDeliverySystem(t, x)
    dxdt    = zeros(8, 1);
    
    % definitions
    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      = 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;
    mu1     = 10^-3;
    mu2     = 10^-3;
    mu3     = 10^-3;
    mu4     = 10^-3;
    mu5     = 10^-3;
    mu6     = 10^-3;
    muo     = 10^-4;
    mup     = 10^-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^2  - K2*M1*M2 - mu2*M2; % - tmacro1 - tmacro2 - dtdab
    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 for your response @Sam Chak

The drug component is here

dtdab = Abi*(dt*(0.14*1.9*10^-4*exp(-1.9*10^-4)))+ (lambdaN*1)+(lambdaA*1)-(dabom*(macrofag1+(teta*macrofag2))+daboM*(micro1+(teta*micro2))*(1+h))*(AoB/(AoB+Kaob));

and the microglia (include pro-inflamation and anti-inflamation (see in the picture before)) is representent by

tmacro1 = (Mg*(lambdaMf*(Fo/(Fo+KFo))+lambdaMa*(Ao/(Ao+KAo)))*(epsilon1*beta/beta*(epsilon1+epsilon2)))-(lambdaM1Tb*(Tb/(Tb+KTb)*macro1))-(dM1*macro1);
tmacro2 = (Mg*(lambdaMf*(Fo/(Fo+KFo))+lambdaMa*(Ao/(Ao+KAo)))*(epsilon2/beta*(epsilon1+epsilon2)))+(lambdaM1Tb*(Tb/(Tb+KTb))*macro1)-(dM2*macro2);

The input of drug and microglia like only the number

[Sam Chak]

Hi @cindyawati cindyawati

Thanks for your clarification.

Now, I understand you want the drug (dtdab) to be administered from Day 1 to 30, then stop the drug from Day 31 to 60, and bring it back from Day 61 to 180. As for the microglia (tmacro1 and tmacro2), they are only activated from Day 61 to 180. Please confirm if I have understood the requirements correctly.

tmacro1 = (Mg*(lambdaMf*(Fo/(Fo+KFo))+lambdaMa*(Ao/(Ao+KAo)))*(epsilon1*beta/beta*(epsilon1+epsilon2)))-(lambdaM1Tb*(Tb/(Tb+KTb)*macro1))-(dM1*macro1)

tmacro1 = 5.0201e-04

tmacro2 = (Mg*(lambdaMf*(Fo/(Fo+KFo))+lambdaMa*(Ao/(Ao+KAo)))*(epsilon2/beta*(epsilon1+epsilon2)))+(lambdaM1Tb*(Tb/(Tb+KTb))*macro1)-(dM2*macro2)

tmacro2 = 1.1544e-05

dtdab = Abi*(dt*(0.14*1.9*10^-4*exp(-1.9*10^-4)))+ (lambdaN*1)+(lambdaA*1)-(dabom*(macrofag1+(teta*macrofag2))+daboM*(micro1+(teta*micro2))*(1+h))*(AoB/(AoB+Kaob))

dtdab = 0.0016

[cindyawati cindyawati]

hi @Sam Chak yes you right

Answer 1

[Sam Chak]

Hi @cindyawati cindyawati

Please review the plots and verify if the microglia and the drugs are correctly injected in your Drug Delivery System. In the code, you are free to set the time to stop the drug and activate the microglia to investigate their effects on the system.

%% call ode45 solver
tspan   = [0, 180];                                 % simulation time
M10     = 10;
M20     = 0;
M30     = 0;
M40     = 0;
M50     = 0;
M60     = 0;
O0      = 0;
P0      = 0;
x0      = [M10; M20; M30; M40; M50; M60; O0; P0];   % initial condition
[t, x]  = ode45(@drugDeliverySystem, tspan, x0);    % solve the system

%% get data
M1      = x(:,1);
M2      = x(:,2);
M3      = x(:,3);
M4      = x(:,4);
M5      = x(:,5);
M6      = x(:,6);
O       = x(:,7);
P       = x(:,8);

tmacro1 = zeros(1, numel(t));
tmacro2 = zeros(1, numel(t));
dtdab   = zeros(1, numel(t));
for j = 1:numel(t)
    [~, tmacro1(j), tmacro2(j), dtdab(j)] = drugDeliverySystem(t(j), x(j,:).');
end

%% plot results
figure(1)
tL1     = tiledlayout(3, 2, 'TileSpacing', 'Compact');
xlabel(tL1, 'days',   'Fontsize', 14, 'FontWeight', 'bold')
ylabel(tL1, '(g/mL)', 'Fontsize', 14, 'FontWeight', 'bold')

nexttile;
plot(t, M1, 'Linewidth', 4, 'color', '#542485'), title ("M1")
grid on, xlim([0 180]), xticks ([30 60 90 120 150 180])

nexttile;
plot(t, M2, 'Linewidth', 4, 'color', '#8B47B5'), title ("M2")
grid on, xlim([0 180]), xticks ([30 60 90 120 150 180])

nexttile;
plot(t, M3, 'Linewidth', 4, 'color', '#5F2993'), title ("M3")
grid on, xlim([0 180]), xticks ([30 60 90 120 150 180])

nexttile;
plot(t, M4, 'Linewidth', 4, 'color', '#A855D4'), title ("M4")
grid on, xlim([0 180]), xticks ([30 60 90 120 150 180])

nexttile;
plot(t, M5, 'Linewidth', 4, 'color', '#6D38A0'), title ("M5")
grid on, xlim([0 180]), xticks ([30 60 90 120 150 180])

nexttile;
plot(t, M6, 'Linewidth', 4, 'color', '#B472E8'), title ("M6")
grid on, xlim([0 180]), xticks ([30 60 90 120 150 180])

figure(2)
tL2     = tiledlayout(2, 2, 'TileSpacing', 'Compact');
xlabel(tL2, 'days',   'Fontsize', 14, 'FontWeight', 'bold')

nexttile;
plot(t, O,       'Linewidth', 4, 'color', '#ECB0B0'), title ("O")
grid on, xlim([0 180]), xticks ([30 60 90 120 150 180])
ylabel('(g/mL)', 'Fontsize', 10)

nexttile;
plot(t, P,       'Linewidth', 4, 'color', '#7B84AA'), title ("P")
grid on, xlim([0 180]), xticks ([30 60 90 120 150 180])
ylabel('(g/mL)', 'Fontsize', 10)

nexttile;
plot(t, tmacro1, 'Linewidth', 2, 'color', '#A49488'), hold on
plot(t, tmacro2, 'Linewidth', 2, 'color', '#E4D295'), hold off, title ("Microglia")
grid on, xlim([0 180]), xticks ([30 60 90 120 150 180]), ylim([-1e-4, 6e-4])
legend('tmacro1', 'tmacro2', 'location', 'east')

nexttile;
plot(t, dtdab,   'Linewidth', 2, 'color', '#71A8A1'), title ("Drugs")
grid on, xlim([0 180]), xticks ([30 60 90 120 150 180]), ylim([-0.5e-3, 2e-3])

%% Drug Delivery System
function [dxdt, tmacro1, tmacro2, dtdab] = drugDeliverySystem(t, x)
    dxdt    = zeros(8, 1);
    
    % definitions
    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          = 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;
    mu1         = 10^-3;
    mu2         = 10^-3;
    mu3         = 10^-3;
    mu4         = 10^-3;
    mu5         = 10^-3;
    mu6         = 10^-3;
    muo         = 10^-4;
    mup         = 10^-5;
    
    % 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;
    macro1      = 0.02;
    macro2      = 0.02;
    dM1         = 0.015;
    dM2         = 0.015;
    tmacro1_amt = (Mg*(lambdaMf*(Fo/(Fo+KFo))+lambdaMa*(Ao/(Ao+KAo)))*(epsilon1*beta/beta*(epsilon1+epsilon2)))-(lambdaM1Tb*(Tb/(Tb+KTb)*macro1))-(dM1*macro1);
    tmacro2_amt = (Mg*(lambdaMf*(Fo/(Fo+KFo))+lambdaMa*(Ao/(Ao+KAo)))*(epsilon2/beta*(epsilon1+epsilon2)))+(lambdaM1Tb*(Tb/(Tb+KTb))*macro1)-(dM2*macro2);
    
    tm_ON       = 60;       % time at which microglia are activated
    tmacro1     = tmacro1_amt*(tanh(46/5*(t - tm_ON)) + 1)/2;
    tmacro2     = tmacro2_amt*(tanh(46/5*(t - tm_ON)) + 1)/2;
    
    % drugs
    Abi         = 10^-6;
    lambdaN     = 8*10^-4; 
    lambdaA     = 8*10^-4; 
    daboM       = 2*10^-3;  % clearance rate by microglia
    dabom       = 10^-2;    % clearance rate by macrophages 
    macrofag1   = 0;
    macrofag2   = 0;
    micro1      = 0.02;
    micro2      = 0.02;
    teta        = 0.9; 
    h           = 10;
    Kaob        = 7*10^-3;
    AoB         = 10^-8;
    dt          = 0.01;     % time interval
    dtdab_dose  = Abi*(dt*(0.14*1.9*10^-4*exp(-1.9*10^-4)))+ (lambdaN*1)+(lambdaA*1)-(dabom*(macrofag1+(teta*macrofag2))+daboM*(micro1+(teta*micro2))*(1+h))*(AoB/(AoB+Kaob));
    
    td_OFF      = 30;       % time at which the drug is stopped
    td_ON       = 60;       % time at which the drug is brought back
    dtdab       = dtdab_dose*(1 - ((tanh(46/5*(t - td_OFF)) + 1)/2 - (tanh(46/5*(t - td_ON)) + 1)/2));
    
    % 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^2  - K2*M1*M2 - mu2*M2 - tmacro1 - tmacro2 - dtdab;
    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 @Sam Chak for your help but I want to confirm the code that you make. Is the drug for conformation M1, M2, O or just only drug? Because i want to compare the gcoderaph into the normal, with microglia, with drug, and drug if we stopped (like my graph in my code). I mean in 1 image have 4 (normaly, microglia, drug, and drug hase been stopped) I think the drug is not only drug but the drug link into the M1, M2 until P so we can compare it and see the difference. Because M1 until P define the concentration of the disease

Thank you

[Sam Chak]

Hi @cindyawati cindyawati

Thanks for your feedback. However, I'm not quite sure I understand the meaning of "the drug for conformation M1, M2, O or just only drug." Your original problem stated that you wanted to do something about the blue curve by stopping the drug and activating the microglia halfway into the simulation of 180 days.

To address this, I have rewritten the dynamics in the drugDeliverySystem() code to follow only the "blue line section" where the values of tmacro1, tmacro2, and dtdab are injected into 'dM2/dt'.

In fact, I didn't change the disease dynamics at all. What I did was merely adding 6 lines of code to produce the desired signals for tmacro1, tmacro2, and dtdab after getting your confirmation in the your comment.

Anyhow, I encourage you to review the code to check if any changes in the dynamics have been incorrectly made by me, and edit the code as necessary to repair it and compare it with the other three cases (1. normal without both, 2. only microglia, 3. only drug) in your original Euler solver-based code.

tm_ON   = 60;       % time at which microglia are activated
tmacro1 = tmacro1_amt*(tanh(46/5*(t - tm_ON)) + 1)/2;
tmacro2 = tmacro2_amt*(tanh(46/5*(t - tm_ON)) + 1)/2;

td_OFF  = 30;       % time at which the drug is stopped
td_ON   = 60;       % time at which the drug is brought back
dtdab   = dtdab_dose*(1 - ((tanh(46/5*(t - td_OFF)) + 1)/2 - (tanh(46/5*(t - td_ON)) + 1)/2));

[cindyawati cindyawati]

Hi @Sam Chak thanks for helping me. M1, M2, until P is the conformation for the plaque of the diseases. I mean the drug is connected to M1, M2 until P so we can see the kinetica when the drug stopped in M1, M2 until P. Because your work is true but I think the drugs not connected to the M1 until P. Anyway thanks and I will try it

[Sam Chak]

Hi @cindyawati cindyawati

It's no problem at all. However, I'm not an expert in this area. From your original dynamics, the drug is injected in dM2/dt. Is this correct, or do you wish to inject the drug from the first state equation dM1/dt all the way to the last state equation dP/dt?

If the time-controlled drug signal (dtdab) is correctly designed, then you can try injecting the drug where you think is most appropriate. Please let me know how you would like to adjust the drug delivery.

Thanks!

[cindyawati cindyawati]

Hi @Sam Chak its oke if you not the expert. Yeah you Right! I want to inject the drug from M1 until P and the result like my pic (see my pic I hope its clear) but I dont know in my code if I injected the drug in dM2/dt the dM1 until dP is connected and give the different result. I think because my equation is coupled. Maybe something wrong in my euler code. So, I give u the freedom to use the ODE45. Thank you so much

[cindyawati cindyawati]

Hi @Sam Chak I already try it in my Euler code only drug interup for M1 until M4 and I get some error. The error says "Unable to perform assignment because the left and right sides have a different number of elements'. I dont know if using ODE45 is give the result or not, i need your help for ODE45 thank you before

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.01; %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_amt = (Mg*(lambdaMf*(Fo/(Fo+KFo))+lambdaMa*(Ao/(Ao+KAo)))*(epsilon1*beta/beta*(epsilon1+epsilon2)))-(lambdaM1Tb*(Tb/(Tb+KTb)*makro1))-(dM1*makro1);
tmakro2_amt = (Mg*(lambdaMf*(Fo/(Fo+KFo))+lambdaMa*(Ao/(Ao+KAo)))*(epsilon2/beta*(epsilon1+epsilon2)))+(lambdaM1Tb*(Tb/(Tb+KTb))*makro1)-(dM2*makro2);

%tm_ON       = 60;       % time at which microglia are activated
%tmacro1     = tmakro1_amt*(tanh(46/5*(t - tm_ON)) + 1)/2;
%tmacro2     = tmakro2_amt*(tanh(46/5*(t - tm_ON)) + 1)/2;

%component drugs
Abi = 10^-6;
lambdaN = 8*10^-4; 
lambdaA = 8*10^-4; 
daboM = 2*10^-3; %clearance rate by microglia
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_dose = 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));

td_OFF      = 30;       % time at which the drug is stopped
td_ON       = 60;       % time at which the drug is brought back
dtdab       = dtdab_dose*(1 - ((tanh(46/5*(t - td_OFF)) + 1)/2 - (tanh(46/5*(t - td_ON)) + 1)/2));

%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;
M25(1)=0;
M35(1)=0;
M45(1)=0;
M55(1)=0;
M65(1)=0;
O5(1)=0;
P5(1)=0;

%initial condition with drug interupt 
M15(1)=10;
M25(1)=0;
M35(1)=0;
M45(1)=0;
M55(1)=0;
M65(1)=0;
O5(1)=0;
P5(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*M25+K4*M35+K4*M45+K5*M55;


%empty array with drugs interupt
M15=zeros(length(t)+1,1); %empty array for M1
M25=zeros(length(t)+1,1); %empty array for M2
M35=zeros(length(t)+1,1); %empty array for M3
M45=zeros(length(t)+1,1); %empty array for M4
M55=zeros(length(t)+1,1); %empty array for M5
M65=zeros(length(t)+1,1); %empty array for M6
O5=zeros(length(t)+1,1); %empty array for O
P5=zeros(length(t)+1,1); %empty array for P

sumter2=K2*M25+K4*M35+K4*M45+K5*M55;


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_amt-tmakro2_amt);
    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_amt-tmakro2_amt-dtdab_dose);
    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)));
    
    
    %with drugs depend interupt
    T(j+1)=T(j)+dt;
    M15(j+1) = M15(j)+1./(1+exp(-T(j)));
    
    M15(j+1) = M15(j)+(dt*(delta*M15(j+1)*(1-(M15(j+1)/gamma))-2*K1*M15(j+1)*M15(j+1)-M15(j+1)*sumter2(j+1))-((Oa-n)*K6*M15(j+1)*M65(j+1))-((Pa-Oa)*Ko*M15(j+1)*O5(j+1))-(mu_1*M15(j+1)));
    M25(j+1) = M25(j)+(dt*(K1*M15(j)*M15(j)-K2*M15(j)*M25(j))-(mu_2*M25(j+1))-tmakro1_amt-tmakro2_amt-dtdab);
    M35(j+1) = M35(j)+(dt*(K2*M15(j)*M25(j)-K3*M15(j)*M35(j))-mu_3*M35(j));
    M45(j+1) = M45(j)+(dt*(K3*M15(j)*M35(j)-K4*M15(j)*M45(j))-mu_4*M45(j));
    M55(j+1) = M55(j)+(dt*(K4*M15(j)*M45(j)-K5*M15(j)*M55(j))-mu_5*M55(j));
    M65(j+1) = M65(j)+(dt*(K5*M15(j)*M55(j)-K6*M15(j)*M65(j))-mu_6*M65(j));
    O5(j+1)  = O5(j)+(dt*(K6*M15(j)*M65(j)-Ko*M15(j)*O5(j)-mu_o*O5(j)));
    P5(j+1)  = P5(j)+(dt*(Ko*M15(j)*O5(j)-mu_p*P5(j)));

end

Unable to perform assignment because the left and right sides have a different number of elements.

subplot (2,2,1)
plot( T,M12,'r',T,M1,'k',T,M14,'b',T,M15,'g','Linewidth',4)
legend ('M1 without drug','M1 with microglia', 'M1 with drug','M1 with drug interupt');
%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('Waktu (hari)','Fontsize',14,'FontWeight','bold')
ylabel('Konsentrasi (g/mL)','Fontsize',14,'FontWeight','bold')
title ('M1')

subplot (2,2,2)
plot(T,M22,'r',T,M2,'k',T,M24,'b',T,M25,'g','Linewidth',4)
legend ('M2 without drug','M2 with microglia', 'M2 with drug','M2 with drug interupt');
xticks ([30 60 90 120 150 180  ])
yticks ([ 0 2 4 6 8 10 12 14 16])
xlabel('Waktu (hari)','Fontsize',14,'FontWeight','bold')
ylabel('Konsentrasi (g/mL)','Fontsize',14,'FontWeight','bold')
title ('M2')

subplot (2,2,3)
plot(T,M32,'r',T,M3,'k',T,M35,'b',T,M35,'g','Linewidth',4)
legend ('M3 without drug','M3 with microglia', 'M3 with drug','M3 with drug interupt');
xticks ([30 60 90 120 150 180  ])
yticks ([0 0.5 1 1.5 2 2.5 3 3.5])
xlabel('Waktu (hari)','Fontsize',14,'FontWeight','bold')
ylabel('Konsentrasi (g/mL)','Fontsize',14,'FontWeight','bold')
title ('M3')

subplot (2,2,4)
plot(T,M42,'r',T,M4,'k',T,M45,'b',T,M45,'g','Linewidth',4)
legend ('M4 without drug','M4 with microglia', 'M4 with drug','M4 with drug interupt');
xticks ([30 60 90 120 150 180  ])
yticks ([0 0.1 0.2 0.3 0.4 0.5])
xlabel('Waktu (hari)','Fontsize',14,'FontWeight','bold')
ylabel('Konsentrasi (g/mL)','Fontsize',14,'FontWeight','bold')
title ('M4')

[Sam Chak]

Hi @cindyawati cindyawati

I revisited the problem and noticed that the state evolution-like law:

M15(j+1) = M15(j) + 1./(1+exp(-T(j)));

is not given in the differential equation set. No wonder the results look all the same.

[cindyawati cindyawati]

hi @Sam Chak Thank you, so I should remove that code?

[Sam Chak]

Hi @cindyawati cindyawati,

Don't get confused, no need to remove anything. I just want to replicate the results of Case 1 (normal without both microglia and drug), but I don't understand these lines which are related to the S-curve sigmoid, known as the Logistic function.

M1(j+1) = M1(j) + 1./(1+exp(-T(j)));

M12(j+1) = M12(j) + 1./(1+exp(-T(j)));

M14(j+1) = M14(j) + 1./(1+exp(-T(j)));

M15(j+1) = M15(j) + 1./(1+exp(-T(j)));

[cindyawati cindyawati]

Hi @Sam Chak which I do my code for normal condition is put the logistic function to get sigmoid graph because my model relates to the growth of the plaque. And I put tmacro1, tmacro2, and dtdab (drugs) for the microglia condition and drugs condition and I compare it so, I get 3 graph in 1 pic. I put tmacro1, tmacro2, and dtadab only in M2 because my eq is couple to each other (M2 to M3, M3 to M4,..). I dont put in M1 because the result give no effect for the microglia and drugs. The problem is I want to add one graph in the same picture when the drug is stopped. The code you are working is correct but I think its only for M2 and isnt couple to each other (M3, M4,...P). Sorry If my code are complicated

[Sam Chak]

Hi @cindyawati cindyawati

I'm actually not good at reading mathematical code with a lot of indexing. That's why I've preferred keeping things simple and intuitive using ode45. Am I correct in understanding that there are two differential equations for ? If not, please provide the full mathematical model so that I can better understand and know how to help you.

[cindyawati cindyawati]

hi @Sam Chak its oke if you using Ode45. Yeah you right for M1 base on my code and than for M2 until P

like a model that I give. Then, for microglia the tmacro1 and macro2 is include to the 1 equation there is "microglia" I want to input in the M2 because the M2 is couple to M3 until P and same for the drug. Thanks