COVID-19 propagation and the usefulness of awareness-based control measures: A mathematical model with delay

1.   Introduction

The FitzHugh-Nagumo equation

where the nonlinearity $f(u) = u(u-a)(1-u)$, $0 < a < \frac{1}{2}$, $\gamma > 0$, and $0 < b\ll 1$, was originally proposed as a simplification of the Hodgkin-Huxley model of nerve axon dynamics ^[8,24]. It has been the subject of a number of papers as a model of nerve conduction in recent years ^{[1,2,3,4,5,6,9,10,11,13,14,15,16,17,20,23,24,25]}. The FitzHugh-Nagumo equation has been used not only as a set of equations for the qualitative description of nerve axon behavior, but also for excitable media in general ^[12,26]. Traveling wave solutions are of special interest.

There are two main trends of research of traveling wave solutions in FitzHugh-Nagumo equation. In the concept based on the geometric singular perturbation theory ^[4,5,19,20], the presence of a small parameter in the right side of the kinetic equation is essentially used. In another approach based on the phase diagram analysis ^[16,17], the phase trajectories of a dynamical system associated with the initial system of PDEs are considered, and their dependence on the parameters of the system are studied.

Recently, Deng ^[7] has derived a mathematical model for neuron by imposing only a principle of symmetry that two modelers must obtain the same model when one models the conductances of ion channels and the other models the channels' resistances. This model is called neuron model with conductance-resistance symmetry, denoted by acronym CRS neuron model. It has many advantages. Firstly, the model gives an explanation to the leak current discovered by Hodgkin-Huxley^[18]. Secondly, the model has a better fit to the experimental data than the Hodgkin-Huxley model does. Thirdly, how such N-shaped nonliearity arises in neuroscience has always been a puzzle ^[8,21] but for this model it is a simple consequence to the underlining bias-free symmetries. In addition, the model can be reduced to a two-dimensional model similar to the FitzHugh-Nagumo equation.

In this paper, we are interested in the existence of solitary wave solution for a two-dimensional reduction model

of the CRS neuron model in a propagated action potential. See Section 2 for the detailed expression of (1.1).

System (1.1) have some apparent similarities to the FitzHugh-Nagumo equation. Both of these traveling wave equations have the N-shaped nonliearity structure, and the origin is their equilibrium point. However, there is much less known about qualitative properties of the system (1.1) than to the FitzHugh-Nagumo equation, due to structural differences between (1.1) and the FitzHugh-Nagumo equation. For instance, for $b = 0$, the fast front and the fast back of the traveling wave equations of the FitzHugh-Nagumo equation can be calculated precisely by taking the appropriate wave velocities $c$ and appropriate $n$ ^[23], while such fast orbits of the Eq (1.1) is very difficult to calculate. Therefore, it is difficult to obtain the existence of solitary wave solution only by using geometric singular perturbation theory and exchange lemma. As a consequence, we use phase diagram analysis to avoid this problem directly.

To the best of our knowledge, the traveling wave solutions of the neural model with conduction-resistance symmetry and its reduction model have not been considered in the literature. Our work may be the first.

In Section 2, we will give the derivation of the conductance-resistance symmetry neuron model from Deng ^[7] to the system (1.1). In Section 3, we investigate the neuron model (1.1). The existence of a traveling pulse solution is proved.

2.   Neuron model with conductance-resistance symmetry

In this section, we shall briefly sketch propagation of the excitation along the nerve axon following Deng ^[7]. Then, we will reduce this neuron model to the 2-dimensional system (1.1).

2.1.   List of symbols

2.2.   Neuron model with CR-symmetry

This is the neuron model with conductance-resistance symmetry ^[7]:

where $I_{L} = \varepsilon_{K}\bar{g}_{K}(V-E_{K})+\varepsilon_{Na}\bar{g}_{Na}(V-E_{Na})+\varepsilon_{G}\bar{g}_{G}(V-E_{G})$, which can be interpreted as the leak current discovered by Hodgkin-Huxley ^[18].

The fact that the local circuit currents have to be provided by the net membrane current leads to the well-known relation ^[18]:

where $i$ is the membrane current per unit length, $r_{1}$ and $r_{2}$ are the external and internal resistances per unit length, and $\bar{x}$ is distance along the fibre. For an axon surrounded by a large volume of conducting fluid, $r_{1}$ is negligible compared with $r_{2}$. Hence, we have

or

where $I$ is the membrane current density, $a$ is the radius of the fibre and $R_{2}$ is the specific resistance of the axoplasm. Inserting this relation in Eq (2.1) and rescale the spacial variable $\bar{x} = x\sqrt{\frac{a}{2R_{2}}}$, we have:

This basic model can be simplified further in three aspects by Deng ^[7]. Firstly, the above model can afford to drop the small leak current $I_{L}$ without affecting the qualitative properties of the system. Secondly, since gated protein current is generated by sodium pore deformation caused by depolarization or hyperpolarization voltage, we can assume that the rate constant $\alpha_{G}$ approaches infinity, that is, the time evolution from $h$ to $\phi_{G}$ instantaneously. Finally, among the rate parameters $\alpha_{K}$ and $\alpha_{Na}$, the former needs to be at least one order of magnitude smaller than the latter. This indicates that sodium kinetics is faster than potassium kinetics. So we can assume that sodium kinetics is instantaneous with $\alpha_{Na} = \infty$. As a result, the 4-dimensional model (2.2) is reduced to the following 2-dimensional system

Using the data show in the list:

Let us normalize this equations, that is, let $u = \frac{V+55}{130}$, then the equations becomes

where

$\psi_{K}(u)$ and $\psi_{Na}(u)$ are increasing functions for $u\geq1/65$ and zero for $u\leq1/65$, $\psi_{G}(u)$ is a decreasing function for $u\leq 1$ and zero for $u\geq1$. $\psi_{X}:\mathbb{R}\rightarrow [0, 1)$, $X = K$, $Na$, $G$. $g$ is a decreasing function for $(-\infty, 0]\cup[1, +\infty)$. $g$ tends to positive infinity as $u$ tends to negative infinity and tends to negative infinity as $u$ tends to positive infinity. Further, $g(0) = g(1) = 0$ and $g'(0) = g'(1) < 0$. However, because of the complexity of the $g$ expression, it is difficult to get the exact monotonicity of $g$ in the interval $(0, 1)$. It is clear that $g$ is continuous, whose figure can be well depicted by numerical simulation, where $g(\frac{1}{65}) = -\frac{1}{65}\psi_{G}(\frac{1}{65}) < 0$ (see Figure 1). Based on the above analysis, $g(u)$ has N-shaped nonliearity which is combined with the sodium and the gating characteristic curves. And $g(u)$ has three zeros: $u = 0$, $a$, and $1$, where $\frac{1}{65} < a < \frac{1}{2}$ is a fixed constant.

For readability, we replace the parameter $\alpha_{K}$ with $b$, and then get the system (1.1). In this paper, we mainly prove the existence of the solitary wave solution for (1.1).

3.   Existence of solitary wave solutions

In this section, we aim to establish the existence of solitary wave solutions for Eq (1.1). As a result, we consider the traveling wave solutions $u(t, x) = u(\xi)$, $n(t, x) = n(\xi)$, $\xi = x+ct$, where $c > 0$ is the wave velocity. Finding such solutions to (1.1) is equivalent to finding solutions of the following system of ODEs:

The origin $O(0, 0, 0)$ is the unique equilibrium point of the system. Note that the homoclinic orbit of system (3.1) corresponds to the solitary wave solution of (1.1). So we only need to prove the existence of homoclinic orbit in system (3.1).

Because the vector field of the system is relatively complex, the nullcline of the system (3.1) on the $(u, n)$-plane is graphically represented for the convenience of readers' subsequent reading and understanding, see Figure 2.

3.1.   Local invariant manifolds of the origin

If the system (3.1) is linearized around the rest point $O(0, 0, 0)$, then linearization matrix takes the form

The characteristic equation of the matrix $A$ is as follows:

It is easy to verify that $\Delta(\lambda)$ has one positive real root $\lambda_{1}$ and two negative real roots $\lambda_{2, 3}$, where

and the eigenvectors corresponding to $\lambda_{i}$ take the form $Y_{i} = \left(1, \lambda_{i}, 0\right)$ where $i = 1, 2$, and the eigenvector corresponding to $\lambda_{3}$ takes the form $Y_{3} = \left(1, -\frac{b}{c}, \eta\right)$, where $\eta = \frac{26}{35}(b+\frac{b^2}{c^2}-2\psi_{G}(0))$. Thus, there exists a one-dimensional local invariant unstable manifold $W_{loc}^{u}$ tangent to the eigenvector $Y_{1}$ at the origin and a two-dimensional local invariant stable manifold $W_{loc}^{s}$ tangent to the plane spanned by the vectors $Y_{2}$, $Y_{3}$.

3.2.   Behavior of saddle separatrices in case $b=0$

In this section, we consider the condition $b = 0$, allowing to separate the first two equations from the third one, having the trivial solution $n = {\rm{const}}$. In particular, let's consider the case $n = 0$. The remaining system then takes the form

For any $c > 0$, it is obvious that this system has three equilibrium points $(0, 0)$, $(a, 0)$ and $(1, 0)$. The points $(0, 0)$ and $(1, 0)$ are saddle points, while $(a, 0)$ is either an unstable node or an unstable focus. The unstable manifolds at $(0, 0)$ and $(1, 0)$ have positive slopes. For $c = 0$, the above system can be presented in the Hamiltonian form with the Hamiltonian function

from which and easy analysis of phase space, we get the following:

Proposition 1. If $c = 0$, then in the $(u, v)$ phase plane, system (3.2) has a nonconstant solution $\pi = (u, v)$ such that $\pi(\pm\infty) = (0, 0)$. Also, the trajectories all lie on curves of the form $H(u, v) = constant$. (see Figure 3).

Proof. Since

$H$ is constant along every solution of system (3.2). This means that the solutions of system lie on the level sets of $H$. All we need to do is figure out the directions of the solutions curves on these level sets, which is easy since we have the vector field.

Let $L_{c}^{u}$ denotes the branch of the unstable manifold at $(0, 0)$ which points into the first quartile, and $L_{c}^{u}$ consists of an orbit. We denote by $\pi_{c} = (u_{c}, v_{c})$, the corresponding solution with $u_{c}(0) = a$, $v_{c}(0) > 0$. Using the comparison principle ^[22], it is easy to arrive at the following:

Proposition 2. There exists a unique speed $\bar{c} > 0$ such that the orbit $\pi_{c}(\xi)$ satisfies the following properties:

${\rm{(i)}}$ For $c = \bar{c}$, $\pi_{c}(\xi)$ is a heteroclinic orbit from $(0, 0)$ to $(1, 0)$.

${\rm{(ii)}}$ For $0 < c < \bar{c}$, $\pi_{c}(\xi)$ intersects the horizontal axis at $u^{*}$ such that $a < u^{*} < 1$, then enters the fourth quadrant and its coordinate $v(\xi)$ remains negative. In particular, $\lim_{\xi\rightarrow +\infty}\pi_{c}(\xi) = (-\infty, -\infty)$.

${\rm{(iii)}}$ For $c > \bar{c}$, $\pi_{c}(\xi)$ always stays in the first quadrant and $\lim_{\xi\rightarrow +\infty}\pi_{c}(\xi) = (+\infty, +\infty)$.

3.3.   Behavior of solutions in case $b > 0$

In this section, we consider the behavior of solutions for system (3.1) in case $b > 0$, and establish the following results.

Lemma 3.1. Suppose that $\pi(\xi) = (u(\xi), v(\xi), n(\xi))$ is a solution of the dynamical system (3.1). For any $c > 0$, the following statements are true:

${\rm{(a)}}\ \pi(\xi)$ cannot meet $n = \psi_{K}(u)$ at the point $(u(\xi_{0}), v(\xi_{0}), n(\xi_{0}))$ with $v(\xi_{0}) > 0$, if there exists some $\delta > 0$, such that $0 < n(\xi) < \psi_{K}(u(\xi))$ for $\xi\in(\xi_{0}-\delta, \xi_{0})$.

${\rm{(b) }}$ The sets

are positively invariant.

Proof. (a) Suppose $0 < n(\xi) < \psi_{K}(u(\xi))$ for $\xi\in(\xi_{0}-\delta, \xi_{0})$, then $u(\xi) > \frac{1}{65}$, $\psi_{K}(u(\xi))-n(\xi) > 0$ for $\xi\in(\xi_{0}-\delta, \xi_{0})$. We also have $u(\xi_{0}) > \frac{1}{65}$ because $u'(\xi_{0}) = v(\xi_{0}) > 0$ from (3.1). If $\pi(\xi)$ meets $n = \psi_{K}(u)$ at $\xi_{0}$, then $\left.(\psi_{K}(u(\xi))-n(\xi))'\right|_{\xi = \xi_{0}} = \frac{\partial\psi_{K}}{\partial u}(u(\xi_{0})) v(\xi_{0})\leq0$ contradicts with $u(\xi_{0}) > \frac{1}{65}$, $v(\xi_{0}) > 0$ and $\psi_{K}(u)$ being increasing in the range $u > \frac{1}{65}$. So $\pi(\xi)$ cannot meet $n = \psi_{K}(u)$ at the point $(u(\xi_{0}), v(\xi_{0}), n(\xi_{0}))$, where $v(\xi_{0}) > 0$ if there exists some $\delta > 0$, such that $0 < n(\xi) < \psi_{K}(u(\xi))$ for $\xi\in(\xi_{0}-\delta, \xi_{0})$.

(b) It is worth underlining that $n$ varies between zero and one. So we can get $v'$ = $cv+35n(u+\frac{1}{26})-g(u)$$> 0$ on the set $Q^{+}$ except at the point $(1, 0, 0)$. $v' = 0$, $v'' = cv'+35n'(u+\frac{1}{26})+35nv-\frac{\partial g(u)}{\partial u}v' > 0$ at the point $(1, 0, 0)$. The relation $v = u'$ imply $u(\xi)$ and $v(\xi)$ are strictly monotonically increasing in $Q^{+}$. It turns out that the solution $\pi(\xi)$ will enter the region $Q^{+} = \{(u, v, n): u\geq1, v\geq0\}$ without leaving, i.e., $Q^{+}$ is positively invariant. The positive invariance of the set $Q^{-}$ can be achieved by $v' < 0$ on $u\leq-\frac{1}{26}$, $v\leq0$.

According to the above proof, we can get the following conclusions:

Corollary 1. The sets

are positively invariant. Moreover, if there exists some $\xi_{0}$ such that $\pi(\xi_{0})\in Q^{\pm}$, then $\pi(\xi)\in\tilde{Q}^{\pm}$ for all $\xi > \xi_{0}$ and $\lim_ {\xi\rightarrow +\infty}(u(\xi), v(\xi)) = (\pm\infty, \pm\infty)$.

Proof. We only prove that if there exists some $\xi_{0}$ such that $\pi(\xi_{0})\in \tilde{Q}^{\pm}$ then

And the other conclusions are obviously followed by the proof of Lemma 3.1. According to the first two equations of system (3.1), $u' = v > 0$ and $v'\geq-g(u) > 0$ on $\tilde{Q}^{+}$. If there exists some $\xi_{0}$ such that $\pi(\xi_{0})\in \tilde{Q}^{+}$, then both $u(\xi)$ and $v(\xi)$ are monotonically increasing and tend to positive infinity. Similarly, the case of $\tilde{Q}^{-}$ can be proved.

Corollary 2. If $\pi(\xi) = (u(\xi), v(\xi), n(\xi))$ is a bounded solution of (3.1) on $[0, +\infty)$, then $\pi(\xi)$ will enter the region $\{(u, v, n): -\frac{1}{26} < u < 1\}$, and stay in this region.

Proof. Assume that the orbit $\pi(\xi)$ is bounded, it is clear that $\pi(\xi)$ will neither enter the region $Q^{+}$ nor enter the region $Q^{-}$ by Lemma 3.1 and Corollary 1. Furthermore, if the orbit $\pi(\xi)$ goes into the region $\{(u, v, n): -\frac{1}{26} < u < 1\}$ at some time, it can't meet the half plane $\{(u, v, n):u = -\frac{1}{26}, v > 0\}$ and $\{(u, v, n):u = 1, v < 0\}$ because $u' = v$. Thus, if the orbit $\pi(\xi)$ goes into the region $\{(u, v, n): -\frac{1}{26} < u < 1\}$ at some time, it will remain in the region thereafter. Now it is necessary to show that the orbit will enter region $\{(u, v, n): -\frac{1}{26} < u < 1\}$ in finite time.

Without loss of generality, assume that $\pi(0)$ belongs to the region

According to the proof of Lemma 3.1, we know that $v' < 0$ on the sets $\{(u, v, n):u\leq-\frac{1}{26}, v = 0\}$, $v' > 0$ on the sets $\{(u, v, n): u > 1, v = 0\}$, $v' = 0$ and $v'' > 0$ at the point $(1, 0, 0)$. Therefore, $\pi(\xi)$ will leave $\Theta$ after a finite time. However, due to the boundedness of $\pi(\xi)$, $\pi(\xi)$ will neither enter $Q^{+}$ nor $Q^{-}$, so it will enter the region $\{(u, v, n): -\frac{1}{26} < u < 1\}$.

Moreover, we can get the following conclusion:

Proposition 3. If $\pi(\xi)$ is a bounded solution on $[0, +\infty)$ such that $u(\xi) < a$ for all large $\xi$, then $\lim_{\xi\rightarrow +\infty}\pi(\xi) = (0, 0, 0)$.

Proof. From Corollary 2 and the above conditions it follows that $-\frac{1}{26} < u(\xi) < a$ for all large $\xi$. Without loss of generality, we assume that $-\frac{1}{26} < u(\xi) < a$ on $[0, +\infty)$. To keep the proof covenient, let's make the following notations

Notice that $v = u'$ and $v'$ is always positive in the set $A$. If the orbit $\pi(\xi)$ enters the region $A$, then it will intersect with $u = a$, which contradicts with $-\frac{1}{26} < u(\xi) < a$ on $[0, +\infty)$. So $\pi(\xi)$ will not enter the region $A$. In addition, $\pi(\xi)$ does not tend to the half-plane $\{(u, v, n):0 < u < a, v = 0\}$ as $\xi$ tends to infinity, because $v' > 0$ in this half-plane. Therefore, there are only three possibilities for the solution $\pi(\xi)$:

${\rm{(i) }} \pi(\xi)\rightarrow (0, 0, 0)$;

${\rm{(ii)}}$ The orbit $\pi(\xi) = (u(\xi), v(\xi), n(\xi))$ tends to the half-plane $\{(u, v, n):-\frac{1}{26} < u < 0, v = 0\}$ as $\xi$ tends to infinity. Since $\psi_{K}(\cdot)$ is equal to zero on the set $(-\infty, \frac{1}{65}]$,

Thus, $n(\xi)$ tends to $0$ and $(u(\xi), v(\xi))$ tends to the point $(u_{0}, 0)$ as $\xi$ tends to infinity, where $u_{0}\in(-\frac{1}{26}, 0)$. But it is impossible because the point $(u_{0}, 0, 0)$ is not the equilibrium point of the system (3.1).

${\rm{(iii)}}$ The orbit $\pi(\xi)$ goes back and forth between the region $C$ and the region $D$. Without loss of generality, let's assume that the orbit goes into the region $D$ and after that it enters the region $C$. In this case, the orbit $\pi(\xi) = (u(\xi), v(\xi), n(\xi))$ must satisfies $-\frac{1}{26} < u(\xi) < 0$, $v(\xi) = 0$ and $v'(\xi)\leq0$ at a certain moment $\xi = \xi_{0}$. In addition, note that

Hence, we have $v'(\xi)\leq0$ on $\xi\geq\xi_{0}$. That is to say, $v(\xi)$ is nonincreasing and $\pi(\xi)\in C$ on the interval $[\xi_{0}, +\infty)$. But this is impossible since the orbit goes back and forth between the region $C$ and the region $D$.

In summary, $\pi(\xi)$ can only tend to the origin as $\xi$ tends to positive infinity.

It follows from Proposition 3 that we need to find the parameters $b$ and $c$ that can make the unstable invariant manifold of the origin conform to the conditions of Proposition 3. In the next section, we study asymptotic behavior of the unstable invariant manifold at the origin.

3.4.   Asymptotic behavior of the unstable invariant manifold in case $b > 0$

Notice that when $u\leq\frac{1}{65}$, $\psi_{K}(u)\equiv0$, then system (3.1) becomes

Now the set $\{(u, v, n):n = 0\}$ is an invariant set. Assume that $\pi_{c, b}^{u}(\xi) = (u_{c, b}^{u}(\xi), v_{c, b}^{u}(\xi), n_{c, b}^{u}(\xi))$ is the unstable invariant manifold of the stationary point $(0, 0, 0)$ corresponding to the given parameters $c$, $b$, pointing into the region $\{(u, v, n):u > 0, v > 0\}$. For simplicity, from now on we will omit the superscript. Since the eigenvector of the positive eigenvalue $\lambda_{1}$ is $(1, \lambda_{1}, 0)$, $\pi_{c, b}(\xi) = (u_{c, b}(\xi), v_{c, b}(\xi), n_{c, b}(\xi))$ satisfies

Clearly, $v' > 0$ at least until $u = a$, since $g(u) < 0$ on $0 < u < a$ and $n\geq0$. Further $n' > 0$ when $u > \frac{1}{65}$ and $n = 0$. So $\pi_{c, b}(\xi)$ will enter the region $\{0 < n < \psi_{K}(u)\}$ at $u = \frac{1}{65}$. Combine with Lemma 3.1(a), we have $0 < n < \psi_{K}(u)$ on $\frac{1}{65} < u\leq a$. Without loss of generality, we will also assume that $u_{c, b}(0) = a, 0 < n_{c, b}(0) < \psi_{K}(a)$ and $v_{c, b}(\xi) > 0$, $v'_{c, b}(\xi) > 0$, $0\leq n_{c, b}(\xi)\leq\psi_{K}(u_{c, b}(\xi))$ when $\xi\leq0$.

Let us define the following subsets of the set $\Omega = \{(c, b): c > 0, b > 0\}$:

We will show that there is no other possible behavior of the trajectory $\pi_{c, b}(\xi)$ different from these presented above.

Lemma 3.2. The following statements are true:

${\rm{(a) }}\ \Omega = \Omega_{1}\cup\Omega_{2}\cup\Omega_{3}$.

${\rm{(b)}}\ \Omega_{2}$ and $\Omega_{3}$ are relatively open in $\Omega$.

Proof. According to the Corollary 1, we have

Notice that $0\leq n < 1$. Assume that the orbit $\pi(\xi)$ is unbounded, then $\pi(\xi)$ will enter either $\tilde{Q}^{+}$ or $\tilde{Q}^{-}$ by the relationship $u' = v$. Thus the statement (a) is true. The statement (b) appears from the fact that the sets $\tilde{Q}^{+}$ and $\tilde{Q}^{-}$ are open and positively invariance, while the solutions of system (3.1) continuously depend on the parameters.

Proposition 4. There exists values $c^{*} > 0$ and $b^{*} > 0$ such that:

Proof. Without loss of generality, we will assume that $u(0) = a$, $0 < n(0) < \psi_{K}(a)$ and $v(\xi) > 0$, $v'(\xi) > 0$, $0\leq n(\xi)\leq\psi_{K}(u(\xi))$ when $\xi\leq0$. As long as $u(\xi)$ is an increasing function, $v$ and $n$ can be expressed as functions of $u$:

Differentiating and using (3.1) gives

Let $g_{max} = \text{sup}_{u\in[0, 1]}g(u) = \text{sup}_{u\in[a, 1]}g(u) > 0$. If $V > g_{max}/c$ and $u > -\frac{1}{26}$, then

Thus, suppose that $V[(u(0)] = V(a) > g_{max}/c$, then $V(u)$ is growing as $u > a$, and $\pi_{c, b}(\xi)$ attains the set $Q^{+}$. According to the Corollary 1, we have $\lim_{\xi\rightarrow +\infty}(u_{c, b}(\xi), v_{c, b}(\xi)) = (+\infty, +\infty)$. Now it is necessary to find the conditions assuring that the inequality $V(a) > g_{max}/c$.

First, let's estimate $b^{*}$ by contradiction. Assume $V(a)\leq g_{max}/c$, then the inequality

holds on the interval $u\in(\frac{1}{65}, a)$. $\psi_{K}(u)$ is identical to zero on the interval $[0, \frac{1}{65}]$ and is a strictly convex function on the interval $(\frac{1}{65}, a)$. Take $u_{0}\in(\frac{1}{65}, a)$, then

From this we get the inequality

Applying the substitution $N(u) = W(u)e^{-\rho u}$, $\rho = \frac{b}{g_{max}} > 0$, we get the inequality

which, after the integration w.r.t. $u$ on the interval $(u_{0}, a)$, takes the form

From this we get the inequality

where $0 < \mu < 1$. So if $b > b^{*}: = \frac{2\psi_{K}(a)u_{0}g_{max}}{\mu\psi_{K}(u_{0})(a-u_{0})^2}$, then $N(a) > \psi_{K}(a)$, which contradicts Lemma 3.1(a).

Next, let us estimate $c^{*}$. The first equation of system (3.4) can be rewritten in the following form

from which appears the inequality

Suppose that $c > g_{max}/H(a)$, then $cV(a) > g_{max}$. To complete the proof, we just need to show that the inequality $c > g_{max}/H(a)$. So take $c^{*} = \frac{g_{max}}{\sqrt{-2\int_{0}^{a}g(u)du}}$, the statement is completely proved.

The analysis made in Section 3.2 about the case $b = 0$ and the continuous dependence of the solution of system (3.1) on parameters show that $\Omega_{3}$ contains an open connected set which takes $I = \left\{(c, b):0 < c < \bar{c}, b = 0\right\}$ as the boundary. It is seen from the geometry of the open sets $\Omega_{2}$ and $\Omega_{3}$ adjacent to the horizontal axis (see Figure 4 (a)), that there should exist one subsets of the set $\Omega_{1}$ lying between them.

Let $h(u): = g(u)/35(u+\frac{1}{26})$. Next define the "$L$-shaped" region consisting of the unoin of two orthogonal half-planes as follows:

where $u_{min}$ is the point in $(0, 1)$ where $h(u)$ has a local minimum. Also, we introduce the following additional subset of the $(c, b)$-plane:

From the definition of $\Lambda$ and $\Omega_{2}$ is open, we can obtain the following conclusion.

Proposition 5. $\Omega_{2}\cap\bar{\Lambda} = \emptyset$, where $\bar{\Lambda}$ is the closure of $\Lambda$ in $\Omega$.

Now we are going to prove the following:

Proposition 6. $\Omega_{3}\cap\partial\Lambda = \emptyset$, where $\partial\Lambda$ is the boundary of $\Lambda$ in $\Omega$.

Proof. The proof of this theorem is based on ideas underlying the proof of the Lemma 9 of the paper ^[17], so we try to adhere to some notations that were adopted in this paper.

Let us suppose the opposite, namely, that there exists a point $(c, b)$ belongs to the set $\Omega_{3}\cap\partial\Lambda$, corresponding the orbit $\pi_{c, b}(\xi)$ with the following properties:

● The orbit $\pi_{b, c}(\xi)$ starts from the origin and points to the first octant as $\xi\ll-1$. Without loss of generality, we can assume that $u'(\xi)$ and $v'(\xi)$ are positive on the interval $(-\infty, 0)$ and $u(0) = a$;

● At some time of the argument, says $\xi = s_{-} > 0$, the orbit crosses the set $\Sigma$ for the first time, intersecting it at a point belonging to the plane $u > u_{min}$, $v = 0$, and next, at $\xi = \check{s} > s_{-}$, intersects $\Sigma$ for the second time at a point belonging to the plane $u = u_{min}$, $v < 0$;

● Before crossing the plane $u = 0$ and going to infinity (suppose that such intersection take place at $\xi = s_{1}$), the orbit must touch the set $\Sigma$, at $\xi = s_{0} > \check{s}$ (otherwise it does not belong to the set $\partial\Lambda$). Analysis of the first equation of system (3.1) tells us that the touch point must be located at the intersection of planes $\{u = u_{min}, v\leq0\}$ and $\{u\geq u_{min}, v = 0\}$.

Projection of an orbit $\pi_{c, b}(\xi)$ on the plane $(u, v)$ is schematically represented in Figure 4(b). It is obviously that, since the orbit $\pi_{b, c}(\xi)$ is tangent to the set $\Sigma$, then there exists a number $\delta$ such that $v(\xi) > 0$ at $(s_{0}-\delta, s_{0})$, $v(\xi) < 0$ at $(s_{0}, s_{0}+\delta)$ and $v(s_{0}) = 0$. Looking at the second equation of system (3.1), we easily conclude that

which contradicts the non-negativity of $n$.

3.5.   Solitary wave solution for CRS neuron model

In this section, we establish the existence of solitary wave solutions for system (1.1). That is, the main result of this paper:

Theorem 3.1. For sufficiently small $b > 0$, there exist a positive number $c(b)$ such that Eq (2.3) has a solitary wave solution in the sense that the equivalent solitary wave Eq (3.1) has a homoclinic orbit. Moreover,

Proof. From Propositions 5 and 6, it follows that the point $(c, b)\in\Omega$ which belongs to $\partial\Lambda$ must be an element of the set $\Omega_{1}$, and $u_{c, b}(\xi)\leq u_{min} < a$ for all large $\xi$. Assume $(c, b)\in\Omega\cap\partial\Lambda$, the solitary wave Eq (3.1) has a homoclinic orbit by Proposition 3. Notice that $\Lambda\cap\Omega_{3}$ is a nonempty open connected subset from Proposition 2 and the continuous dependence of the solution of system (3.1) on parameters. For sufficiently small positive $b$ the half-line $L_{b}$ intersects $\Omega\cap\partial\Lambda$ at least once. Thus the above theorem is proved.

Remark 1. For sufficiently small $b > 0$, the homoclinic orbit above falls in the region $\{(u, v, n):u\geq0\}$ from the definition of $\Lambda$ and the proof of Propositions 3.

Conflict of interest

The authors declare that there are no conflicts of interest.