restart
:

Digits:= 15: #best value for most purposes, IMO

local gamma: #Avoid potential conflict with Maple's predefined gamma constant.
#

#Numeric values of IVP parameters from the paper:
#
params:= [
# symbol    value    description                                         units
#------------------------------------------------------------------------------------
   r      = 1e-2,  #cancer growth rate                                   /day
   k      = 1e12,  #max possible cancer cells                            cell   
   c__1   = 5e-11, #interaction coefficient affecting cancer cells       /cell/day
   c__2   = 1e-13, #interaction coefficient affecting immune cells       /cell/day
   mu     = 8,     #death rate of cancer cells from chemo                /day
   lambda = 4.16,  #drug elimination rate                                /day             
   a      = 2e3,   #half-max drug to kill cancer                         mg
   s__0   = 3e5,   #normal immune cell influx                            cell/day
   d      = 1e-3,  #normal immune cell death rate                        /day
   rho    = 1e-12, #cancer-stimulated immune cell production rate        /day
   gamma  = 1e2,   #number cancer cells for half-max immune response     cell
   delta  = 1e4,   #death rate of immune cells from chemo                /day
   b      = 5e6,   #half-max drug to kill immune cells                   mg
   N__0   = 2e10,  #initial cancer cells                                 cell
   P__0   = 5e7,   #initial immune cells                                 cell
   #substitute symbolic parameters to represent constant input functions:
   s(t)   = s__inf,#immunotherapty rate                                  cell/day
   q(t)   = q__inf #chemotherapy rate                                    mg/day
]:
t__max   := 3e4:   #max integration time                                 day
N__inf   := 0.5:   #number of cancer cells remaining for patient to be    
                   #considered cured:                                     cell
#

#Complete ODE-IVP system:
ODEs:=
   #(Eq. 1) N(t) = number of cancer cells at time t (day):
   diff(N(t),t) = r*N(t)*(1 - N(t)/k) - c__1*N(t)*P(t) - mu*N(t)*Q(t)/(a + Q(t)),

   #(Eq. 2) P(t) = number of immume cells at time t (day):
   diff(P(t),t) =
      s(t) + s__0 - d*P(t) + rho*N(t)*P(t)/(gamma + N(t)) - c__2*N(t)*P(t) -
      delta*P(t)*Q(t)/(b + Q(t)),
   
   #(Eq. 3) = amount of chemo drug (mg) in body at time t (day):
   diff(Q(t),t) = q(t) - lambda*Q(t)
:
#Initial conditions: Q(0) = 0, always, because the total amount of chemo in the patient
#starts at 0:
ICs:= N(0) = N__0, P(0) = P__0, Q(0) = 0:
#

<ODEs>; #This line is just to display the ODEs prettyprinted.

dsol:= dsolve(
   eval({ODEs, ICs}, params), 'numeric',
   'parameters'= [s__inf, q__inf],
   #Halt at time of cure, i.e., number of cancer cells = N__inf.
   'events'= [[N(t)-N__inf, 'halt']],
   'abserr'= 1e-1, 'relerr'= 1e-13,
   #Using 'events', available methods are rkf45, ck45, and rosenbrock. Rosenbrock works
   #over a wider variety of parameter values for this problem.
   'method'= 'rosenbrock'
);

t_cure:= proc(q_inf::realcons)
local
   t__cure,
   wl:= interface('warnlevel'= 0) #Suppress 'events'-related warnings.
;
   dsol('parameters'= [q__inf= q_inf]);
   dsol('eventenable'= {1}); #Reinitialize integration (i.e., at t=0).
   t__cure:= eval(t, dsol(t__max));
   interface('warnlevel'= wl); #Restore its original value.
   t__cure
end proc
:  

#Set values of s__infinity (immunotherapy infusion rate) to use:
S_infs:= [0.0, 4e5, 7e5] #cell/day
:

#Using Fig. 5(a) from the paper as a guide, set the correspondence between
#s__infinity and the minimum value of q__infinity used.
q_min:= table(S_infs=~ [9e0, 9e0, 2e0]): #mg/day
q_max:= 1e5 #mg/day #It's the same for all plots.
:

(q_inf, time_cure, D_chemo):= (1,2,3): #columns in the Data matrix
#

st:= time(): #Track computation time (just for my curiousity).
for s_inf in S_infs do
   dsol('parameters'= [s__inf= s_inf]);
   #Semilogplot is used to get the correct logarithmic spacing for smaller q values:
   t_vs_q:= op([1,1], plots:-semilogplot(t_cure, q_min[s_inf]..q_max));
   #Store (as 3rd column) computation D__chemo = q__infinity * t__cure
   Data[s_inf]:= <t_vs_q | t_vs_q[.., q_inf] *~ t_vs_q[.., time_cure]>
od:
time()-st;

#formatting options common to both plots:
common_opts:=
   'axis'[1]= ['mode'= 'log', 'tickmarks'= [6, 'subticks'= 8]],
   'axis'[2]= ['mode'= 'log', 'tickmarks'= [4, 'subticks'= 8]],
   'labeldirections'= ['horizontal', 'vertical'],
   'view'= [1e0 .. 1e5, 1e3 .. 1e6],
   'legend'= [seq(s__infinity= s_inf, s_inf= S_infs)],
   'titlefont'= ["TIMES", 16],
   'axes'= 'boxed',
   'size'= [500, 500],
   'gridlines'= false
:

#Figure 5(a) (D__chemo vs q__infinity):
plot(
   [seq(Data[s_inf][..,[q_inf, D_chemo]], s_inf= S_infs)],
   common_opts,
   'labels'= [q__infinity*Unit('mg'/'day'), 'D__chemo'*Unit('mg')],   
   'title'= "Total chemo dose vs. Chemo infusion rate",
   'caption'=
      "\nFigure 5(a):"
      " Total chemo dose as a function of the chemo infusion rate\n"
);

#Figure 5(b) (D__chemo vs t__cure):
plot(
   [seq(Data[s_inf][..,[time_cure, D_chemo]], s_inf= S_infs)],
   common_opts,
   'labels'= [t__cure*Unit('day'), 'D__chemo'*Unit('mg')],
   'title'= "Chemo dose vs. Time for cure",
   'caption'=
      "\nFigure 5(b):"
      " Total chemo dose as an implicit function of the time required for cure\n"
);