Global stability of a novel nonlinear diffusion online game addiction model with unsustainable control

1.   Introduction

The World Health Organization ^[1] regards the online game addiction (OGA) as a new disease. The general defining feature of this disease is "persistent and recurrent use of the internet to engage in games, often with other players, leading to clinically significant impairment or distress" (see ^[2]). OGA is often an addictive mental disease. People with OGA tend to be impulsive, violent, misanthropic and withdrawn. As the number of internet users continues to increase, the number of people addicted to online games, especially children, is also growing. The research on OGA has attracted wide attention. The diagnostic and statistical manual of mental disorders ^[3,4] proposed some classification criteria for OGA. To better cure and reduce the population of OGA, scholars from all walks of life began to study it from various perspectives. One of the important research methods to reveal the dynamic behavior of OGA is using mathematical models. In recent years, some mathematical workers have made use of various mathematical theories and methods to study the OGA problem. For example, Li and Guo ^[5,6] studied the stability and optimal control of two kinds of OGA mathematical models. Viriyapong and Sookpiam ^[7] discussed the stability and qualitative behavior of an OGA model for children and youth in Thailand. Therefore, we try to apply a differential equation model to study the OGA problem in this manuscript.

Inspired by the previous works ^[5,6,7,8] in terms of population division of online gamers, based on the ideas and methods of the papers ^{[9,10,11,12,13]}, we put forward a mathematical model to describe the changes and interactions of population density with OGA under unsustainable treatment as follows:

where $\mathcal{U}(\tau)$ and $\mathcal{V}(\tau)$ stand for the population densities of moderate gamers and addictive gamers at time $\tau$, respectively. $r$ stands for the proportion coefficient of population density of new gamers. $a\mathcal{U}(\tau)$ is the population density of stopping gaming for moderate gamers through self-control. $b\mathcal{U}(\tau)\mathcal{V}(\tau)$ indicates the number of population density converted from moderate gamers to addictive gamers. $\beta\Theta(\tau)\mathcal{V}(\tau)$ measures the unsustainable treatment for an addictive gamer. $r, a, b, \beta > 0$ are some constants.

As we all know, the number of moderate and addictive gamers not only depends on time variables $\tau$, but also often shows a certain diffusion in space. Therefore, we modify the system (1.1) to the following diffusion partial differential equation model:

where the practical meanings of $\mathcal{U}(\tau, \cdot)$, $\mathcal{V}(\tau, \cdot)$, $r$, $a\mathcal{U}(\tau, \cdot)$, $b\mathcal{U}(\tau, \cdot)\mathcal{V}(\tau, \cdot)$ and $\beta\Theta(\tau, \cdot)\mathcal{V}(\tau, \cdot)$ are the same as in system (1.1). $d_i > 0 \; (i = 1, 2)$ is the diffusion coefficient. $\Delta$ is the Laplacian. $\Omega\subset\mathbb{R}$ is bounded, and its boundary $\partial\Omega$ is smooth. The unsustainable treatment function $\Theta(\tau, \cdot)$ has the following characteristics:

$(A_1)$ $\Theta(\tau, \cdot)$ is piecewise continuous with respect to $\tau\in[0, \infty)$, that is, $\Theta(\tau, \cdot)$ is continuous in the interval $I\in[0, \infty)\backslash\mathcal{B}$, where $\mathcal{B} = \{\tau_k\in[0, \infty):\, \Theta(\tau_k^-, \cdot) < \Theta(\tau_k^+, \cdot)\}$ is countable. For any compact subset of $[0, \infty)$, $\Theta(\tau, \cdot)$ has a finite number of discontinuities. For all $\tau_1, \tau_2\in[0, \infty)$, $\Theta(\tau_1, \cdot)\leq\Theta(\tau_2, \cdot)$ provided that $\tau_1 < \tau_2$. $0 < \gamma\leq\Theta(\tau, \cdot)$, and $\Theta(0, \cdot) = \Theta(0^+, \cdot) = \gamma$.

To study the dynamic behavior of model (1.2), the following initial and boundary conditions are necessary:

where $\overrightarrow{n}$ is the outward normal vector of $\partial\Omega$. The initial values $\mathcal{U}_0(\cdot)$ and $\mathcal{V}_0(\cdot)$ satisfy the below assumption:

$(A_2)$ For $\forall\, x\in\Omega$, $\mathcal{U}_0(\cdot), \mathcal{V}_0(\cdot) > 0$, and $\mathcal{U}_0(\cdot), \mathcal{V}_0(\cdot)\in L^{\infty}(\Omega)$.

It is similar to ^[14,15] that the model (1.2) with initial and boundary conditions (1.3) is rewritten of the form

where $\overline{[\Theta(\tau, \cdot)]} = [\Theta(\tau^-, \cdot), \Theta(\tau^+, \cdot)]$, $\Theta(\tau^-, \cdot) = \lim_{\tau\to \tau^-}\Theta(\tau, \cdot)$, $\Theta(\tau^+, \cdot) = \lim_{\tau\to \tau^+}\Theta(\tau, \cdot)$. Obviously, $\overline{[\Theta(\tau, \cdot)]} = \Theta(\tau_k, \cdot)$ provided that $\Theta(\tau, \cdot)$ is continuous at $\tau = \tau_k$.

This paper mainly studies the well-posedness and stability of the solution of the model (1.2) with initial and boundary conditions (1.3) with the help of (1.4). The principal illuminating points of our research work are as below:

$(\mathrm{i})$ By analyzing the mechanism of OGA and considering the unsustainable treatment, we establish a new ordinary differential equation model (1.1). Considering that the population of gamers has a certain diffusion in space, we propose a novel diffusion partial differential equation model (1.2).

$(\mathrm{ii})$ We first convert the model (1.2) with initial and boundary conditions (1.3) into differential inclusion (1.4). Next, we study the boundness and existence of solutions to the model (1.2) with initial and boundary conditions (1.3). Meanwhile, we build some Lyapunov type functionals and prove the global stability of the model. Finally, an example is numerically simulated and analyzed using MATLAB. The numerical simulation is in good agreement with the theoretical results.

$(\mathrm{iii})$ The major outcomes of the paper provide theoretical help for the research and treatment of online game addiction. Our work provides an example for applying mathematical theories and methods to solve social problems such as online game addiction, which makes the study of this kind of problem transform from the qualitative to the quantitative.

The paper's remaining structure is as follows. In Section 2, we prove that the model has bounded solutions. In Section 3, the global stability of the steady states of the model is studied. A numerical example and simulation are given in Section 4. Finally, a brief summary is made in Section 5.

2.   Existence and boundedness

Definition 2.1. Let $\phi_1(\mathcal{U}, \mathcal{V}) = r-a\mathcal{U}-b\mathcal{U}\mathcal{V}$, $\phi_2(\mathcal{U}, \mathcal{V}) = b\mathcal{U}\mathcal{V}-\beta\overline{[\Theta(\tau, \cdot)]}\mathcal{V}$. $(\mathcal{U}, \mathcal{V})$ is a solution of (1.4) if and only if, for $\mathcal{U}, \mathcal{V}\in C([0, \sigma], \mathscr{H})$, $\exists\, \Lambda_1(\mathcal{U}, \mathcal{V})\in L^1([0, \sigma], \mathscr{H})$, $\Lambda_2(\mathcal{U}, \mathcal{V})\in \phi_2(\mathcal{U}, \mathcal{V})$ within $(0, \sigma)$ satisfying that $(\mathcal{U}, \mathcal{V})$ is a solution in $(0, \sigma)$ of the following equation:

Set $\mathscr{Y} = C([0, \sigma], \mathscr{H})\times C([0, \sigma], \mathscr{H})$ and define $\mathcal{W} = (\mathcal{U}, \mathcal{V})^T\in\mathscr{Y}$, and then Eq (2.1) is rewritten by

where the operators $\mathcal{L}$ and $\mathcal{F}$ are defined as

The solution of (2.2) is formulated by

where $\mathcal{T}(\tau) = (\mathcal{T}_1(\tau), \mathcal{T}_2(\tau))^T$, and $\mathcal{T}_1(\tau)$ and $\mathcal{T}_2(\tau)$ are compact and strong $C_0$-semigroup induced by $d_1\Delta$ and $d_2\Delta$ with $\frac{\partial \mathcal{U}(\tau, \cdot)}{\partial \overrightarrow{n}} = \frac{\partial \mathcal{V}(\tau, \cdot)}{\partial \overrightarrow{n}} = 0$, respectively.

Theorem 2.1. Assume that $(A_1)$ and $(A_2)$ are true. Then, for any solution $\mathcal{W}(\tau, \cdot) = (\mathcal{U}(\tau, \cdot), \mathcal{V}(\tau, \cdot))$ of system (1.2) with initial and boundary conditions (1.3), there exist the constants $N_1, N_2 > 0$, such that $N_1\leq \mathcal{U}(\tau, \cdot), \mathcal{V}(\tau, \cdot)\leq N_2$, $\forall\, (\tau, x)\in[0, \infty)\times\Omega$.

Proof. Taking $L^-_1 = \min\{\frac{\beta}{b},$ $\min_{x\in\Omega}\mathcal{U}_0(\cdot)\}$, $L^+_1 = \max\{\frac{\beta(\gamma+1)}{b},$ $\max_{x\in\Omega}\mathcal{U}_0(\cdot)\}$, $L^-_2 = \min\{\frac{r-aL^-_1}{bL^-_1},$ $\min_{x\in\Omega}\mathcal{V}_0(\cdot)\}$ and $L^+_2 = \max\{\frac{r}{bL^-_1},$ $\max_{x\in\Omega}\mathcal{V}_0(\cdot)\}$, we define $\mathcal{R}$ by

Obviously, $\mathcal{W}_0(\cdot) = (\mathcal{U}_0(\cdot), \mathcal{V}_0(\cdot))^T\in\mathcal{R}$. It is similar to ^[16] that $\mathcal{T}(t)\mathcal{W}_0(\cdot)\in\mathcal{R}$. Define a vector field

On the left of $\mathcal{R}$, choosing $\mathcal{U} = L^-_1$, $L^-_2\leq \mathcal{V}\leq L^+_2$, one derives from (2.3) and (2.4) that

On the right of $\mathcal{R}$, taking $\mathcal{U} = L^+_1$, $L^-_2\leq \mathcal{V}\leq L^+_2$, it is similar to (2.5) that

On the below of $\mathcal{R}$, associating with (2.3), (2.4), $(A_1)$ and $(A_2)$, letting $\mathcal{V} = L^-_2$, $L^-_1\leq \mathcal{U}\leq L^+_1$, $0 < \gamma\leq\Theta(\tau, \cdot)$, we have

On the above of $\mathcal{R}$, choosing $\mathcal{V} = L^+_2$, $L^-_1\leq \mathcal{U}\leq L^+_1$, $0 < \gamma\leq\Theta(\tau, \cdot)$, we similarly obtain

From (2.5)–(2.8), similar to ^[16,17], one knows that $\mathcal{R}$ is the invariant set of $\mathscr{F}$. Thereby, taking $N_1 = \min\{L^-_1, L^-_2\}$, $N_2 = \max\{L^+_1, L^+_2\}$, one has $N_1\leq \mathcal{U}(\tau, \cdot), \mathcal{V}(\tau, \cdot)\leq N_2$, $\forall\, (\tau, x)\in[0, \infty)\times\Omega$. The proof is completed. To prove the existence of solutions, we need the below fixed point theorem.

Lemma 2.1. ^[18,19] Let $\mathbb{K}$ be a nonempty and weakly compact subset in a real Hilbertspace $\mathscr{H}$ and let $\mathcal{E}:\mathbb{K}\to P(\mathbb{K})$ be such that for each $u\in\mathbb{K}$, $\mathcal{E}(u)$ is closed and convex. If the graph $Graf(\mathcal{E})$ of $\mathcal{E}$ is weakly $\times$ weakly sequentially closed, then $\mathcal{E}$ has at least one fixedpoint, i.e., there exists at least one element $u\in\mathbb{K}$ such that $u\in\mathcal{E}(u)$.

Definition 2.2. ^[20] Let $\mathbb{X}$ be a Lebesgue measurable subset in $\mathbb{R}^q$, $q\geq1$, and $\mathcal{E}:\mathbb{X}\to P(\mathscr{H})$ be a mapping. A function $f:\mathbb{X}\to\mathscr{H}$ is called a selection of $\mathcal{E}$ if $f(z)\in\mathcal{E}(z)$ a.e. $z\in\mathbb{X}$. The set of all the measurable selections of $\mathcal{E}$ is denoted by $Sel(\mathcal{E})$.

Theorem 2.2. The model (1.2) with initial and boundary conditions (1.3) has at least one solution $\mathcal{W}^*(\tau, \cdot) = (\mathcal{U}^*(\tau, \cdot), \mathcal{V}^*(\tau, \cdot))\in C([0, \infty), \mathscr{H})\times C([0, \infty), \mathscr{H})$, for $(\tau, x)\in[0, \infty)\times\Omega$, provided that $(A_1)$ and $(A_2)$ are true.

Proof. Our proof is thanks to ^[20,21]. We first prove that the model (1.2) with initial and boundary conditions (1.3) has at least a solution $\mathcal{W}^*(\tau, \cdot) = (\mathcal{U}^*(\tau, \cdot), \mathcal{V}^*(\tau, \cdot))\in\mathscr{Y}$. In fact, from Theorem 2.1 together with $(A_1)$ and $(A_2)$, one knows that there exist some constants $M, R > 0$ such that $\max\{\|\mathcal{U}^*\|_{\mathscr{H}}, \, \|\mathcal{V}^*\|_{\mathscr{H}}\}\leq M$, and for all $(x, y)\in(\phi_1(\mathcal{U}, \mathcal{V}), \, \phi_2(\mathcal{U}, \mathcal{V}))$, $\max\{\|\mathcal{U}\|_{\mathscr{H}}, \, \|\mathcal{V}\|_{\mathscr{H}}\}\leq M+1$ implies that $\max\{\|x\|_{\mathscr{H}}, \, \|y\|_{\mathscr{H}}\}\leq R+1$. Let

$\mathbb{K}$ is nonempty and weakly compact in $L^2([0, \sigma], \mathscr{H})\times L^2([0, \sigma], \mathscr{H})$. Define $P_{\sigma}:\mathbb{K}\to\mathscr{Y}$ as

where $(\mathcal{U}^*, \mathcal{V}^*)$ is the unique solution on $[0, \sigma]$ of (2.1). Adopting fallacy reduction, one easily verifies that $\max\{\|\mathcal{U}\|_{\mathscr{H}}\, \|\mathcal{V}\|_{\mathscr{H}}\}\leq M+1$. Define the operator $\mathcal{E}:\mathbb{K}\to P(\mathbb{K})$ as

where $P_{\sigma}(\Lambda_1, \Lambda_2) = (\mathcal{U}^*, \mathcal{V}^*)$. From Theorem 3.4 of ^[22], Theorem 2.3.3 of ^[23], Proposition 3.6 of ^[22] and Theorem 3.3 of ^[21], we know that the graph $Graf(\mathcal{E})$ of $\mathcal{E}$ is weakly $\times$ weakly sequentially closed in $\mathbb{K}$. It follows from Lemma 1 that there is $(\Lambda_1, \Lambda_2)\in\mathbb{K}$ such that $(\Lambda_1, \Lambda_2)\in\mathcal{E}((\Lambda_1, \Lambda_2))$, and consequently $P_{\sigma}(\Lambda_1, \Lambda_2) = (\mathcal{U}^*, \mathcal{V}^*)$ is a weak solution of (2.1). Since $(\Lambda_1, \Lambda_2)\in L^2([0, \sigma], \mathscr{H})\times L^2([0, \sigma], \mathscr{H})$, Theorem 3.6 of ^[22] ensures that $P_{\sigma}(\Lambda_1, \Lambda_2) = (\mathcal{U}^*, \mathcal{V}^*)$ is in fact a strong solution of (2.1).

Next, we need to show that this strong solution $\mathcal{W}^*(\tau, \cdot) = (\mathcal{U}^*(\tau, \cdot), \mathcal{V}^*(\tau, \cdot))$ of (2.1) defined on $[0, \infty)$. Indeed, according to the continuation theorem, the maximum existence interval of solution of the strong solution $\mathcal{W}^*(\tau, \cdot) = (\mathcal{U}^*(\tau, \cdot), \mathcal{V}^*(\tau, \cdot))$ is a form of $[0, \sigma^*)$. Assume that $\sigma^* < \infty$, we derive from ^[24] that

which contradicts he boundedness of $\mathcal{W}^*(\tau, \cdot) = (\mathcal{U}^*(\tau, \cdot), \mathcal{V}^*(\tau, \cdot))$. Thus, the model (1.2) with initial and boundary conditions (1.3) has a strong solution $\mathcal{W}^*(\tau, \cdot) = (\mathcal{U}^*(\tau, \cdot), \mathcal{V}^*(\tau, \cdot))$ on $[0, \infty)$. The proof is completed.

3.   Global asymptotic stability

It is easy to see that the model (1.2) with initial and boundary conditions (1.3) has two steady states $S = (\mathcal{U}_0^*, \mathcal{V}_0^*) = (\frac{r}{a}, 0)$ and $S_+ = (\mathcal{U}^*_+, \mathcal{V}^*_+)$ determined by

In light of (3.1), we obtain $h(\mathcal{V})\triangleq\frac{rb}{\beta(a+b\mathcal{V})} = \overline{[\Theta(\tau, \cdot)]} = [\Theta(\tau^-, \cdot), \Theta(\tau^+, \cdot)]$, and the steady state $S_+ = (\mathcal{U}^*_+, \mathcal{V}^*_+)$ is represented by $\mathcal{U}^*_+ = \frac{\beta\Theta^*}{b}$, $\mathcal{V}^*_+ = \frac{rb-a\beta\Theta^*}{b\beta\Theta^*}$, where $\Theta^*\in\overline{[\Theta(\tau, \cdot)]}$. Next, we shall discuss the global asymptotic stability of steady states $S$ and $S_+$ in $\mathcal{R}$.

Theorem 3.1. The semi-positive steady state $S$ of the model (1.2) with initial and boundary conditions (1.3) possesses global asymptotic stability in $\mathcal{R}$, provided that $rb < a\beta\gamma$.

Proof. Clearly, the steady state $S$ is semi-positive. Build a functional

Obviously, $V(\tau)$ is smooth. Let $g(\mathcal{U}) = \mathcal{U}-\frac{r}{a}-\frac{r}{a}\ln\frac{a\mathcal{U}}{r}$, and then $g'(\mathcal{U}) = 1-\frac{r}{a\mathcal{U}}$, which implies that $g(\mathcal{U})\geq g(\frac{r}{a}) = 0$. Thus, we have $V(\tau) > 0$ for all $\tau\neq0$ and $V(0) = 0$. Moreover, according to ^[15], one knows that the set $\{\tau\in\mathbb{R}:V(\tau)\leq\alpha\}$ is bounded for $\alpha\geq0$. Define $\mathcal{G}(\mathcal{U}, \mathcal{V})$ as

From (1.2), $(A_1)$ and $(A_2)$, one knows that the set-valued map $\mathcal{G}(\mathcal{U}, \mathcal{V})$ is upper semi-continuous and non-empty compact convex. For any $\zeta = (\zeta_1, \zeta_2)^T\in\mathcal{G}(\mathcal{U}, \mathcal{V})$, there is a $\Theta^*\in\overline{[\Theta(\tau, \cdot)]}$ such that

From (3.2), one calculates $(\frac{\partial V}{\partial\mathcal{U}}, \frac{\partial V}{\partial\mathcal{V}})\cdot\zeta$ as

After simplification, (3.3) gives

By Green's first identity and (1.3), one has

Moreover, from $rb < a\beta\gamma$, we obtain the following estimation:

Noticing that $V(0) = 0$, provided that $(\mathcal{U}, \mathcal{V}) = (\frac{r}{a}, 0)$, based on ^[15] and (3.6), one knows that $V(\tau)$ is a Lyapunov function of (1.2). According to the boundedness of (1.2) in Theorem 2.1, the parabolic $L^p$-theory, the Sobolev Embedding Theorem and the standard compactness argument (see ^[25]), one knows that there exist some constants $M, \tau_0 > 0$ such that $\|\mathcal{U}(\tau, \cdot)\|_{C^2(\overline{\Omega})}+\|\mathcal{V}(\tau, \cdot)\|_{C^2(\overline{\Omega})}\leq M$, for all $\tau\geq\tau_0$. Thereby, it follows from the Sobolev Embedding Theorem ^[20] that $(\mathcal{U}, \mathcal{V})\to(\frac{r}{a}, 0)$ in $L^2(\Omega)\times L^2(\Omega)$, as $\tau\to\infty$. Also, $\frac{dV}{d\tau} = 0$ iff $(\mathcal{U}, \mathcal{V}) = (\frac{r}{a}, 0)$, which implies that $\{(\mathcal{U}, \mathcal{V}):\frac{dV}{d\tau} = 0\} = \{(\frac{r}{a}, 0)\}$. Thus, the steady state $S$ attracts all the solutions of system (1.2). The proof is completed.

Theorem 3.2. The positive steady state $S_+$ of the model (1.2) with initial and boundary conditions (1.3) possesses global asymptotic stability in $\mathcal{R}$, provided that $rb\geq a\beta\gamma$.

Proof. When the condition $rb\geq a\beta\gamma$ is true, $\mathcal{V}^*_+ = \frac{rb-a\beta\Theta^*}{b\beta\Theta^*} > 0$. Thus, one knows that $S_+$ is a unique positive steady state of system (1.2). The key of the proof is the construction of Lyapunov functional $W(\tau)$ and the estimation of $\frac{dW}{d\tau}$. In fact, one can construct the Lyapunov functional $W(\tau)$ as

By (3.7), one has

By simplification, similar to (3.5), one derives from the condition $rb\geq a\beta\gamma$ and (3.8) that

In view of (3.7), (3.9) and the remaining discussion being similar to Theorem 3.1, one concludes that the positive steady state $S_+$ attracts all the solutions of system (1.2). The proof is completed.

4.   Illustrative example and simulation

This section considers the following reaction-diffusion OGA model:

where $\Omega = (0, 10)$, $d_1 = 0.6$, $d_2 = 0.3$, $\Theta(\tau, \cdot) = 0.8+\frac{0.2}{\pi}\arctan([\tau])$, $\mathcal{U}_0(x) = 1+|\sin(x)|$, $\mathcal{V}_0(x) = 1+|\cos(x)|$, and $[\tau]$ is a function such that only integral numbers are taken for the variable $\tau$. Obviously, the set of unsustainable points of unsustainable treatment $\Theta(\tau, \cdot)$ is the natural number set $\mathbb{N}^+$, which is infinitely countable.

4.1.   Theoretical analysis

Case 1: Take $r = 1$, $a = 0.2$, $b = 0.6$, $\beta = 4$. From (4.1), one has $\gamma = 0.8$, $\min_{x\in\Omega}\mathcal{U}_0(x) = \min_{x\in\Omega}\mathcal{V}_0(x) = 1$, $\max_{x\in\Omega}\mathcal{U}_0(x) = \max_{x\in\Omega}\mathcal{V}_0(x) = 2$, $L^-_1 = \frac{20}{3}$, $L^+_1 = 12$, $L^-_2 = 1$, $L^+_2 = 2$. By Theorems 2.1 and 2.2, one knows that Eq (4.1) has at least a strong solution $\mathcal{W}^*(\tau, \cdot) = (\mathcal{U}^*(\tau, \cdot), \mathcal{V}^*(\tau, \cdot))$, for $(\tau, x)\in[0, +\infty)\times(0, 10)$, satisfying $1\leq \mathcal{U}^*(\tau, \cdot), \mathcal{V}^*(\tau, \cdot)\leq12$. Moreover, the condition $0.6 = rb < a\beta\gamma = 0.64$ is true, and it follows from Theorem 3.1 that the semi-positive steady state $S = (\mathcal{U}_0^*, \mathcal{V}_0^*) = (\frac{r}{a}, 0) = (5, 0)$ of (4.1) has global asymptotic stability.

Case 2: Take $r = 1$, $a = 0.2$, $b = 0.6$, $\beta = 3$. Similar to Case 1, $\gamma = 0.8$, $L_1^- = 5$, $L_1^+ = 9$, $L_2^- = 1$, $L_2^+ = 2$. We conclude from Theorems 2.1 and 2.2 that Eq (4.1) has at least a strong solution $\mathcal{W}^*(\tau, \cdot) = (\mathcal{U}^*(\tau, \cdot), \mathcal{V}^*(\tau, \cdot))$, for $(\tau, x)\in[0, +\infty)\times(0, 10)$, satisfying $1\leq \mathcal{U}^*(\tau, \cdot), \mathcal{V}^*(\tau, \cdot)\leq9$. By the condition $0.6 = rb\geq a\beta\gamma = 0.48$, we know that the positive steady state $S_+ = (\mathcal{U}^*_+, \mathcal{V}^*_+)$ is given by $\mathcal{U}^* = \frac{\beta\Theta^*}{b} = 5\Theta^*$, $\mathcal{V}^* = \frac{rb-a\beta\Theta^*}{b\beta\Theta^*} = \frac{1-\Theta^*}{3\Theta^*}$ of (4.1), and has global asymptotic stability, where

Remark 4.1. Indeed, in Case 1, when the cure ratio of addictive gamers $\beta > \frac{rb}{a\gamma} = 3.75$, mathematically speaking, as long as the cure is powerful enough, the addictive gamers will vanish, and the population density of moderate gamers will tend to $\frac{r}{a} = 5$. Unfortunately, this case is difficult to achieve in practice. Generally speaking, the moderate gamers and addictive gamers coexist. People can only reduce the number of addictive gamers through treatment, and this treatment is intermittent and unsustainable. Therefore, Case 2 is more realistic.

Remark 4.2. In Case 2, due to the influence of unsustainable treatment $\Theta(\tau, \cdot) = 0.8+\frac{0.2}{\pi}\arctan([\tau])$, the positive steady state $S_+ = (\mathcal{U}^*_+, \mathcal{V}^*_+)$ varies with time variable $\tau$. By formula (4.2), when $\tau\in(0, \infty)\setminus\mathbb{N}^+$, $\Theta^* = 0.8+\frac{0.2}{\pi}\arctan([\tau])$, the positive steady state $S_+ = (\mathcal{U}^*_+, \mathcal{V}^*_+)$ is uniquely determined by variable $\tau$. For example, taking $\tau = 1.5$, $\Theta^* = 0.8626$, and the positive steady state $S_+ = (\mathcal{U}^*_+, \mathcal{V}^*_+)\approx(4.3128, 0.0531)$. However, when $\tau = n\in\mathbb{N}^+$, $\Theta^*$ belongs to the closed interval $\big[0.8+\frac{0.2}{\pi}\arctan(n-1), \, \, 0.8+\frac{0.2}{\pi}\arctan(n)\big]$, and then the positive steady state $S_+ = (\mathcal{U}^*_+, \mathcal{V}^*_+)$ is still not unique. For example, taking $\tau = 2$, $\Theta^*\in[0.8500, \, 0.8705]$, and there exist infinite positive steady states $S_+ = (\mathcal{U}^*_+, \mathcal{V}^*_+)$.

4.2.   Numerical simulation

(1) Figures 1–4 are some numerical simulations of Case 1. These simulations show that the equilibrium $S = (\mathcal{U}_0^*, \mathcal{V}_0^*) = (5, 0)$ of example (4.1) has global asymptotic stability. For any initial population density of moderate gamers and addicted gamers, after a long time of discontinuous treatment, the population density of the moderate gamers will be stable at $5$, and the addicted gamers will vanish.

(2) Figures 5–8 are some numerical simulations of Case 2 when $\tau = 1.5\in\mathbb{R}^+\setminus\mathbb{N}^+$. These simulations show that the positive steady state $S_+ = (\mathcal{U}^*_+, \mathcal{V}^*_+)\approx(4.3128, 0.0531)$ of example (4.1) has global asymptotic stability. For any initial population density of moderate gamers and addicted gamers, the population densities of the moderate gamers and addicted gamers will be stable at $S_+ = (4.3128, 0.0531)$.

(3) Figures 9–12 are some numerical simulations of Case 2 when $\tau = 2\in\mathbb{N}^+$. In this case, we use the random number generator of MATLAB to get $\Theta^* = 0.8605\in[0.8500, \, 0.8705]$. A theoretical calculation gives the positive steady state $S_+ = (\mathcal{U}^*_+, \mathcal{V}^*_+)\approx(4.3025, 0.0540)$ of example (4.1). For any initial population density of moderate gamers and addicted gamers, the population densities of the moderate gamers and addicted gamers will be stable at $S_+ = (4.3025, 0.0540)$.

5.   Conclusions

In the last decade, with the popularity of the internet, the number of internet users has continued to increase. While people enjoy the convenience and benefits brought by the internet, some disadvantages brought by the internet also begin to appear gradually. For example, online game addiction endangers the physical and mental health of players. In particular, many young addictive gamers are trapped in it. Many scholars, including mathematicians, began to pay attention to and study this phenomenon. Through the analysis of the dynamic change process of internet gamers, we put forward a new nonlinear reaction-diffusion model (1.2) of internet gamers having unsustainable treatment in this paper. Utilizing differential inclusion theory and stability theory, we study the existence and boundedness of solutions to the model as well as the global asymptotic stability of the steady states. Based on the MATLAB toolbox, an example is numerically simulated in detail. It is easy to see that the condition $rb < a\beta\gamma$ in Theorem 3.1 and the condition $rb\geq a\beta\gamma$ in Theorem 3.2 are completely opposite. For any initial population density $(\mathcal{U}_0(x), \mathcal{V}_0(x))$ of moderate gamers and addicted gamers, after a long time of discontinuous treatment, the population density of the moderate gamers will be stable at $\frac{r}{a}$, and the addicted gamers will vanish in the case of Theorem 3.1. Otherwise, the population densities of the moderate gamers and addicted gamers will be stable at $S_+$ in the case of Theorem 3.2. Inspired by some recent research works ^[26,27,28], we can further apply fractional differential equation theory to study the model (1.2) in the future.

Acknowledgments

The work was funded by research start-up funds for high-level talents of Taizhou University.

Conflict of interest

The author declares no competing interest.