\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{amssymb}

\AtBeginDocument{{\noindent\small
Seventh Mississippi State Conference on Differential Equations and
Computational Simulations,
{\em Electronic Journal of Differential Equations},
Conf. 17 (2009),  pp. 197--206.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2009 Texas State University - San Marcos.}
\vspace{9mm}}

\begin{document} \setcounter{page}{197}
\title[\hfilneg EJDE-2009/Conf/17\hfil A model for gene activation]
{A model for gene activation}

\author[S. F. Oppenheimer, R. Fan, S. Pruett\hfil EJDE/Conf/17 \hfilneg]
{Seth F. Oppenheimer, Ruping Fan, Stephan Pruett}  % in alphabetical order

\address{Seth F. Oppenheimer \newline
Department of Mathematics and Statistics\\
Mississippi State University\\
Drawer MA, Mississippi State, MS 39762, USA}
\email{seth@math.msstate.edu}

\address{Ruping Fan \newline
Department of Cellular Biology and Anatomy\\
Louisiana State University Health Science Center \\
1501 Kings Highway\\
Shreveport, LA 71130, USA}

\address{Stephan Pruett \newline
 CVM Basic Science Department\\
Mississippi State University\\
PO box 6100\\
Mississippi State, MS 39762, USA}
\email{sbp7@msstate.edu}

\thanks{Published April 15, 2009.}
\thanks{SFO was supported by National Institutes of Health grants
DHHS/NIH 5 P20 RR17661 \hfill\break\indent
and NIH ES 11278. Center for Environmental Health Sciences Publication \#119}
\subjclass[2000]{92C17, 34A12}
\keywords{Gene activation; neurotoxons}

\begin{abstract}
 The purpose of this paper is to develop a model for the activation of the
 gene responsible for the production of the cytokine interleukin 6, IL-6.
 This is motivated by experimental work that indicates that exposure to
 certain exogenous chemicals results in changes in cytokine production.  In
 particular, exposure to both the widely used pesticide atrazine and the
 legacy pesticide dieldrin, still very much present in the environment,
 resulted in the reduction of the production of IL-6.  We develop of model
 of twelve ordinary differential equations to model the effect of changes in
 transcription factor levels on IL-6 production rates and establish basic
 qualitative properties of solutions.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]

\section{Introduction}

The purpose of this paper is to develop a model for the activation of the
gene responsible for the production of the cytokine interleukin 6, IL-6.
This is motivated by experimental work that indicates that exposure to
certain exogenous chemicals results in changes in cytokine production.  In
particular, exposure to both the widely used pesticide atrazine and the
legacy pesticide dieldrin, still very much present in the environment,
resulted in the reduction of the production of IL-6 \cite{p1}.  We start with a
description of the biological problem.

We recall that a cytokine is a signaling protein or glycoprotein used in
cellular communication. In particular, IL-6 is important in dealing with
immune response to trauma, the bone formation and breakdown cycle, and,
perhaps, response to some bacterial attacks \cite{a1,p1}.  The following
surprising result motivated the search for a mathematical model.  Doses of
dieldrin and Atrazine, which separately produced roughly  10\% drops in
IL-6 levels, when combined, produced an 80\% drop in IL-6 levels, far
greater than the expected additive or subadditive effect \cite{p1}:
\begin{equation*}
\begin{tabular}{ll}
Treatment & IL-6 level post treatment \\
Dieldrin & 0.92 \\
Atrazine & 0.88 \\
Both & 0.19
\end{tabular}
\end{equation*}
Given that public policy on acceptable exposure levels are based on single
exposure measurements, this is bad news.

Further measurements showed that the exposure to the dieldrin and atrazine
resulted in the reduction of certain transcription factors within the
nucleus of the cell. A transcription factor is a protein that binds to a
specific part of the DNA.  called a response element. The transcription
factors we consider are named $AP-1$, $NF-\kappa B$, and $AP-1-NF-\kappa B$.

The experimental results now look like
\begin{equation*}
\begin{tabular}{llll}
Treatment & $AP-1$ & $NF-\kappa B$ & $AP-1-NF-\kappa B$ \\
Dieldrin & 0.7 & 0.7 & 0.49 \\
Atrazine & 0.65 & 0.7 & 0.455 \\
Both & 0.35 & 0.4 & 0.14
\end{tabular}
\end{equation*}
\begin{equation*}
\begin{tabular}{ll}
Treatment & IL-6 level post treatment \\
Dieldrin & 0.92 \\
Atrazine & 0.88 \\
Both & 0.19
\end{tabular}
\end{equation*}

The desire now was to model IL-6 production in terms of transcription
factors. Thus, we will develop a model that will give the rate of IL-6
production as a function of transcription factor concentration within the
nuclei of the cells.  To do this we will need to also know the what
response elements are important. In the each cell of a mouse there is a gene
with two response elements, $ARE$ and $NRE$.  If both response elements are
activated by having the appropriate proteins bind to each of them, the gene
will start the cell producing IL-6.  We have $AP-1$ binding reversibly to
the $ARE$ response element. We expect competition for the $NRE$ response
element between the $NF-\kappa B$ and the $AP-1-NF-\kappa B$ \cite{f1}.

In section 2 of this paper, we will produce a model that gives the fraction
of response elements filled with a particular transcription factor as a
function of time.  Section 3 will develop a model showing what fraction of
genes are activated and producing IL-6 and how they are activated. Finally,
section 4 will give some basic mathematical results for the model developed
in section 2.

\section{The main physical model}

Our initial assumptions are as follows.
\begin{itemize}
\item[1.]  All genes will be treated as being in a single well mixed
solution of cytoplasm.

\item[2.] The likelihood of binding to a response element is independent
of the status of the other response element.

\item[3.] We will only consider three chemical species $AP-1$,
$NF-\kappa B$, and $ AP-1-NF-\kappa B$.

\end{itemize}
In the solution, we assume the following reaction is occurring:
\begin{equation*}
AP-1+NF-\kappa B\rightleftarrows _{k_2}^{k_1}AP-1-NF-\kappa B
\end{equation*}
and that it is governed by mass action. Using the following variable names
\begin{equation*}
\begin{array}{cc}
\text{Chemical} & \text{Variable} \\
AP-1\text{ molarity} & c_1 \\
NF-\kappa B\text{ molarity} & c_2 \\
AP-1-NF-\kappa B\text{ molarity} & c_3
\end{array}
\end{equation*}
we obtain the three nonlinear differential equations below\
\begin{gather*}
\frac{dc_1}{dt} = -k_1c_1c_2+k_2c_3, \\
\frac{dc_2}{dt} = -k_1c_1c_2+k_2c_3, \\
\frac{dc_3}{dt} = k_1c_1c_2-k_2c_3.
\end{gather*}
These equations describe the chemical reactions in solution assuming no
other mechanisms are in play.

We now consider the other physical processes in isolation as well.  These
are $AP-1$ binding reversibly to the $ARE$ response element and competition
for binding to the $NRE$ response element between the $NF-\kappa B$ and the
$AP-1-NF-\kappa B$.

We assume that we have a total volume $V$ of cytoplasm, our fluid, and that
there are $M$ moles of cells. The units for the solution concentrations,
$c_{j}$, will be molarities, and the sorbed concentrations, $q_{j}$, will be
in moles per mole.

We model this as we would a sorption process \cite{o1}.
\begin{equation*}
\begin{array}{cc}
\text{Chemical} & \text{Variable} \\
AP-1\text{ molarity} & c_1 \\
NF-\kappa B\text{ molarity} & c_2 \\
AP-1-NF-\kappa B\text{ molarity} & c_3 \\
\begin{array}{c}
\text{Moles of }AP-1 \\
\text{per mole of }ARE
\end{array}
& q_1 \\
\begin{array}{c}
\text{Moles of }NF-\kappa B \\
\text{per mole of }NRE
\end{array}
& q_2 \\
\begin{array}{c}
\text{Moles of }AP-1-NF-\kappa B \\
\text{per mole of }NRE
\end{array}
& q_3
\end{array}
\end{equation*}
We make the assumption that at equilibrium for a fixed set of solution
concentrations we have a fixed fraction of response elements filled:
\begin{gather*}
q_1 = f_1\left( c_1\right) , \\
q_2 = f_2\left( c_1,c_3\right) , \\
q_3 = f_3\left( c_1,c_3\right) .
\end{gather*}
Competition for the $NRE$ response element between the $NF-\kappa B$ and the
$AP-1-NF-\kappa B$ will mean that $f_2$ will be increasing in $c_2$ and
decreasing in $c_3$ and $f_3$ will be increasing in $c_3$ and
decreasing in $c_2$. We will assume $f_1$ is an increasing function and
that $f_{j}=0$ when $c_{j}\leq 0$ and positive otherwise.

Our initial hope was to use the Langmuir model for a single layer sorption
process
\begin{gather*}
q_1 = f_1\left( c_1\right) =\frac{c_1}{\beta _1+c_1}, \\
q_2 = f_2\left( c_2,c_3\right) =\frac{c_2}{\beta
_2+c_2+\gamma _3c_3}, \\
q_3 = f_3\left( c_2,c_3\right) =\frac{c_3}{\beta
_3+c_3+\gamma _2c_2}.
\end{gather*}
Unfortunately, experimental measurements of equilibrium data \cite{p1}
did not conform to our expectations. We will thus use a
Langmuir-Freundlich type model that will allow for competition.
\begin{gather*}
q_1  = f_1\left( c_1\right) =\frac{c_1^{\pi _1}}{\beta
_1+c_1^{\pi _1}}, \\
q_2 = f_2\left( c_2,c_3\right) =\frac{c_2^{\pi _2}}{\beta
_2+c_2^{\pi _2}+\gamma _3c_3^{\pi _{23}}}, \\
q_3 = f_3\left( c_2,c_3\right) =\frac{c_3^{\pi _3}}{\beta
_3+c_3^{\pi _3}+\gamma _2c_2^{\pi _{32}}}.
\end{gather*}
However, our analysis will not depend on the functional form. That
functional form was used to fit data in \cite{p1}.

We assume that rate of change between sorbed concentration and solution
concentration is proportional to the difference from equilibrium:
\begin{gather*}
\frac{dc_1}{dt} = r_{c1}\left( q_1-f_1\left( c_1\right) \right) , \\
\frac{dc_2}{dt} = r_{c2}\left( q_2-f_2\left( c_2,c_3\right)
\right) , \\
\frac{dc_3}{dt} = r_{c3}\left( q_3-f_3\left( c_2,c_3\right)
\right) , \\
\frac{dq_1}{dt} = r_{q1}\left( f_1\left( c_1\right) -q_1\right) , \\
\frac{dq_2}{dt} = r_{q2}\left( f_2\left( c_2,c_3\right)
-q_2\right) , \\
\frac{dq_3}{dt} = r_{q3}\left( f_3\left( c_2,c_3\right)
-q_3\right) .
\end{gather*}
We observe that the total number of moles of a given transcription factor is
\begin{equation*}
q_{j}M+c_{j}V.
\end{equation*}
To maintain conservation of mass, we require
\begin{equation*}
r_{qj}=\frac{V}{M}r_{cj}.
\end{equation*}
We will rename
\begin{equation*}
r_{j}=r_{cj}
\end{equation*}
Our full model (so far) is now
\begin{equation} \label{CEQ}
\begin{gathered}
\frac{dc_1}{dt} = r_1\left( q_1-f_1\left( c_1\right) \right)
-k_1c_1c_2+k_2c_3,   \\
\frac{dc_2}{dt} = r_2\left( q_2-f_2\left( c_2,c_3\right)
\right) -k_1c_1c_2+k_2c_3,   \\
\frac{dc_3}{dt} = r_3\left( q_3-f_3\left( c_2,c_3\right)
\right) +k_1c_1c_2-k_2c_3,   \\
\frac{dq_1}{dt} = \frac{V}{M}r_1\left( f_1\left( c_1\right)
-q_1\right) ,   \\
\frac{dq_2}{dt} = \frac{V}{M}r_2\left( f_2\left( c_2,c_3\right)
-q_2\right) ,   \\
\frac{dq_3}{dt} = \frac{V}{M}r_3\left( f_3\left( c_2,c_3\right)
-q_3\right).
\end{gathered}
\end{equation}

\section{The book keeping equations}

We will build a probabilistic model on top of the differential equations
model which will follow the fraction of cells in each state of activation.
We have 6 classes of cells:
\begin{equation*}
\begin{array}{cc}
\text{Cell type} &
\begin{array}{c}
\text{Moles of cells in class} \\
\text{ per mole of cells}
\end{array}
\\
\begin{array}{c}
\text{No bound } \\
\text{response sites}
\end{array}
& p_{0} \\
\begin{array}{c}
ARE\text{ occupied } \\
\text{by }AP-1 \\
\text{but }NRE\text{ unoccupied}
\end{array}
& p_1 \\
\begin{array}{c}
NRE\text{ occupied} \\
\text{ by }NF-\kappa B \\
\text{but }ARE\text{ unoccupied}
\end{array}
& p_2 \\
\begin{array}{c}
NRE\text{ occupied } \\
\text{by }NF-\kappa B-AP-1 \\
\text{but }ARE\text{ unoccupied}
\end{array}
& p_3 \\
\begin{array}{c}
ARE\text{ occupied } \\
\text{by }AP-1 \\
\text{and }NRE\text{ occupied } \\
\text{by }NF-\kappa B
\end{array}
& p_{4} \\
\begin{array}{c}
ARE\text{ occupied } \\
\text{by }AP-1 \\
\text{and }NRE\text{ occupied } \\
\text{by }NF-\kappa B-AP-1
\end{array}
& p_{5}
\end{array}
\end{equation*}
We now must use this information to find the values of the $p_{j}$ as
functions of time. We will assume that only one change occurs at a time.
That is a single site may become occupied or unoccupied at a time. The
diagram below indicates how the states connect,
where the edges indicate reversible transformations.

\begin{figure}[ht] \label{fig1}
\begin{center}
\setlength{\unitlength}{0.8mm}
\begin{picture}(80,60)(0,10)
\put(38.5,69.5){$p_4$}
\put(38,69){\line(-5,-2){33}}
\put(43,69){\line(4,-3){36}}
\put(0,54.5){$p_2$}
\put(38,41){\line(-5,2){33}}
\put(0,24.5){$p_3$}
\put(38,39){\line(-5,-2){33}}
\put(38.5,39.5){$p_0$}
\put(43,40){\line(1,0){36}}
\put(38.5,9.5){$p_5$}
\put(38,11){\line(-5,2){33}}
\put(43,11){\line(4,3){36}}
\put(79.5,39.5){$p_1$}
\end{picture}
\end{center}
\caption{Graph of the activation states of genes}
\end{figure}


Our model is probabilistic in the sense that the rate of change from one
state to another depends on the proportion of cells in the given state at a
given time. We will use the notation
\[
t^{+} = \max (0,t) , \quad
t^{-} = (-t) ^{+}.
\]
The charts below will indicate how each state is changing, what state it is
changing to, and at what rate.
\begin{equation*}
\begin{tabular}{ll}
& Transition rate for $p_{0}$ \\
Into $p_1$ & $\big( \frac{dq_1}{dt}\big) ^{+}\frac{p_{0}}{
p_{0}+p_2+p_3}$ \\
Into $p_2$ & $\big( \frac{dq_2}{dt}\big) ^{+}\frac{p_{0}}{p_{0}+p_1}
$ \\
Into $p_3$ & $\big( \frac{dq_3}{dt}\big) ^{+}\frac{p_{0}}{p_{0}+p_1}
$
\end{tabular}
\end{equation*}
\begin{equation*}
\begin{tabular}{ll}
& Transition rate for $p_1$ \\
Into $p_{0}$ & $\big( \frac{dq_1}{dt}\big) ^{-}\frac{p_1}{
p_1+p_{4}+p_{5}}$ \\
Into $p_{4}$ & $\big( \frac{dq_2}{dt}\big) ^{+}\frac{p_1}{p_{0}+p_1}
$ \\
Into $p_{5}$ & $\big( \frac{dq_3}{dt}\big) ^{+}\frac{p_1}{p_{0}+p_1}
$
\end{tabular}
\end{equation*}
\begin{equation*}
\begin{tabular}{ll}
& Transition rate for $p_2$ \\
Into $p_{0}$ & $\big( \frac{dq_2}{dt}\big) ^{-}\frac{p_2}{p_2+p_{4}}
$ \\
Into $p_{4}$ & $\big( \frac{dq_1}{dt}\big) ^{+}\frac{p_2}{
p_{0}+p_2+p_3}$
\end{tabular}
\end{equation*}
\begin{equation*}
\begin{tabular}{ll}
& Transition rate for $p_3$ \\
Into $p_{0}$ & $\big( \frac{dq_3}{dt}\big) ^{-}\frac{p_3}{p_3+p_{5}}
$ \\
Into $p_{5}$ & $\big( \frac{dq_1}{dt}\big) ^{+}\frac{p_3}{
p_{0}+p_2+p_3}$
\end{tabular}
\end{equation*}
\begin{equation*}
\begin{tabular}{ll}
& Transition rate for $p_{4}$ \\
Into $p_1$ & $\big( \frac{dq_2}{dt}\big) ^{-}\frac{p_{4}}{p_2+p_{4}}
$ \\
Into $p_2$ & $\big( \frac{dq_1}{dt}\big) ^{-}\frac{p_{4}}{
p_1+p_{4}p_{5}}$
\end{tabular}
\end{equation*}
\begin{equation*}
\begin{tabular}{ll}
& Transition rate for $p_{5}$ \\
Into $p_1$ & $\big( \frac{dq_3}{dt}\big) ^{-}\frac{p_{5}}{p_3+p_{5}}
$ \\
Into $p_3$ & $\big( \frac{dq_1}{dt}\big) ^{-}\frac{p_{5}}{
p_1+p_{4}+p_{5}}$
\end{tabular}
\end{equation*}

We can now write out the six additional equations for the fraction of cells
in each state. For brevity, instead of writing
\begin{equation*}
\frac{V}{M}r_1\left( f_1\left( c_1\right) -q_1\right) ,
\end{equation*}
we will simply write
$\frac{dq_1}{dt}$.
The equations are:
\begin{align*}
\frac{dp_{0}}{dt} &= -\big( \frac{dq_1}{dt}\big) ^{+}\frac{p_{0}}{
p_{0}+p_2+p_3}-\big( \frac{dq_2}{dt}\big) ^{+}\frac{p_{0}}{
p_{0}+p_1}-\big( \frac{dq_3}{dt}\big) ^{+}\frac{p_{0}}{p_{0}+p_1} \\
&\quad+\big( \frac{dq_1}{dt}\big) ^{-}\frac{p_1}{p_1+p_{4}+p_{5}}
+\big( \frac{dq_2}{dt}\big) ^{-}\frac{p_2}{p_2+p_{4}}+\big( \frac{dq_3}{dt}\big) ^{-}\frac{p_3}{p_3+p_{5}},
\end{align*}
\begin{align*}
\frac{dp_1}{dt} &= \big( \frac{dq_1}{dt}\big) ^{+}\frac{p_{0}}{
p_{0}+p_2+p_3}-\big( \frac{dq_1}{dt}\big) ^{-}\frac{p_1}{
p_1+p_{4}+p_{5}}-\big( \frac{dq_2}{dt}\big) ^{+}\frac{p_1}{
p_{0}+p_1} \\
&\quad -\big( \frac{dq_3}{dt}\big) ^{+}\frac{p_{0}}{p_{0}+p_1}+\big( \frac{dq_2}{dt}\big) ^{-}\frac{p_{4}}{p_2+p_{4}}+\big( \frac{dq_3}{dt}\big) ^{-}\frac{p_{5}}{p_3+p_{5}},
\end{align*}
\begin{align*}
\frac{dp_2}{dt} &= \big( \frac{dq_2}{dt}\big) ^{+}\frac{p_{0}}{
p_{0}+p_1}-\big( \frac{dq_2}{dt}\big) ^{-}\frac{p_2}{p_2+p_{4}}
-\big( \frac{dq_1}{dt}\big) ^{+}\frac{p_2}{p_{0}+p_2+p_3} \\
&\quad +\big( \frac{dq_1}{dt}\big) ^{-}\frac{p_{5}}{p_1+p_{4}+p_{5}}
\end{align*}
\begin{align}
\frac{dp_3}{dt}
&= \big( \frac{dq_3}{dt}\big) ^{+}\frac{p_{0}}{
p_{0}+p_1}-\big( \frac{dq_3}{dt}\big) ^{-}\frac{p_3}{p_3+p_{5}}
-\big( \frac{dq_1}{dt}\big) ^{+}\frac{p_3}{p_{0}+p_2+p_3} \notag \\
&\quad +\big( \frac{dq_1}{dt}\big) ^{-}\frac{p_{5}}{p_1+p_{4}+p_{5}},
\label{BEQ}
\end{align}
\begin{align*}
\frac{dp_{4}}{dt} &= \big( \frac{dq_2}{dt}\big) ^{+}\frac{p_1}{
p_{0}+p_1}+\big( \frac{dq_1}{dt}\big) ^{+}\frac{p_2}{
p_{0}+p_2+p_3}-\big( \frac{dq_2}{dt}\big) ^{-}\frac{p_{4}}{
p_2+p_{4}} \\
&\quad -\big( \frac{dq_1}{dt}\big) ^{-}\frac{p_{4}}{p_1+p_{4}p_{5}}
\end{align*}
\begin{align*}
\frac{dp_{5}}{dt}
&= \big( \frac{dq_3}{dt}\big) ^{+}\frac{p_1}{
p_{0}+p_1}+\big( \frac{dq_1}{dt}\big) ^{+}\frac{p_3}{
p_{0}+p_2+p_3}-\big( \frac{dq_3}{dt}\big) ^{-}\frac{p_{5}}{
p_3+p_{5}} \\
&\quad -\big( \frac{dq_1}{dt}\big) ^{-}\frac{p_{5}}{p_1+p_{4}+p_{5}}.
\end{align*}

The IL-6 production only occurs in cells in class $p_{4}$ and $p_{5}$ with
the production rate for class $p_{5}$ about 5 times as high as for class
$p_{4}$ \cite{f1}. So IL-6 production is proportional to
\begin{equation*}
p_{4}+5p_{5}.
\end{equation*}

The book keeping equations are purely for the short term. After the
equilibrium values for $q_1,\;q_2,$ and $q_3$ are reached, we expect
the cell classes to equilibrate according to the expected probability
distribution. That is, after a longer period of time, with the rates of
desorption and resorption balancing out, we will expect
\begin{gather*}
p_{4} = q_1q_2, \\
p_{5} = q_1q_3.
\end{gather*}

\section{Mathematical results for the model}

We start by considering existence and uniqueness for \eqref{CEQ}.
If we assume the
$f_{j}$ are Lipschitz, the Picard theorem \cite{b1} yields the following.

\begin{theorem} \label{thm1}
The system \eqref{CEQ} has a unique local solution for each set of initial
conditions. Furthermore, if the solutions stay bounded on each interval of
existence, the system \eqref{CEQ} has a unique global solution.
\end{theorem}

There will be a difficulty for the Langmuir-Freundlich functions proposed
above at $c_{j}=0$ if any of the $\pi _{j}<1$. However, we will show that
when the initial conditions are strictly positive, then the solutions remain
bounded and strictly positive, so the above theorem holds for every choice
of positive initial conditions.

If we take the following masses per mole
\begin{equation*}
\begin{tabular}{ll}
Transcription factor & mass per mole \\
$AP-1$ & $m_1$ \\
$NF-\kappa B$ & $m_2$ \\
$AP-1-NF-\kappa B$ & $m_1+m_2$
\end{tabular}
\end{equation*}
we obtain the following result.

\begin{theorem} \label{thm2}
For any solution to \eqref{CEQ} the quantity
\begin{equation*}
m_1Vc_1+m_2Vc_2+( m_1+m_2)
Vc_3+m_1Mq_1+m_2Mq_2+( m_1+m_2) Mq_3
\end{equation*}
is constant in time
\end{theorem}

\begin{proof}
Multiply the first equation in \eqref{CEQ} by $m_1V$, the second equation by
$m_2V,$ the third equation by $(m_1+m_2) V$, the fourth
equation by $m_1M$, the fifth equation by $m_2M$, and the sixth by
$(m_1+m_2) M$ and then add all of the equations together we
obtain
\begin{equation*}
\frac{d}{dt}\left( m_1Vc_1+m_2Vc_2+(m_1+m_2)
Vc_3+m_1Mq_1 + m_2Mq_2+(m_1+m_2) Mq_3\right) =0,
\end{equation*}
and we are done. We will call that constant value $H_1$.
\end{proof}

We observe that if all of the components of the solution are nonnegative,
than the above result implies each component will be bounded by some
positive constant, $G$. With this observation, we now obtain positivity of
solution for positive initial conditions.

\begin{theorem} \label{thm3}
Assume that there is are positive constant $K$ and $p\in (0,1]$ so that for
positive $c_{j}$ in a neighborhood of $0$, $f_{j}\leq Kc_{j}^{p}$. For any
solution to \eqref{CEQ}, if all of the initial conditions are positive,
then the solutions will remain positive. Given the previous result, any such
solution will be confined to a bounded region in the positive cone of
$\mathbb{R}^{6}$.
\end{theorem}

\begin{proof}
We start with the observation that for each $j$,
\begin{equation*}
\frac{dq_{j}}{dt}\geq -\frac{V}{M}r_{j}q_{j},
\end{equation*}
Thus
\begin{equation*}
q_{j}(t) \geq \exp \big( -\frac{V}{M}r_{j}t\big) q_{j}(0) >0.
\end{equation*}
Because any solution will be continuous on its interval of existence, there
will be a interval on which all components are nonnegative. Let $[0,T]$ be
the maximum such interval. Let $t\in \lbrack 0,T]$; then
\begin{align*}
\frac{dc_1}{dt} &= r_1\left( q_1-f_1\left( c_1\right) \right)
-k_1c_1c_2+k_2c_3 \\
&\geq  r_1\exp \big( -\frac{V}{M}rt\big) q_1(0)
-r_1Kc_1^{p}-k_1Gc_1.
\end{align*}
Thus, if $r_1\exp \left( -\frac{V}{M}rt\right) q_1(0)>r_1Kc_1^{p}+k_1Gc_1$,
\begin{equation*}
\frac{dc_1}{dt}>0
\end{equation*}
and $c_1$ has a positive lower bound on the entire interval. A similar
argument provides positive lower bounds for $c_2$ and $c_3$ on the whole
interval $[0,T] $. Therefore, the interval is actually
infinite and the solutions components are strictly positive for all positive
$t$.
\end{proof}

We will now obtain two more conservation relationships that will allow us to
state what the unique equilibrium will be for a given set of positive
initial values.

\begin{theorem} \label{thm4}
For any solution to \eqref{CEQ} the quantities
\begin{equation*}
Vc_1+Mq_1+Vc_3+Mq_3
\end{equation*}
and
\begin{equation*}
Vc_2+Mq_2+Vc_3+Mq_3
\end{equation*}
are constant. We will call these constants $H_2$ and $H_3$ respectively
\end{theorem}

\begin{proof}
To obtain the first, multiply the first equation in \eqref{CEQ} by $V$
the third equation by $V$, the fourth equation by $M$,and the sixth
 by $M$ and then add all of the equations together we obtain
\begin{equation*}
\frac{d}{dt}\left( Vc_1+Vc_3+Mq_1+Mq_3\right) =0.
\end{equation*}
The second sum is analyzed the same way.
\end{proof}

We now see that we have an equilibrium if the following equations are
satisfied
\begin{gather*}
\begin{aligned}
&m_1Vc_1+m_2Vc_2+(m_1+m_2)Vc_3\\
&+m_1Mq_1+m_2Mq_2+(m_1+m_2) Mq_3 = H_1,
\end{aligned} \\
Vc_1+Mq_1+Vc_3+Mq_3 = H_2, \\
Vc_2+Mq_2+Vc_3+Mq_3 = H_3, \\
q_1 = f\left( c_1\right) , \\
q_2 = f_2\left( c_2,c_3\right) , \\
q_3 = f_3\left( c_2,c_3\right) , \\
c_3 = \frac{k_1}{k_2}c_1c_2,
\end{gather*}
The monotonicity properties of the $f_{j}$ should guarantee a unique
positive equilibrium point.

Turning to the book keeping equations, \eqref{BEQ} the observation that the
forcing functions are Lipschitz in the $p_{j}$ and continuous in $t$ will
allow the application of the Picard theorem again to guarantee unique local
solutions. Simply adding the equations up will yield
\begin{equation*}
\frac{d}{dt}\left( p_{0}+p_1+p_2+p_3+p_{4}+p_{5}\right) =0;
\end{equation*}
i.e., the total number of cells is conserved by the model. Similar
arguments as the ones used for \eqref{CEQ} will give us that
the cell classes stay nonnegative.

\begin{thebibliography}{0}

\bibitem{a1} Abul K. Abbas and Andrew H. Lichtman;
\emph{Basic Immunology}, W. B. Saunders Company, Philadelphia, 2001.

\bibitem{b1} M. Braun;
 \emph{Differential Equations and Their Applications, Second
Edition}, Springer-Verlag, New York, 1975.

\bibitem{f1} L. Faggioli, C. Costanzo, M. Donadelli, and M. Palmieri;
\emph{Activation of the Interleukin-6 promoter by a dominant negative
 mutant of c-Jun}, Biochim Biophys Acta Vol. 1692(1) pages 17-24, 2004

\bibitem{o1} Seth F. Oppenheimer, William L. Kingery and FengXiang Han;
\emph{Phase Plane Analysis and Dynamical Systems Approaches to the
Study of Metal Sorption in Soils.} ``Sorption/desorption of heavy metals
in soils'', H. M. Salim and D. L. Sparks Editors. CRC Press (2001),
pp 109--130.

\bibitem{p1} Stephen Pruett, Ruping Fan, and Seth Oppenheimer;
\emph{Greater Than Additive Suppression of TLR3-Induced IL-6 Responses
by Administration of Dieldrin and Atrazine, }.
J. Immunotoxicology Volume 3, Issue 4 December
2006, pages 253 - 262..

\bibitem{s1} Lauren Sompayrac;
 \emph{How the Immune System Works, Second Edition},
Blackwell Publishing, Malden, MA, 2003.

\end{thebibliography}

\end{document}