A numerical method using Legendre polynomials for solving two-point interface problems

1.   Introduction

Molecular communication, which is dramatically different from conventional communication, is attracting the attention of a number of researchers these days. Numerous interdisciplinary discussions are taking place among researchers in a variety of fields such as biology, mathematics, statistics, information theory, and so forth. Researchers in molecular communication are also closely interacting with those in another active research area, namely, nanonetworks and body area networks.

Molecular communication is also expected to possess various applications for engineering fields, such as bio-, medical, industrial and communications engineering [21]. Among them, the most important application in the medical field is drug delivery and intercellular therapy [21]. In those processes, a bio-nanomachine detects some target area, such as tumor cells in the human body, and releases drug molecules in an appropriate location, with a suitable volume and timing. Nowadays, many theoretical frameworks concerning the method of target detection in molecular communication networks are being proposed, and implementation and experiments are taking place in laboratories.

One of the key factors for carrying signals and materials in molecular communication is molecule diffusion in the medium. The early research in this direction was carried out by Einolghozati et al. [4][5] and is widely applied in the arguments for the coding and encoding methods in molecular communication networks these days. Another key factor in molecular communication is that it utilizes cooperation between bio-nanomachines to achieve a certain objective. For instance, Nakano et al. [22] recently proposed a model of target detection that imposes two different functionals on bio-nanomachines: leader and follower. The leader nanomachines search for a target in the human body and release attractant molecules upon detecting it. Follower nanomachines move according to the concentration gradient of the attractant toward the source of it, and then release drug molecules.

Nakano's model utilizes the chemotaxis of bio-nanomachines. Investigations in this direction can also be found in several papers [24][25], and these activated the modeling of the behavior of bio-nanomachines based on the Keller-Segel model [17]. The discussions by Nakano et al. [22] were limited to those of agent-based simulations and experiments. Recently, Iwasaki, Yang and Nakano [15] proposed a theoretical study based on arguments by Okaie and other authors [24][25]. In these models, both the attractant and repellent exist to effectively deliver drug molecules. The attractant plays the role of absorbing bio-nanomachines closer to the target cells, while the repellent plays the role of diffusing bio-nanomachines over an area and enabling them to search for the target over a broader area. Iwasaki, Yang and Nakano focused on the temporal behaviors of the concentration of bio-nanomachines, attractant and repellent.

The proposed model was a coupling of a one-dimensional reaction-diffusion-type partial differential equation and two ordinary differential equations.

It is also to be noted that the proposed model in that paper [15] was based on the variant model of the Keller-Segel equation, but the diffusion terms in attractant and repellent equations are neglected. After showing the existence of a positive stationary solution, Iwasaki numerically verified its stability. However, few theoretical analyses concerning the proposed model were conducted in that research. Besides, their models were limited to only the one-dimensional case.

In this paper, we extend their formulation to the two-dimensional. Although higher dimensions (such as the third dimension) can be discussed in a similar framework, two dimensions seem appropriate since bio-nanomachines move along the tissue surface. We rigorously discuss the well-posedness of the model. More concretely, we discuss the following issues:

(ⅰ) the unique existence of a positive stationary solution,

(ⅱ) the global-in-time solvability of the non-stationary problem under the smallness of data,

(ⅲ) and the non-negativeness of the non-stationary solution.

Our argument in this paper directly applies to the one-dimensional case. The remainder of this paper is organized as follows: In the next section, we give an overview of the existing mathematical arguments. In Section 3, we formulate the problem and, in Section 4, we introduce the notations used throughout this paper. The main results of this paper are stated in Section 5, followed by their proofs in Sections 6 and 7. The final section briefly concludes this paper.

2.   Existing arguments

There exist a number of works concerning molecular communication networks these days, and the following arguments are limited to those concerning the theoretical modeling of target detection in the human body. Recently, Nakano et al. [22] proposed a mathematical model that describes the temporal behavior of attractant molecules. In it, they proposed that two functionals are imposed on bio-nanomachines, which they call the leaders and followers. They showed the effectiveness of the proposed method through numerical simulations. They also clarified the situations in which their proposed model works better than the conventional method, but few mathematically rigorous discussions have been obtained.

Following these works, Iwasaki, Yang and Nakano [15] proposed a mathematical model that describes non-diffusion-based mobile molecular communication networks. They focused only on the temporal behavior of the concentration of the attractant, repellent, and bio-nanomachines under the assumption that the concentration of the target is time invariant. A similar model was discussed in a previous paper [14], including the existence and uniqueness of the solution, and the stability of the stationary solution. It reads

where and hereafter $I\equiv(0, 1)$ and ${\bf R}_+ \equiv (0, \infty)$. It was proposed as a simple version of a model proposed by Okaie et al. [24][25] that was similar to the Keller-Segel model [17] but with a variant target concentration:

Iwasaki [14] proved the global-in-time solvability of (2.1), and argued the stability of the stationary solution by constructing the Lyapunov function. However, he did not discuss the convergence rate as time tends to infinity. We also point out that, although the model in that research admits constructing a global-in-time solution without the smallness of the initial data, the method does not apply to the model studied in this paper. For other arguments concerning the model by Okaie et al. [24][25], see the review by Iwasaki [14] and the references therein.

The model discussed in this paper is a coupling of the reaction-diffusion equation and ordinary differential equations. Recently, Marciniak-Czochra et al. [31] studied the nonstability of such systems under certain conditions. That seems meaningful since the reaction-diffusion equation reflects the denovo patterns or the Turing instability.

On the other hand, Iwasaki, Yang and Nakano [15] numerically indicated the stability of a stationary solution of their model. It seems correct under some assumptions through the analysis of the corresponding eigenvalue problem.

Since the models presented so far arise from the Keller-Segel model, we will give a brief overview of the mathematical arguments concerning the Keller-Segel equations. A huge number of contributions concerning the mathematical arguments of Keller-Segel equations and its variations exists, and therefore, we limit ourselves to the arguments that closely relate ours.

Schaaf [27] studied the stationary solution to the Keller-Segel equation under the general non-linearity, and reduced the problem to a scalar equation by using the bifurcation technique. She also provided a criterion for bifurcation of the solution. Osaki and Yagi [26] provided a global-in-time solution of the classical one-dimensional Keller–Segel equation.

Kang [16] investigated the existence and stability of a spike solution in the asymptotic limit of a large mass. They also studied the global-in-time existence of a solution to a reduced version of the Keller–Segel equation. The latter part is conducted by using the energy method.

Thorough surveys are provided by Horstmann [13] and the references therein.

3.   Formulation

In this section, we formulate the problem to be discussed in this paper. From Iwasaki, Yang and Nakano [15], the temporal behavior of the concentrations of bio-nanomachines, attractant, and repellent in one-dimensional space, denoted as $C_b(x, t)$, $C_a(x, t)$ and $C_r(x, t)$, respectively, are represented as follows.

with boundary and initial conditions

For the sake of simplicity, we introduced the notation $\mathit{\boldsymbol{u}} = \bigl( C_b, C_a, C_r \bigr)^{\rm T}$ above. Here, $t$ is time; $x$ is the location of the materials on the tissue surface; $V_a$ and $V_r$ are positive constants denoting attractant and repulsive coefficients, respectively; $a_1(x) = g_ah(x)/(h(x)+K_a)$, $a_2(x) = g_rK_r/(h(x)+K_r)$, with $h(x)$ being the target concentration; $\nabla$ is the two-dimensional gradient; $\mathit{\boldsymbol{n}}$, is the outer unit normal to $\Omega$; and $K_a$ (or $K_r$) is the positive constant standing for the target concentration leading to the half maximum attractant (or repellent) production rate, respectively. The notation $g_a$ (or $g_r$) is also a positive constant representing the maximum attractant (or repellent) production rate, respectively; and the positive constant $k_a$(or $k_r$) is the decay rate of the attractant (or repellent), respectively.

Note that in this type of formulation, $C_b$ is a probability density so that $C_b \geq 0$ and $\int_\Omega C_b \; {\rm d}x = 1$ must always be satisfied, while $C_a$ and $C_r$ are the concentrations, and satisfy only $C_\tau \geq 0 \; (\tau = a, r)$. For instance, we can take $C_r\bigr|_{t = 0} = 0$ as was the case in [24].

With the term of classical Keller–Segel equation, this corresponds to the case when the sensitivity function $\chi(u) = u$. This type of model is used to formulate the angiogenesis model (see, e.g., [8]). Friedman and Tello [8] studied a system of equations close to ours:

There, the authors imposed the assumption $g(p, w) = \phi(p, w) h(p, w)$, which satisfies $\phi >0$ and $h(p_1, w_1) = h(p_2, w_2)$ with some constants $(p_i, w_i) \; (i = 1, 2)$ satisfying $0 \leq p_1 < p_2$ and $w_1 < w_2$. Together with some other assumptions, they showed the unique existence of a global-in-time solution in $C^{2+\beta, \frac{2+\beta}{2}}(\Omega_\infty)$. They also showed the asymptotic stability of the solution when the initial data is sufficiently close to the stationary solution. However, their assumptions are not satisfied in general in our case, except for a special case (which will be stated later).

The arguments in this direction are expounded on by Guarguaglini and Natalini [9] [10] with a more general setting:

With some assumptions on $F$, $G$, $\phi$ and $\psi$, Guarguaglini and Natalini [9] obtained the global-in-time weak solution to this system with small data and the spatial dimension of $2$ or more. They also discussed the case when the first equation contains the gradient of the first unknown variable:

By the change of variable, this case is reduced to (3.3):

with $\mu\phi^\prime(c) = \phi(c)\chi(c)$. Guarguaglini and Natalini [9] considered an example, which is close to ours:

Here, if

holds with some other assumptions, the a-priori estimate is obtained on which the existence of a weak global-in-time solution follows. Since $\phi(c) = \phi_0{\rm e}^{c/\mu}$ with some constant $\phi_0$ in our case, the conditions above are not satisfied in our case since $\chi(c) \equiv 1$. In the paper by Iwasaki, Yang and Nakano [15], the authors provided the existence of a stationary solution ${\mathit{\boldsymbol{\bar u}}}(x) \equiv \bigl( \overline{C}_b(x), \overline{C}_a(x), \overline{C}_r(x) \bigr)^{\rm T}$ to (3.1)-(3.2):

They also showed the stability of ${\mathit{\boldsymbol{\bar u}}}(x)$, that is, the convergence of the solution of (3.1)–(3.2) to that of (3.4) through numerical simulations under specific parameter values.

To extract the mathematical essence, we introduce the notation $\mathit{\boldsymbol{\Phi}} = D_b \nabla \overline{C}_b - \overline{C}_b \bigl( V_a \nabla \overline{C}_a - V_r \nabla \overline{C}_r \bigr)$ for now. Then, we have

We will study (3.5) later.

As for the stationary solution, we again mention the discussion by Friedman and Tello [8], in which they obtained only a constant stationary solution under some assumptions. This is a special case of ours.

Schaaf [27] first discussed the one-dimensional stationary solution to the Keller–Segel system:

They also discussed the bifurcation from the stationary solution. Therein, they showed that the system (3.6) is reduced to an equation:

with $u(x) = \phi(v(x), \lambda)$. (See also, [13].) In our case, we can reduce our problem to this form by taking $v(x) \equiv V_a \overline{C}_a - V_r \overline{C}_r$, $k_c = 0$, and $g(u, v) \equiv F(x)u-v$ with $F(x) \equiv \frac{V_aa_1(x)}{k_a} - \frac{V_ra_2(x)}{k_r}$, and the result partially gives us an insight. But we will have more general results, which will be stated in Section 4. For the equation similar to (3.5), there have been some arguments in the past. Struwe and Tarantello [29] discussed the problem

where $\Omega$ is a two-dimensional torus. By using variational techniques and mountain-pass solution, they showed the existence of a nontrivial solution to (3.8) for $\lambda \in (8\pi, 4\pi^2).$ The studies by Wang and Wei [33] and Senba and Suzuki [28] extended these results to more general cases. They studied the following problem:

As in Schaaf [27], they translated (3.9) by using $u = \sigma {\rm e}^v$ into the form:

However, the studies in this direction concern the stationary solution to the elliptic problem, and are not applicable to our case.

Consider the coupled system of chemotaxis equation and ODE in the Keller–Segel literature. Corrias et al. [3] thoroughly studied the system

They found the global-in-time solutions to (3.10) for $d \geq 2$ with very general settings on $\chi(c)$. They showed that for $d \geq 3$, there exists a global-in-time solution under the smallness of initial data, whereas for $d = 2$, without smallness of initial data. Their method is based on changing the form of unknown variables, and then the energy estimate of the $p$-norm of new variables. Although it seems simple and widely applicable, it does not apply to our case. The reason is in our case, $C_a$ is represented as:

which will be shown later as (7.3). Similar representation holds for $C_r$ also. Therefore, we cannot assume the essential boundedness of them as in case of [3]. In their estimate, they made use of the boundedness of $\sup_{0 \leq c \leq \|c_0\|_\infty} \phi^2(c) \chi(c) c^m$, which no longer holds in our case. Sugiyama et al. [30] studied the system

They obtained the global-in-time weak solution to (3.11) in Besov spaces under the smallness of initial data for $n \geq 3$, and local-in-time solution for $n \geq 1$. Ahn and Kang [1] studied the system

They classified the approaches to (3.12) according to the range of $\lambda$, among which the case of $\lambda = 1$ is close to ours. Indeed, in that case, they transformed the problem into the form:

in case $\Omega = {\bf R}^d$, with decay at infinity of unknown variables. They implied that the blow-up of the solution to (3.13) in the Sobolev space occurs. Fontelos et al. [7] studied the following system in one-dimensional case.

They transformed (3.14) into a simpler one:

After further change of variables, they considered

Then, they adopted the change of variables $c = {\rm e}^{-V}$ and $u = P{\rm e}^V$. However, this approach is not applicable to our case, for in our case, $k_a \ne k_r$ in general. Even if $k_a = k_r$ holds, since $C_a$ and $C_r$ take the form of (7.3) again, which prevents us from assuming the boundedness of $c$.

The previous paper of the author [12] provides the local-in-time solvability of (3.1)–(3.2) in one-dimensional space but here we will go further. In this paper, we study the well-posedness of (3.1)–(3.2). Since we wish to consider the global-in-time solvability around the stationary solution, we first subtract ${\mathit{\boldsymbol{\bar u}}}(x)$ from the original problem, and consider the problem concerning

It reads

where ${\mathit{\boldsymbol{\tilde u}}}_0 = \bigl( \widetilde{C}_{b0}, \widetilde{C}_{a0}, \widetilde{C}_{r0} \bigr)^{\rm T}$.

4.   Notations

In the following, let ${\mathcal G}$ be an arbitrary open set in ${\bf R}^2$, and let $T>0$, and ${\mathcal G}_T \equiv {\mathcal G} \times (0, T).$ For a set $S$ in general, $\overline{S}$ denotes its closure and $\partial S$ its boundary.

Hereafter, $C^l({\mathcal G}) \; (l \in {\bf N}$) stands for a set of functions defined on ${\mathcal G} \subset {\bf R}^2$ that have $l$-th order continuous derivatives.

By $C^{r+\alpha}( {\mathcal G} )$ with a non-negative integer $r$ and $\alpha \in (0, 1)$, we mean the Banach space of functions from $C^r(\overline{ \mathcal G})$, whose $r$th derivatives satisfy the Hölder condition with exponent $\alpha$, i.e. the space of functions with the finite norm

where

By the notation $C^{l, 1}({\mathcal G}) \; (l \in {\bf N}$), we denote a set of functions defined on ${\mathcal G}$ that have $l$-th order Lipschitz continuous derivatives. We also use similar notations to represent the regularity of domain boundaries. $L_2({\mathcal G})$ means a set of square-integrable functions defined on ${\mathcal G}$, equipped with the norm

The inner product in $L_2({\mathcal G})$ is defined by

where $\overline{z}$ stands for the complex conjugate of $z \in {\bf C}$.

Analogously, we define

By $\|\cdot\|_{L_p({\mathcal G})}$, we denote the usual $L_p$-norm with $1<p\leq +\infty$ on ${\mathcal G}$:

For simplicity, we hereafter denote the $L_p$-norm of a function $f$ over $\Omega$ ($\Omega$ is the region where we consider the problem) merely by $|f|_p$. Especially, we denote the $L_2$-norm of a function $f$ over $\Omega$ merely as $|f|$ when it is obvious.

By $W_2^r({\mathcal G}) \; (r>0)$, we mean a space of functions $f(x), \; x\in {\mathcal G}$ equipped with the norm $\|f\|_{W_2^r({\mathcal G})}^2 \equiv \mathop{\sum}_{|\alpha|<r} \|D^{\alpha}f \|_{L_2({\mathcal G})}^2 + \|f\|_{\dot{W}_2^r({\mathcal G})}^2$, where

Next, for arbitrary $T>0$, we introduce anisotropic Sobolev–Slobodetskiĭ spaces [34]:

whose norms are defined by

The set of functions with vanishing initial data, $\overset{\circ}{W}_2^{r, \frac{r}{2}}({\mathcal G}_T)$, is defined as [18]:

We also introduce

We also define a function space

Its norm is defined as

For simplicity, we shall use the following notations later.

The norms of these spaces are denoted as $\|\cdot\|_{{W}^{(l)}}$ and so forth.

For a Banach space ${\mathcal B}$ with the norm $\|\cdot\|_{\mathcal B}$, we denote the space of ${\mathcal B}$-valued measurable functions $f(t)$ on the interval $(a, b)$ by $L_p(a, b;{\mathcal B})$, whose norm is defined by

The norms of vector and product spaces are defined in the usual manner.

5.   Main results

In this section, we state the main results of this paper. Detailed proofs are provided in Sections 6 and 7. Before discussing the solvability of the non-stationary problem, we first argue the solvability of the stationary problem.

Theorem 5.1. Assume that $\Omega$ is bounded, $l \in (1/2, 1)$, $\Gamma = \partial \Omega \in C^{2, 1}$. In addition, we consider either of the two cases below.

(i) $F(x) \equiv V_aa_1(x)/k_a - V_ra_2(x)/k_r$ is a constant;

(ii) $F(x) > 0$ and $\nabla^2 F <0$ on $\overline{\Omega}$, $\mathit{\boldsymbol{n}} \cdot \nabla F >0$ on $\Gamma$, $a_j \in C^{3+l}(\overline{\Omega}) \; (j = 1, 2)$, and $a_j>0 \; (j = 1, 2)$.

Then, there exists the unique solution ${\mathit{\boldsymbol{\bar u}}} = \Bigl( \overline{C}_b, \overline{C}_a, \overline{C}_r \Bigr) \in C^{3+l}(\Omega)$ to $\rm(3.4)$. It also satisfies $\overline{C}_\tau(x, t) > 0 \; (\tau = b, a, r)$ on $\overline{\Omega}$.

Proof. The first case falls into the case of Friedman and Tello. [8], and their result guarantees that there exists only a constant solution to (3.4).

Our argument basically follows the one by Iwasaki [14] applied to the 1- dimensional case. There is little modification due to the difference of dimension. In his argument, however, he derived the linear differential equation of $\log \overline{C}_b$ from the beginning. Since it is not clear whether $\overline{C}_b>0$ holds or not in advance, we need to somewhat modify the proof.

Recall (3.4) is equivalent to (3.5), and further introduce $\mathit{\boldsymbol{\Phi}} = \nabla \psi$ there. Then, we have

and, consequently,

It is clear that $\psi = c_0$ with a constant $c_0 \in {\bf R}$ by virtue of the maximum principle, and thus $\mathit{\boldsymbol{\Phi}} = {\bf{0}}.$

Now, since $\overline{C}_b$ is subjected to a linear ordinary differential equation, $\overline{C}_b>0$ holds everywhere in $\overline{\Omega}$, and thus on $\overline{\Omega}$, if there exists a point $x_0 \in \overline{\Omega}$ on which $\overline{C}_b(x_0) >0$ holds.

Such a point does exist by virtue of Lemma 5.2 below. Thus, if there exists a non-zero solution to (3.4), it has to satisfy $\overline{C}_b>0$.

This and the fact $\mathit{\boldsymbol{\Phi}} = {\bf{0}}$ lead to

Thus, following the arguments by Iwasaki [14], we arrive at the unique existence of a solution to (3.4), and the positiveness of $\overline{C}_b$. The positiveness of $\overline{C}_\tau \; (\tau = a, r)$ follows directly from (3.4)$_2$ and (3.4)$_3$.

Lemma 5.2. Let us assume $l \in (1/2, 1)$, $F(x)>0$ and $\nabla^2 F < 0$ on $\overline{\Omega}$, and $\mathit{\boldsymbol{n}} \cdot \nabla F >0$ on $\Gamma$. Then, if there exists the non-zero solution ${\mathit{\boldsymbol{\bar u}}} = \Bigl( \overline{C}_b, \overline{C}_a, \overline{C}_r \Bigr) \in C^{3+l}(\Omega)$ to $\rm(3.4)$, it satisfies $\overline{C}_b \geq 0$ everywhere in $\overline{\Omega}$.

Proof. Substitute (3.4)$_2$ and (3.4)$_3$ into (3.4)$_1$, and we obtain

We first assume that there exists a point $x_1 \in \Omega$ such that

holds.

Then, it satisfies $\nabla \overline{C}_b (x_1) = 0$ and $\nabla^2 \overline{C}_b (x_1)>0$. Otherwise, there exists a neighborhood $U_1 \subset \Omega$ of $x_1$ and $x_1^\prime \in U_1$, such that the same situation holds at $x_1^\prime$. This comes from (3.4)$_5$.

Under the assumption of the lemma, the first issue and (5.3) yield

a contradiction.

Next, assume that there exists a point $x_2 \in \Gamma$, such that

holds. We note that the case $\mathit{\boldsymbol{n}} \cdot \nabla \overline{C}_b(x_2) = 0$ has been discussed in the argument above.

This time $\mathit{\boldsymbol{n}} \cdot \nabla \overline{C}_b <0$ on $x = x_2$. The boundary condition (3.4)$_2$ reads

From the assumption, this yields $\mathit{\boldsymbol{n}} \cdot \nabla \overline{C}_b >0$ on $x = x_2$, a contradiction. Thus, we can conclude that $\overline{C}_b(x) \geq 0$ on $\overline{\Omega}$.

Next, we state the existence theorem of the solution to (3.17).

Theorem 5.3. In addition to the assumptions in ${\rm Theorem}~~~ 5.1$, let $l \in (1/2, 1)$. We also assume:

(i) ${\mathit{\boldsymbol{\tilde u}}}_0 \in W_2^{2+l}(\Omega)$, $\bigl\| {\mathit{\boldsymbol{\tilde u}}}_0 \bigr\|_{W_2^{2+l}(\Omega)} < \delta_0$ with a sufficiently small $\delta_0>0$,

(ii) $h(x) \in C^{3+l}(\overline{\Omega})$, $C_{b0}(x) \geq 0 \; \forall x \in \overline{\Omega}$, $\int_\Omega \widetilde{C}_{b0}(x) \; {\rm d}x = 0$,

(iii) $C_{a0}(x), \; C_{r0}(x) \geq 0 \; \forall x \in \overline{\Omega}$,

(iv) the following inequality is satisfied:

where $c_\Omega$ is a constant dependent on the size of $\Omega$,

(v) and the following inequality is satisfied with a certain $c_{52}>0$:

In addition, let the compatibility condition up to order $1$ [18] be satisfied.

Then, problem $\rm(3.17)$ has a unique solution,

that satisfies $\widetilde{C}_b \in \widetilde{W}_2^{3+l, \frac{3+l}{2}}(\Omega_\infty)$ also.

In order to prove ${\rm Theorem}~~ 5.3$, we first consider the linear problem. Next, we make use of the multiplicative inequalities and fixed point theorem to verify the existence of the global-in-time solution. These discussions are held in the next section. We also show

Theorem 5.4. Let $T>0$ be an arbitrary number. Under the assumptions in ${\rm Theorem}~~ 5.3$, the solution of $\rm(3.1)$–$\rm(3.2)$, if it exists, satisfies $C_\tau(x, t) \geq 0 \; (\tau = b, a, r)$ for $(x, t) \in \Omega_T.$

The proof of Theorem 5.4 is provided in Section 7.

6.   Proof of solvability

In this section, we discuss the solvability of (3.17). In the following, we use general positive constant $c$.

6.1.   Linear problem

In this subsection, we first prove the solvability of the linear problem associated with (3.17), and then Theorem 5.3. We first consider the following linearized problem with vanishing initial data.

We have

Theorem 6.1. Let $\Omega$ be bounded, $\Gamma \in C^{2, 1}$, $l \in (1/2, 1)$, and assume conditions (ⅰ)–(ⅴ) in ${\rm Theorem}~~ 5.3$. We also assume $F_1 \in \overset{\circ}{W}_2^{1+l, \frac{1+l}{2}}(\Omega_{\infty})$, $F_2 \in \overset{\circ}{W}_2^{\frac32+l, \frac34+\frac{l}{2}}(\Gamma_\infty)$, and

In addition, let the compatibility conditions up to order $1$ be satisfied. That is,

Then, there exists a unique solution ${\mathit{\boldsymbol{\tilde u}}} = (\widetilde{C}_b, \widetilde{C}_a, \widetilde{C}_r)^{\rm T} \in\overset{\circ}{W}_2^{3+l, \frac{3+l}{2}}(\Omega_{\infty})$ to $\rm(6.1)$ satisfying $\widetilde{C}_b \in \widetilde{W}_2^{3+l, \frac{3+l}{2}}(\Omega_{\infty})$ and

In order to prove ${\rm Theorem}~~ 6.1$, we apply the Fourier transform with respect to $t$ to unknown functions extended into the region $t<0$ by zero [2]:

for a function $f$ in general. After applying this transform to (6.1), we further substitute $\tau = -{\rm i}\lambda$ to obtain

We substitute (6.4)$_2$ and (6.4)$_3$ into (6.4)$_1$, and then reduce the problem to one concerning $\hat{C}_b(x, -{\rm i}\lambda)$.

where

Note that $\hat{F}_1(x, -{\rm i}\lambda) = \int_{\bf R} {\rm e}^{-\lambda t} F_1(x, t) \; {\rm d}t$ is the Fourier transform of ${\rm e}^{-\sigma_0 t} F_1$ with respect to $\sigma_1$. Therefore, from the assumption $F_1 \in W_2^{1+l, \frac{1+l}{2}}(\Omega_\infty)$,

is finite for $\sigma_0>0$. Thus, $\hat{F}_1(x, -{\rm i}\lambda)$ has an analytic extension onto the region ${\rm Re} \; \lambda>0$. Similar facts hold for other functions in (6.5).

Below, we write $\lambda = \sigma_0 + {\rm i}\sigma_1$, and $D_{(+)} \equiv \bigl\{ z \in {\bf C} \bigr| {\rm Re} \; z > 0 \bigr\}$. In order to define the term generalized solution to (6.5), we first introduce the following operator.

for $\forall \eta(x) \in W_2^1(\Omega).$

Definition 6.2. A function $u(x) \in {W}_2^1(\Omega)$ is called a generalized solution to (6.5) if it satisfies

First, we state

Lemma 6.3. Under the assumptions of ${\rm Theorem}~~ 6.1$, we have

Proof. Since the assumption (6.2) in Theorem 6.1 implies that $\int_\Omega \widetilde{C}_b(x, t) \; {\rm d}x = 0 \; \forall t>0$ if it exists in $W_2^1(\Omega)$, we seek for the solution satisfying it. Observe that

thanks to the Poincaré inequality. Next, note that when $\sigma_0 \geq 0$,

Elementary calculations make us to obtain (6.8), so we omit the detail here. Now, (6.8) yields

Then, the assumption (ⅳ) in Theorem 5.3 implies

This is sufficient for the following estimate to hold:

which is equivalent to (6.7).

Next, we state the following lemma.

Lemma 6.4. Under the assumptions of ${\rm Theorem}~~ 6.1$, the problem $\rm(6.5)$ has no more than one generalized solution in $\widetilde{W}_2^1(\Omega)$.

Proof. From the definition of $L[u, \eta]$ and Lemma 6.3, there exists a constant $c_{62}$ and

where we have applied the trace and Sobolev embedding theorems. Hereafter, $\varepsilon$ is an arbitrary small positive number and $C_\varepsilon$, a positive constant non-decreasingly dependent on $\varepsilon$. Note that the Poincaré inequality holds for $u \in \widetilde{W}_2^1(\Omega)$. This yields

Thus, it is clear that $u$ is equivalently zero if we replace $\hat{F}_1$ and $\hat{F}_2$ with zero.

Next, we show the existence of a generalized solution to (6.5). Before that, we prepare an inner product:

As is easily seen, ${\rm Re}(G_1(x, \lambda)) \geq \exists c_{63}$ holds under the assumption (v) of Theorem 5.3. Thus, we have

By virtue of the Cauchy-Schwartz inequality, it is also seen that

Thus, $[u, u]$ is equivalent to the norm of $W_2^1(\Omega)$.

Next, we define

Then, it is easily observed that

Thus, the Riesz representation theorem enables us to represent $I_1$ in the form of the scalar product:

This operator ${\mathcal K}$ is a bounded map on $W_2^1(\Omega)$ in virtue of the following.

Similarly, since

it is represented in the form of the scalar product:

We also state the following lemma.

Lemma 6.5. The operator ${\mathcal K}$ defined above is a compact operator in $W_2^1(\Omega)$.

Proof. Let $\{v_m\}_m$ denote a sequence of elements in $W_2^1(\Omega)$, satisfying $\|v_m\|_{W_2^1(\Omega)} \leq c_{67}$. Then, the sequence $\{ {\mathcal K} v_m\}_m$ are also uniformly bounded. Thanks to the Rellich-Kondrachov theorem, $W_2^1(\Omega)$ is compactly embedded into $L_2(\Omega)$. Therefore, there are subsequences of $\{v_m\}$ and $\{ {\mathcal K} v_m\}$ that converge strongly in $L_2(\Omega)$. By virtue of the equality

it is clear that

and therefore, $\{ {\mathcal K} v_m\}$ converges strongly in $W_2^1(\Omega)$.

Now, (6.6) enables us to re-formulate (6.5) in the following form.

Since (6.12) should be satisfied for all $\eta \in W_2^1(\Omega)$, it is equivalent to the abstract equation in $W_2^1(\Omega)$:

Since ${\mathcal K}$ is a linear and compact operator, we can apply Fredholm's theorems. Especially, the first statement of it guarantees that (6.13) has a solution if the homogeneous equation

has only a trivial solution $w = 0$. But such $w$ is a solution to (6.5) with $\hat{F}_1 = \hat{F}_2 = 0$, since (6.14) is equivalent to the integral identity

which is simply the identity $L[w, \eta] = 0$. Lemma 6.4 states that $w = 0$, and therefore (6.14) certainly has a solution. Since (6.13) is equivalent to (6.12), and consequently to (6.6), this is the desired solution.

Now we are in a position to state

Proposition 6.1. Let us assume the same assumptions as in ${\rm Theorem}~~ 6.1$. Then, for each $\lambda \in D_{(+)}$, there exists a unique solution $\hat{C}_b(x, -{\rm i} \lambda) \in \widetilde{W}_2^{3+l}(\Omega)$ to $\rm(6.5)$ satisfying

Proof. Since the existence of the generalized solution in $W_2^1(\Omega)$ has been established before, we limit ourselves to the regularity discussions.

Under the assumptions $F_1 \in L_2(\Omega)$ and $F_2 \in W_2^\frac12(\Gamma)$, the generalized solution constructed above actually satisfies $\hat{C}_b \in W_2^2(\Omega)$. This is derived from the argument in [19]. You may also refer to [23].

From the relationship $L[u, \bar{u}] = l(u)$, we also have

Thus, together with (6.9), we have the desired estimate. Similarly, if we assume $F_1 \in W_2^2(\Omega)$ and $F_2 \in W_2^\frac52(\Gamma)$, we have $\hat{C}_b \in W_2^4(\Omega)$.

Finally, applying the molifier and interpolation argument [20] leads us to the desired results.

Now we prove Theorem 6.1. In virtue of (6.15), $\hat{C}_b(x, -{\rm i} \lambda)$ can be defined and is holomorphic with respect to $\lambda$ on $D_{(+)}$, since $\hat{F}_j(x, -{\rm i}\lambda) \; (j = 1, 2)$ are as well.

If we fix ${\rm Re}(\lambda) = \sigma_0>0$, $\hat{C}_b(x, -{\rm i}\lambda)$ on the line ${\rm Im}(\lambda)$ is the Fourier transform of ${\rm e}^{-\sigma_0 t}C_b(x, t)$.

Thanks to (6.15), we thus have $C_b(x, t){\rm e}^{-\sigma_0 t} \in W_2^{3+l}(\Omega)$ for each $t>0$ and $\sigma_0 > 0.$

Due to the Paley-Wiener theorem, ${\rm e}^{-\sigma_0 t} F_1 = 0$ for $t<0$. The Plancherel theorem yields

Since the right-hand side is finite due to the assumption, if we make $\sigma_0$ tend to zero, it tends to $\|F_1\|_{L_2\bigl( {\bf R}_+; W_2^{1+l}(\Omega)\bigr)}^2$. Similarly, the $L_2$-norm of $\bigl\| \lambda^\frac{1+l}{2} \hat{F}_2(\cdot, -{\rm i}\lambda) \bigr\|_{W_2^\frac12(\Gamma)}$ converges as $\sigma_0 \rightarrow 0$, whose limit is estimated from above by $\|F_2\|_{W_2^{\frac32+l, \frac34+\frac{l}{2}}(\Gamma_\infty)}$.

Thus, the right-hand side of (6.15) is finite for each $\lambda \in D_{(+)}$, and we apply the Lebesgue convergence theorem, to obtain

and ${\rm e}^{-\sigma_0 t} C_b \rightarrow C_b$ as $\sigma_0 \rightarrow 0$ in this space.

Similarly, the inverse Fourier transform of $\lambda^\frac{3+l}{2} \hat{C}_b(x, -{\rm i}\lambda)$ converges in $L_2({\bf R}; L_2(\Omega))$ as $\sigma_0 \rightarrow 0$. From these, we have

in $W_2^{3+l, \frac{3+l}{2}}(\Omega_\infty).$

The uniqueness of the solution is obvious thanks again to the coerciveness of $B_\lambda[\cdot, \cdot]$. This completes the proof of Theorem 6.1.

6.2.   Nonlinear problem

Next, we consider the nonlinear problem (3.17). Before proceeding to the detailed arguments, we prepare some lemmas. The following lemma is well known (see, for instance, [32]).

Lemma 6.6. Let $l \in (1/2, 1)$ and $m \geq 2$ and $m \geq k$. Then, for $f \in W_2^{m+l, \frac{m+l}{2}}(\Omega_\infty)$ and $g \in W_2^{k+l, \frac{k+l}{2}}(\Omega_\infty)$ in general, $fg \in W_2^{k+l, \frac{k+l}{2}}(\Omega_\infty)$ and

Our problem (3.17) is described as

where ${\mathcal A}$ is a linear operator defined on $W_2^{3+l, \frac{3+l}{2}}(\Omega_\infty)$, that associates ${\mathit{\boldsymbol{\tilde u}}}$ with $\Bigl( \frac{\partial {\mathit{\boldsymbol{\tilde u}}}}{\partial t} - {\mathcal L} {\mathit{\boldsymbol{\tilde u}}}, B {\mathit{\boldsymbol{\tilde u}}} \Bigr)^{\rm T}$, and ${\mathcal F} [ {\mathit{\boldsymbol{\tilde u}}}] = ( {\mathcal F}_1[ {\mathit{\boldsymbol{\tilde u}}}], 0, 0, {\mathcal F}_2[ {\mathit{\boldsymbol{\tilde u}}}] )^{\rm T},$ with

From Lemma 6.6, we have

for ${\mathit{\boldsymbol{\tilde u}}}$ in a bounded set in $W_2^{3+l, \frac{3+l}{2}}(\Omega_\infty)$. Similarly, it is seen that

for ${\mathit{\boldsymbol{\tilde u}}}^{(j)} \in W_2^{3+l, \frac{3+l}{2}}(\Omega_\infty) \; (j = 1, 2).$ To solve (6.16) iteratively, we first find ${\mathit{\boldsymbol{\tilde u}}}^{(0)} \in W_2^{3+l, \frac{3+l}{2}}(\Omega_\infty)$, which satisfies (3.17) at $t = 0$, and

This is achieved by a well-known method as follows (see, for instance, Theorem IV.4.3 in [18], p.298). Now, define $\mathit{\boldsymbol{\phi}}_1(x) \equiv {\mathcal L}{\mathit{\boldsymbol{\tilde u}}}\bigr|_{t = 0}$, and extend it, together with ${\mathit{\boldsymbol{\tilde u}}}_0$ onto the whole space ${\bf R}^2$ preserving the regularity [18]. Then, consider the solution $\mathit{\boldsymbol{u}}^*$ of the following problem.

Next, define ${\mathit{\boldsymbol{\tilde u}}}^{(0)}$ as a solution to

This ${\mathit{\boldsymbol{\tilde u}}}^{(0)}$ actually satisfies the desired features. Now, we rewrite the problem (6.16) for new variable ${\mathit{\boldsymbol{\tilde u}}}^{(1)} \equiv {\mathit{\boldsymbol{\tilde u}}} -{\mathit{\boldsymbol{\tilde u}}}^{(0)}$:

If ${\mathit{\boldsymbol{\tilde u}}}^{(1)} \in \overset{\circ}{W}_2^{3+l, \frac{3+l}{2}}(\Omega_\infty)$, then, using our method of constructing ${\mathit{\boldsymbol{\tilde u}}}^{(0)}$, ${\mathcal A} {\mathit{\boldsymbol{\tilde u}}}^{(0)} = {\mathcal F}[{\mathit{\boldsymbol{\tilde u}}}^{(0)}] = {\mathcal F}[{\mathit{\boldsymbol{\tilde u}}}^{(0)} +{\mathit{\boldsymbol{\tilde u}}}^{(1)} ]$ at $t = 0$. Thus, the right-hand side of (6.19) belongs to $\overset{\circ}{W}^{(1+l)}$. Let ${\mathcal A}_0$ be a solution operator of Theorem 6.1 for the linear problem with zero initial data. Then, by virtue of Theorem 6.1, if

then the above assumption is satisfied, for ${\mathcal A}_0^{-1}$ is a bounded operator from $\overset{\circ}{W}^{(1+l)}$ to $\overset{\circ}{W}_2^{3+l, \frac{3+l}{2}}(\Omega_\infty)$. In order to show the solvability of (6.20), let us define a map

and show that it has a fixed point, assuming that

with a sufficiently small $\delta_0>0$. By our choice of ${\mathit{\boldsymbol{\tilde u}}}^{(0)}$, we have

and $\| {\mathcal A} {\mathit{\boldsymbol{\tilde u}}}^{(0)} \|_{W_2^{1+l, \frac{1+l}{2}}(\Omega_\infty)} \leq c_{614} \delta_0.$ By (6.17), we have

By combining this with the boundedness of ${\mathcal A}_0^{-1}$, we have

Thus, if we take $\overline{B} \equiv \Bigl\{ {\mathit{\boldsymbol{\tilde v}}} \in\overset{\circ}{W}_2^{3+l, \frac{3+l}{2}}(\Omega_\infty) \Bigr| \|{\mathit{\boldsymbol{\tilde v}}} \|_{W_2^{3+l, \frac{3+l}{2}}(\Omega_\infty)} \leq 2c_{616} \delta_0 \Bigr\}$, ${\mathcal M}$ maps $\overline{B}$ to itself provided $\delta_0$ is sufficiently small, satisfying $(4c_{616}^2+1)\delta_0 \leq 1$. In a similar manner, thanks to (6.18), we have

Thus, if we take $\delta_0 < 1/c_{617},$ then ${\mathcal M}$ is a contraction map on $\overline{B}$, and has a unique fixed point there. Thus, the existence of a solution ${\mathit{\boldsymbol{\tilde u}}} \in W_2^{3+l, \frac{3+l}{2}}(\Omega_\infty)$ has been established.

Next, we argue the uniqueness of the solution in the similar line with Beale [2]. Let us assume that there exists a number $T_1>0$ and another solution ${\mathit{\boldsymbol{\tilde u}}}^{(2)} \in W_2^{3+l, \frac{3+l}{2}}(\Omega_{T_1})$ to (3.17), and define

If $T_2 < T_1$, we may set $T_2 = 0$ and, thereby, we assume ${\mathit{\boldsymbol{\tilde u}}}^{(2)}(t) \ne {\mathit{\boldsymbol{\tilde u}}} (t)$ near $t = 0$. Take $T_0>0$ so that $\|{\mathit{\boldsymbol{\tilde u}}}^{(2)} \|_{W_2^{3+l, \frac{3+l}{2}}(\Omega_{T_0})} < c_{618} \delta_0,$ where $c_{618}$ and $\delta_0$ are as before.

Introduce ${\mathit{\boldsymbol{\tilde w}}} = {\mathit{\boldsymbol{\tilde u}}} - {\mathit{\boldsymbol{\tilde u}}}^{(2)}$ for $0<t<T_0$. Obviously, ${\mathit{\boldsymbol{\tilde w}}} \in W_2^{3+l, \frac{3+l}{2}}(\Omega_{T_0})$, and

From the assumption, the second term of the right-hand side in (6.24) is estimated by $c_{615} \delta_0$, while the first term is also estimated as

By using the extension by reflection [20], we are able to extend ${\mathit{\boldsymbol{\tilde w}}}$ onto the region $t \in (0, +\infty)$ so that

with a certain $c_{620}>0.$

On the other hand, if we introduce $\breve{\mathit{\boldsymbol{w}}} \equiv {\mathcal A}_0^{-1} \bigl( {\mathcal F}[{\mathit{\boldsymbol{\tilde u}}}] - {\mathcal F}[{\mathit{\boldsymbol{\tilde u}}}^{(2)}] \bigr)$, we have

By virtue of (6.25)–(6.27), we then have

with a certain $c_{622}>0$. Now, from the definition, ${\mathcal A} \breve{\mathit{\boldsymbol{w}}} = {\mathcal F}[ {\mathit{\boldsymbol{\tilde u}}} ] - {\mathcal F}[ {\mathit{\boldsymbol{\tilde u}}}^{(2)} ]$, and also ${\mathcal A} {\mathit{\boldsymbol{\tilde w}}} = {\mathcal F}[ {\mathit{\boldsymbol{\tilde u}}} ] - {\mathcal F}[ {\mathit{\boldsymbol{\tilde u}}}^{(2)} ]$.

Thus, we have ${\mathcal A} \bigl( \breve{\mathit{\boldsymbol{w}}} - {\mathit{\boldsymbol{\tilde w}}} \bigr) = {\bf{0}}$ on $t \in (0, T_0)$. From the uniqueness of the solution of the linear problem, we then observe $\breve{\mathit{\boldsymbol{w}}} = {\mathit{\boldsymbol{\tilde w}}}$ on $t \in (0, T_0)$. Replacing $\breve{\mathit{\boldsymbol{w}}}$ by ${\mathit{\boldsymbol{\tilde w}}}$ in (6.28) and using (6.26), we then have

Take $\delta_0$ so small that $c_{620}c_{622} \delta_0 < 1$, and then, we have ${\mathit{\boldsymbol{\tilde w}}} = {\bf{0}}$ on $t \in (0, T_0)$. This means ${\mathit{\boldsymbol{\tilde u}}}$ and ${\mathit{\boldsymbol{\tilde u}}}^{(2)}$ coincide on $t \in (0, T_0)$, which contradicts the original assumption. Thus, we obtain the uniqueness of the solution to (3.17).

Finally, it is easily seen that $\int_\Omega \widetilde{C}_b(x, t) \; {\rm d}x = 0$, by

by virtue of the boundary condition. This completes the proof of Theorem 5.3.

7.   Non-negativeness of the solution

In this section, we argue the non- negativeness of solution ${C}_b$ to (3.1)–(3.2) to prove Theorem 5.4. We use the so-called Stampacchia truncation method (see, for instance, [6] and [11]). To do that, we divide $C_b$ into its positive and negative parts:

where $C_b^{(+)} = (|C_b| + C_b)/2 \geq 0$ and $C_b^{(-)} = (|C_b| - C_b)/2 \geq 0$. We also define

for each $t>0$. Let us multiply (3.1)$_1$ by $C_b^{(-)}$, and integrate with respect to $x$ over $\Omega_{(-)}(t)$. We show

Indeed, by Green's theorem and the boundary condition (3.2)$_1$, we observe

and therefore, we have

Noting that $C_b^{(-)} (x, t) = 0$ on $\Omega \backslash \Omega _{(-)}(t)$, we arrive at (7.1). Next, taking into mind $C_b^{(+)}(x, t) C_b^{(-)}(x, t) = 0$ on each $(x, t) \in \Omega_\infty$, we note

Thereby, we have

Finally, we further modify the first term of the right-hand side in (7.1). Note that since

we have

From these, we have

Then, by virtue of the Young's inequality, the right-hand side is estimated as

By noting that $W_2^{1+l}(\Omega) \subset L_\infty(\Omega)$, and taking this into mind, we modify (7.2) and apply the Gronwall inequality. We then have an estimate of the form

By noting that $C_b^{(-)} \bigr|_{t = 0} = 0$, we arrive at $C_b^{(-)} (x, t) = 0$, i.e., $C_b(x, t) \geq 0$ for each $x \in\Omega$, $t>0$. By virtue of the Sobolev embedding theorem, this holds in the pointwise sense. Finally, from (3.1)$_2$, we have

This and assumption (ⅲ) of ${\rm Theorem}~~ 5.3$ yield $C_a(x, t) \geq 0$ on $\Omega_\infty$. The fact $C_r(x, t) \geq 0$ is derived in a similar manner. This completes the proof of Theorem 5.4.

8.   Conclusion

In this paper, we provided the global-in-time solvability of the two-dimensional non-stationary problem of a target detection model in a molecular communication network in Sobolev–Slobodetskiĭ $\quad$ space. We also showed the non-negativeness of the non-stationary solution. We will tackle the stability arguments in the near future.