\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}