How i use break in continues data in couple differential? it is possible?
3 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
cindyawati cindyawati
il 17 Mag 2024
Modificato: cindyawati cindyawati
il 20 Mag 2024
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")
9 Commenti
Sam Chak
il 18 Mag 2024
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)
tmacro2 = (Mg*(lambdaMf*(Fo/(Fo+KFo))+lambdaMa*(Ao/(Ao+KAo)))*(epsilon2/beta*(epsilon1+epsilon2)))+(lambdaM1Tb*(Tb/(Tb+KTb))*macro1)-(dM2*macro2)
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))
Risposta accettata
Sam Chak
il 19 Mag 2024
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
12 Commenti
Sam Chak
il 20 Mag 2024
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.
Più risposte (0)
Vedere anche
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!