Chemotherapy Scheduling

Mathematical Modeling

Richard Hitt

Introduction

In many cancer tumors, the cells can be divided into two groups:  proliferating cells and quiescent cells.  Cells can move from either class to the other.  The cells in the quiescent class do not reproduce wheras those in the proliferating class do.  Treatment (in the form of chemotherapy, say) can be administered periodically, but only affects proliferating cells.  We assume a single drug treatment destroys a constant fraction of the proliferating cells rather than a given number of cells (called the log cell kill hypothesis).

Acknowledgement

The material in this notebook was first explained to me by G. F.  Webb of Vanderbilt University.  He presented a very stimulating colloquium talk on this subject here a couple of years ago.  This material is based on my notes from the talk.

Linear Quiescence

We will begin with a simplified model of tumor cell population growth and treatment.  Let

p(t)  = population of proliferating cells at time t
q(t)  = polulation of quiescent cells at time t
β      = natural growth rate of tumor
σ     = transition rate from p class to q class
τ      = transition rate from q class to p class
μ      = treatment function with period c

The treatment function will be a pulse function  μ  of period  c  where the drug is administered for c/2 units time and withheld for the next c/2 units of time.  The drug will be therapeutic only while administered.

With the above assumptions, we have the following system of differential equations.

    p'[t] == (β - σ - μ[t,c])*p[t] + τ*q[t],
    q'[t] == σ*p[t] - τ*q[t],
     p[0] == (Subscript[p, 0]),
     q[0] == (Subscript[q, 0])

Can you see where all the terms come from?  What assumptions are built in to the system beyond what we have already mentioned?

An Example

We consider a specific example using the following constants.

\(β\ \  = 0.9; 
σ\ \  = 0.5; 
τ\ \  = 0.5; 

    p\_0 = 1.0; 
q\_0 = 4.0; \)

We build a pulse function variable period  c  for 31 periods (enough to get us through the following examples).

m[t_, c_] := 3 /; 0 ≤ t ≤ (c/2); m[t_, c_] := 0 /; t < 0; \(m[t_, c_] := 0 /; t > c\/2; 

    μ[t_, c_] := ∑\+\(n = 0\)\%30 m[t - n\ c, c]\)

Now for the differential equatons.  We will still leave  c  variable.

equations[c_]:={p'[t]==(β - σ - μ[t, c]) p[t] + τ q[t],q'[t]==σ p[t] - τ q[t],p[0] == (Subscript[p, 0]),q[0] == (Subscript[q, 0])}

Now we can let Mathematica solve the system.  Since we only want a graph in the end, we can relax the accuracy requirements so that the calculations don't go on forever.

sol1 = 
      NDSolve[equations[2], {p[t], q[t]}, {t, 0, 30}, 
        PrecisionGoal -> 2, WorkingPrecision -> 4]

Now we can plot the solutions for  p  and  q .  We are really interested in their sum, though, since that is the total population of the tumor.

Plot[Evaluate[Flatten[{p[t] + q[t], p[t], q[t]} /. sol1]], {t, 0, 30}, 
      PlotStyle 
        -> {RGBColor[0, 0, 0], RGBColor[1, 0, 0], RGBColor[0, 0, 1]}, 
      PlotRange -> {0, 5}]

So, under all the assumptions made, this looks like an effective treatment.  Unfortunately, we cannot predict what will happen once the chemothreapy is stopped.

Another Example

This is identical to the first example, except the constant c is 6 instead of 2.

sol2 = 
      NDSolve[equations[6], {p[t], q[t]}, {t, 0, 30}, 
        PrecisionGoal -> 2, WorkingPrecision -> 4]

Plot[Evaluate[Flatten[{p[t] + q[t], p[t], q[t]} /. sol2]], {t, 0, 30}, 
      PlotStyle 
        -> {RGBColor[0, 0, 0], RGBColor[1, 0, 0], RGBColor[0, 0, 1]}]

Another Example

Now for c=8.

sol3 = 
      NDSolve[equations[8], {p[t], q[t]}, {t, 0, 30}, 
        PrecisionGoal -> 2, WorkingPrecision -> 4]

Plot[Evaluate[Flatten[{p[t] + q[t], p[t], q[t]} /. sol3]], {t, 0, 30}, 
      PlotStyle 
        -> {RGBColor[0, 0, 0], RGBColor[1, 0, 0], RGBColor[0, 0, 1]}]

Yes, Another Example

Here is c=10.

sol4 = 
      NDSolve[equations[10], {p[t], q[t]}, {t, 0, 30}, 
        PrecisionGoal -> 2, WorkingPrecision -> 4]

Plot[Evaluate[Flatten[{p[t] + q[t], p[t], q[t]} /. sol4]], {t, 0, 30}, 
      PlotStyle 
        -> {RGBColor[0, 0, 0], RGBColor[1, 0, 0], RGBColor[0, 0, 1]}, 
      PlotRange -> {0, 15}]

Here we see tumor growth in spite of treatment.

Putting It All Together

Well, this could go on for a long time.  What we need is to assemble all of these curves into a unified whole to get the big picture.  The idea here is to compute a whole bunch of them and stack them together side-by-side.  This takes a little effort since we are dealing with interpolating functions instead of real formulae.  Here goes.

f[c_] := 
      NDSolve[equations[c], {p[t], q[t]}, {t, 0, 30}, 
        PrecisionGoal -> 2, WorkingPrecision -> 4, 
        MaxSteps -> 2000]

garray = {}; Do[{h[t_] = p[t] + q[t] /. f[c]; 
        AppendTo[garray, 
          Table[h[t][[1]], {t, 0, 
              30, .5}]]}, {c, 1, 10, .2}]

ListPlot3D[garray, 
      AxesLabel -> {"time", "period", "n(t)"}, 
      Axes -> True, MeshRange -> {{0, 30}, {1, 10}}, 
      ViewPoint -> {-1.2, -1., 1.}, PlotRange -> {0, 13}, 
      LightSources 
        -> {{{0, 1, 1}, RGBColor[1, 0, 0]}, {{-1, 1, 1}, 
            RGBColor[0, 1, 0]}, {{-1, 0, 1}, RGBColor[0, 0, 1]}}]

A Non-Linear Model of Quiescence

In the previous analysis, we assumed that the transition rates (tau and sigma) from one the proliferating class to the quiescent class were constant.  In particular, we assumed they did not depend on the size of the tumor.  In reality, those rates often do depend on the size of the tumor.  Proliferating cells have an increasing tendency to become quiescent as tumor size increases and quiescent cells have an increasing tendency to remain quiescent as tumor size increases.  This occurs because cells in the interior of the tumor have less access to nutrients and oxygen.

An Example

In this example, we will let σ and τ change over time to reflect the increasing tendency of the cells to remain or become quiescent.

\(Clear[σ, τ]; 
n[t_] := p[t] + q[t]\)σ[t_] := n[t]/(1 + n[t])τ[t_] := (1/(1 + n[t]))\(β = 0.9; 
p\_0 = 1.0; 
q\_0 = 2.0; \)

We set up a pulse function as before so that the period is  c .

m[t_, c_] := 5 /; 0 ≤ t ≤ (c/2); m[t_, c_] := 0 /; t < 0; \(m[t_, c_] := 0 /; t > c\/2; 

    μ[t_, c_] := ∑\+\(n = 0\)\%60 m[t - n\ c, c]\)

The system of differential equations is as before except that sigma and tau are now functions of time.

equations2[c_]:=
	  {p'[t]==(β - σ[t] - μ[t, c]) p[t] + 
              τ[t] q[t],
		  q'[t]==σ[t] p[t] - τ[t] q[t],
		  p[0] == (Subscript[p, 0]),
		  q[0] == (Subscript[q, 0])}

We can now construct the graph of the surface  n(t,c) as before.

f[c_] := 
      NDSolve[equations2[c], {p[t], q[t]}, {t, 0, 60}, 
        PrecisionGoal -> 2, WorkingPrecision -> 4, 
        MaxSteps -> 2000]

garray2 = {}; Do[{h[t_] = p[t] + q[t] /. f[c]; 
        AppendTo[garray2, 
          Table[h[t][[1]], {t, 0, 
              60, 2}]]}, {c, 1, 31, 1}]

ListPlot3D[garray2, 
      AxesLabel -> {"time", "period", "n(t)"}, 
      Axes -> True, MeshRange -> {{0, 60}, {1, 31}}, 
      ViewPoint -> {-1.2, -1., 1.}, PlotRange -> {0, 25}, 
      LightSources 
        -> {{{0, 1, 1}, RGBColor[1, 0, 0]}, {{-1, 1, 1}, 
            RGBColor[0, 1, 0]}, {{-1, 0, 1}, RGBColor[0, 0, 1]}}]

Exercise

In the linear model of quiescence, the shorter the period was, the better the outcome of the treatment.  In the non-linear model, this is not the case.  Graph the curves  p(t), q(t), and n(t)  for the non-linear model on the same set of axes for various values of  t .  Choose one c-value below 8.1, choose 8.1, choose one between 8.1 and 25.5, choose 25.5, and choose one above 25.5 .  Explain what is happening.