\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