Plotting differential equations using ODEs with multiple initial conditions, trying to recreate a plot

I'm very stuck and very new to ODEs. Trying to recreate the plot below and I'm not completely sure what I'm doing wrong. So far I only have the plot with the singular blue line that's labeled "untreated".

Attached are the equations given and parameters.

Article Link with more information: https://iji.sums.ac.ir/article_48257_93803012f585b515b2ad0e32a601db7c.pdf

a = 4.3*10^-1; b = 2.17*10^-8; p = 2*10^2; g = 1*10^7;
m = 2*10^-2; q = 3.4*10^-10; r1 = 7.25*10^-15; r2 = 6.9*10^-15; r3L = 1.95*10^-12; r3C = 1.95*10^-12; k1 = 5*10^-7; j = 1.245*10^-2; k2 = 2.019*10^7;
              
Yinit = [1.5*10^7];
[t1a,y1a]=ode45(@TumorCells,[0 70],Yinit,[],a,b,p,g);
plot(t1a,9*y1a)
[t1b, y1b] = ode45(@TumorCells,[0 70],Yinit,[],m,q,r1,r2,r3L,r3C,k1,j,k2);
plot(t1b,9*y1b)
ylabel('Tumor (Cells)');
xlabel('Time(day)');
legend('untreated', 'CpG 1', 'CpG 2', 'CpG 3')

function dTdt = TumorCells(t,T,a,b,p,g,m,q,r1,r2,r3L,r3C,k1,j,k2)
a = 4.3*10^-1; b = 2.17*10^-8; p = 2*10^2; g = 1*10^7;
m = 2*10^-2; q = 3.4*10^-10; r1 = 7.25*10^-15; r2 = 6.9*10^-15; r3L = 1.95*10^-12; r3C = 1.95*10^-12; k1 = 5*10^-7; j = 1.245*10^-2; k2 = 2.019*10^7;

E=1;Th1=1;Treg=1;DCC=1;DCL=1;
dTdt = a*T*(1-b*T)-[(p*E*T)/(g+T)];
dEdt = (-m*E) - (q*E*T) + (r1*Th1*E*T) - (r2*Treg*E*T) + (r3L*DCL+r3C*DCC)*[(E*T)/(1+k1*T)]+[j*(T*E)/(k2+T)];
y1a = [dTdt;dEdt];
y1b = [dTdt;dEdt];
end

Answer 1

[Davide Masiello]

I could not run the code, because several constants are missing (e.g. alpha1,beta1 etc.)

The error in your code is that you have not written the entire system of ODEs.

Something like the code below should work when completed.

If you share all the necessary info, I could try to run it and see if it works.

clear,clc

% Constants
a   = 4.3*10^-1; 
b   = 2.17*10^-8; 
p   = 2*10^2; 
g   = 1*10^7;
m   = 2*10^-2; 
q   = 3.4*10^-10; 
r1  = 7.25*10^-15; 
r2  = 6.9*10^-15; 
r3L = 1.95*10^-12; 
r3C = 1.95*10^-12; 
k1  = 5*10^-7; 
j   = 1.245*10^-2; 
k2  = 2.019*10^7;

              
Yinit = ones(1,5)*1.5*10^7;
[t,y]=ode45(@(t,y)TumorCells(t,y,a,b,p,g,m,q,r1,r2,r3L,r3C,k1,j,k2),[0 70],Yinit);
plot(t,y)
ylabel('Tumor (Cells)');
xlabel('Time(day)');
legend('untreated', 'CpG 1', 'CpG 2', 'CpG 3')

function dydt = TumorCells(t,y,a,b,p,g,m,q,r1,r2,r3L,r3C,k1,j,k2)
T       = y(1);
E       = y(2);
Th1     = y(3);
Treg    = y(4);
DCL     = y(5);

dTdt    = a*T*(1-b*T)-(p*E*T)/(g+T);
dEdt    = -m*E-q*E*T+r1*Th1*E*T-r2*Treg*E*T+(r3L*DCL+r3C*DCC)*((E*T)/(1+k1*T))+j*(T*E)/(k2+T);
dThqdt  = % expression
dTregdt = % expression
dDCLdt  = % expression
end

[Madeleine Yee]

Here are the rest of the parameteres:

[Davide Masiello]

Still it is not clear what the functions v_DCC(t) and v_DCL(t) appearing in the last two ODEs are.

[Madeleine Yee]

I couldn't find any information about those functions... unless I didn't look hard enouigh :/ article with the information: https://iji.sums.ac.ir/article_48257_93803012f585b515b2ad0e32a601db7c.pdf