\documentclass[reqno]{amsart} \usepackage{hyperref} \AtBeginDocument{{\noindent\small {\em Electronic Journal of Differential Equations}, Vol. 2005(2005), No. 62, pp. 1--9.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu \newline ftp ejde.math.txstate.edu (login: ftp)} \thanks{\copyright 2005 Texas State University - San Marcos.} \vspace{9mm}} \begin{document} \title[\hfilneg EJDE-2005/62\hfil A steady state of morphogen gradients] {A steady state of morphogen gradients for semilinear elliptic systems} \author[E. H. Kim\hfil EJDE-2005/62\hfilneg] {Eun Heui Kim} \address{Eun Heui Kim \hfill\break Department of Mathematics, California State University, Long Beach, CA 90840-1001, USA} \email{ekim4@csulb.edu \; tel 562-985-5338 \; fax 562-985-8227} \date{} \thanks{Submitted September 9, 2004. Published June 15, 2005.} \thanks{Supported by grant DMS-0103823 from the National Science Foundation, and by \hfill\break\indent grant DE-FG02-03ER25571 from the Department of Energy} \subjclass[2000]{35J55, 35J45} \keywords{Elliptic systems; nonquasimonotone; morphogen gradients} \begin{abstract} In this paper we establish the existence of positive solutions to a system of steady-state Neumann boundary problems. This system has been observed in some biological experiments, morphogen gradients; effects of {\it Decapentaplegic} (Dpp) and {\it short gastrulation} (Sog). Mathematical difficulties arise from this system being nonquasimonotone and semilinear. We overcome such difficulties by using the fixed point iteration via upper-lower solution methods. We also discuss an example, the Dpp-Sog system, of such problems. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \section{Introduction} We present a mathematical analysis on the model problem of morphogen gradients. Morphogens are molecules that diffuse from a local source to form a concentration. The concentration determines the reactions of all cells and activates genes to form patterns of cell differentiation \cite{L}. Experiments in morphogen diffusions to understand the pattern formations of tissue in developing animals give rise to a broad range of systems of reaction-diffusion equations see for example \cite{GM,GB,KNF,LNW,r}. This paper focuses on a mathematical analysis on a general system including the model problem of morphogen gradients. More precisely, we establish the existence of positive solutions to the nonlinear steady state system arising morphogen gradients. Interesting features of the system we study in this paper are that it is nonquasimonotone, where few theories can be found for nonquasimonotone system \cite{CM,gm,gm1,KS} (note that although the first three references are systems with the zero Dirichlet boundary conditions, the techniques therein can be employed in the Neumann boundary conditions) and it has certain nonlinearities on the source terms. As an example of such system we also discuss model equations arising in experiments of morphogen. In particular it is of our interest to study the Dpp-Sog system \cite{KNF,LNW}, see Section 3, which has certain growth rates and nonlinearities, see \emph{structure conditions} for the general system in the following section. We note that to understand the stabilities of the reaction-diffusion system, in this paper, we focus on the steady state solutions. This paper is organized in two parts; first we state structure conditions and establish an existence result for a class of elliptic systems. In the second part of paper we discuss the model problem and outline a proof by constructing upper and lower solutions to the model system. \section{A class of systems}\label{exsec} In this section we show the existence of positive solutions to the following elliptic system \begin{gather*} \Delta u + f(x,u,v,w)=0 \quad \mbox{in } \Omega, \\ \Delta v + g(x,u,v,w) =0 \quad \mbox{in } \Omega, \\ \Delta w + h(x,u,v)=0\quad \mbox{in } \Omega, \\ \frac{\partial u}{\partial n} = \frac{\partial v}{\partial n} = \frac{\partial w}{\partial n} = 0 \quad \mbox{on } \partial \Omega , \end{gather*} where the source terms $f$, $g$ and $h$ satisfy the following \textbf{structure conditions}: \begin{itemize} \item[(F)] $f(x,p,q,r)$ is $L^P$ in $x\in\Omega\subset \mathbb{R}^N$ with $P>N$ and $f(x,p,q,r)$ is locally Lipschitz in $p,q,r$. For $p,q,r \ge 0$, \begin{equation}\label{f0} f(x, p,q,r) = f^1(p,q,r) + f^2(x,p,q,r), \quad f(x,0,q,r)\ge f_0(x)\ge 0, \end{equation} \begin{equation}\label{f1} f^1(p,q,r) \le - \lambda p^\gamma \end{equation} where $\lambda >0$ and $\gamma >1$, and \begin{equation}\label{f2} |f^2(x,p,q,r)|\le |\tilde f(q,r)|(|p|+1) +C_f,\quad {\rm and} \quad 0 \le f_0(x) \le C_f, \end{equation} with a positive constant $C_f$. For $(x,p,q,r)\in \overline\Omega\times [0,\infty)^3$, $f(x,p,q,r)$ is nondecreasing in $q$ and non-increasing in $r$. \item[(G)] $g(x,p,q,r)$ is $L^P$ in $x\in\Omega\subset \mathbb{R}^N$ with $P>N$ and $g(x,p,q,r)$ is locally Lipschitz in $p,q,r$. \begin{equation}\label{g} g(x,u,v,w)= g(u,v,w)+ g_0(x). \end{equation} For $(x,p,q,r)\in \overline\Omega\times [0,\infty)^3$, there exist positive constants $C_g$ and $g_1$ such that \begin{equation}\label{g1} 0\ge g(p,q,r)\ge -C_g \quad \mbox{and} \quad 0\le g_0(x)\le g_1. \end{equation} \item[(H)] $h(x,p,q)$ is $L^P$ in $x\in\Omega\subset \mathbb{R}^N$ with $P>N$ and $h(x,p,q)$ is locally Lipschitz in $p,q$. For $(x,p,q)\in \overline\Omega\times [0,\infty)^2$, $h(x,p,q)$ is nondecreasing in $p$ and non-increasing in $q$. \end{itemize} The condition (F) implies the source term can be separated in two parts based on the growth rates in terms of $u$. We note that the conditions \eqref{f1} and \eqref{f2} can be generalized further, namely the result still holds if we replace the assumptions \eqref{f1} and \eqref{f2} to \begin{equation}\label{flim} \lim_{p\rightarrow +\infty} f(x,p,\cdot,\cdot)p^{-\gamma}=-\lambda \end{equation} with $\lambda >0$ and $\gamma >1$. In fact the conditions (F)(G) and (H) implies that the existence result of the system is governed by the first equation and the other two equations act almost like ``shadow" equations where such terminology was introduced by \cite{Ni}. We use the upper-lower solutions method and the Schauder fixed point theorem to establish the existence result \cite{S, GT}. We point out that the upper-lower solutions are not necessary in a classical sense, namely we allow them to be distributional solutions for the corresponding inequalities. That is, we say $u$ is a lower (upper) solution if $u\in C(\overline\Omega)$ satisfies \[ \Delta u +F(x) \ge (\le)0 \quad\mbox{in } \mathcal{D}'(\Omega),\quad \frac{\partial u}{\partial n}\big|_{\partial\Omega} \le(\ge)0\mbox{ a.e.}. \] For $x_0\in \partial\Omega$ where the normal derivative is undefined we impose \[ \frac{\partial u}{\partial n}(x_0)\equiv \lim_{x\rightarrow x_0} \sup (\inf) \frac{u(x_0)-u(x)}{|x_0-x|} \le (\ge) 0, \] where $x_0-x$ and the normal at $x_0$ is less than $\pi/2 -\delta$ for some fixed $\delta>0$. We now establish the following existence theorem. \begin{theorem}\label{ethm} For a given bounded open domain $\Omega\subset \mathbb{R}^N$ and if $\Omega$ is a region of the class $W^{2,P}$ then there exist $u_i, v_i, w_i \in C(\overline\Omega)$ ,$i=1,2$, upper-lower solutions in a distributional sense, and exist $u,v,w\in W^{2,P}(\Omega)\cap C^{1+\alpha}(\overline\Omega)$-solutions such that $0=u_1 < u < u_2$, $0\le v_1 \le v\le v_2$ and $0\le w_1 \le w\le w_2$ in $\Omega$. \end{theorem} \begin{proof} We first construct a set of upper-lower solutions in a distributional sense to the system. We find $v_i\in C(\overline\Omega)$, $i=1,2$, lower and upper solutions respectively. Since $0 \ge g(u,v,w) \ge -C_g$ provided $u,v,w\ge 0$ we find $v_1$ (the lower solution to $v$) to be a positive solution to $$ \Delta v_1 - C_g + g_0(x) \ge 0 \ {\rm in}\ \ \mathcal{D}'(\Omega), \quad \frac{\partial v_1}{\partial n}\big|_{\partial\Omega} \le 0 \mbox{ a.e.}. $$ We find $v_2 >v_1$ (the upper solution to $v$) to be a positive solution to $$ \Delta v_2 + g_0(x)\le 0 \ {\rm in }\ \mathcal{D}'(\Omega), \quad \frac{\partial v_2}{\partial n}\big|_{\partial\Omega} \ge 0\mbox{ a.e.}. $$ This can be done by choosing integration constants for $v_1$ and $v_2$ correspondingly. We now fix $v_i$, $i=1,2$. Now with $v_2$ we find a positive solution $w_1$ (the lower solution to $w$) which satisfies $$ \Delta w_1 + h(x,0,v_2) \ge 0\quad\mbox{in }\mathcal{D}'(\Omega), \quad \frac{\partial w_1}{\partial n}\big|_{\partial\Omega} \le 0 \mbox{ a.e.}. $$ With $v_2$ and $w_1$ we now find a positive solution $u_2$ (the upper solution to $u$) to $$ \Delta u_2 +f(x,u_2, v_2,w_1)\le 0 \quad\mbox{in }\mathcal{D}'(\Omega), \quad \frac{\partial u_2}{\partial n}\big|_{\partial\Omega} \ge 0 \mbox{ a.e.}. $$ Let $\phi$ be the first eigenfunction satisfying $\Delta \phi+ \lambda_1\phi =0$ with $\partial \phi/\partial n = 0$ on $\partial \Omega$ $\lambda_1>0$, and $\|\phi\|_\infty =1$. Then by letting $u_2 = K(\phi+2)\ge 1$ with a constant $K>1$ to be determined we show $u_2$ satisfies the last inequalities. More precisely since we have $f^1(p,\cdot,\cdot)\le -\lambda p^\gamma $ with $\gamma>1$ and thus \begin{align*} &\Delta u_2 +f(x,u_2,v_2,w_1) \\ &\le -K\lambda_1\phi -\lambda (K(\phi+2))^\gamma + \sup \tilde f(v_2, w_1)(K\phi +2K+1)+ C_f \\ &\le K\lambda_1 -\lambda K^\gamma + \sup \tilde f(v_2, w_1)(3K+1)+ C_f < 0 \end{align*} by taking $K\ge K_1(\lambda_1,\gamma,\lambda,\sup \tilde f(v_2, w_1),C_f)>0$. In fact $u_2$ is the upper solution in a classical sense. We now find $w_2$ (the upper solution to $w$) such that $w_2 >w_1$ and it is a positive solution to $$ \Delta w_2 +h(x,u_2,v_1)\le 0,\ \mathcal{D}'(\Omega), \quad \frac{\partial w_2}{\partial n}\big|_{\partial\Omega} \ge 0 \mbox{ a.e.}. $$ Finally since $w_2>0$ and $v_1>0$, we let $u_1\equiv 0$ so that $u_1$ (the lower solution to $u$) satisfies $$ \Delta u_1 +f(x,u_1,v_1,w_2)\ge f_0(x) \ge 0. $$ We now define a set $S\subset C^1(\overline\Omega)\times C^1(\overline\Omega)\times C^1(\overline\Omega)$; \begin{align*} S=\big\{&(u,v,w): u_1\le u\le u_2, \ v_1\le v \le v_2,\ w_1\le w \le w_2 \mbox{ in }\Omega,\\ &\|u\|_{C^1(\overline\Omega)}\le A,\; \|v\|_{C^1(\overline\Omega)}\le B,\; \|w\|_{C^1(\overline\Omega)}\le C\\ &\frac{\partial u}{\partial n}=\frac{\partial v}{\partial n} =\frac{\partial w}{\partial n}=0 \; \mbox{on } {\partial \Omega} \big\}. \end{align*} The set $S$ is clearly closed, bounded and convex. Define a map $T$ on $S$ such that \begin{gather} \label{mapu} \Delta Tu - M Tu +f(x,u,v,w) = -M u\\ \label{mapv} \Delta Tv + g(x,u,v,w)=0\\ \label{mapw} \Delta Tw +h(x,u,v)=0\\ \label{BCs} \frac{\partial Tu}{\partial n}=\frac{\partial Tv}{\partial n} =\frac{\partial Tw}{\partial n}=0 \quad\mbox{on } {\partial \Omega}\,, \end{gather} where $M$ is a positive constant so that $f_u +M \ge 0$ for $u,v,w,\in S$ and $f_u$ is the Lipschitz constant of $f$ in $u$. Such $M$ can be found independently to $u$ since $u_1 \le u\le u_2$ and $f\in C^{0,1}$. Since $f$, $g$ and $h$ are uniformly bounded in $L^\infty$ with respect to $u,v$ and $w$, there exist unique solutions $Tu$, $Tv$ and $Tw$ in $W^{2,P}\cap C^{1,\beta}(\overline\Omega)$ with $\beta =1-N/P$, see for the existence of the unique solution in \cite[Proposition 7.18]{Gary} and the solution space \cite[Theorem 7.26]{GT} to \eqref{mapu}, \eqref{mapv} and \eqref{mapw} correspondingly, and thus the map $T$ is well-defined. We first show the map $T$ satisfies the first inequalities. First since $v,w \in S$, evaluate $u_1=0$ we get \[ \Delta u_1 +f(x, u_1,v,w) \ge f_0(x) \ge 0 \] and \begin{align*} 0&\le \Delta (u_1-Tu) -M(u_1-Tu) - Mu + f(x,u_1,v,w)-f(x,u,v,w)\\ &\le \Delta (u_1-Tu)-M(u_1-Tu) + (f_u+M)(u_1-u)\\ &\le \Delta (u_1-Tu)-M(u_1-Tu), \end{align*} since $u\in S$ and $f_u +M \ge 0$ by the choice of $M>0$. In fact since $u_1$ satisfies the last inequality point-wise and also holds the zero Neumann boundary condition point-wise as well, we can apply the strong maximum principle \cite[Theorem 3.5]{GT} to get $u_1 -Tu <0$ or $u_1 -Tu\equiv c$ for some constant $c$. Since $M>0$ the constant $c$ must be zero and since $\Delta u_1 \not\equiv \Delta Tu$ thus we get $Tu >u_1$ in $\overline\Omega$. To show $u_2$ be an upper solution, we evaluate \begin{align*} 0& \le \Delta (Tu-u_2) +f(x,u,v,w)-f(x,u_2,v_2,w_1) +Mu -MTu\\ &\le \Delta (Tu-u_2) - M(Tu-u_2) + (f_u + M)(u-u_2)+f_v(v-v_2) +f_w(w-w_1)\\ &\le \Delta (Tu-u_2) -M (Tu-u_2), \end{align*} since $u,v,w \in S$, $f_u +M\ge 0$, and $f_v \ge 0$ and $f_w \le 0$ where $f_v$ and $f_w$ are some bounded functions. Now since $\partial (u_2 -Tu)/\partial n \ge 0$ on $\partial \Omega$ we apply the same argument as before to get $u_2 > Tu$ in $\overline\Omega$. Therefore $u_1 < Tu < u_2$ in $\overline\Omega$. We also show $v_1\le Tv\le v_2$ and $w_1 \le Tw \le Tw_2$ by using $u,v,w\in S$ and the differential inequalities of $v_i$ and $w_i$. More precisely as before we have $\Delta (v_1 -Tv)\ge 0$ in $\mathcal{D}'(\Omega)$. By the weak Maximum principle \cite[Corollary 3.2]{GT} we get $\sup_\Omega (v_1-Tv) \le \sup_{\partial \Omega} (v_1 -Tv)^+ .$ Now suppose there exist a constant $k>0$ and a point $x_0\in\partial\Omega$ such that $$ \sup_{\partial\Omega}(v_1 -Tv)^+=(v_1-Tv)(x_0)=k. $$ Since $v_1$ and $Tv$ are continuous, we can find a set $\Omega_k$ such that $\Omega_k=\{x\in\overline\Omega: (v_1-Tv) \frac{\partial Tv}{\partial n}(x_0)=0. \] The contradiction is apparent. Thus there is no such $k>0$,(hence $\sup_{\partial\Omega}(v_1 -Tv)^+=0$) and thus $Tv \ge v_1$. Similarly we get the rest of inequalities, and for $i=1,2$, $Tv\not\equiv v_i$ and $Tw \not\equiv w_i$. This shows the map $T$ satisfies the first inequalities in $S$. (In the case if we find $u_i$,$i=1,2$ as lower/upper solutions in a distributional sense we apply the same arguments as we did for $v_1-Tv$ to get the desired inequalities.) We show the map $T$ is compact, satisfies the second inequalities (this leads $T$ being into), and continuous in $S$ to get a fixed point in $S$. The map is compact since the source term is uniformly bounded in $u$, $v$ and $w$. To be precise, since $u,v,w \in S$ and the source terms $f$, $g$ and $h$ are uniformly bounded in $L^\infty$, thus we apply the $L^P$-theory \cite[Proposition 7.18]{Gary} to obtain the uniform bounds of the solutions $Tu$, $Tv$, and $Tw$ in $W^{2,P}(\Omega)$ for given $P>N$, namely, $\|Tu\|_{W^{2,P}}\le C(N,P,M)\|f\|_{L^P}$, $\|Tv\|_{W^{2,p}}\le C(N,P)\|g\|_{L^P},$ and $\|Tw\|_{W^{2,p}}\le C(N,P)\|h\|_{L^P}$. Apply imbedding theory \cite[Theorem 7.26]{GT} to get the uniform bounds of the solutions $Tu$, $Tv$, and $Tw$ in $C^{1,\alpha}(\overline\Omega)$ where $0<\alpha < 1- N/P$ is independent of solutions. Since we now have $C^{1,\alpha}(\overline\Omega)$ bounds uniformly in $Tu$, $Tv$ and $Tw$ we simply let their uniform bounds to $A$, $B$ and $C$ respectively. Thus $T$ maps $S$ into. Furthermore, $T(S)\subset C^{1,\alpha}(\overline\Omega)$ which is precompact in $C^1(\overline\Omega)$ with $0<\alpha<1-N/P$. This leads $T$ is compact in $S$. Finally to show the map $T$ is continuous, we take convergence sequences $u_i$, $v_i$ and $w_i$ and show that the sequences $Tu_i$, $Tv_i$ and $Tw_i$ have limits in $S$. Calculate \begin{align*} &\Delta (Tu_i -Tu_j) -M(Tu_i-Tu_j) + f(x,u_i,v_i,w_i)- f(x,u_j,v_j,w_j)+M(u_i-u_j) \\ &=\Delta (Tu_i -Tu_j) -M (Tu_i-Tu_j)+ (f_u+M)(u_i-u_j)\\ &\quad + f_v(v_i-v_j) + f_w (w_i-w_j) \end{align*} Since $Tu_i$, $v_i$, $w_j$ are in $S$ and $M>0$ the $L^P$-theory \cite{Gary} leads to $$ |Tu_i-Tu_j|_{W^{2,P}} \le C(|u_i-u_j|_{L^\infty} + |v_i - v_j|_{L^\infty}+|w_i-w_j|_{L^\infty}) $$ and this implies $Tu_i$ is a Cauchy sequence in $W^{2,P}(\Omega)$ and thus (by imbedding) the sequence $Tu_i$ has a limit in $S$. Similarly $v_i$ and $w_i$ have limits in $S$ as well. Therefore there exists a fixed point $Tu=u$, $Tv=v$ and $Tw=w$ in $S$. Now apply regularity arguments \cite{GT} to obtain $W^{2,P}(\Omega)\cap C^{1,\alpha}(\overline\Omega)$-solutions. This completes the proof. \end{proof} We note that the regularity of the solutions can be improved to $C^{2,\alpha}(\overline\Omega)$ if we allow \\ $f(x,\dots), g(x,\dots), h(x,\dots) \in C^{0,\beta}(\overline\Omega)$ and $\partial \Omega\in C^{2,\gamma}$ for some $0<\beta, \gamma<1$. \section{An example: Morphogen gradients}\label{egsec} In this section we present a biological example of the system where the details of modeling viewpoints can be found in \cite{KNF,LNW,L,MSC,PA,r}. In multicellular systems, experiments suggest that the {\it Drosophila} wing disc of fly depends on the {\it decapentaplegic} (Dpp) gene \cite{L}. In development of the Dpp concentration, proteins including {\it short gastrulation} (Sog) activate to inhibit Dpp. As a result of the activation of Sog, Dpp decreases its activity \cite{MSC,PA,r}. It is of our interest to understand mathematical structures, in particular the steady state solutions, on the model problem of the Dpp-Sog system. In fact, the model problem also includes effects of the co-inhibitor {\it twisted gastrulation} (Tsg), extracellular protease {\it tolloid} (Tld) and a second ligand {\it screw} (Scw) \cite{MSC,PA}. For a mathematical simplification, this paper focuses only on the effect of Dpp and Sog. It is noted by \cite{KNF} that although the mathematical simplification may reduce the biological factors the numerical simulation suggests that the simplified system preserves the fundamental features of the model system. More precisely we consider the following system on the interval $\Omega\equiv (0,1)\subset R$ \begin{gather*} \frac{\partial A}{\partial t}=\Delta A - h_L A(1-B) -h_{LS} AD +f_L B+(f_{LS} +g_{LS})C +v_{OL}\\ \frac{\partial B}{\partial t}= h_L A(1-B)-(f_L +g_L)B \\ \frac{\partial C}{\partial t}= \Delta C +h_{LS} AD -(f_{LS} +g_{LS})C\\ \frac{\partial D}{\partial t}=\Delta D -h_{LS} AD +f_{LS} C + v_{OS} \end{gather*} with zero Neumann boundary conditions. Here $A$ and $B$ are the concentrations of free ligand and of receptor-bound ligand of Dpp, respectively, and $C$ and $D$ are the concentrations of the degradation of the bound complex and of the destruction of the inhibitors Sog, respectively. Coefficients $h_L$ {\it etc} are positive constants (biological factors) and $v_{OL} =v_{L} H(1/2-x)$ and $v_{OS}=v_{S} H(x-1/2)$ where $v_{L}$ and $v_{S}$ are positive constants and $H$ is the Heaviside function, where the Heaviside functions incorporate that ligand $A$ is produced on the half of the domain, {\sl i.e.}, in the dorsal half of the embryo, and the new factor $D$ is produced on the other half of the domain, {\sl i.e.}, the ventral half of the embryo, see \cite{KNF,MSC,PA} for details. To understand the stabilities of the governing system, we consider the steady state system: \begin{gather}\label{A} \Delta A - h_L A(1-B) -h_{LS} AD +f_L B+(f_{LS} +g_{LS})C +v_{OL}=0\\ \label{B} h_L A(1-B)-(f_L +g_L)B=0 \\ \label{C} \Delta C +h_{LS} AD -(f_{LS} +g_{LS})C=0\\ \label{D}\Delta D -h_{LS} AD +f_{LS} C + v_{OS}=0 \end{gather} with zero Neumann boundary conditions. Now notice that using equation \eqref{B} we can write \eqref{A} to \begin{equation}\label{A1} \Delta A -h_{LS} AD -g_L B +(f_{LS} +g_{LS})C +v_{OL}=0. \end{equation} Define $v=A+C$ and denote $u=A$ so that $C=v-u$ and from \eqref{A1} and \eqref{C} to get \begin{equation}\label{P} \Delta v - g_L B +v_{OL}=0. \end{equation} Also define $w=C+D$ such that $D=w-v+u$ and from \eqref{A1} and \eqref{D} to get \begin{equation}\label{Q} \Delta w - g_{LS} (v-u) +v_{OS} =0. \end{equation} From \eqref{B} we get \begin{eqnarray*} B=\frac{h_L A}{h_L A + f_L +g_L}=\frac{u}{u+\alpha}\equiv B(u) \end{eqnarray*} where $\alpha = (f_L +g_L)/h_L>0$, and thus $B$ is increasing in $u$. Finally we obtain an equivalent system in the following; \begin{gather} \label{eqA}\Delta u - h_{LS} u(w-v +u) -g_L B(u) +(f_{LS} +g_{LS})(v-u) +v_{OL}=0\\ \label{eqP} \Delta v - g_L B(u) +v_{OL}=0\\ \label{eqQ} \Delta w - g_{LS} (v-u) +v_{OS}=0\\ \label{bcs} \frac{\partial u}{\partial n} =\frac{\partial v}{\partial n} =\frac{\partial w}{\partial n} =0 \end{gather} Note that in \eqref{eqA} for $u,v,w\ge 0$ we have \begin{gather*} f(x,0,v,w)= (f_{LS}+g_{LS})v +v_{OL}\ge v_{OL} \ge 0\,,\\ f^1(u,v,w)=-h_{LS} u^2 -h_{LS} uw -g_L B(u) -(f_{LS}+g_{LS})u < -h_{LS} u^2\,,\\ f^2(x,u,v,w)= h_{LS} uv + (f_{LS}+g_{LS})v +v_{OL}, \end{gather*} and clearly all the conditions in (F) hold. Also in \eqref{eqP} and in \eqref{eqQ}, we have $g(x,u,v,w) =-g_L B(u) +v_{OS}$ which satisfies the conditions in (G) provided $u\ge 0$, and $h(x,u,v)=- g_{LS} (v-u) +v_{OS}$ which holds the conditions (H) as well. Thus we establish the existence of positive solutions to the system \eqref{eqA}-\eqref{eqQ} in the following theorem. Since the proof follows exactly as in the Theorem ~\ref{ethm} we only construct upper-lower solutions explicitly. \begin{theorem}\label{eg} There exist positive solutions to the steady state system \eqref{eqA}--\eqref{eqQ} and \eqref{bcs}. \end{theorem} \begin{proof} Since the proof follows exactly as in Theorem ~\ref{ethm} we only construct a set of upper-lower solutions to the system. We first find $v_i$, $i=1,2$, distributional lower and upper solutions respectively. Since $0 \le B(u)<1$ provided $u\ge 0$ we find $v_1$ (the lower solution to $v$) to be a positive solution to $$ v_1'' - g_L + v_{L} H(1/2-x) \ge 0\,, $$ where $H$ is the Heaviside function. Set $$ v_1=\begin{cases} g_L\frac{x^2}{2} & x\in[0,1/2]\\ g_L \frac{(x-1)^2}{2} & x\in (1/2,1] \end{cases} $$ so that $v_1$ is positive and continuous on $(0,1)$, and $\partial v_1/\partial n =0$ at $x=0$ and $x=1$. Since $v_{L}>0$ and $v_1'' = g_L$ {\sl a.e} in $(0,1)$, the inequality holds. We find $v_2$ (the upper solution to $v$) to be a positive solution to $$v_2'' + v_{L} H(1/2-x)\le 0.$$ Similar calculation as before we set $$ v_2=\begin{cases} -v_{L}\frac{x^2}{2} + c_1 & x\in[0,1/2]\\ -v_{L} \frac{x^2}{2} +v_{L} x+ c_2 & x\in (1/2,1] \end{cases} $$ so that $\partial v_2/\partial n =0$ at $x=0$ and $x=1$ and $v_2'' =-v_{L}$ {\sl a.e.} in $(0,1)$. We can find integration constants $c_1$ and $c_2$ so that $v_2(1/2)= -\frac{v_{L}}{8} + c_1= -\frac{v_{L}}{8} + \frac{v_{L}}{2} + c_2$ to get $v_2$ is continuous on $[0,1]$ and $v_2 > v_1$ that is $c_1$ and $c_2$ satisfy $$ c_1 > \frac{v_{L}}{8} + \frac{g_L}{8}, \quad \mbox{and}\quad c_2 > \frac{3v_{L}}{8}+\frac{g_L}{8}. $$ We now fix the integration constants $c_i$, $i=1,2$. With the $v_2$ we just found, we look for a positive solution $w_1$ (the lower solution to $w$) to $$ w_1'' -g_{LS} v_2 +v_{S} H(x-1/2) \ge 0. $$ Set $$ w_1=\begin{cases} g_{LS} \max v_2\frac{x^2}{2} & x\in[0,1/2]\\ g_{LS} \max v_2 \frac{(x-1)^2}{2} & x\in (1/2,1] \end{cases} $$ so that $\partial w_1/\partial n =0$ at $x=0$ and $x=1$ and $w_1'' =g_{LS} \max v_2$ {\sl a.e} and $w_1$ is positive and continuous on $(0,1)$. With $v_2$ and $w_1$, we now find a positive solution $u_2$ (the upper solution to $u$). By letting $u_2 = K(\cos(\pi x)+2)\ge 1$ with some large $K>1$ we can see that $u_2$ satisfies \begin{align*} &u_2'' +f^1(u_2, v_2, w_1) + f^2(x,u_2, v_2, w_1)\\ & \le u_2'' - (h_{LS} \min w_1 + f_{LS} +g_{LS})u_2 - h_{LS} u_2^2 \\ &\quad + h_{LS}\max v_2 u_2 +(f_{LS} +g_{LS})\max v_2 +v_{L} H(1/2-x) \\ &\le -K\pi^2 \cos(\pi x) - h_{LS} K^2 (\cos(\pi x)+2)^2\\ &\quad + h_{LS} \max v_2 K(\cos(\pi x)+2) +(f_{LS} +g_{LS})\max v_2 +v_{L}\\ &\le - h_{LS} K^2 + (3h_{LS} \max v_2 +\pi^2) K +(f_{LS} +g_{LS})\max v_2 +v_{L} < 0 \end{align*} by taking $K=\max\{ 6h_{LS} \max v_2 + 2\pi^2, [h_{LS}^{-1}2((f_{LS} +g_{LS})\max v_2 +v_{L})]^{1/2}\}$. We now find $w_2$ (the upper solution to $w$) to be a positive solution to $$ w''_2 -g_{LS} v_1+ g_{LS} u_2 +v_{S} H(x-1/2) \le 0. $$ Since $\min v_1 =0$ and $\max u_2 =2K$ we can let $w_2$ be a positive solution to $$w''_2 + g_{LS} 2K +v_{S} H(x-1/2) \le 0.$$ Again we let $$ w_2=\begin{cases} -(g_{LS} 2K +v_{S})\frac{x^2}{2} +d_1 & x\in[0,1/2]\\ -(g_{LS} 2K +v_{S})(\frac{1}{2} x^2 - x) + d_2 & x\in (1/2,1] \end{cases} $$ so that $\partial w_2/\partial n =0$ at $x=0$ and $x=1$ and $w_2'' =- (g_{LS} 2K +v_{S})$ {\sl a.e.} in $(0,1)$. This again brings two integration constants and so we choose them such that $w_2 >w_1$ and $w_2$ is continuous on $[0,1]$. Namely, we find $d_1 > (g_{LS} 2K +v_{S})/8 + (g_{LS} \max v_2)/8$ and $d_1 = (g_{LS} 2K +v_{S})/2 +d_2$ where $d_2$ holds $d_2 > 3(g_{LS} 2K +v_{S})/8 + (g_{LS} \max v_2)/8$. Finally since now we have $v_1 \ge 0$ it is easy to see that $f(x,0,v_1,w_2)=(f_{LS} +g_{LS})v_1 + v_{L} H(1/2-x) \ge v_{L} H(1/2-1) \ge 0$ and thus we let $u_1 =0$ so that $u_1'' +f(x,0,v_1,w_2) \ge 0$. Clearly $0=u_1