\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{graphicx}

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2013 (2013), No. 272, pp. 1--14.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2013 Texas State University - San Marcos.}
\vspace{9mm}}

\begin{document}
\title[\hfilneg EJDE-2013/272\hfil Existence of periodic solutions]
{Existence of periodic solutions in the modified Wheldon model of CML}

\author[P. Amster, R. Balderrama, L. Idels \hfil EJDE-2013/272\hfilneg]
{Pablo Amster, Roc\'io Balderrama, Lev Idels}  % in alphabetical order

\address{Pablo Amster \newline
Departamento de Matem\'atica\\
FCEyN - Universidad de Buenos Aires  \& IMAS-CONICET \\
Ciudad Universitaria, Pab. I, 1428 Buenos Aires, Argentina}
\email{pamster@dm.uba.ar}

\address{Roc\'io Balderrama \newline
Departamento de Matem\'atica\\
FCEyN - Universidad de Buenos Aires  \& IMAS-CONICET \\
Ciudad Universitaria, Pab. I, 1428 Buenos Aires, Argentina}
\email{rbalde@dm.uba.ar}

\address{Lev Idels \newline
Department of Mathematics\\
Vancouver Island University (VIU) \\
 900 Fith St. Nanaimo BC Canada}
\email{lev.idels@viu.ca}

\thanks{Submitted  October 11, 2013. Published December 17, 2013.}
\subjclass[2000]{34K20, 92D25, 34K45, 34K12, 34K25}
\keywords{Nonlinear nonautonomous delay differential equation; \hfill\break\indent
 positive periodic solution; Leray-Schauder degree;
 chronic myelogenous leukemia; \hfill\break\indent model with pharmacokinetics}

\begin{abstract}
 The Wheldon model (1975) of a chronic myelogenous leukemia (CML) dynamics is
 modified and enriched by introduction of a time-varying microenvironment
 and time-dependent drug efficacies.  The resulting model is a special
 class of nonautonomous nonlinear system of differential equations with delays.
 Via topological methods, the existence of positive periodic solutions is proven.
 We introduce our main insight and formulate some relevant open problems
 and conjectures.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{remark}[theorem]{Remark}
\allowdisplaybreaks

\section{Modified Wheldon Model of CML}

\subsection{Background}
Chronic myelogenous leukemia (CML) is cancer of the blood in which too 
many granulocytes, a type of white blood cell, are produced in the marrow,
and it makes up about 10 to 15 percent of all leukemias 
(see, for example, \cite{Gold, Ho, Maa, Pe,Ro}).
In 1974  Wheldon in the paper \cite{We0} (see also \cite{Wek0})
introduced the following model of granulopoiesis (granulocyte production)
\begin{equation}
\label{Wheldon}
\begin{gathered}
\frac{dM}{dt}=\frac{\alpha}{1+\beta M^n(t-\tau)}
 -\frac{\lambda M(t)}{1+\mu B^m (t)},\\
\frac{dB}{dt}=-\omega B(t)+\frac{\lambda M(t)}{1+\mu B^m (t)},
\end{gathered}
\end{equation}
where all  parameters are positive constants.
In model \eqref{Wheldon},  $M(t)$ is the number of cells in the marrow;
$B(t)$ is the number of white blood cells;
$\beta$ is the coupling constant for cell production loop;
 $\alpha$ is the maximum rate of cell production;
$\lambda$ is the maximum rate of release of mature cells from marrow;
 $\mu$ is the coupling constant for  release loop;
 $\omega$ is the constant rate for loss of granulocytes from blood to tissue;
 $\tau$ represents mean time for  stem cell maturity;
 $n$ controls gain of cell production loop and
 $m$ controls gain of release loop.

A different mechanism of CML was modeled and studied by  Mackey (see, for example, \cite{Mac2}).
\begin{equation}\label{Macc}
\begin{gathered}
\frac{dN}{dt}=-\delta (t)N(t) - \beta( N(t)) +2e^{-\delta \tau}\beta( N(t-\tau)),\\
\frac{dP}{dt}=- \gamma P(t) +\beta ( N(t)) -e^{-\delta \tau}\beta( N(t-\tau)).
\end{gathered}
\end{equation}
This model consists of a proliferating phase cellular population $P(t)$ 
and a $G_0$ resting phase with a population of cells $N(t)$, where 
$$
\beta(N)=\frac{\beta_{0}\theta^{n}N}{\theta^{n}+N^{n}} ~ (n >0) .
$$
This is simply a model of stem cells dynamics - daughter cells either 
differentiate or return to the stem cell compartment follows by another 
division cycle. There is only the implicit suggestion above that there are 
positive and negative feedback signals regulating the rates at which cells
 will move through these ``decisions".


However,   model \eqref{Wheldon} has a major drawback, i.e., 
it describes a wrong mechanism. At the (unique)
nontrivial equilibrium point $(M_*,B_*)$
of system \eqref{Wheldon}, we have:
\begin{equation}\label{Fail}
\omega B_*=\frac{\alpha}{1+\beta M_*^n}.
\end{equation}
Thus, the $B$-population in the Wheldon model is
inversely proportional to the $M$-population; the latter does not
 have any biological explanation.

To reanimate the Wheldon model, we used Wheldon's remarks in his 
later work \cite{W1}  to introduce a new mechanism:
\begin{equation}\label{WheldonMod}
\begin{gathered}
\frac{dM}{dt}=\frac{\alpha M(t)}{1+\beta M^n(t-\tau_{1})}
 -\frac{\lambda M(t)}{1+\mu B^m (t-\tau_{2})},\\
\frac{dB}{dt}=-\omega B(t)+\frac{\lambda M(t)}{1+\mu B^m (t-\tau_{2})}.
\end{gathered}
\end{equation}

This model creates a time-delay loop triggering stem cell production and 
a fast loop regulating release of mature cells in the blood. Studies of 
the  model  imply that the oscillatory pattern in  leukemia  may be bring 
forth in two principal ways, either by an increased cell production rate 
or by an increased maturation time.
Note also that model \eqref{WheldonMod} assumes that there is a direct 
negative feedback from mature to the precursors of those cells. 
Time delay $\tau_{1}$ ($\tau$ in model \eqref{Wheldon}) represents a mean 
time for $M-$ cell maturity.
A stimulator/inhibitor mechanism is presented by the second term in both 
equations, where a time delay $\tau_{2}$
is a lag between when $B-$cells are initiated and when an apparent 
tumor progressed (the latency time) since each cell cycle phase is 
dependent on the completion of the previous ones.

\begin{remark} \label{rmk1.1} \rm
Note that the first term in \eqref{Wheldon} is a decreasing function of $M$
$$\frac{\alpha }{1+\beta M^n},$$
 whereas in model \eqref{WheldonMod}
$$
\frac{\alpha M }{1+\beta M^n}
$$ 
is a one-hump function, resulting in a
relationship between stem cells and white blood cells more realistic than
in \eqref{Fail}:
\begin{equation}\label{New}
\omega B_*=\frac{\alpha M_*}{1+\beta M^n_*}.
\end{equation}
\end{remark}

Exposure to chemoradiation therapy will kill not only cancer cells, but 
other rapidly dividing cells in the body as well (e.g. the cells in the 
bone marrow that go on to become white blood cells),
and will therefore suppress immune system  \cite{Be}--\cite{ Ca} 
\cite{ Gold,  Ma, Pe, Vap, We2}. Note that for a new model the 
complete recovery is possible for sufficiently high drug dosage 
(see Figure \ref{fig1}).

\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.8\textwidth, viewport=110 520 510 720,clip=true]{fig1a}
\\
\includegraphics[width=0.8\textwidth, viewport=110 520 510 725,clip=true]{fig1b}
\end{center}
\caption{Dynamics before therapy and after therapy}
\label{fig1}
\end{figure}

It is well recognized that tumor microenvironment changes with time and 
in response to treatment.
These fluctuations can modulate tumor progression and acquired treatment 
resistance (see, for example,  \cite{Fe,Ma, Maa, Vap}). Henceforth, 
to mimic changes of the tumor microenvironment, we incorporate time-dependent 
parameters.
\begin{equation}\label{Wheldon0}
\begin{gathered}
\frac{dM}{dt}=\frac{\alpha(t) M(t)}{1+\beta(t) M^n(t-\tau_{1})}
 -\frac{\lambda(t) M(t)}{1+\mu(t) B^m (t-\tau_{2})}-\delta p(t)M(t),\\
\frac{dB}{dt}=-\omega(t) B(t)+\frac{\lambda(t) M(t)}{1+\mu(t) B^m (t-\tau_{2})}
 - \delta q(t)B(t),
\end{gathered}
\end{equation}
where $p(t)=p(c)$ and $q(t)=q(c)$ are the varying effectiveness of the drug, 
and $c=c(t)$ is the drug concentration at time $t$.
Traditionally, this pharmokinetic is modeled by linear functions, namely
$p(c)=\alpha c(t)$ and $ g(c)=\beta c(t)$ where  $\alpha$ and $\beta$ are 
the appropriate drug sensitivity parameters. Clearly, $\alpha=\beta$ 
if the drugs are cycle-non-specific, i.e.,
they will be equally toxic to all types of cells.
Some types of chemotherapy can be modeled  based on a non-monotone one-humped 
functions- $p(c)=\alpha c(t) e^{-a c(t)}$ and  $q(c)=\beta c(t) e^{-b c(t)}$.
Throughout the paper, it shall be assumed that
$\alpha(t), \beta(t), \omega(t), \lambda(t), \mu(t), p(t)$ and $q(t)$ 
are continuous, positive and $T$-periodic functions and $\tau_{1,2} >0$ 
are fixed delays.
The parameter $\delta$ is assumed to be $1$ or $0$ according the presence 
or absence of pharmacokinetics.
Different and interesting models of CML were recently examined 
in \cite{Ba, Ef, Ho, Sw}.

It is worth noticing that, given set of nonnegative initial conditions,
the solution of problem \eqref{Wheldon0} is globally defined and positive 
over $[0,+\infty)$. Indeed,

\begin{theorem} \label{thm1.1}
Let $\varphi_{i}:[-\tau_i,0]\to [0,+\infty)$ be continuous functions such 
that $\varphi_i(0)>0$. Then there exists
a unique global positive solution of problem \eqref{Wheldon0}
under initial conditions
\begin{gather*}
M(t)=\varphi_{1}(t) \quad -\tau_{1}\leq t\leq 0,\\
B(t)=\varphi_{2}(t)\quad -\tau_{2}\leq t\leq 0.
\end{gather*}
\end{theorem}

\begin{proof}
Set $R(t):= \ln M(t)$, then the system
becomes
\begin{equation}\label{generalc2}
 \begin{gathered}
  R'(t)= \frac {\alpha(t) }{1+\beta(t) e^{n R(t-\tau_{1})}} 
- \frac{\lambda(t)}{1+\mu(t) B^m(t-\tau_{2})}-\delta p(t),\\
 B'(t)= -\omega(t) B(t)+ \frac{\lambda(t) e^{R(t)}}{1+\mu(t) B^m(t-\tau_{2})}
-\delta q(t)B(t).
\end{gathered}
\end{equation}
Suppose that $M(t)$ and $B(t)$ are defined and positive for $t<t_0$,
then from the inequalities $-\lambda(t)-\delta p(t)<R'(t)<\alpha(t)$
it is clear that $R(t)$ is defined up to $t_0$. Moreover, 
$B'(t)< \lambda e^{R(t)}$ and hence $B(t)$ is defined in $t_0$.
Finally, if $B(t_0)=0$ then $B'(t_0)>0$, a contradiction.
\end{proof}

In next section we shall prove, under appropriate conditions, 
the existence of at least one positive $T$-periodic solution: namely, a pair
$(M,B)$ of $C^1$ functions satisfying
$$
M(t+T)=M(t)>0, \quad B(t+T)=B(t)>0
$$
for all $t\in \mathbb{R}$.
In view of the preceding result, one might attempt to define a Poincar\'e-like
operator in order to apply some fixed point theorem. However, 
the conditions for such a procedure seem to be very
restrictive; thus we apply, instead, the Leray-Schauder degree theory  
\cite{Lo,Maw} over an appropriate open
subset of $C_T\times C_T$, where $C_T$ denotes the
space of continuous and $T$-periodic real functions.

For the reader's convenience, we make a short account of the main properties 
of the degree that shall be used in this work.
Let $X$ be a Banach space, let $\Omega\subset X$ be open and bounded and 
denote by $cl(\Omega)$ the closure of $\Omega$.
If $\mathcal K:cl(\Omega)\to X$ is compact with $\mathcal Ku\neq u$ for 
all $u\in \partial \Omega$, then the
Leray-Schauder degree of the Fredholm operator $\mathcal F=Id-\mathcal K$ 
at $0$ shall be denoted by
$\deg(\mathcal F,\Omega,0)$. Roughly speaking, this (whole) number can be 
regarded as an algebraic count of the zeros of $\mathcal F$.

(1) (Solution) If
 $\deg(\mathcal F,\Omega,0)\neq 0$, then $\mathcal F$ has at least one zero in $\Omega$.

(2) (Homotopy invariance)
 If  $\mathcal F_\sigma=Id-\mathcal K_\sigma$ with 
$\mathcal K_\sigma:cl(\Omega) \to X$ compact such that
 $\mathcal K_\sigma u\neq u$ for all $u\in \partial \Omega$, $\sigma\in [0,1]$ 
and $\mathcal K:cl(\Omega)\times [0,1]\to X$ given by
 $\mathcal K(u,\sigma):= \mathcal K_\sigma(u)$ continuous, then 
$\deg(\mathcal F_\sigma,\Omega,0)$ is independent on $\sigma$.

(3) If $\mathcal K(cl(\Omega))\subset V$, with
 $V\subset X$ a finite dimensional subspace, then
 $$
\deg(\mathcal F, \Omega,0)=  \deg(\mathcal F|_{cl(\Omega)\cap V}, \Omega\cap V,0).
$$
Identifying $V$ with $\mathbb{R}^n$, the latter term is simply the so-called 
{\sl Brouwer degree}.
In this paper, we only need to know that if $\Omega_0\subset \mathbb{R}^n$ 
is open and bounded with $0\in \Omega_0$, then
$\deg(-Id, \Omega_0, 0)=(-1)^n$.

\section{Existence of periodic solutions}
\subsection{Case 1: No pharmokinetic}

\begin{theorem} \label{thm2.1}
Assume that $\alpha(t), \beta(t), \lambda(t),\mu(t)$
and $\omega(t)$ are continuous, positive and $T$-periodic. Furthermore, 
assume that
\begin{enumerate}
 \item $n>\frac m{m+1}$.
\item $\alpha(t)>\lambda(t)>\omega(t)$ for all $t$.
\end{enumerate}
Then system \eqref{Wheldon0} with $\delta=0$
admits at least one positive $T$-periodic solution.
\end{theorem}

\begin{proof}
Set $u(t)=\ln M(t)$ and $v(t)=\ln B(t)$, then
\eqref{Wheldon0} with $\delta=0$ reads
\begin{gather*}
u'(t)= \frac {\alpha(t)}{1+\beta(t)e^{nu(t-\tau_1)}} -
 \frac {\lambda(t) }{1+\mu(t)e^{mv(t-\tau_2)}}:=\psi_1(u,v)(t), \\
v'(t)= -\omega(t)+
 \frac {\lambda(t)e^{u(t)-v(t)} }{1+\mu(t)e^{mv(t-\tau_2)}}:=\psi_2(u,v)(t).
\end{gather*}
To prove the existence of $T$-periodic solutions of this system,
we shall apply the
continuation method \cite{Maw}. Adapted to this case, the method guarantees 
the existence of solutions,
provided there exists an open bounded set $\Omega\subset C_T\times C_T$ such that
\begin{enumerate}
\item For $\sigma\in (0,1]$, the system
\begin{gather*}
u'(t)= \sigma\psi_1(u,v)(t), \\
v'(t)= \sigma\psi_2(u,v)(t)
\end{gather*}
has no $T$-periodic solutions on $\partial \Omega$.
\item $\deg(F, \Omega\cap \mathbb{R}^2,0)$ is well defined and different from $0$,
where the function $F:\mathbb{R}^2\to \mathbb{R}^2$ is defined by
$$
F(u,v):=\frac 1T \int_0^T\Big(
\frac {\alpha(t)}{1+\beta(t)e^{nu}} -
 \frac {\lambda(t) }{1+\mu(t)e^{mv}}, 
\frac {\lambda(t)e^{u-v}}{1+\mu(t)e^{mv}}-\omega(t)  \Big)\, dt.
$$
\end{enumerate}

For simplicity, we divide the proof in two steps.
\medskip

\noindent \textbf{First step:}
Let $\Omega_0:=(-R,R)\times (-R, c R)\subset \mathbb{R}^2$, where
$c$ is a fixed constant such that
$\frac 1{m+1} < c < \frac nm.$
We claim that $\deg(F,\Omega_0 ,0) =1$ for $R>0$ large enough.

Indeed, let us firstly assume that $-R\le v\le cR$, then
$$
F_1(R,v)=\frac 1{Te^{nR}}\int_0^T
\frac {\alpha(t)e^{nR}}{1+\beta(t)e^{nR}} -
 \frac {\lambda(t)e^{nR}}{1+\mu(t)e^{mv}}\, dt .
$$
As $nR > mcR$, it follows that $F_1(R,v)\le F_1(R,cR)
<0$ for $R\gg 0$. On the other hand,
$$
F_1(-R,v) = \frac 1{T} \int_0^T
\frac {\alpha(t)}{1+\beta(t)e^{-nR}} -
 \frac {\lambda(t)}{1+\mu(t)e^{mv}}\, dt \ge
\frac 1{T} \int_0^T
\frac {\alpha(t)}{1+\beta(t)e^{-nR}} \, dt -\overline \lambda.
$$
The right-hand side term tends to $\overline \alpha-\overline \lambda$ as 
$R\to +\infty$; thus, as $\alpha(t)>\lambda(t)$ for all $t$, we deduce that
$F_1(-R,v)>0$ for $R\gg 0$.

Next, assume that $|u|\le R$, and compute
$$
F_2(u, c R ) =-\overline \omega + \frac 1T
\int_0^T\frac {\lambda(t)e^{u-c R}}{1+\mu(t)e^{mcR}}\, dt\le
-\overline \omega + \int_0^T\frac {\lambda(t)e^{(1-c)R}}{1+\mu(t)e^{{mcR}}}\, dt\to -\overline\omega
$$
as $R\to+\infty$ since $c(m+1)>1$, and
$$
F_2(u,-R)=
-\overline \omega + \frac 1T
\int_0^T\frac {\lambda(t)e^{u+R}}{1+\mu(t)e^{-mR}}\, dt\ge
-\overline \omega + \frac 1T
\int_0^T\frac {\lambda(t)}{1+\mu(t)e^{-mR}}\, dt.
$$
Here, the right-hand side term tends to
$\overline \lambda-\overline  \omega$ as $R\to +\infty$. 
This quantity is positive since $\lambda(t)>\omega(t)$ for all $t$,
so we conclude that
$F_2(u,cR) < 0 < F_2(u,-R)$ for $R\gg 0$.
Thus, we may define the homotopy
$$
H(u,v,\sigma):= \sigma F(u,v) - (1-\sigma)(u,v),
$$
which does not vanish on $\partial \Omega_0$.
It follows that $\deg(F,\Omega_0,0)=\deg(-Id,\Omega_0,0)=(-1)^2=1$.

\begin{remark} \label{rmk2.1} \rm
As a consequence, it is deduced that $F$ vanishes in $\Omega_0$.
In particular, when $\alpha$, $\beta$, $\lambda$ and $\mu$  are positive 
constants we deduce that the system has a positive equilibrium, as it 
shall be proven in section \ref{equil} by direct computation.
\end{remark}


\noindent\textbf{Second step:}
Let
$$
\Omega:= \{ (u,v)\in C_T\times C_T: \| u\|_\infty<R, -R< v(t)< cR 
\text{ for all $t$}\}.
$$
We claim that if $R$ is large enough then the $T$-periodic solutions 
of the system
\begin{gather*}
u'(t)= \sigma\psi_1(u,v)(t), \\
v'(t)= \sigma\psi_2(u,v)(t)
\end{gather*}
with $0<\sigma\le 1$ do not belong to $\partial \Omega$.

Indeed, suppose firstly that $u_{\rm max}=R > \frac{v_{\rm max}}c$
and take $\xi\in [0,T]$
is such that $u_{\rm max}=u(\xi)$. From the first equation of the system
 we obtain
$$
\frac {\alpha(\xi)}{1+\beta(\xi)e^{nu(\xi-\tau_1)}} =
 \frac {\lambda(\xi)}{1+\mu(\xi)e^{mv(\xi-\tau_2)}} >
 \frac {\lambda(\xi)}{1+\mu(\xi)e^{{mcR}}}.
 $$
Moreover, observe that $u'(t) > -\lambda(t)$ for all $t$,
so by periodicity we deduce that
$$
u(\xi-\tau_1) -R \ge -\int_{\xi}^{kT+\xi-\tau_1}
\lambda(t)\, dt \ge -\int_{0}^{T} \lambda(t)\, dt:=-C_1
$$
where $k$ is the first natural number such that $kT>\tau_1$. It follows that
$$
\alpha(\xi) > \lambda(\xi)
\frac{1+\beta(\xi)e^{nu(\xi-\tau_1)}} {1+\mu(\xi)e^{{mcR}}} >
\lambda(\xi) \frac{1+\beta(\xi)e^{n(R-C_1)}} {1+\mu(\xi)e^{{mcR}}}.
$$
The right-hand side of this inequality tends uniformly to $+\infty$ as 
$R\to+\infty$.
Now assume that $v_{\rm max}=cR \ge cu_{\rm max}$, then take $\eta\in [0,T]$ 
such that $v(\eta)=v_{\rm max}$ and deduce,
from the second equation of the system:
$$
\omega(\eta)=
 \frac {\lambda(\eta)e^{u(\eta)-v(\eta)} }{1+\mu(\eta)e^{mv(\eta-\tau_2)}}\le
 \frac {\lambda(\eta)e^{(1-c)R}}{1+\mu(\eta)e^{mv(\eta-\tau_2)}}.
$$
As before, from the inequality $v'(t)\ge -\omega(t)$ it is seen that
$$
v(\eta-\tau_2) - cR  \ge -\int_{\eta}^{lT+\eta-\tau_2} \omega(t)\, dt 
\ge -\int_{0}^{T} \omega(t)\, dt:=-C_2,
$$
where $l$ is the first natural number such that $lT>\tau_2$. This implies
$$
\omega(\eta) \le \frac {\lambda(\eta) e^{(1-c)R}}{1+\mu(\eta)e^{m(cR-C_2)}}\to 0
$$
uniformly as $R\to +\infty$.
We conclude that $u_{\rm max}$ and $v_{\rm max}$ cannot be arbitrarily large.

Next, suppose that $u_{\rm min}=-R<v_{\rm min}$ and
$\xi\in [0,T]$ be such that $u_{\rm min}=u(\xi)$. As before,
$$
\frac {\alpha(\xi)}{1+\beta(\xi)e^{nu(\xi-\tau_1)}} =
 \frac {\lambda(\xi)}{1+\mu(\xi)e^{mv(\xi-\tau_2)}}<
 \frac {\lambda(\xi)}{1+\mu(\xi)e^{-mR}}
 $$
and hence
$$
\alpha(\xi) < \lambda(\xi) \frac{1+\beta(\xi)e^{nu(\xi-\tau_1)}} {1+\mu(\xi)e^{-mR}}.
$$
As $u(\xi-\tau_1) \le -R + \int_{\xi-\tau_1}^\xi \lambda(t)\, dt$, 
the right-hand side of the last inequality tends uniformly to $\lambda(\xi)$ 
as $R\to +\infty$.
In the same way, if
$v(\eta)=v_{\rm min}=-R\le u_{\rm min}$, then it is seen that
$$
\omega(\eta)\ge  \frac {\lambda(\eta)}{1+\mu(\eta)
e^{mv(\eta-\tau_2)}}\to \lambda(\eta)
$$
uniformly as $R\to +\infty$.
As
$\alpha(t)>\lambda(t)>\omega(t)$ for all $t$, we deduce that $R$ cannot 
be arbitrarily large and the claim is proven.
\end{proof}

\subsection{Case 2: With pharmokinetic}

\begin{theorem} \label{thm2.2}
Assume that $\alpha(t), \beta(t), \lambda(t),\mu(t),\omega(t), p(t)$ and $q(t)$ 
are positive and $T$-periodic. Furthermore, assume that:
$$
\alpha(t)-p(t) >\lambda(t) > \omega(t)+q(t)
$$
for all $t$.
Then system \eqref{Wheldon0} with $\delta=1$
admits at least one positive $T$-periodic solution.
\end{theorem}

\begin{proof}
We shall follow the general outline of the previous proof.
As before,
set $u(t)=\ln M(t)$ and $v(t)=\ln B(t)$,
then the model with $\delta =1$
reads
\begin{gather*}
u'(t)= \frac {\alpha(t)}{1+\beta(t)e^{nu(t-\tau_1)}} -
 \frac {\lambda(t) }{1+\mu(t)e^{mv(t-\tau_2)}} -p(t):=\psi_1^{p,q}(u,v)(t), \\
v'(t)= -\omega(t)+
 \frac {\lambda(t)e^{u(t)-v(t)} }{1+\mu(t)e^{mv(t-\tau_2)}} - q(t):=\psi_2^{p,q}(u,v)(t).
\end{gather*}
For the first step, let us consider now
$F^{p,q}:\mathbb{R}^2\to \mathbb{R}^2$ given by
$$
F^{p,q}(u,v):= F(u,v) - (\overline p,\overline q)
$$
with $F$ as in the previous proof.
First, assume that $|v|\le R$. Then
$$
F_1^{p,q}(R,v)=\frac 1{T}\int_0^T
\frac {\alpha(t)}{1+\beta(t)e^{nR}} -
 \frac {\lambda(t)}{1+\mu(t)e^{mv}}\, dt -\overline p <0
$$
for $R\gg 0$. On the other hand,
\begin{align*}
F_1^{p,q}(-R,v) 
&= \frac 1{T} \int_0^T \frac {\alpha(t)}{1+\beta(t)e^{-nR}} -
 \frac {\lambda(t)}{1+\mu(t)e^{mv}}\, dt -\overline p \\
&\ge \frac 1{T} \int_0^T
\frac {\alpha(t)}{1+\beta(t)e^{-nR}} \, dt -\overline {\lambda} - \overline {p}.
\end{align*}
The last term tends to $\overline \alpha-\overline \lambda-\overline {p}$ as 
$R\to +\infty$; thus, as $\alpha(t)>\lambda(t)+p(t)$ for all $t$, we deduce that
$F_1^{p,q}(-R,v)<0$ for $R\gg 0$.

Next, assume that $|u|\le R$ and compute
$$
F_2^{p,q}(u,R) 
\le  \int_0^T\frac {\lambda(t)}{1+\mu(t)e^{mR}}\, dt -\overline 
\omega -\overline q < 0
$$
for $R\gg 0$ and
$$
F_2^{p,q}(u,-R)\ge  \frac 1T
\int_0^T\frac {\lambda(t)}{1+\mu(t)e^{-mR}}\, dt-\overline \omega -\overline q .
$$
Here, the right-hand side term tends to
$\overline \lambda-\overline \omega -\overline q$ as $R\to +\infty$. 
This quantity is positive since $\lambda(t)>\omega(t)+q(t)$ for all $t$; 
so we conclude that $F_2^{p,q}(u,R) < 0 < F_2^{p,q}(u,-R)$ for $R\gg 0$. 
As in the previous proof, we have $\deg(F^{p,q},(-R,R)^2,0)=1$.

For the second step, set
$$
\Omega:= \{ (u(t),v(t))\in C_T\times C_T: \| u\|_\infty<R,\| v\|_\infty<R\}.
$$
As before, we claim that if $R$ is large enough then the $T$-periodic 
solutions of the system
\begin{gather*}
u'(t)= \sigma\psi^{p,q}_1(u,v)(t), \\
v'(t)= \sigma\psi^{p,q}_2(u,v)(t)
\end{gather*}
with $0<\sigma\le 1$ do not belong to $\partial \Omega$.
Indeed, suppose firstly that $u_{\rm max}=R > v_{\rm max}$, then take 
$\xi\in [0,T]$ is such that $u_{\rm max}=u(\xi)$ and
from the first equation we obtain
$$
\frac {\alpha(\xi)}{1+\beta(\xi)e^{nu(\xi-\tau_1)}} >
 \frac {\lambda(\xi)}{1+\mu(\xi)e^{mR}} + p(\xi).
 $$
As before, using now the fact that
$u'(t) > -\lambda(t) - p(t)$ for all $t$ we deduce that
$$
u(\xi-\tau_1) -R \ge -\int_{0}^{T} [\lambda(t)+p(t)]\, dt:=-C_1^{p,q}.
$$
It follows that
$$
\frac{\alpha(\xi)}{p(\xi)} >  1+\beta(\xi) e^{nu(\xi-\tau_1)} 
\ge 1+\beta(\xi) e^{n(R-C_1^{p,q})}
$$
and hence $R$ cannot be arbitrarily large.
On the other hand, assume that $u_{\rm max}\le v_{\rm max}=R$, 
then take $\eta\in [0,T]$ such that $v(\eta)=v_{\rm max}$ and deduce, 
from the second equation of the system, that
$$
\omega(\eta) + q(\eta)\le
 \frac {\lambda(\eta)}{1+\mu(\eta)e^{mv(\eta-\tau_2)}}
$$
and, from the inequality $v'(t)\ge -\omega(t)-q(t)$, that
$$
v(\eta-\tau_2) - R \ge -\int_{0}^{T} [\omega(t)+q(t)]\, dt:=-C_2^{p,q}.
$$
This implies
$$
\omega(\eta) +q(\eta)\le \frac {\lambda(\eta)}{1+\mu(\eta)
e^{m(R-C_2^{p,q})}}\to 0
$$
uniformly as $R\to +\infty$.
We conclude that $u_{\rm max}$ and $v_{\rm max}$ cannot be arbitrarily large.

Next, suppose that $u_{\rm min}=-R<v_{\rm min}$ and
$\xi\in [0,T]$ be such that $u_{\rm min}=u(\xi)$. As before, it follows that
$$
\frac {\alpha(\xi)}{1+\beta(\xi)e^{nu(\xi-\tau_1)}}  <
 \frac {\lambda(\xi)}{1+\mu(\xi)e^{-mR}}+p(\xi)
 $$
and hence
$$
\alpha(\xi) < \Big(\frac{\lambda(\xi)} {1+\mu(\xi)e^{-mR}}+p(\xi)\Big)
\Big(1+\beta(\xi)e^{nu(\xi-\tau_1)}\Big).
$$
Thus, the right-hand side of the last inequality tends uniformly to
 $\lambda(\xi) + p(\xi)$ as $R\to +\infty$.
In the same way, if
$v(\eta)=v_{\rm min}=-R\le u_{\rm min}$, then
$$
\omega(\eta)+q(\eta)\ge
 \frac {\lambda(\eta)}{1+\mu(\eta)e^{mv(\eta-\tau_2)}}\to \lambda(\eta)
$$
uniformly as $R\to +\infty$.
As $\alpha(t)-p(t) >\lambda(t) > \omega(t) +q(t)$ for all $t$, we 
deduce that $R$ cannot be arbitrarily large and the proof is complete.
\end{proof}

\section{Remarks about equilibrium points}
\label{equil}

In this section, we briefly discuss the uniqueness or multiplicity of 
positive equilibrium points for the autonomous case and make some comments 
on possible oscillation properties of the solutions.

With this aim, assume that all the parameters of \eqref{Wheldon0}
are constant, then the existence of at least one positive
equilibrium $(M_*,B_*)$ is easily shown, provided that
$$
n > (1-\delta)\frac m{m+1},\quad \alpha>\lambda-\delta p.
$$
Indeed, consider the system
\begin{equation} \label{Eq}
\begin{gathered}
\frac {\alpha }{1+\beta M^{n}} = \frac{\lambda }{1+\mu B^m}+\delta p,\\
(\omega + \delta q) B=\frac{\lambda M}{1+\mu B^m }
\end{gathered}
\end{equation}
and let 
$$
c(B) :=  \frac{B(1+ \mu B^m)(\omega +\delta q)}\lambda.
$$
Then \eqref{Eq} has at least a positive solution if and only if the
function $\varphi:[0,+\infty)\to \mathbb{R}$ given by
$$
\varphi(B):= \frac \alpha{1+\beta c(B)^n} -\frac{\lambda}{1+\mu B^m}-\delta p
$$
has at least a positive root. This is easily verified, since
$$
\varphi(0)=\alpha-\lambda-\delta p>0
$$
and
$$
\lim_{B\to +\infty}\varphi(B)=-\delta p.
$$
Thus, the result follows for $\delta=1$. When $\delta =0$, 
condition $n >\frac m{m+1}$ implies
$\varphi(B) <0$ for $B\gg 0$ and so completes the proof.
\smallskip

It is worth noticing that the number of equilibria depends on the parameters 
of the system. Although more
precise computations are possible, we shall not pursue a detailed analysis
 here and restrict ourselves to some elementary comments.
Consider, for instance, the case $\delta =0$, then
$$
B_* = \frac {\alpha M}{\omega (1+\beta M^n)}.
$$
Calling $z=1+\beta M^n$, we obtain the following equation for $z$:
$$
z = \frac \alpha \lambda + r \big[
\frac{\sqrt[n]{z-1}}{z}\big]^m:=\psi(z),
$$
where
$r=\frac{\alpha^{m+1}\mu}{\omega^{m} \beta^{m/n}\lambda}$.
The function $z-\psi(z)$ is negative for $z=1$ and, as $n > \frac m{m+1}$, 
tends to $+\infty$ as $z\to +\infty$.
Next, we compute
\begin{gather*}
\psi'(z)= \frac {rm{(z-1)^{\frac{m-n}n}}} {nz^{m+1}} [n - (n-1)z],\\
\psi''(z)=\frac {rm {(z-1)^{\frac{m-2n}n}}} {nz^{m+2}}[a z^2 + b z + c],
\end{gather*}
where
$$
a= \frac {n-1}n[n + m(n-1)],\quad
b = -2[n + m(n-1)],\quad
c = (m+1)n.
$$
In particular, $\psi$ vanishes at most twice in $(1,+\infty)$, which
implies that the system cannot have more than $3$ positive equilibrium points.

When $n\neq 1$,
the quadratic $a z^2 + bz + c$ has two different real roots, namely
$$
R_\pm = \frac n{n-1}\Big(1 \pm \frac 1{\sqrt{n+m(n-1)}}\Big).
$$

Let us prove, in the first place, that the positive equilibrium is unique
 when $m\le n$. This is immediate for $m<n$, since the function $z-\psi(z)$
is strictly decreasing near $1$, and $\psi''$ vanish at most once in 
$(1,+\infty)$. When $m=n$, there are two cases:
\begin{itemize}
 \item
If $n\le 1$, then $\psi''$ does not vanish in $(1,+\infty)$.
\item
If $n>1$, then direct
computation shows that the equation $\psi'(z)=1$ has at most one solution 
in $(1,+\infty)$.
\end{itemize}
In both cases, the function $z-\psi(z)$ has at most one critical point 
in $(1, +\infty)$ and the claim follows.

The situation is different when $m>n$: for instance, if $r$ is large enough 
then there are $3$ positive equilibria, provided that
$\frac \alpha\lambda$ is sufficiently close to $1$.
Indeed, we may set, for example, $R>1$ as the largest root of the quadratic 
function $a z^2 + b z +c$, namely
$$
R= \begin{cases}
 \frac {m+1}2 &\text{if } n=1,\\
R_- &\text{if } n<1,\\
R_+  &\text{if }  n>1,
\end{cases}
$$
with $R_\pm$ as before.
Next, consider the function $g(z)= z - \psi(z) + \frac\alpha\lambda -1$ and
fix $r$ such that $r > \frac{R^m}{(R-1)^{\frac{m-n}n}}$. Then
$g(R)<0$ and, as $g(1)=0$ and $g'(1)=1$, it is seen that $g$
has exactly one zero in $(1,R)$ and another one in $(R,+\infty)$. Now let
$$
\varepsilon = \max_{1\le z\le R} g(z),
$$
then the function $z -\psi(z)$ has $3$ zeros when 
$\frac\alpha\lambda < 1+\varepsilon$.

In view of the previous example, a natural question arises:
is it possible to find a sharp set of
sufficient conditions for the uniqueness of the positive equilibrium when $m>n$?
For example, a sufficient condition when $n\le 1$ is
$$
\frac\alpha\lambda \ge R
$$
with $R$ as before: indeed, in this case $\psi'(z) >0$ in $(1,+\infty)$, 
so $\psi(z)>z$ in $[1,R]$ and $\psi''$ does not vanish after $R$, so
the equation $\psi'(z) =1$ has at most one solution in 
$(R, +\infty)$.

When $n>1$, a sufficient condition for uniqueness of the positive equilibrium is:
$$
\frac \alpha\lambda \ge \frac{n}{n-1}.
$$
Indeed, in this case $\psi$ strictly increases
up to $z = \frac n{n-1}$ and strictly decreases after that point. 
As $\psi(z)>z$ on $(1,\frac n{n-1})$ it follows that the equation 
$\psi(z)=z$ has exactly one solution.
Observe that $R > \frac n{n-1}$, so the previous condition is sharper 
than the condition $\alpha/\lambda \ge R$.

Also, it is worth noticing that, in all cases, if $r$ is small then the 
equilibrium is unique. More precisely,
for $n\le 1$ the function $\psi'$ is positive and achieves its absolute 
maximum at $z=R$; thus, a sufficient condition for uniqueness is:
\begin{equation} \label{n le 1}
\psi'(R) <1.
\end{equation}
For $n>1$, the function $\psi'$ achieves its absolute maximum at $z=R_->1$. 
This yields the sufficient condition
\begin{equation}\label{n > 1}
\psi'(R_-)<1.
\end{equation}
Conditions \eqref{n le 1} and \eqref{n > 1} are obviously satisfied when 
$r$ is small.
\medskip

The presence of delays yields also an interesting matter about the
oscillation properties of the autonomous model.
This is an interesting field of research that can be the object of a 
future work; here, we shall
only prove some behavior that might indicate the presence of oscillation.

In more precise terms, we set a positive equilibrium $(M_{*},B_{*})$  
as the center of coordinates and denote
by $Q_j$ the $j$-th quadrant, namely
\begin{gather*}
Q_1:= \{ (M,B): M > M_{*}, B > B_{*}\},\\
Q_2:= \{ (M,B): M < M_{*}, B > B_{*}\},\\
Q_3:= \{ (M,B): M <M_{*}, B < B_{*}\},\\
Q_4:= \{ (M,B): M > M_{*}, B < B_{*}\}.
\end{gather*}
We shall prove that, under appropriate conditions,
if a non-constant positive solution starts in $\overline Q_2$ or 
$\overline Q_4$ then it cannot remain there for all $t$.

\begin{proposition} \label{prop3.1}
Let $\tau_{1} > (1+\beta M_{*}^{n})^2/(n\alpha\beta M_{*}^{n})$
and assume that there are no equilibrium points in $Q_4$.
Then there exists
a sequence $t_n\to +\infty$ such that, for all $n$, $M(t_n)<M_{*}$ or 
$B(t_n) > B_{*}$.
\end{proposition}

\begin{proposition} \label{prop3.2}
Let $\tau_{1} > \frac {(1+\beta M_{*}^n)^2}{n\alpha\beta M_{*}^{n}}$ and assume that
there are no equilibrium points in $Q_2$.
Then there exists a sequence $t_n\to +\infty$ such that, for all $n$, 
$M(t_n)>M_{*}$ or $B(t_n) < B_{*}$.
\end{proposition}
In other words, a non-constant positive solution starting at $\overline Q_2$ or
$\overline Q_4$
might abandon the respective quadrant and never return,
or it might eventually come back but then it leaves the quadrant again and so on.
A proof of Proposition \ref{prop3.1} is given below; 
the proof of Proposition \ref{prop3.2} is similar so we omit it.

\begin{lemma} \label{lem1}
Assume that $R(t_1-\tau_{1})\ge R(t_2-\tau_{1})$ and 
$B(t_1-\tau_{2})\le B(t_2-\tau_{2})$.
If $R(t_1)\ge R(t_2)$ and $B(t_1)\le B(t_2)$, at least one of the 
inequalities being strict,
then $R'(t_1) < R'(t_2)$ and $B'(t_1) > B'(t_2)$.
\end{lemma}

\begin{proof}
 It suffices to observe that the right hand side of the first equation 
of \eqref{generalc2} is strictly decreasing in the variables $R(t-\tau_{1})$ 
and strictly increasing in the variable $B(t-\tau_{2})$, and the
right hand side of the second equation of \eqref{generalc2}
is strictly increasing in the variable $R(t)$ and strictly decreasing 
in the variables $B(t)$ and $B(t-\tau_{2})$.
Then $R'(t_{1})<R'(t_{2})$ and $B'(t_{1})>B'(t_{2})$.
\end{proof}

 \begin{remark} \label{rmk3.1} \rm
As in  Lemma \ref{lem1}, it is easily seen that
if $R(t) > R_{*}:=\ln(M_{*})$
for $t \in [t_0-\tau_{1},t_1)$ and $B(t)<B_{*}$ for all $t\in [t_0-\tau_{2},t_1)$
then $R'(t)<0<B'(t)$ for all $t \in [t_0,t_1]$.
If
$R(t_1)=R_{*}$ or $B(t_1)=B_{*}$, then there exists $\eta >0$ such that
$(R(t),B(t))\notin \overline Q_4$ for $t\in (t_1,t_1+\eta)$.
On the other hand, if
$R(t) > R_{*}$ for all $t\ge t_0-\tau_{1}$ and $B(t)<B_{*}$ for all $t\ge t_0-\tau_{2}$ then $R'(t)<0< B'(t)$ for all $t\ge t_0$ and,
if there are no equilibrium points in $Q_4$, then $R(t)\to R_{*}$ and $B(t)\to B_{*}$.
\end{remark}



\begin{proof}[Proof of Proposition \ref{prop3.1}]
Suppose that $M(t)>M_{*}$ for all $t\ge t_0-\tau_{1}$ and $B(t)< B_{*}$ 
for all $t\ge t_0-\tau_{2}$. A simple computation shows that
$$
R'(t) = -A (R(t-\tau_{1})-R_{*}) -C(B_{*}-B(t-\tau_{2})) ,$$
with
\begin{gather*}
A = A(R(t),R(t-\tau_{1})):= \frac{\alpha\beta (e^{nR_{*}}-e^{nR(t-\tau_{1})})}
{(1+\beta e^{nR(t-\tau_{1})})(1+\beta e^{nR_*})(R_{*}-R(t-\tau_{1}))}>0,
\\
C= C(B(t),B(t-\tau_{2})):= \frac {\lambda\mu(B_{*}^{m}-B^{m}
 (t-\tau_{2}))}{(1+\mu B_{*}^m)(1+\mu B^{m}(t-\tau_{2}))(B_{*}-B(t-\tau_{2}))}>0,
\\
A(R(t),R(t-\tau_{1}))\to \frac{n\alpha\beta e^{nR_{*}}}{(1+\beta e^{nR_{*}})^2}\quad 
\text{ as } \, t\to +\infty,
\\
C(B(t),B(t-\tau_{2}))\to \frac{\lambda \mu m B_{*}^{m-1}}{(1+\mu B_{*}^m)^2}\quad 
\text{as }  t\to +\infty.
\end{gather*}
Moreover,
$$
R(t-\tau_{1})-R_{*} =
R(t-\tau_{1})- R(t) + R(t) -
R_{*} = -\tau_{1} R'(\theta) + R(t) -
R_{*}
$$
for some mean value $\theta \in (t-\tau_{1},t)$. From Lemma \ref{lem1} with 
$t_1=\theta$ and $t_2=t$, it follows that $R'(\theta)< R'(t)$.
Thus,
$$
R'(t) < -A (R(t)-R_{*})- C(B_{*} - B(t-\tau_{2})) + \tau_{1} A R'(t).
$$
Observe that the hypothesis says that
$\tau_{1} > \frac {(1+\beta e^{n R_{*}})^2}{n\alpha\beta e^{n R_{*}}}$.
Without loss of generality, we may assume that $t_0$ is large enough so that
$\tau_{1} A(R(t),R(t-\tau_{1})) > 1$, then
$$
(\tau_{1} A-1)R'(t) > A (R(t)-R_{*})+C(B_{*}-B(t-\tau_{2}))>0,
$$
a contradiction.
\end{proof}

\subsection*{Open Problems}

We outline some problems that might be of interest for scientists who plan
to start future research in this field.

(1) Use Lyapunov-like functionals to find sufficient conditions for the global 
stability of a non-trivial equilibrium of the autonomous model.

(2) Prove or disprove that for a new model the complete recovery is possible 
for sufficiently high drug dosage; examine permanence, persistence and 
extinction of the solutions.
 
(3) Define what is the required type, frequency and  intensity of the cancer 
treatment that switch unfavorable oscillatory dynamics of a system  
to a non-oscillatory state.


\subsection*{Acknowledgments}
 We thank the anonymous referee for insightful comments that led to an 
improvement of this manuscript.


\begin{thebibliography}{99}

\bibitem{Ba} J. Batzel, F. Kappel;
\emph{Time delay in physiological systems: Analyzing and modeling its impact}, 
Mathematical Biosciences,  234 (2011) 61-74.

\bibitem{Be} N. Bellomo, A. Bellouquid, J. Nieto, J. Soler;
\emph{Complexity and mathematical tools toward the modelling of multicellular 
growing systems}, Math. Comput. Modelling, 51 (2010) 441-551.

\bibitem{BKS} N. Bellomo, D. Knopoff,  J. Soler;
\emph{On the difficult interplay between life, ``complexity'', 
and mathematical sciences}, Math. Models Methods Appl. Sci., 23(10) (2013), 
1861-1913.

\bibitem{BDK} A. Bellouquid, E. De Angelis, D. Knopoff;
\emph{From the modeling of the immune hallmarks of cancer to a black swan 
in biology}, Math. Models Methods Appl. Sci., 23 (2013), 949-978.

\bibitem{Ca} P. Castorina, T. Deisboeck, P. Gabriele, C. Guiot;
\emph{Growth Laws in Cancer: Implications for Radiotherapy}. 
Radiat. Res. 168 (2007) 349–-356.

\bibitem{Mac2} C. Colijn, M. Mackey;
\emph{A mathematical model of hematopoiesis I. Periodic chronic myelogenous 
leukemia}, Journal of Theoretical Biology,  237 (2005) 117-132.

\bibitem{Ef} R. Eftimie, J. Bramson, D. Earn;
\emph{Interactions between the immune system and cancer: a brief review 
of non-spatial mathematical models}, Bull. Math. Biol. 73 (2011) 2-32.

\bibitem{Fe} E. Fessler, F. Dijkgraaf, F. De Sousa, E. Melo, J. Medema;
\emph{Cancer stem cell dynamics in tumor progression and metastasis: 
Is the microenvironment to blame?}
Cancer Letters, In Press, Corrected Proof, Available online 22 October 2012.

\bibitem{Gold} J. Goldman, J. Melo;
\emph{Chronic Myeloid Leukemia}, Advances in Biology and New Approaches 
to Treatment, N Engl J Med, 349 (2003) 1451-1464.

\bibitem{Ho} M. Horn, M. Loeffler, I. Roeder;
\emph{Mathematical modeling of genesis and treatment of chronic myeloid 
leukemia}, Cells Tissues Organs, 188 (2008) 236-247.

\bibitem{Lo} N. Lloyd;
\emph{Degree Theory}, Cambridge University. Press, Cambridge 1978.

\bibitem{Ma} P. Macklina, J. Lowengrub;
\emph{Nonlinear simulation of the effect of microenvironment on tumor growth}, 
Journal of Theoretical Biology, 245 (2007) 677-704.

\bibitem{Maa} H. Mayani, E. Flores-Figueroa, A. Chavez-Gonzalez;
\emph{In vitro biology of human myeloid leukemia}, 
Leukemia Research, 33 (2009)  624-637.

\bibitem{Maw} J. Mawhin;
\emph{Topological degree methods in nonlinear boundary value
  problems, volume 40 of  CBMS Regional Conference Series in Mathematics}.
American Mathematical Society, Providence, RI, 1979.
Expository lectures from the CBMS Regional Conference held at Harvey
  Mudd College, Claremont, Calif., June 9--15, 1977.

\bibitem{Pe} D. Perrotti, C. Jamieson, J. Goldman, T. Skorski;
\emph{Chronic myeloid leukemia: mechanisms of blastic transformation}, 
J. Clin Invest. 120 (2010) 2254-2264.

\bibitem{Ro} I. Roeder, M. d'Inverno, et al.;
\emph{New experimental and theoretical investigations of hematopoietic 
stem cells and chronic myeloid leukemia}, Blood Cells, Molecules, and Diseases,  
43  (2009) 88-97.

\bibitem{Sw} A. Swierniak, M. Kimmel, J. Smieja;
\emph{Mathematical modeling as a tool for planning anticancer therapy}, 
European Journal of Pharmacology, 625  (2009)  108-121.

\bibitem{Vap} P. Vaupel;
\emph{Tumor microenvironmental physiology and its implications for radiation 
oncology}, Seminars in Radiation Oncology,  14 (2004) 198-206.


\bibitem{We2} E. Wheldon, K. Lindsay, T.  Wheldon;
\emph{The dose-response relationship for cancer incidence in a two-stage 
radiation carcinogenesis model incorporating cellular repopulation}, 
International Journal of Radiation Biology, 76 (2000) 699-710.

\bibitem{W1} T. Wheldon;
\emph{Mathematical Models in Cancer Research}, 
Bristol and Philadelphia, PA: Adam Hilger 1988.

\bibitem{Wek0} T. Wheldon;
\emph{Mathematical models of oscillatory blood cell production}, 
Mathematical Biosciences, 24 (1975) 289-305.

\bibitem{We0}T. Wheldon, J. Kirk, H. Finlay;
\emph{Cyclical granulopoiesis in chronic granulocytic leukemia: 
a simulation study},  Blood,  43 (1974) 379-387.

\end{thebibliography}

\end{document}
