Stochastic simulations of the Schnakenberg model with spatial inhomogeneities using reactive multiparticle collision dynamics

1.   Introduction

1.1.   Background

This paper focuses on the boundary control of large amplitude classical solutions to an initial-boundary value problem of the system of hyperbolic balance laws:

which originates from the chemotaxis model with logarithmic sensitivity:

Chemotaxis is the biochemical process through which the movement of an organism or entity is not only regulated by random diffusion, but also controlled by the concentration gradient of a chemical stimulus in the local environment. Biological processes undergoing chemotaxis include bacterial foraging, immune response, embryonic development and tissue homeostasis ^[1,2], blood vessel formation ^[3], fish pigmentation patterning ^[4], tumor angiogenesis ^[5], primitive streak formation ^[6], slime mould formation ^[7], and wound healing ^[8], just to name a few.

Interpretation of the underlying mechanism of chemotaxis using continuum PDE models dates bak to the 1950s by Patlak ^[9] and the 1970s by Keller and Segel ^[10,11,12], the former from a probabilistic perspective, and the latter based on a phenomenological approach. In a general form, the Patlak-Keller-Segel model can be expressed as

which is a system of reaction-diffusion-advection equations. Here, the unknown functions $u(x, t)$ and $c(x, t)$ denote the density of the organic population and concentration of the chemical signal at position $x$ at time $t$, respectively. The parameters: $D > 0$ stands for the diffusion coefficient of the organic density; $\chi\neq0$ is the coefficient of chemotactic sensitivity, the sign of $\chi$ dictates whether the chemotaxis is attractive ($\chi > 0$) or repulsive ($\chi < 0$), with $|\chi|$ measuring the strength of chemotactic response; and $\varepsilon\ge0$ is the diffusion coefficient of the chemical signal. The function $\Phi(c)$ denotes the chemotactic sensitivity, which underlines the main character of such a model, whose spatial derivative depicts the mechanistic feature of chemotactic movement – advection of the organic population induced by the spatial gradient of the chemical signal in the local environment. Moreover, the function $f(u, c)$ accounts for consumption/production/degradation of the chemical signal. Formally, the model portrays the evolution (biased movement) of the organic population in response to the chemical stimulus in the local environment. Mathematically, the synergy of random diffusion, nonlinear advection and nonlinear reaction makes the dynamics of the model an intriguing problem to investigate. Depending on the specific forms of $\Phi(c)$ and $f(u, c)$, the model can be utilized to explain the underlying mechanisms of various biological processes experiencing chemotaxis. We refer the readers to the survey papers ^[13,14,15] for a comprehensive list and the mathematical development of such models.

To introduce the specific problem, we first present the following model, which is a special version of (1.3) when $\Phi$ is the logarithmic function of $c$ and $f(u, c)$ is a polynomial of $u$ and $c$:

where the parameter $\mu\neq0$ denotes the density-dependent production/degradation rate of the chemical signal and $\sigma\ge0$ stands for the natural degradation rate of the chemical signal. When $\sigma = 0$, this model can be viewed as a borderline case of one of the original models developed by Keller and Segel:

which appeared in ^[12] to understand J. Adler's experimental result ^[16] on the formation of traveling bands in nutrient-enticed E. Coli population. The logarithmic sensitivity entails that the chemotactic movement obeys the Weber-Fechner's law, which is a fundamental principle in psychophysics and has prominent applications in biology (^{[12,17,18,19,20]}).

System (1.4) appeared in a variety of contexts to elucidate the fundamental principles of different chemotactic processes. For example, when $\chi < 0$, $\mu < 0$ and $\sigma > 0$, the model was designed in ^[21,22] to illustrate the chemotactic movement of reinforced random walkers, such as surface or matrix-bound adhesive molecules, which deposit non-diffusive ($\varepsilon = 0$) or slowly moving ($0 < \varepsilon\ll 1$) chemical signals that modify the local environment for succeeding passages. In light of the biological contexts of the parameters $\chi$ and $\mu$, the model characterizes the process of repulsive chemotaxis in response to a chemical stimulus that is being produced by the organism. On the other hand, when $\chi > 0$, $\mu > 0$ and $\sigma = 0$, the model was employed by Levine { et al.} ^[23] to interpret the dynamical interactions between vascular endothelial cells (VECs) and signaling molecules vascular endothelial growth factor (VEGF) in the onset of tumor angiogenesis. In this case, the model portrays the process of attractive chemotaxis in response to a chemical stimulus that is being consumed by the organism.

It should be stressed that although the logarithmic sensitivity plays an important role in biological modeling, its singular nature brings out an obstruction for the qualitative (numerical and rigorous) analysis of the mathematical equations. The existing results concerning the qualitative behavior of the models with logarithmic sensitivity is much less than those with other types of sensitivity, such as linear sensitivity, saturating sensitivity, {et al}, see e.g., ^[13,14,15]. Fortunately, it has been recognized that the singular nature of the logarithmic sensitivity can be removed by applying the Cole-Hopf type transformation: $v = \big[\log (e^{\sigma t}c)\big]_x = \frac{c_x}{c}$. Combining such a transformation with the temporal-spatial rescaling:

one can convert the original chemotaxis model into

where the equations are written in the conservative form to facilitate further discussions in the context of balance laws. A direct calculation shows that the characteristics associated with the flux in (1.6) are

from which we see that the principle part of (1.6) is hyperbolic in biologically relevant regimes (where the cellular density $u > 0$), provided $\chi\mu > 0$. Hence, the analytical tools in hyperbolic balance laws become feasible for studying the qualitative behavior of the model. We focus on this case throughout the paper, since otherwise the characteristic fields may change the type, which could alter the dynamics of the model drastically. This is supported by the finite-time blowup of the explicit and numerical solutions constructed in ^[21].

Rigorous mathematical study of (1.6) has been carried out for nearly two decades. Most of the results are established for the case when $\chi\mu > 0$. The global well-posedness of large amplitude classical solutions was first established in ^[24,25], and upgraded in a series of recent works ^{[26,27,28,29,30,31,32,33,34,35]} under various initial and/or boundary conditions. The existence and local stability of large-strength traveling waves and boundary spike-layer solutions have been settled in ^{[36,37,38,39,40,41,42,43,44]} and ^[45], respectively. There are also works investigating the qualitative behavior of the appended version of (1.6) with logistic growth ^{[46,47,48,49,50,51]}. Collectively, these results suggest that large-time homogenization is generic in the processes of logarithmically sensitive chemo-attraction with chemical consumption and chemo-repulsion with chemical production, while finite-time singularities are not expected in such processes.

Inspired by the biological background and the mathematical research recently conducted on (1.6), the following model was proposed and analyzed in ^[52,53]:

which is an amended version of (1.4) by replacing the linear dependence of the chemical production/consumption rate on the organic density by a power-like one. Since monomials are building blocks of general nonlinear functions, this model paves a way for the understanding of the dynamical behavior of the Keller-Segel type model with genuinely nonlinear production/consumption rates. By applying the Cole-Hopf transformation and the rescaling applied to (1.4), we obtain

which is the focal point of this paper. From the perspective of rigorous mathematical analysis, (1.9) is considerably more complicated than (1.6), due to the power-like nonlinearity. So far, the global stability of constant ground states and convergence rates in the regime of classical solutions to the Cauchy problem has been proven in ^[53,54], and the global stabilization of classical solutions on finite intervals subject to dynamic boundary conditions has been demonstrated in ^[52].

1.2.   Motivation

Mathematically, the initial-boundary value problem of a PDE model subject to dynamic (time-dependent) boundary conditions can be equivalently viewed as a control problem with the boundary data as external inputs. A typical question in this type of problem is under what conditions on the boundary data the solution stabilizes in the long run. Such a problem has been studied for (1.9) in ^[52], where the global stabilization of large data classical solutions under dynamic boundary conditions is established. However, we note that the result of ^[52] required the boundary values of the unknown functions to match at the endpoints, i.e., $u(a, t) = u(b, t)$ and $v(a, t) = v(b, t)$. Apparently, such a restriction is not desirable from the physical/biological point of view. Rigorous mathematical study of the models subject to unmatched boundary conditions, i.e., $u(a, t)\neq u(b, t)$ and $v(a, t)\neq v(b, t)$, thus becomes physically/biologically relevant. Nevertheless, the long-time dynamics of large data solutions to the model subject to unmatched boundary conditions still remains to be investigated. This has been highlighted in Remark 1.1 of ^[52]. To enrich the contemporary body of knowledge, we dedicate this paper to the investigation of such a problem. Since $\gamma$ serves as a tuning parameter and the dynamic boundary data are regarded as external inputs, our major task is to identify the conditions on the parameter value and the boundary data, under which solutions to the initial-boundary value problem of the model stabilize when time getting large. Comparing with the case of matched boundary data, the analysis under unmatched boundary conditions is much more involved, since the reference profile depends on the spatial variable. Hence, the combined effect of the spatially dependent reference profile coupled with power-like/quadratic nonlinearities creates an intriguing mathematical problem to pursue.

1.3.   Statement of results

To simplify the presentation, we take $\chi = D = \varepsilon = 1$, since the specific values of the parameters do not affect the underlying analysis. Also, recall that we are focusing on the case when $\chi\mu > 0$. Under these circumstances, we obtain the cleaner version of (1.9):

The system is supplemented with the initial conditions

and the boundary conditions:

The main result of this paper addresses the global nonlinear stability of large data classical solutions to (1.10) when $\gamma = 2$ and subject to unmatched dynamic boundary conditions.

Theorem 1.1. Consider the initial-boundary value problems (1.10)–(1.12) with $\gamma = 2$. Suppose the initial data satisfy $u_0 > 0$, $(u_0, v_0) \in [H^2((0, 1))]^2$, and are compatible with the boundary conditions. Assume $\alpha_1$, $\alpha_2$, $\beta_1$, and $\beta_2$ are smooth functions on $[0, \infty)$, and satisfy the following conditions:

where $\underline{\alpha}_i$ are constants. Then there exists a unique solution to the IBVP, such that

where

and the constant $C > 0$ is independent of $t$. Moreover, the solution obeys the long-time behavior

We have the following remarks concerning Theorem 1.1.

Remark 1.1. Since we did not specify any decay rate for the boundary data, the rate of convergence of the perturbed functions can not be identified under the assumptions of Theorem 1.1. The integrability conditions of the boundary data present general criteria that guarantee the convergence of the perturbations. On the other hand, if we specify certain decay rates for $\alpha_i$ and $\beta_i$, then the convergence rate of ${\tilde u}$ and ${\tilde v}$ could be determined by using ODE arguments. We refer the reader to the example in the Appendix of ^[55], whose proof can be adapted to fit into the situation encountered in this paper.

Remark 1.2. The technical assumptions in Theorem 1.1 imply the difference between the corresponding boundary data will converge to zero as $t\to\infty$, i.e., the unmatched boundary data will eventually match. However, based on the assumptions we see that the boundary functions do not necessarily equal to each other at any finite time. In addition, there is no smallness restriction on the amplitude of the initial perturbation around the reference profile.

Remark 1.3. Since the boundary data are smooth functions on $[0, \infty)$ and their first order derivatives belong to $W^{1, 1}((0, \infty))$, it follows from the Fundamental Theorem of Calculus that the functions themselves and their first order derivatives are uniformly bounded with respect to $t$. Such information will be frequently utilized in the proof of the theorem.

Remark 1.4. Our result only covers a single value of $\gamma$, i.e., $\gamma = 2$. For other values of $\gamma$, it is not clear whether similar results hold true or not. We provide numerical results towards the end of the paper to present some positive evidences, while leave the rigorous investigation for future work. On the other hand, we would like to mention that by combining the analysis of this paper with the techniques in ^[33], especially the Orlicz-type free energy functional, one can establish a similar result for the original model (1.6) under the boundary conditions: $u(0, t) = u(1, t) = \alpha(t)$, $v(0, t) = \beta_1(t)$, $v(1, t) = \beta_2(t)$, where $\beta_1(t)\neq \beta_2(t)$, and $\alpha$, $\beta_1$, $\beta_2$ satisfy similar conditions as those in Theorem 1.1. However, the problem with unmatched boundary data for $u(x, t)$ can not be settled at this point, due to technical obstructions caused by the Orlicz-type functional. We also leave the investigation for future work.

We prove Theorem 1.1 in Section 2. The proof is self-contained, utilizing delicate energy methods. We apply various elementary inequalities, such as Cauchy-Schwarz, Young, Hölder and Grönwall. The energy estimates are carefully crafted to unveil the assumptions about the growth/decaying properties of the boundary data. In Section 3, we present numerical results to illustrate some problems that are not covered by Theorem 1.1, including dynamic boundary data with non-trivial final profiles, time periodic boundary conditions, and the case when $\gamma$ is a fraction. For notational convenience, throughout the rest of the paper, we use $\|\cdot\|$ and $\|\cdot\|_\infty$ to denote the standard norms $\|\cdot\|_{L^2}$ and $\|\cdot\|_{L^\infty}$, respectively. The functional space $H^s((0, 1))$ is denoted by $H^s$. Moreover, we use $C$ to denote a generic constant which is independent of time, but may depend on the initial and/or boundary data. The value of the constant may vary line by line according to the context.

2.   Proof of Theorem 1.1

We now prove Theorem 1.1. The proof is divided into several steps which are contained in a series of subsections. We define the reference profiles:

which interpolate the boundary data in the linear fashion. To facilitate the asymptotic analysis, we let

In terms of ${\tilde u}$ and ${\tilde v}$, system (1.10) with $\gamma = 2$ reads as

with the initial and boundary conditions

First of all, under the assumptions of Theorem 1.1, the local well-posedness of solutions to (2.1), (2.2) and (2.3) can be established by using standard approaches, such as Galerkin approximation and fixed point argument, along with the energy estimates derived in this section. We present the result without going through the technical details in order to simplify the presentation.

Lemma 2.1. Under the assumptions of Theorem 1.1, there exist a finite $T_0 > 0$ and a unique solution $({\tilde u}, {\tilde v})$ to (2.1), (2.2) and (2.3), such that $({\tilde u}, {\tilde v}) \in \left[C([0, T_0];H^2) \cap L^2(0, T_0;H^3)\right]^2$.

Next, we shall derive the a priori estimates of the local solution, in order to extend it to a global one. The a priori estimates alongside their proofs are recorded in the following subsections.

2.1.   $L^\infty_tL^2_x$–$L^2_tH^1_x$–Estimates

Lemma 2.2. Under the assumptions of Theorem 1.1, there exists a constant $C > 0$ which is independent of $t$, such that

Proof. Taking $L^2$ inner product of (2.1a) with $2 {\tilde u}$ and integrating by parts, we can show that

Taking $L^2$ inner product of (2.1b) with ${\tilde v}$ and integrating by parts, we can show that

Taking the sum of (2.4) and (2.5), we obtain

Since $\alpha_x = \alpha_2(t)-\alpha_1(t)$ and $\beta_x = \beta_2(t)-\beta_1(t)$, we infer that

Using the uniform boundedness of $\alpha$ and $\beta$, see Remark 1.3, we can show that

Similarly, it can be shown that

Substituting the above estimates into (2.6), we can show that

Applying Grönwall's inequality and using the assumptions of Theorem 1.1, we can show that

for some constant which is independent of $t$. This completes the proof of the lemma.

2.2.   $L^\infty_tH^1_x$–$L^2_tH^2_x$–Estimates

Lemma 2.3. Under the assumptions of Theorem 1.1, there exists a constant $C > 0$ which is independent of $t$, such that

Proof. Step 1. Taking $L^2$ inner product of (2.1a) with $- {\tilde u}_{xx}$, we obtain

Applying the Sobolev embedding theorem and Poincaré inequality, we can show that

Using the uniform boundedness of $\alpha_x = \alpha_2(t)-\alpha_1(t)$ and $\alpha$, we deduce

Similarly, we can show that

For $R_{2d}$ and $R_{2e}$, we can show that

and

Feeding the above estimates to (2.9), applying {Cauchy's inequality}, and invoking the uniform boundedness of $\alpha$, $\beta$ and their derivatives (see Remark 1.3), we can show that

Applying Grönwall's inequality to (2.10) and utilizing Lemma 2.2, we can show that

where the constant is independent of $t$.

Step 2. Taking $L^2$ inner product of (2.1b) with $- {\tilde v}_{xx}$, we deduce

Similar to the estimate of $R_{2a}$, we can show that

where we applied (2.11) for $\| {\tilde u}_x\|^2$. Similar to the estimates of $R_{2b}$ and $R_{2c}$, we can show that

Similar to the estimate of $R_{2d}$, we can show that

Similar to the estimate of $R_{2e}$, we can show that

Substituting the above estimates into (2.12), we can show that

Applying Grönwall's inequality to (2.13) and utilizing Lemma 2.2, we can show that

for some constant which is independent of $t$. This completes the proof of the lemma.

2.3.   $L^\infty_tH^2_x$–$L^2_tH^3_x$–Estimates

Lemma 2.4. Under the assumptions of Theorem 1.1, there exists a constant $C > 0$ which is independent of $t$, such that

Proof. Since the information about the higher order spatial derivatives of the solution is unknown at the boundary points, the usual procedure (differentiating with respect to $x$) for estimating the $L^\infty_tH^2_x$ and $L^2_tH^3_x$ norms of the solution can not be directly implemented here. To circumvent such a technical issue, we turn to the estimation of the temporal derivatives of the solution, then utilize the model equations to recover the estimates of the spatial derivatives.

Step 1. Taking $\partial_t$ of (2.1a) and (2.1b), we have

Taking $L^2$ inner product of (2.15b) with $\tilde{u}_t$ and integrating by parts, we can show that

Using the Sobolev embedding theorem and Lemmas 2.2–2.3, we can show that

Since $\alpha_t = [\alpha_2'(t)-\alpha_1'(t)]x+\alpha_1'(t)$ is uniformly bounded (see Remark 1.3), we infer that

Since $\alpha$ is uniformly bounded, we can show that

Since $\alpha_{xt} = \alpha_2'(t)-\alpha_1'(t)$ and $\alpha_x = \alpha_2(t)-\alpha_1(t)$ are uniformly bounded, we can show that

where we applied Poincaré inequality to ${\tilde v}$, and

Similarly, using the uniform boundedness of $\beta_t$ and $\beta_{xt}$, we can show that

where we applied Poincaré inequality to ${\tilde u}$. Lastly, $R_{4h}$ is estimated as

Substituting the above estimates into (2.16), we have

Step 2. Taking $L^2$ inner product of (2.15b) with $\tilde{v}_t$ and integrating by parts, we can show that

Using similar arguments as those in deriving the estimates of $R_{4a}$–$R_{4h}$, we can show that

Step 3. Using (2.1a) and similar arguments as in Step 1, we can show that

In the same spirit, by using (2.1b), we can show that

Using (2.20) and (2.21), we update (2.17) and (2.19) as

and

Applying Grönwall's inequality to (2.22) and (2.23), using Lemmas 2.2 and 2.3 and the assumptions in Theorem 1.1, we can show that

As a consequence of the above estimates, we can show by using (2.1a) and (2.1b) that

This completes the proof of Lemma 2.4.

Lemmas 2.2–2.4 established the desired energy estimates of the solution as stated in Theorem 1.1. The global well-posedness of the initial-boundary value problem thus follows from the local well-posedness in Lemma 2.1, a priori estimates in Lemmas 2.2–2.4 and standard continuation argument. To complete the proof of Theorem 1.1, it remains to derive the temporal decaying of the perturbed variables, which is carried out in the next subsection.

2.4.   Decay estimate

We first recall a technical lemma which is essential in deriving the temporal decaying of the perturbations. The proof of the lemma can be found in ^[53].

Lemma 2.5. Let $f \in W^{1, 1}(0, \infty)$ be a nonnegative function. Then $f(t) \to 0$ as $t \to \infty$.

Using Lemma 2.5, we can show the following:

Lemma 2.6. Under the assumptions of Theorem 1.1, $\|\tilde{u}(t)\|_{H^2}^2+\|\tilde{v}(t)\|_{H^2}^2 \to 0$, as $t\to\infty$.

Proof. Step 1. Based on Lemma 2.2, we know that

which, together with Poincaré inequality, implies

Using the assumptions of Theorem 1.1 and the energy estimates established in Lemmas 2.2–2.4, we can deduce from (2.6) and (2.7) that

where the constant is independent of $t$. Integrating the above inequality, we obtain

Collectively, (2.26) and (2.27) imply

Hence, $\|\tilde{u}(t)\|^2+\|\tilde{v}(t)\|^2 \to 0$, as $t\to\infty$, according to Lemma 2.5. Moreover, since $\| {\tilde u}_{xx}(t)\|$ and $\| {\tilde v}_{xx}(t)\|$ are uniformly bounded with respect to time, the decaying of $\|\tilde{u}_x(t)\|^2+\|\tilde{v}_x(t)\|^2$ follows from the decaying of $\|\tilde{u}(t)\|^2+\|\tilde{v}(t)\|^2$ and the interpolation inequality $\|f_x\|^2 \le \|f\|\cdot\|f_{xx}\|$, for any $f\in H_0^2((0, 1))$.

Step 2. For the second order spatial derivatives, we know from (2.24) and Poincaré inequality that

Based on (2.17), (2.19), and previous estimates, it can be shown that

which implies

Hence, $\|\tilde{u}_t(t)\|^2+\|\tilde{v}_t(t)\|^2 \in W^{1, 1}((0, \infty))$, and thus $\|\tilde{u}_t(t)\|^2+\|\tilde{v}_t(t)\|^2 \to 0$, as $t\to\infty$. Using (2.1), we can show that

Since $\alpha_1-\alpha_2$, $\alpha_1$, $\alpha_2$, $\beta_1-\beta_2$, $\beta_1'$ and $\beta_2'$ belong to $W^{1, 1}((0, \infty))$, they tend to zero as $t\to\infty$. Therefore, the decaying of $\|\tilde{u}_{xx}(t)\|^2+\|\tilde{v}_{xx}(t)\|^2$ follows from (2.28) and the decaying of the first order derivatives of the perturbed solution. This completes the proof of Lemma 2.6.

Collectively, Theorem 1.1 follows from Lemmas 2.1–2.4 and 2.6.

3.   Numerical experiments

We have shown that large data classical solutions to an initial boundary value problem of the nonlinear PDE system (1.10), subject to unmatched Dirichlet type dynamic boundary conditions, stabilize in the long run under appropriate assumptions on the boundary data. However, it was mentioned in Remark 2 that the technical assumptions in Theorem 1.1 imply the unmatched boundary data will eventually match at the endpoints, i.e., the final states are constants. This indeed rules out the more physically/biologically relevant situations in which the final equilibria are non-trivial. Another situation that is not covered by Theorem 1.1 is when the boundary data are time periodic. Because of the ubiquitous presence of time periodic phenomena in natural sciences, it is worth investigating whether time periodic boundary data generate time periodic solutions of the nonlinear PDE model. Moreover, our theoretical result only included the special case $\gamma = 2$. It is not clear whether the results hold true for other values of $\gamma$, especially when $\gamma$ is a fraction. Inspired by these observations, we devote the last part of this paper to numerically showcasing two problems that are not covered by the analytical results. These include dynamic boundary data with unmatched final states and time periodic boundary conditions. We use a finite difference scheme to simulate the problems. We hope that our numerical studies will shed some light on, and provide guidance for further analytical studies in this specific area of research.

Recall the initial-boundary value problem

where the initial and boundary data are compatible, i.e.,

3.1.   Dynamic boundary data with unmatched final states

In Example 1, we set $u_0 = 5\sin(\pi x/2)+10$, $v_0 = 2+\cos(\pi x/2)$, $\alpha_1(t) = 9+\frac{1}{1+t}$, $\alpha_2(t) = 14+e^{-t}$, $\beta_1(t) = 4-e^{-t}$, $\beta_2(t) = 2+\frac{t}{1+t}$. The value of $\gamma$ is set as $\gamma$ = 1.5, 2 and 2.5. The numerical solutions of $u$ and $v$ evolve to non-trivial steady states as shown in Figure 1.

Meanwhile, it is interesting to observe from Figure 2 that the graphs of the steady state solutions steepen near the boundary points as $\gamma$ increases. We thus expect boundary layers will develop near the endpoints of the interval as $\gamma$ grows. Mathematical analysis of this problem requires more sophisticated tools and we leave the investigation in a future paper.

3.2.   Time-Periodic boundary conditions

In Examples 2A, 2B and 2C, we simulate the problem by using time-periodic boundary data. The simulations are run for the same values of $\gamma$ as in Example 1. It turns out, however, the simulation results are almost the same for different values of $\gamma$ and the plots are visibly indistinguishable. The plots presented below are for $\gamma = 1.5$.

In Example 1A, the initial values are $u_0 = 2+x(1-x)$ and $v_0 = 4+\sin(2\pi x)$. The boundary conditions are $\alpha_1(t) = \alpha_2(t) = 2+\sin(2\pi t)$ and $\beta_1(t) = \beta_2(t) = 4+\sin(2\pi t)$. Thus the boundary conditions of $u$ and $v$ have the same period and same phase angle. The solutions are shown in Figure 3, where both $u$ and $v$ evolve to periodic solutions at later time, although initially they are not. Note the maximum and minimum values of $u$ and $v$ occur at the same time.

In Example 2B, the initial values of $u$ and $v$ are the same as those in Example 2A. The boundary condition of $u$ is still $\alpha_1(t) = \alpha_2(t) = 2+\sin(2\pi t)$, but that of $v$ is changed to $\beta_1(t) = \beta_2(t) = 3+\cos(2\pi t)$. That is, they share the same period but have different phase angles. The solutions are shown in Figure 4, where both $u$ and $v$ evolve to periodic solutions. Note the maximum and minimum values of $u$ and $v$ occur at a gap of $1/4$ time units, which originates from the difference of the phase angles of $\alpha$ and $\beta$.

In Example 2C, the initial values of $u$ and $v$ are the same as in Examples 2A and 2B. The boundary condition of $u$ is changed to $\alpha_1(t) = 2+\sin(2\pi t)$ and $\alpha_2(t) = 1+\cos(2\pi t)$, and that of $v$ is changed to $\beta_1(t) = 4+\sin(2\pi t)$ and $\beta_2(t) = 3+\cos(2\pi t)$. That is, the left boundary conditions for $u$ and $v$ have the same phase angle and the right boundary conditions have another phase angle. The solutions are shown in Figure 5, where both $u$ and $v$ evolve to periodic solutions.

4.   Conclusions

In this paper, we studied the global dynamics of classical solutions to the nonlinear PDE system (1.1) subject to dynamic boundary conditions. When $\gamma = 2$, we proved that if the time-dependent Dirichlet-type boundary conditions satisfy (1.13)–(1.15), there exists a unique global-in-time solution to the initial-boundary value problem (IBVP) for any initial data in $H^2((0, 1))$, and the solution is shown to converge asymptotically to the steady state determined by the boundary data as time goes to infinity. In our analytical results, the boundary values of the solution are not required to match at any time. This generalized the previous work ^[52], where the boundary values must match for all time. To investigate some problems that are not covered by our analytical results, we numerically simulated the IBVP for dynamic boundary conditions with unmatched end-states or periodic in time. Our numerical results suggested that in the former case, the solution converges to a non-trivial steady state in the long run, and in the latter case, time periodicity was observed in the solution after certain duration of time. Moreover, we simulated the problem for several values of $\gamma$ different from 2 and discovered the phenomenon of boundary steepening of the steady state as the value of $\gamma$ increases. Rigorous mathematical studies of these phenomena will be carried out in future works.

Acknowledgments

The authors would like to thank the anonymous referees for their valuable comments and suggestions which help substantially improve the quality of the paper. Support for this work came in part from a National Natural Science Foundation of China Award 12171116 (LX), a Fundamental Research Funds for Central Universities of China Award 3072020CFT2402 (LX), and from Simons Foundation Collaboration Grant for Mathematicians Award 413028 (KZ).

Conflict of interest

The authors declare that there is no conflict of interest.