Stability and Optimal Control Analysis for a Fractional-Order Industrial Virus Propagation Model Based on SCADA System

1. Introduction

The supervisory control and date acquisition (SCADA) system is widely used in industrial processes such as power system, the petroleum industry, medicine, metallurgy, chemical industry, transportation and other industries[1,2]. By the development of SCADA systems over time, modern SCADA systems have adopted standard network protocols that allow networked SCADA systems to be accessed via Ethernet. Unfortunately, this has brought about an increase in the threat of viruses in computer networks. It would bring an unimaginable disaster if SCADA systems are attacked by hackers since they are mainly used to monitor critical national infrastructures[3,4]. For example, the Stuxnet virus that spread across Iran in 2010 was able to target vulnerabilities in Microsoft systems and Siemens industrial systems, particularly Siemens SCADA systems. The Programmable Logical Controller (PLC) code can be modified by Stuxnet according to the characteristics of the target system, causing the centrifuge to accelerate abnormally until it is scrapped[5,6]. Therefore, there is an essential need to establish mathematical models to study the behavior of industrial viruses and to control their propagation. As noted in 1991 by Kephart and White, there are many similarities between biological virus propagation and computer virus propagation[7]. Since then, many scholars have also conducted in-depth research on computer viruses using biological models and analytical techniques[6,7,19,27,28].

Similar to biological population models, computer virus models usually exhibit an aspect of ecological memory: the ability of a community to respond in the present or future can be influenced by its past state or experiences[8]. This characteristic coincides with that of fractional calculus. Fractional differentiation can better express the biological memory of a system than integer differentiation[9]. The fractional-order model was initially proposed in the study of basic mathematical theory, but it has not been utilized extensively for a long time due to the difficulty of its calculation[10]. Thanks to the development of computer technology, researchers can do complex calculation using various software packages and fractional models re-entered the attention of many scholars[6,12,14,15,16,17,18,22,23,24,25,31,34,35,37,38]. The most used fractional operators are Gr$\ddot{u}$nwald-Letnikov, Riemann-Liouville and Caputo derivates. When talking about real problem, the Caputo derivative is highly useful since it allows traditional starting and boundary conditions be included in the derivation, and the derivative of a constant is zero, that is not the case with the Riemann-Liouville fractional derivative. The present work aims to analyze the corresponding model via the Caputo fractional derivative.

The optimal control tool is quite effective in the development of models, which refers to the search for control to optimize system performance under certain constrains[11]. These constraints are described by a dynamic system of differential equations. An integer or classic optimal control problem (OCPs) is the dynamic system governed by ordinary differential equations (ODEs) where academics have found some gaps and limitations in the formulation of these OCPs. As a generalization of classic OCPs, the fractional optimal control problems (FOCPs) appeared to address some limitations and gaps in the classical OCPs. The FOCPs consist of either the objective functional or the dynamic system or both include at least one fractional-order term[12,13]. In recent years, various optimization problems associated with FOCPs arising in biology, ecology and epidemiology were investigated by many authors[14,15,16,17,18]. In the present study, to use the limited resources to control the computer virus propagation, optimal control is introduced in the model.

In light of the discussion above, we consider a fractional industrial virus model of Caputo type. The model is reposed on an integer-order industrial virus model without memory effects established in Zhu et al[19]. One purpose in this work is to study dynamics of the proposed model. While the positivity and boundedness of solutions are established with the help of standard comparison results, the establishment of stability properties is not a simple and trivial task. As is known to all, the establishment of asymptotic stability of dynamics systems has an important role in both theory and practice but not a simple work in general. One of the most effective and powerful techniques to this problem is the Lyapunov stability theory[20,21]. However, it is not easy to construct appropriate Lyapunov functions for a given dynamical system. Another purpose in this work is to formulate the FOCP depending on the suggested fractional-order model. The proposed optimal control problem aims not only to control the spread of viruses but to achieve this economically. The objectives of this work are:

(i). The existence conditions and the locally asymptotic stability criterion are established for virus-free and virus-present equilibrium points in the proposed model.

(ii). By constructing a suitable Lyapunov function, we study the global stability of virus-free and virus-present equilibrium points of the system.

(iii). The model system is subjected to an optimal control analysis by incorporating two control efforts.

(iv). We employ the Adam-Bashforth-Moulton predictor-corrector technique to obtain numerical solution.

The rest of the paper is organized as follows: in Section 2, we present some basic properties of the Caputo fractional derivative operator and Mittag-Leffler function. In Section 3, we propose a fractional-order industrial virus model. In Section 4, we show the well-posedness of solutions. In Section 5, we study the stability of the virus-free and virus-present equilibrium points. In Section 6, we investigate FOCP. In Section 7, we illustrate the numerical simulations of the obtained results. The conclusion of this work is given in Section 8.

2. Preliminaries

In this section, we present some definitions, notions, and theorems that will be used in the next sections.

Definition 1. 

([9]) The Riemann-Liouville fractional integral of order $\alpha >0$ is given by

$$\left(Left\right){}_a{I}_t^{\alpha }f\left(t\right)={\displaystyle \frac{1}{\mathsf{\Gamma }\left(\alpha \right)}}{\int }_a^t{(t-\xi )}^{\alpha -1}f\left(\xi \right)d\xi ,$$

$$\left(Right\right){}_t{I}_b^{\alpha }f\left(t\right)={\displaystyle \frac{1}{\mathsf{\Gamma }\left(\alpha \right)}}{\int }_t^b{(\xi -t)}^{\alpha -1}f\left(\xi \right)d\xi .$$

where Γ (.) is the well-known gamma function.

Definition 2. 

([9]) The Riemann-Liouville fractional derivative of order $\alpha \in (n-1,n];n\in \mathbb{N},$ is defined by

$$\left(Left\right){}_a{D}_t^{\alpha }f\left(t\right)={\displaystyle \frac{d^n}{d{t}^n}}{(}_a{I}_t^{n-\alpha }f\left(t\right))={\displaystyle \frac{1}{\mathsf{\Gamma }(n-\alpha )}}{\left({\displaystyle \frac{d}{dt}}\right)}^n{\int }_a^t{(t-\xi )}^{n-\alpha -1}f\left(\xi \right)d\xi ,$$

$$\left(Right\right){}_t{D}_b^{\alpha }f\left(t\right)={\left(-{\displaystyle \frac{d}{dt}}\right)}^n{(}_t{I}_b^{n-\alpha }f\left(t\right))={\displaystyle \frac{1}{\mathsf{\Gamma }(n-\alpha )}}{\left(-{\displaystyle \frac{d}{dt}}\right)}^n{\int }_t^b(\xi -t)f\left(\xi \right)d\xi .$$

Definition 3. 

([9]) The Caputo fractional derivative of order $\alpha \in (n-1,n];n\in \mathbb{N}$ is defined as follows:

$$\left(Left\right){}_a^c{D}_t^{\alpha }f\left(t\right){=}_a{I}_t^{n-\alpha }{\displaystyle \frac{d^n}{d{t}^n}}f\left(t\right)={\displaystyle \frac{1}{\mathsf{\Gamma }(n-\alpha )}}{\int }_a^t{(t-\xi )}^{n-\alpha -1}{f}^{\left(n\right)}\left(\xi \right)d\xi ,$$

$$\left(Right\right){}_t^c{D}_b^{\alpha }f\left(t\right){=}_t{I}_b^{n-\alpha }{\left(-{\displaystyle \frac{d}{dt}}\right)}^{n}f\left(t\right)={\displaystyle \frac{{(-1)}^n}{\mathsf{\Gamma }(n-\alpha )}}{\int }_t^b{(\xi -t)}^{n-\alpha -1}{f}^{\left(n\right)}\left(\xi \right)d\xi ,$$

where $f(.)$ is a given function in the interval $[a,b]$.

Theorem 1. 

([9]) Suppose that $t>0$ and α∈$(n-1,n]$; n$\in \mathbb{N}.$ Then, we have

$${}_a{D}_t^{\alpha }f\left(t\right){=}_a^c{D}_t^{\alpha }f\left(t\right)+\sum _{j=0}^{n-1}{\displaystyle \frac{f^{\left(j\right)}\left(a\right)}{\mathsf{\Gamma }(j-\alpha +1)}}{(t-a)}^{(j-\alpha )},$$

$${}_t{D}_b^{\alpha }f\left(t\right){=}_t^c{D}_b^{\alpha }f\left(t\right)+\sum _{j=0}^{n-1}{\displaystyle \frac{f^{\left(j\right)}\left(b\right)}{\mathsf{\Gamma }(j-\alpha +1)}}{(b-t)}^{(j-\alpha )}.$$

Therefore,

$$iff\left(a\right)=f^{\prime }\left(a\right)=\dots =f^{(n-1)}\left(a\right)=0,then{}_a{D}_t^{\alpha }f\left(t\right){=}_a^c{D}_t^{\alpha }f\left(t\right),$$

$$iff\left(b\right)=f^{\prime }\left(b\right)=\dots =f^{(n-1)}\left(b\right)=0,then{}_t{D}_b^{\alpha }f\left(t\right){=}_t^c{D}_b^{\alpha }f\left(t\right).$$

Lemma 1. 

For $0<\alpha \le 1,t\in [a,b].$ Then, the left and right Caputo fractional derivative satisfies

$${}_t^c{D}_b^{\alpha }f\left(t\right){=}_a^c{D}_t^{\alpha }f(b-t).$$

The proof of this Lemma can be found in [18,22].

Definition 4.

Let $F\left(s\right)$ is the Laplace transform of the $f\left(t\right).$ Then,

$$\mathcal{L}\left\{{}^c{D}^{\alpha }f\left(t\right),s\right\}=s^{\alpha }F\left(s\right)-\sum _{i=0}^{n-1}{s}^{\alpha -i-1}{f}^{\left(i\right)}\left(0\right),\alpha \in (n-1,n];n\in \mathbb{N}.$$

Definition 5.

The Mittag-Leffler function $E_{l,m}\left(x\right)$ is given by

$$E_{l,m}\left(x\right)=\sum _{n=0}^{\infty }{\displaystyle \frac{x^n}{\mathsf{\Gamma }(ln+m)}},x\in \mathbb{R},l>0,m>0,$$

and satisfies:

(1) the property (see[23])

$$E_{l,m}\left(x\right)=x{E}_{l,l+m}\left(x\right)+{\displaystyle \frac{1}{\mathsf{\Gamma }\left(m\right)}},$$

(2) the Laplace transform of $t^{m-1}{E}_{l,m}(\pm \lambda t^l),$ as follows:

$$\mathcal{L}\left[t^{m-1}{E}_{l,m}(\pm \lambda t^l)\right]={\displaystyle \frac{s^{l-m}}{s^l\mp \lambda }}.$$

To prove the existence and uniqueness of the solution for fractional system, we need the following Lemma.

Lemma 2.([24]) Consider the system

$$\begin{array}{c}\hfill {}_{t_0}^c{D}_t^{q}x\left(t\right)=f(t,x),t>t_0\end{array}$$

(1)

with initial conditions $x_{t_0},$ where $q\in (0,1],f:[t_0,\infty )\times \Omega \to {\mathbb{R}}^n,\Omega \in {\mathbb{R}}^n,$ if $f(t,x)$ satisfies the locally lipschitz condition with respect to $x,$ then there exists a unique solution of (1) on $[t_0,\infty )\times \Omega .$

The following lemma proved in Vargas-De-León, which describes the Volterra-type Lyapunov function for the fractional order epidemic systems.

Lemma 3.([25]) Let $0<\alpha <1,$ and $\varsigma \in [0,T]$ be positive valued function. Then, for all $t\in [0,T)$, one has ${}_0^c{D}_t^{\alpha }\left(\varsigma \left(t\right)-{\varsigma }^{*}-{\varsigma }^{*}ln{\displaystyle \frac{\varsigma \left(t\right)}{{\varsigma }^{*}}}\right)\le \left(1-{\displaystyle \frac{{\varsigma }^{*}}{\varsigma \left(t\right)}}\right){}_0^c{D}_t^{\alpha }\varsigma \left(t\right),$ for all ${\varsigma }^{*}\in {\mathbb{R}}_{+}.$

In order to prove the global stability of all possible equilibrium problems, we now state the following Lemma.

Lemma 4.([26]) (Fractional LaSalle’s invariance principle) Let Π be a positively invariant subset of D where $D\subseteq {\mathbb{R}}^n,U:D\to \mathbb{R}$ be continuously differentiable function such that $U\left(x\right)>0$ and ${}^c{D}_t^{\alpha }U\left(x\left(t\right)\right)\le 0$ in Π for the solutions $x\left(t\right).$ Let ϝ be the set contains all points in Π where ${}^c{D}_t^{\alpha }U\left(x\left(t\right)\right)=0$ and Υ the largest invariant set in $ϝ,$ then every bounded solution starting in Π approaches Υ as $t\to \infty .$

3. Model Formulation

The SCADA system is mainly composed of communication devices, the remote terminal unit ( RTU ), the main terminal unit ( MTU ) and human-machine interface ( HMI ). The RTU is mainly responsible for collecting real-time data from the site and send it to the master station for analysis, and then the master station will send specific commands to the RTU according to the analysis results. In this process, a virus may make the RTU send data that contradicts the field data to the master station, thus causing the MTU to send wrong commands to the RTU. Because of the self-replication ability of the virus, an infected RTU may infect other vulnerable RTUs, eventually leading to unimaginable consequences [27]. To accurately describe the propagation characteristics of industrial viruses in SCADA systems, in [19] , the following integer-order model has been given to document the propagation of industrial viruses

$$\left\{\begin{array}{cc}\hfill & {\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}=b-\beta S(L+B)-\mu S,\hfill \\ \hfill & {\displaystyle \frac{\mathrm{d}L}{\mathrm{d}t}}=(1-p)\beta (L+B)S-\gamma L+\epsilon B-\mu \mathrm{L},\hfill \\ \hfill & {\displaystyle \frac{\mathrm{d}B}{\mathrm{d}t}}=p\beta (L+B)S+\gamma L-\epsilon B-aB-\mu B,\hfill \\ \hfill & {\displaystyle \frac{\mathrm{d}R}{\mathrm{d}t}}=aB-\mu R,\hfill \end{array}\right.$$

(2)

where $S\left(t\right),L\left(t\right),B\left(t\right)$ and $R\left(t\right)$ denote the susceptible compartments, latent compartments, breaking compartments and recovered compartments, respectively. This model is based on the following assumptions:

(A1) All newly connected RTUs to the SCADA system with a constant access rate of b.

(A2) Due to non-functional use, the RTU leaves the SCADA system with a probability of $\mu$.

(A3) The probability of each susceptible RTU being infected is $\beta$, and the probability of the susceptible node being infected at time t is $\beta S(L+B)$.

(A4) Each susceptible RTU either latent with probability $(1-p)\beta (L+B)$ or breaks out directly with probability $p\beta (L+B)$.

(A5) RTUs outbreak in the latent state with probability $\gamma$.

(A6) The virus may be removed with a probability of a due to installing and updating the latest version of antivirus software.

(A7) The spread of the virus will be suppressed with probability of inhibition $\epsilon$ when the virus outbreak reaches a certain level.

To the best of our knowledge, the model (2) has been studied and developed at some levels(see[28,29]and references therein), but its fractional-order versions have not been considered. This motives us to study the model (2) in the context of the Caputo fractional derivative. More precisely, we consider the following SLBR model

$$\left\{\begin{array}{cc}\hfill & {}_0^c{D}_t^{\alpha }S=b-\beta S(L+B)-\mu S,\hfill \\ \hfill & {}_0^c{D}_t^{\alpha }L=(1-p)\beta (L+B)S-\gamma L+\epsilon B-\mu L,\hfill \\ \hfill & {}_0^c{D}_t^{\alpha }B=p\beta (L+B)S+\gamma L-\epsilon B-aB-\mu B,\hfill \\ \hfill & {}_0^c{D}_t^{\alpha }R=aB-\mu R,\hfill \end{array}\right.$$

(3)

with initial conditions

$$S\left(0\right)=S^0⩾0,L\left(0\right)=L^0⩾0,B\left(0\right)=B^0⩾0,R\left(0\right)=R^0⩾0,$$

where ${}_0^c{D}_t^{\alpha }\omega \left(t\right)$ with $0<\alpha <1$ stands for the Caputo derivative of a given function $\omega \left(t\right)$. Apparently, model (3) degenerates to system (2) when $\alpha =1$.

4. Well-Posedness

In this section, we investigate the dynamics of the existence, uniqueness, non-negativity and uniform boundedness of the solutions for the proposed SLBR model (3).

4.1. Existence and uniqueness of solutions

In this subsection, we show the existence and uniqueness for the solutions of SLBR model (3) in the region:

$$\Omega =\left\{(S,L,B,R)\in R^4:max\left(\right|S|,|L|,|B|,|R\left|\right)\le M\right\}$$

Theorem 2.

There always exists a unique solution $X\left(t\right)\in \Omega$ for each initial condition $X\left(0\right)=(S^0,L^0,B^0,R^0)$, and $\forall t\ge 0$, in the SLBR model (3).

Proof. 

We use the approach in [24] to prove the existence and uniqueness. Consider $F\left(X\right)=(F_1\left(X\right),F_2\left(X\right),F_3\left(X\right),F_4\left(X\right))$, with

$$\begin{array}{ccc}\hfill F_1\left(X\right)& =& b-\beta S(L+B)-\mu S,\hfill \end{array}$$

(4)

$$\begin{array}{ccc}\hfill F_2\left(X\right)& =& (1-p)\beta (L+B)S-\gamma L+\epsilon B-\mu L,\hfill \end{array}$$

(5)

$$\begin{array}{ccc}\hfill F_3\left(X\right)& =& p\beta (L+B)S+\gamma L-\epsilon B-aB-\mu B,\hfill \end{array}$$

(6)

$$\begin{array}{ccc}\hfill F_4\left(X\right)& =& aB-\mu R.\hfill \end{array}$$

(7)

For any $X,\widehat{X}\in \Omega$, it follows from (4.1)-(4.4) that

$$\begin{array}{ccc}\hfill \parallel F\left(X\right)-F\left(\widehat{X}\right)\parallel & =& \left|F_1\left(X\right)-F_1\left(\widehat{X}\right)\right|+\left|F_2\left(X\right)-F_2\left(\widehat{X}\right)\right|+\left|F_3\left(X\right)-F_3\left(\widehat{X}\right)\right|\hfill \\ & & +\left|F_4\left(X\right)-F_4\left(\widehat{X}\right)\right|\hfill \\ & =& |b-\beta S(L+B)-\mu S-b+\beta \widehat{S}(\widehat{L}+\widehat{B})+\mu \widehat{S}|\hfill \\ & & +|(1-p)\beta (L+B)S-\gamma L+\epsilon B-\mu L-(1-p)\beta (\widehat{L}+\widehat{B})\widehat{S}+\gamma \widehat{L}\hfill \\ & & -\epsilon \widehat{B}+\mu \widehat{L}|+|p\beta (L+B)S+\gamma L-\epsilon B-aB-\mu B-p\beta (\widehat{L}+\widehat{B})\widehat{S}\hfill \\ & & -\gamma \widehat{L}+\epsilon \widehat{B}+a\widehat{B}+\mu \widehat{B}|+|aB-\mu R-a\widehat{B}+\mu \widehat{R}|\hfill \\ & \le & l_1|S-\widehat{S}|+l_2|L-\widehat{L}|+l_3|B-\widehat{B}|+l_4|R-\widehat{R}|\hfill \\ & \le & K\parallel (S,L,B,R)-(\widehat{S},\widehat{L},\widehat{B},\widehat{R})\parallel \hfill \\ & \le & K\parallel X-\widehat{X}\parallel ,\hfill \\ \hfill where,K& =& max(l_1,l_2,l_3,l_4),\hfill \\ \hfill l_1& =& 4\beta M+\mu ,\hfill \\ \hfill l_2& =& 2\beta M+\mu +2\gamma ,\hfill \\ \hfill l_3& =& 2\beta M+\mu +2\epsilon +2a,\hfill \\ \hfill l_4& =& \mu .\hfill \end{array}$$

Hence, $F\left(X\right)$ satisfies the Lipschitz condition with respect to X. Therefore, Lemma 2.2 confirms that there exists a unique solution $X\left(t\right)$ of system (3) with initial condition $X\left(0\right)=(S^0,L^0,B^0,R^0)$. And that is what we proved. □

4.2. Non-Negativity and Uniform Boundedness

Consider the set $\widehat{G}$=$\{X\left(t\right):X\left(t\right)=(S\left(t\right),L\left(t\right),B\left(t\right),R\left(t\right))\in R^4,X\left(t\right)\ge 0\}.$ To prove that each solution of SLBR model (3) is non-negative and belong to $\widehat{G}$, we need the following generalized mean value theorem [30] and corollary.

Lemma 5.(Generalized Mean Value Theorem) Let $f\left(x\right)\in C[a,b]$ and ${}_a^c{D}_x^{\alpha }f\left(x\right)\in C(a,b]$ for $0<\alpha \le 1$, then

$$f\left(x\right)=f\left(a\right)+{\displaystyle \frac{1}{\mathsf{\Gamma }\left(\alpha \right)}}\left({}_a^c{D}_x^{\alpha }f\right)\left(\xi \right){(x-a)}^{\alpha }$$

with $a\le \xi \le x,\forall x\in (a,b].$

Corollary 1.

Suppose that $f\left(x\right)\in C[a,b]$ and ${}_a^c{D}_x^{\alpha }f\left(x\right)\in C(a,b]$, for $0<\alpha \le 1$. If ${}_a^c{D}_x^{\alpha }f\left(x\right)\ge 0,\forall x\in (a,b)$, then $f\left(x\right)$ is non-decreasing for each $x\in [a,b]$. If ${}_a^c{D}_x^{\alpha }f\left(x\right)\le 0,\forall x\in (a,b)$, then $f\left(x\right)$ is non-increasing for each $x\in [a,b].$

Proof. 

This is clear from Lemma 4.1. □

Theorem 3.

The solution of SLBR model (3) is a positive invariant set and belongs to $\widehat{G}$.

Proof. 

In order to demonstrate our conclusion, we use the following steps to prove the non-negativity of the solution of SLBR model (3).

$$\begin{array}{ccc}\hfill {\left.{}_0^c{D}_t^{\alpha }S\right|}_{S=0}& =& b>0,\hfill \\ \hfill {\left.{}_0^c{D}_t^{\alpha }L\right|}_{L=0}& =& (1-p)\beta BS+\epsilon B,\hfill \\ \hfill {\left.{}_0^c{D}_t^{\alpha }B\right|}_{B=0}& =& p\beta LS+\gamma L,\hfill \\ \hfill {\left.{}_0^c{D}_t^{\alpha }R\right|}_{R=0}& =& aB.\hfill \end{array}$$

Using Corollary 4.1, we can get $S\left(t\right)$ is non-decreasing at the neighbourhood of time $t=t^{*}(>0)$ where $S\left(t^{*}\right)=0$ and $S\left(t\right)$ cannot cross the axis $S\left(t\right)=0$ . Hence, $S\left(t\right)\ge 0$ for all $t>0.$

Now,we claim that the solution of $L\left(t\right)$ starts from $\widehat{G}$ and remains non-negative. If not, then there exists a ${\tau }_c<\infty ,0\le t<{\tau }_c$ such that

$$\left\{\begin{array}{c}L\left(t\right)>0,\mathrm{for}0\le t<{\tau }_c,\hfill \\ L\left({\tau }_c\right)=0,\hfill \\ L\left({\tau }_c^{+}\right)<0.\hfill \end{array}\right.$$

If $B\left({\tau }_c\right)\ge 0$, then from ${\left.{}_0^c{D}_t^{\alpha }L\right|}_{L=0}=(1-p)\beta BS+\epsilon B$ we have ${\left.{}_0^c{D}_t^{\alpha }L\left(t\right)\right|}_{L\left({\tau }_c\right)=0}$=$(1-p)\beta B\left({\tau }_c\right)S+\epsilon B\left({\tau }_c\right)\ge 0$. From Corollary 4.1, it is evident that $L\left(t\right)$ is non-decreasing at the neighbourhood of $t={\tau }_c$ and which concludes $L\left({\tau }_c^{+}\right)=0.$ Hence, we arrive at a contradiction.

If $B\left({\tau }_c\right)<0$, then there exists a $\tau$ such that $0<\tau <{\tau }_c$.

$$\left\{\begin{array}{c}B\left(t\right)>0,\mathrm{for}0\le t<\tau ,\hfill \\ B\left(\tau \right)=0,\hfill \\ B\left({\tau }^{+}\right)<0.\hfill \end{array}\right.$$

Thanks to ${\left.{}_0^c{D}_t^{\alpha }B\right|}_{B=0}=p\beta LS+\gamma L,$ we have ${\left.{}_0^c{D}_t^{\alpha }B\left(t\right)\right|}_{B\left(\tau \right)=0}\ge 0$, which indicates $B\left({\tau }^{+}\right)\nless 0$ and it opposes our assumption. Therefore, we have $L\left(t\right)\ge 0,\forall t\in [0,\infty )$. Again from ${\left.{}_0^c{D}_t^{\alpha }B\right|}_{B=0}=p\beta LS+\gamma L$ we have ${\left.{}_0^c{D}_t^{\alpha }B\right|}_{B\left(t\right)=0}=p\beta LS+\gamma L\ge 0$, which means $B\left(t\right)$ is non-decreasing in the neighbourhood of time $t=t^{**}(>0)$ where $B(t^{**}=0)$ and $B\left(t\right)$ cannot cross the axis $B\left(t\right)=0$. Hence, $B\left(t\right)\ge 0,\forall t\in [0,\infty ).$ Similarly, from ${\left.{}_0^c{D}_t^{\alpha }R\right|}_{R=0}=aB,$ it is easy to see that $R\left(t\right)\ge 0,\forall t\in [0,\infty ).$ Therefore, $\widehat{G}$ is a positive invariant set for the proposed fractional order SLBR model (3). □

Now, we are going to prove that each solution of the SLBR model (3) is uniformly bounded.

Theorem 4.

Each solution of SLBR model (3) starting in $\widehat{G}$ is uniformly bounded.

Proof. 

Adding all equations in system (3), we get

$${}_0^c{D}_t^{\alpha }(S+L+B+R)=b-\mu (S+L+B+R).$$

Since $N\left(t\right)=S\left(t\right)+L\left(t\right)+B\left(t\right)+R\left(t\right)$, then ${}_0^c{D}_t^{\alpha }N\left(t\right)+\mu N\left(t\right)\le b$. By Lemma 3 in [31], it follows that

$$N\left(t\right)\le (N\left(0\right)-{\displaystyle \frac{b}{\mu }})E_{\alpha }\left(-\mu t^{\alpha }\right)+{\displaystyle \frac{b}{\mu }},$$

where ${\displaystyle E_{\alpha }\left(z\right)=\sum _{k=0}^{\infty }\frac{z^k}{\mathsf{\Gamma }(\alpha k+1)}}$ is the Mittag-Leffler function of parameter $\alpha$ [23]. Thus, ${\displaystyle \underset{t\to \infty }{lim\; sup}N\left(t\right)\le \frac{b}{\mu },}$ which implies that $S\left(t\right),L\left(t\right),B\left(t\right)$ and $R\left(t\right)$ are uniformly bounded. □

The last equation of recovered compartments $R\left(t\right)$ in system (3) is independent of other equations. With this in mind, we focus on the reduced system as follows:

$$\left\{\begin{array}{cc}\hfill & {}_0^c{D}_t^{\alpha }S=b-\beta S(L+B)-\mu S,\hfill \\ \hfill & {}_0^c{D}_t^{\alpha }L=(1-p)\beta (L+B)S-\gamma L+\epsilon B-\mu L,\hfill \\ \hfill & {}_0^c{D}_t^{\alpha }B=p\beta (L+B)S+\gamma L-\epsilon B-aB-\mu B,\hfill \end{array}\right.$$

(8)

with the following non-negative initial conditions:

$$S\left(0\right)=S^0⩾0,L\left(0\right)=L^0⩾0,B\left(0\right)=B^0⩾0.$$

5. Stability of the Equilibrium Points

The stability theory of the FDEs was introduced by Petras [32], which can be summarized as:

Lemma 6.

Consider the fractional order system

$$D^{\alpha }x\left(t\right)=h\left(x\left(t\right)\right),x\left(0\right)=x_0,$$

where $\alpha \in (0,1],x\left(t\right)\in {\mathbb{R}}^n$ and $h:{\mathbb{R}}^n\to {\mathbb{R}}^n$ is a $C^1$ function. The equilibrium points of the above system are solutions to equation $h\left(x\right)=0.$ An equilibrium point is locally asymptocially stable if all the eigenvalues ${\xi }_j(j=1,2,...,n)$ of the Jacobian matrix $J={\displaystyle \frac{\partial h}{\partial x}}$ evaluated at the equilibrium satisfy $|arg\left({\xi }_j\right)|>{\displaystyle \frac{\alpha \pi }{2}},$ and unstable if there exist an eigenvalue ${\xi }_j$ such that $|arg\left({\xi }_j\right)|<{\displaystyle \frac{\alpha \pi }{2}}.$

We now discuss the existence of equilibria of system (8). To evaluate equilibrium points, we set,

$$\left\{\begin{array}{c}{}_0^c{D}_t^{\alpha }S=0\hfill \\ {}_0^c{D}_t^{\alpha }L=0\hfill \\ {}_0^c{D}_t^{\alpha }B=0\hfill \end{array}\right.$$

we observe that system (8) has only two equilibrium points, a virus-free equilibrium point and a virus-present equilibrium point.

5.1. Equilibria and Basic Reproduction Number

The virus-free equilibrium point is given as follows: $E_1^{*}\left(S_0,0,0\right),$ where $S_0={\displaystyle \frac{b}{\mu }}.$

To consider the existence and uniqueness of the virus-present equilibrium $E_2^{*}=(S^{*},L^{*},B^{*}),$ we first calculate the basic reproduction number ${\mathcal{R}}_0^{\alpha }$ applying the next generation matrix approach[33]. The right-hand system (8) can be written as $\mathcal{F}-\mathcal{V}.$ Here, $\mathcal{F}$ calculates new infections in the model, and $\mathcal{V}$ calculates other transmissions in the infected parts. $\mathcal{F}$ and $\mathcal{V},$ associated with system (8), are given, respectively,by

$$\begin{array}{c}\hfill \mathcal{F}=\left[\begin{array}{c}(1-p)\beta (L+B)\\ p\beta (L+B)\end{array}\right],\mathcal{V}=\left[\begin{array}{c}\gamma L-\epsilon B+\mu L\\ -\gamma L+\epsilon B+aB+\mu B\end{array}\right].\end{array}$$

The Jacobian matrices of $\mathcal{F}$ and $\mathcal{V}$ at the virus-free equilibrium point $E_1^{*}$ are given by

$$\begin{array}{c}\hfill F=\left[\begin{array}{cc}{\displaystyle \frac{(1-p)\beta b}{\mu }}& {\displaystyle \frac{(1-p)\beta b}{\mu }}\\ {\displaystyle \frac{p\beta b}{\mu }}& {\displaystyle \frac{p\beta b}{\mu }}\end{array}\right],\widehat{V}=\left[\begin{array}{cc}\gamma +\mu & -\epsilon \\ -\gamma & \epsilon +a+\mu \end{array}\right].\end{array}$$

The basic reproduction ratio ${\mathcal{R}}_0^{\alpha },$ defined as the spectral radius of the matrix $F{\widehat{V}}^{-1}$, is obtained as

$${\mathcal{R}}_0^{\alpha }=\rho \left(F{\widehat{V}}^{-1}\right)={\displaystyle \frac{\beta S_0[\epsilon +\mu +\gamma +(1-p)a]}{(\mu +\gamma )(\mu +a)+\mu \epsilon }}.$$

In the case when $(L\ne 0)$ and $(B\ne 0)$, system (8) admits $E_2^{*}$ as a unique virus-present equilibrium point, where

$$S^{*}={\displaystyle \frac{(\gamma +\mu )(\mu +a)+\mu \epsilon }{\beta [\epsilon +\mu +\gamma +(1-p\left)a\right]}},$$

$$L^{*}={\displaystyle \frac{(b-\mu S^{*})[(1-p)\beta S^{*}+\epsilon ]}{\beta S^{*}(\epsilon +\gamma +\mu )}}={\displaystyle \frac{b(1-\frac{1}{{\mathcal{R}}_0^{\alpha }})[(\mu +a)(1-p)+\epsilon ]}{(\mu +\gamma )(a+\mu )+\epsilon \mu }},$$

$$B^{*}={\displaystyle \frac{(b-\mu S^{*})[\gamma +\mu -(1-p)\beta S^{*}]}{\beta S^{*}(\epsilon +\gamma +\mu )}}={\displaystyle \frac{(b-\mu S^{*})[\epsilon \gamma +{\gamma }^2+(1+p)\mu \gamma +p\mu \epsilon +p{\mu }^2]}{\beta S^{*}(\epsilon +\gamma +\mu )[\epsilon +\mu +\gamma +(1-p)a]}}.$$

Clearly, it is evident that if ${\mathcal{R}}_0^{\alpha }<1,$ then system (8) does not admit any positive virus-present equilibrium. Thus, we require ${\mathcal{R}}_0^{\alpha }>1,$ to assure the existence and positivity of the virus-present equilibrium point.

5.2. Local Stability

In this subsection, we investigate the local stability of the equilibria $E_1^{*}$ and $E_2^{*}$.

Theorem 5.

If ${\mathcal{R}}_0^{\alpha }<1,$ then the virus-free equilibrium point $E_1^{*}\left({\displaystyle \frac{b}{\mu }},0,0\right)$ is locally asymptotically stable. If ${\mathcal{R}}_0^{\alpha }>1,$ then $E_1^{*}\left({\displaystyle \frac{b}{\mu }},0,0\right)$ is unstable.

Proof. 

The following matrix gives the Jacobian matrix of system (8) at the virus-free equilibrium point $E_1^{*}\left({\displaystyle \frac{b}{\mu }},0,0\right)$

$$J\left(E_1^{*}\right)=\left[\begin{array}{ccc}-\mu & {\displaystyle \frac{\beta b}{\mu }}& {\displaystyle \frac{\beta b}{\mu }}\\ 0& {\displaystyle \frac{(1-p)\beta b}{\mu }}-(\gamma +\mu )& {\displaystyle \frac{(1-p)\beta b}{\mu }}+\epsilon \\ 0& {\displaystyle \frac{p\beta b}{\mu }}+\gamma & {\displaystyle \frac{(1-p)\beta b}{\mu }}-(\epsilon +a+\mu )\end{array}\right]$$

Obviously, $-\mu$ is one eigenvalue of $J\left(E_1^{*}\right)$ and the other two eigenvalues of $J\left(E_1^{*}\right)$ are given by the roots of the following equation

$$P\left(\lambda \right)={\lambda }^2+a_1\lambda +a_2,$$

where,

$$\begin{array}{cc}\hfill a_1& =-{\displaystyle \frac{\beta b}{\mu }}+\gamma +\epsilon +a+2\mu \hfill \\ \hfill & >\epsilon +a+2\mu +\gamma -{\displaystyle \frac{(\mu +\gamma )(\mu +a)+\mu \epsilon }{\epsilon +\mu +\gamma +(1-p)a}}\hfill \\ \hfill & ={\displaystyle \frac{(\epsilon +2\mu +\gamma +a)[(1-p)a+\epsilon ]+{\mu }^2+{\gamma }^2+\epsilon \gamma +2\mu \gamma }{\epsilon +\mu +\gamma +(1-p)a}}>0,\hfill \\ \hfill a_2& =[(\gamma +\mu )(\epsilon +a+\mu )-\epsilon \gamma ][1-{\mathcal{R}}_0^{\alpha }].\hfill \end{array}$$

If ${\mathcal{R}}_0^{\alpha }<1$, one immediately gets $a_1,a_2>0$. Thus, by the Routh-Hurwitz criterion, the eigenvalue ${\xi }_j(j=2,3)$ of $J\left(E_1^{*}\right)$ have negative real part if ${\mathcal{R}}_0^{\alpha }<1$, so that $|arg\left({\xi }_j\right)|>{\displaystyle \frac{\pi }{2}}>{\displaystyle \frac{\alpha \pi }{2}}$ for all $\alpha \in (0,1)$ if ${\mathcal{R}}_0^{\alpha }<1.$ If ${\mathcal{R}}_0^{\alpha }>1,$ then $a_2<0,$ and this suggests that $J\left(E_1^{*}\right)$ admits a positive real eigenvalue ${\xi }^{*}$, then $arg|\left({\xi }^{*}\right)|=0<{\displaystyle \frac{\alpha \pi }{2}}$ for all $\alpha \in (0,1)$ if ${\mathcal{R}}_0^{\alpha }>1.$ Consequently, by Lemma 5.1, one can obtain Theorem 5.1. □

Theorem 6.

If ${\mathcal{R}}_0^{\alpha }>1$, then the virus-present equilibrium point $E_2^{*}$ is locally asymptotically stable for all $\alpha \in (0,1)$.

Proof. 

Now, we investigate the local stability of the virus-present equilibrium of system (8) by assuming that ${\mathcal{R}}_0^{\alpha }>1.$ We compute and define the jacobian matrix at the virus-present equilibrium point as follows:

$$J\left(E_2^{*}\right)=\left[\begin{array}{ccc}-\beta (L^{*}+B^{*})-\mu & -\beta S^{*}& -\beta S^{*}\\ (1-p)\beta (L^{*}+B^{*})& (1-p)\beta S^{*}-(\gamma +\mu )& (1-p)\beta S^{*}+\epsilon \\ p\beta (L^{*}+B^{*})& p\beta S^{*}+\gamma & p\beta S^{*}-(\epsilon +a+\mu )\end{array}\right].$$

The characteristic equation of $J\left(E_2^{*}\right)$ is given by

$${\xi }^3+M_2{\xi }^2+M_1\xi +M_0=0,$$

where the coefficients utilizing $S^{*}={\displaystyle \frac{(\mu +\gamma )(\mu +a)+\mu \epsilon }{\mu +\epsilon +\gamma +(1-p)a}}$ are given by

$$\begin{array}{ccc}\hfill M_2& =& \epsilon +a+3\mu +\gamma +\beta (L^{*}+B^{*}-S^{*})\hfill \\ & >& \epsilon +a+3\mu +\gamma -\beta S^{*}\hfill \\ & =& \epsilon +a+3\mu +\gamma -{\displaystyle \frac{(\mu +\gamma )(\mu +a)+\mu \epsilon }{\mu +\epsilon +\gamma +(1-p)a}}\hfill \\ & =& {\displaystyle \frac{(\epsilon +a+3\mu +\gamma )[\mu +\epsilon +\gamma +(1-p)a]-({\mu }^2+a\mu +\gamma \mu +a\gamma +\mu \epsilon )}{\mu +\epsilon +\gamma +(1-p)a}}\hfill \\ & >& {\displaystyle \frac{(\epsilon +a+3\mu +\gamma )(\mu +\epsilon +\gamma )-({\mu }^2+a\mu +\gamma \mu +a\gamma +\mu \epsilon )}{\mu +\epsilon +\gamma +(1-p)a}}\hfill \\ & =& {\displaystyle \frac{{(\epsilon +\gamma )}^2+\mu (2\mu +3\epsilon +3\gamma )+a\epsilon }{\mu +\epsilon +\gamma +(1-p)a}}>0,\hfill \\ \hfill M_1& =& (\mu -\beta S^{*})(\epsilon +a+\mu +\gamma )+\mu +\beta (L^{*}+B^{*})+{\beta }^2{S}^{*}(L^{*}+B^{*})+a(p\beta S^{*}+\gamma )\hfill \\ & =& [(\epsilon +a+\mu +\gamma )\mu +a\gamma ]-\beta S^{*}[\mu +\epsilon +\gamma +(1-p)a]+\mu +\beta (L^{*}+B^{*})\hfill \\ & & +{\beta }^2{S}^{*}(L^{*}+B^{*})\hfill \\ & =& \mu +\beta (L^{*}+B^{*})+{\beta }^2{S}^{*}(L^{*}+B^{*})>0,\hfill \\ \hfill M_0& =& [{\mu }^2+\mu \beta (L^{*}+B^{*}-S^{*})](\epsilon +a+\mu +\gamma )+a[\mu p\beta S^{*}+\mu \gamma +\gamma \beta (L^{*}+B^{*})]\hfill \\ & =& \mu \left\{(\epsilon +a+\mu +\gamma )\mu +a\gamma -\left[(\epsilon +\mu +\gamma +(1-p)a)\right]\beta S^{*}\right\}+a\gamma \beta (L^{*}+B^{*})\hfill \\ & =& a\gamma \beta (L^{*}+B^{*})>0.\hfill \end{array}$$

Furtherly,

$$\begin{array}{ccc}\hfill M_1{M}_2-M_0& =& \left[a+\epsilon +3\mu +\gamma +\beta (L^{*}+B^{*})-\beta S^{*}\right]\left[\mu +(\beta +{\beta }^2{S}^{*})(L^{*}+B^{*})\right]\hfill \\ & & -a\gamma \beta (L^{*}+B^{*})\hfill \\ & =& \mu (a+\epsilon +3\mu +\gamma )+\mu \beta (L^{*}+B^{*})-\mu \beta S^{*}+(a+\epsilon +3\mu +\gamma )(\beta +{\beta }^2{S}^{*})\hfill \\ & & (L^{*}+B^{*})+\beta (\beta +{\beta }^2{S}^{*}){(L^{*}+B^{*})}^2-\beta S^{*}(\beta +{\beta }^2{S}^{*})(L^{*}+B^{*})\hfill \\ & & -a\gamma \beta (L^{*}+B^{*})\hfill \\ & =& \mu (a+\epsilon +3\mu +\gamma )+\mu \beta (L^{*}+B^{*}-S^{*})+\beta (\beta +{\beta }^2{S}^{*})(L^{*}+B^{*}-S^{*})\hfill \\ & & (L^{*}+B^{*})+(L^{*}+B^{*})\left[(a+\epsilon +3\mu +\gamma )(\beta +{\beta }^2{S}^{*})-a\gamma \beta \right]\hfill \\ & =& \mu (a+\epsilon +3\mu +\gamma )+\mu \beta (L^{*}+B^{*}-S^{*})+\beta (\beta +{\beta }^2{S}^{*})(L^{*}+B^{*}-S^{*})(L^{*}\hfill \\ & & +B^{*})+(L^{*}+B^{*})\left[(a+\epsilon +3\mu +\gamma )\left(\beta +\beta {\displaystyle \frac{(\mu +\gamma )(\mu +a)+\mu \epsilon }{\mu +\epsilon +\gamma +(1-p)a}}\right)-a\gamma \beta \right]\hfill \\ & =& \mu (a+\epsilon +3\mu +\gamma )+\mu \beta (L^{*}+B^{*}-S^{*})+\beta (\beta +{\beta }^2{S}^{*})(L^{*}+B^{*}-S^{*})\hfill \\ & & (L^{*}+B^{*})+(L^{*}+B^{*}){\displaystyle \frac{1}{\mu +\epsilon +\gamma +(1-p)a}}[a(a+1)\beta \mu +a\beta \epsilon \hfill \\ & & +(a+3\mu +\epsilon )(1-p)a\beta +(4a+\gamma +\epsilon +3)\beta {\mu }^2+(2a+3\mu +\epsilon +4)\beta \mu \epsilon \hfill \\ & & +\beta {\epsilon }^2+(4+3\mu +2\epsilon +\gamma )\mu \beta \gamma +3{\mu }^3\beta +(1+a+\gamma +4\mu )a\beta \gamma +2\gamma \beta \epsilon \hfill \\ & & +(1-a){\gamma }^2\beta +(1-a)\gamma \beta (1-p)a]\hfill \\ & =& \mu \left[a+\epsilon +3\mu +\gamma +\beta (L^{*}+B^{*}-S^{*})\right]+\beta (L^{*}+B^{*})\{(\beta +{\beta }^2{S}^{*})(L^{*}+B^{*}\hfill \\ & & -S^{*})+{\displaystyle \frac{1}{\mu +\epsilon +\gamma +(1-p)a}}[a(a+1)\mu +a\epsilon +(a+3\mu +\epsilon )(1-p)a\hfill \\ & & +(4a+\gamma +\epsilon +3){\mu }^2+(2a+3\mu +\epsilon +4)\mu \epsilon +{\epsilon }^2+(4+3\mu +2\epsilon +\gamma )\mu \gamma \hfill \\ & & +3{\mu }^3+(1+a+\gamma +4\mu )a\gamma +2\gamma \epsilon +(1-a){\gamma }^2+(1-a)(1-p)a\gamma \left]\right\}.\hfill \end{array}$$

Since $S^{*}={\displaystyle \frac{(\mu +\gamma )(\mu +a)+\mu \epsilon }{\beta [\mu +\epsilon +\gamma +(1-p\left)a\right]}},$ we deduce that

$$\begin{array}{ccc}\hfill M_1{M}_2-M_0& =& \mu \left[a+\epsilon +3\mu +\gamma +\beta (L^{*}+B^{*}-S^{*})\right]+\beta (L^{*}+B^{*})\{\beta (L^{*}+B^{*})+{\beta }^2{S}^{*}\hfill \\ & & (L^{*}+B^{*})-{\left[{\displaystyle \frac{(\mu +\gamma )(\mu +a)+\mu \epsilon }{\mu +\epsilon +\gamma +(1-p)a}}\right]}^2-{\displaystyle \frac{(\mu +\gamma )(\mu +a)+\mu \epsilon }{\mu +\epsilon +\gamma +(1-p)a}}\hfill \\ & & +{\displaystyle \frac{{\mu }^2+\mu \gamma +a\gamma +a\mu +\mu \epsilon }{\mu +\epsilon +\gamma +(1-p)a}}+{\displaystyle \frac{1}{\mu +\epsilon +\gamma +(1-p)a}}[a^2\mu +a\epsilon +(a+3\mu +\epsilon )\hfill \\ & & (1-p)a+(4a+\gamma +\epsilon +2){\mu }^2+(2a+3\mu +\epsilon +3)\mu \epsilon +{\epsilon }^2+(3+3\mu +2\epsilon \hfill \\ & & +\gamma )\mu \gamma +3{\mu }^3+(a+\gamma +4\mu )a\gamma +2\gamma \epsilon +(1-a){\gamma }^2+(1-a)(1-p)a\gamma \left]\right\}\hfill \\ & =& \mu \left[a+\epsilon +3\mu +\gamma +\beta (L^{*}+B^{*}-S^{*})\right]+\beta (L^{*}+B^{*})\{\beta (L^{*}+B^{*})+{\beta }^2{S}^{*}\hfill \\ & & (L^{*}+B^{*})-{\displaystyle \frac{{\mu }^4+2a{\mu }^3+{\mu }^2{a}^2+2{\mu }^3\gamma +4{\mu }^{2}a\gamma +2a^2\mu \gamma +{\gamma }^2{\mu }^2+2a\mu {\gamma }^2}{{[\mu +\epsilon +\gamma +(1-p)a]}^2}}\hfill \\ & & -{\displaystyle \frac{a^2{\gamma }^2+2{\mu }^2\epsilon \gamma +2a\mu \epsilon \gamma +{\mu }^2{\epsilon }^2+2{\mu }^3\epsilon +2{\mu }^2\epsilon a}{{[\mu +\epsilon +\gamma +(1-p)a]}^2}}+{\displaystyle \frac{1}{{[\mu +\epsilon +\gamma +(1-p)a]}^2}}\hfill \\ & & \left[\right({\mu }^3+2a{\mu }^2+\mu a^2+2{\mu }^2\gamma +4\mu a\gamma +a^2\gamma +\epsilon {\mu }^2+{\mu }^2\gamma +2a\mu \epsilon +{\mu }^2\epsilon )\hfill \\ & & (\mu +\epsilon +\gamma +(1-p)a)]+{\displaystyle \frac{1}{\mu +\epsilon +\gamma +(1-p)a}}[a\epsilon +(a+3\mu +\epsilon )(1-p)a\hfill \\ & & +(2a+2+2\epsilon ){\mu }^2+(\epsilon +3)\mu \epsilon +{\epsilon }^2+(3+2\epsilon +\gamma +\mu )\mu \gamma +2{\mu }^3+2\gamma \epsilon \hfill \\ & & +{\gamma }^2+(1-a)(1-p)a\gamma \left]\right\}\hfill \\ & >& \mu \left[a+\epsilon +3\mu +\gamma +\beta (L^{*}+B^{*}-S^{*})\right]+\beta (L^{*}+B^{*})\{\beta (L^{*}+B^{*})+{\beta }^2{S}^{*}\hfill \\ & & (L^{*}+B^{*})-{\displaystyle \frac{{\mu }^4+2a{\mu }^3+{\mu }^2{a}^2+2{\mu }^3\gamma +4{\mu }^{2}a\gamma +2a^2\mu \gamma +{\gamma }^2{\mu }^2+2a\mu {\gamma }^2}{{[\mu +\epsilon +\gamma +(1-p)a]}^2}}\hfill \\ & & -{\displaystyle \frac{a^2{\gamma }^2+2{\mu }^3\epsilon +2{\mu }^2\epsilon a+2{\mu }^2\epsilon \gamma +2a\mu \epsilon \gamma +{\mu }^2{\epsilon }^2}{{[\mu +\epsilon +\gamma +(1-p)a]}^2}}+{\displaystyle \frac{1}{{[\mu +\epsilon +\gamma +(1-p)a]}^2}}\hfill \\ & & ({\mu }^4+2a{\mu }^3+{\mu }^2{a}^2+2{\mu }^3\gamma +4{\mu }^{2}a\gamma +2a^2\mu \gamma +{\gamma }^2{\mu }^2+4a\mu {\gamma }^2+a^2{\gamma }^2+2{\mu }^3\epsilon \hfill \\ & & +2{\mu }^2\epsilon a+2{\mu }^2\epsilon \gamma +2a\mu \epsilon \gamma +{\mu }^2{\epsilon }^2)+{\displaystyle \frac{1}{\mu +\epsilon +\gamma +(1-p)a}}[a\epsilon +(a+3\mu +\epsilon )\hfill \\ & & (1-p)a+(2a+2+2\epsilon ){\mu }^2+(\epsilon +3)\mu \epsilon +{\epsilon }^2+(3+2\epsilon +\gamma +\mu )\mu \gamma \hfill \\ & & +2{\mu }^3+2\gamma \epsilon +{\gamma }^2+(1-a)(1-p)a\gamma \left]\right\}\hfill \\ & >& 0.\hfill \end{array}$$

Based on the above inequalities, we obtain $M_0>0,M_1>0,M_2>0$ and $M_1{M}_2>M_0$. Thus, according to the Routh-Hurwitz criterion, all roots ${\xi }_j(j=1,2,3)$ of system (8) have negative real part, so that $|arg\left({\xi }_j\right)|>{\displaystyle \frac{\pi }{2}}>{\displaystyle \frac{\alpha \pi }{2}}$ for all $\alpha \in (0,1)$ if ${\mathcal{R}}_0^{\alpha }>1.$ By Lemma 5.1, the virus-present equilibrium point $E_2^{*}$ is locally asymptotically stable. □

5.3. Global Stability

This subsection focus on the global stability analysis of two equilibria $E_1^{*}$ and $E_2^{*}$ based on the direct Lyapunov method, which aims to construct an appropriate Lyapunov functionals, and the fractional LaSalle’s invariance principle. Firstly, we demonstrate the following global stability result for virus-free equilibrium point $E_1^{*}.$

Theorem 7.

If ${\mathcal{R}}_0^{\alpha }<1,$ then the virus-free equilibrium point $E_1^{*}$ is globally asymptotically stable for all $\alpha \in (0,1).$

Proof. 

Let ${\mathcal{E}}_0$ be the Lyapunov functional defined as

$${\mathcal{E}}_0\left(t\right)=S\left(t\right)-S_0-S_{0}ln{\displaystyle \frac{S\left(t\right)}{S_0}}+{\displaystyle \frac{\epsilon +a+\mu +\gamma }{\gamma +\mu +\epsilon +(1-p)a}}L\left(t\right)+{\displaystyle \frac{\gamma +\mu +\epsilon }{\gamma +\mu +\epsilon +(1-p)a}}B\left(t\right).$$

By differentiating ${\mathcal{E}}_0$ with respect to time, we have

$${}_0^c{D}_t^{\alpha }{\mathcal{E}}_0=\left(1-{\displaystyle \frac{S_0}{S}}\right){}_0^c{D}_t^{\alpha }S+{\displaystyle \frac{\epsilon +a+\mu +\gamma }{\gamma +\mu +\epsilon +(1-p)a}}{}_0^c{D}_t^{\alpha }L+{\displaystyle \frac{\gamma +\mu +\epsilon }{\gamma +\mu +\epsilon +(1-p)a}}{}_0^c{D}_t^{\alpha }B.$$

Invoking system (8), we get

$$\begin{array}{cc}\hfill {}_0^c{D}_t^{\alpha }{\mathcal{E}}_0=& \left(1-{\displaystyle \frac{S_0}{S}}\right)\left[b-\beta S(L+B)-\mu S\right]+{\displaystyle \frac{\epsilon +a+\mu +\gamma }{\gamma +\mu +\epsilon +(1-p)a}}[(1-p)\beta (L+B)S-\gamma L\hfill \\ & +\epsilon B-\mu L]+{\displaystyle \frac{\gamma +\mu +\epsilon }{\gamma +\mu +\epsilon +(1-p)a}}\left[p\beta (L+B)S+\gamma L-\epsilon B-aB-\mu B\right],\hfill \end{array}$$

which implies

$$\begin{array}{cc}\hfill {}_0^c{D}_t^{\alpha }{\mathcal{E}}_0& =-{\displaystyle \frac{\mu }{S}}{(S-S_0)}^2-\beta S(L+B)+\beta S_0(L+B)+{\displaystyle \frac{(\epsilon +a+\mu +\gamma )(1-p)+(\gamma +\mu +\epsilon )p}{\gamma +\mu +\epsilon +(1-p)a}}\hfill \\ & \beta S(L+B)-{\displaystyle \frac{(\epsilon +a+\mu +\gamma )(\gamma +\mu )-(\gamma +\mu +\epsilon )\gamma }{\gamma +\mu +\epsilon +(1-p)a}}L-{\displaystyle \frac{(\gamma +\mu +\epsilon )(\epsilon +a+\mu )}{\gamma +\mu +\epsilon +(1-p)a}}B\hfill \\ & +{\displaystyle \frac{(\epsilon +a+\mu +\gamma )\epsilon }{\gamma +\mu +\epsilon +(1-p)a}}B\hfill \\ \hfill & =-{\displaystyle \frac{\mu }{S}}{(S-S_0)}^2-\beta S(L+B)+\beta S_0(L+B)+\beta S(L+B)-{\displaystyle \frac{(\mu +\gamma )(\mu +a)+\mu \epsilon }{\gamma +\mu +\epsilon +(1-p)a}}L\hfill \\ \hfill & -{\displaystyle \frac{(\mu +\gamma )(\mu +a)+\mu \epsilon }{\gamma +\mu +\epsilon +(1-p)a}}B\hfill \\ \hfill & =-{\displaystyle \frac{\mu }{S}}{(S-S_0)}^2+\beta S_0(L+B)-{\displaystyle \frac{\beta S_0}{{\mathcal{R}}_0^{\alpha }}}L-{\displaystyle \frac{\beta S_0}{{\mathcal{R}}_0^{\alpha }}}B\hfill \\ \hfill & =-{\displaystyle \frac{\mu }{S}}{(S-S_0)}^2+\beta S_{0}L\left(1-{\displaystyle \frac{1}{{\mathcal{R}}_0^{\alpha }}}\right)+\beta S_{0}B\left(1-{\displaystyle \frac{1}{{\mathcal{R}}_0^{\alpha }}}\right).\hfill \end{array}$$

Since ${\mathcal{R}}_0^{\alpha }<1,$ we have ${}_0^c{D}_t^{\alpha }{\mathcal{E}}_0\le 0,$ for all $t\ge 0.$ In addition, it is easy to verify that ${}_0^c{D}_t^{\alpha }{\mathcal{E}}_0=0$ if and only if $S\left(t\right)=S_0,L\left(t\right)=0$ and $B\left(t\right)=0.$ Hence, the largest compact invariant set in ${\mathsf{\Lambda }}_1=\left\{(S,L,B)\in {\mathbb{R}}_{+}^3:{}_0^c{D}_t^{\alpha }{\mathcal{E}}_0=0\right\}$ is the singleton $\left\{E_1^{*}\right\}.$ By Lemma 4.6 in [26], which generalized the integer order LaSalle’s invariance principle to fractional order system, we gain that $E_1^{*}$ is globally asymptotically stable if ${\mathcal{R}}_0^{\alpha }<1.$

Furthermore, under a mild condition on the parameters, we prove in the following theorems the global asymptotically stability of $E_2^{*}$ by using Lyapunov method[34]. □

Theorem 8.

If ${\mathcal{R}}_0^{\alpha }>1$ and $(1-p)\gamma =p\epsilon$, then the virus-present equilibrium $E_2^{*}\left(S^{*},L^{*},B^{*}\right)$ of system (8) is globally asymptotically stable.

Proof. 

Let the appropriate Lyapunov function $W:\mathcal{L}\to {\mathbb{R}}_{+}$ as:

$$\begin{array}{c}\hfill W=k_1\left(S-S^{*}-S^{*}ln{\displaystyle \frac{S}{S^{*}}}\right)+k_2\left(L-L^{*}-L^{*}ln{\displaystyle \frac{L}{L^{*}}}\right)+k_3\left(B-B^{*}-B^{*}ln{\displaystyle \frac{B}{B^{*}}}\right)\end{array}$$

where $k_i,(i=1,2,3)$ are positive constants specified later.

It is easy to see that $W\left(E_2^{*}\right)=0$. For $b=\beta S^{*}\left(L^{*}+B^{*}\right)+\mu S^{*}$, we have by the fractional derivative of W that

$$\begin{array}{cc}\hfill {}_0^c{D}_t^{\alpha }W& =k_1\left(1-{\displaystyle \frac{S^{*}}{S}}\right){}_0^c{D}_t^{\alpha }S+k_2\left(1-{\displaystyle \frac{L^{*}}{L}}\right){}_0^c{D}_t^{\alpha }L+k_3\left(1-{\displaystyle \frac{B^{*}}{B}}\right){}_0^c{D}_t^{\alpha }B\hfill \\ \hfill & =k_1\left(1-{\displaystyle \frac{S^{*}}{S}}\right)\left[\beta S^{*}(L^{*}+B^{*})+\mu S^{*}-\beta S(L+B)-\mu S\right]\hfill \\ \hfill & +k_2\left(1-{\displaystyle \frac{L^{*}}{L}}\right)\left[(1-p)\beta (L+B)S-(\gamma +\mu )L+\epsilon B\right]\hfill \\ \hfill & +k_3\left(1-{\displaystyle \frac{B^{*}}{B}}\right)\left[p\beta (L+B)S+\gamma L-(\epsilon +a+\mu )B\right]\hfill \\ \hfill & =k_1\mu S^{*}\left(2-{\displaystyle \frac{S}{S^{*}}}-{\displaystyle \frac{S^{*}}{S}}\right)+\left[k_1\beta (L^{*}+B^{*})S^{*}+k_2(\gamma +\mu )L^{*}+k_3(\epsilon +a+\mu )B^{*}\right]\hfill \\ \hfill & +S\beta (L+B)\left[k_2(1-p)+k_{3}p-k_1]+L[k_3\gamma +k_1\beta S^{*}-k_2(\gamma +\mu )\right]\hfill \\ \hfill & +B\left[k_1\beta S^{*}+k_2\epsilon -k_3(\epsilon +a+\mu )\right]-\left[k_1\beta (L^{*}+B^{*})S^{*}{\displaystyle \frac{S^{*}}{S}}+k_2(1-p)\beta (L+B)S{\displaystyle \frac{L^{*}}{L}}\right]\hfill \\ \hfill & -\left[k_2\epsilon B{\displaystyle \frac{L^{*}}{L}}+k_{3}p\beta (L+B)S{\displaystyle \frac{B^{*}}{B}}+k_3\gamma L{\displaystyle \frac{B^{*}}{B}}\right].\hfill \end{array}$$

Let us choose $k_1=1$ and $k_2,k_3$ in such a way that $k_i,i=1,2,3$ satisfy:

$$\left\{\begin{array}{cc}\hfill & -k_1+k_2(1-p)+k_{3}p=0,\hfill \\ \hfill & -k_2(\gamma +\mu )+k_3\gamma +k_1\beta S^{*}=0,\hfill \\ \hfill & k_1\beta S^{*}+k_2\epsilon -k_3(\epsilon +a+\mu )=0.\hfill \end{array}\right.$$

By solving above equations, we can obtain that

$$\begin{array}{cc}\hfill & k_3={\displaystyle \frac{\epsilon +\gamma +\mu }{(1-p)(\epsilon +a+\mu +\gamma )+p(\epsilon +\gamma +\mu )}},\hfill \\ \hfill & k_2={\displaystyle \frac{\epsilon +a+\gamma +\mu }{(1-p)(\epsilon +a+\gamma +\mu )+p(\epsilon +\gamma +\mu )}},\hfill \\ \hfill & k_1=1.\hfill \end{array}$$

Let us rearrange the terms of ${}_0^c{D}_t^{\alpha }W$ in such a way that ${}_0^c{D}_t^{\alpha }W=W_1+W_2+W_3,$ where

$W_1=k_1\mu S^{*}\left(2-{\displaystyle \frac{S}{S^{*}}}-{\displaystyle \frac{S^{*}}{S}}\right),W_2=k_1\beta \left(L^{*}+B^{*}\right)S^{*}+k_2\left(\gamma +\mu \right)L^{*}+k_3\left(\epsilon +a+\mu \right)B^{*},W_3=-\left[k_1\beta \left(L^{*}+B^{*}\right)S^{*}{\displaystyle \frac{S^{*}}{S}}+k_2(1-p)\beta \left(L+B\right)S{\displaystyle \frac{L^{*}}{L}}+k_2\epsilon B{\displaystyle \frac{L^{*}}{L}}+k_{3}p\beta \left(L+B\right)S{\displaystyle \frac{B^{*}}{B}}+k_3\gamma L{\displaystyle \frac{B^{*}}{B}}\right].$

Since the arithmetic mean is greater than or equal to the geometric mean, we gain that $W_1⩽0.$ Next, we are going to prove $W_2+W_3⩽0.$

Since $\gamma +\mu ={\displaystyle \frac{(1-p)\beta (L^{*}+B^{*})S^{*}+\epsilon B^{*}}{L^{*}}},$$\epsilon +a+\mu ={\displaystyle \frac{p\beta (L^{*}+B^{*})S^{*}+\gamma L^{*}}{B^{*}}},$ and $(1-p)\gamma =p\epsilon ,$ we conclude that

$$\begin{array}{ccc}\hfill W_2& =& 2k_2(1-p)\beta S^{*}{L}^{*}+2k_{3}p\beta S^{*}{B}^{*}+[k_2(1-p)\beta S^{*}{B}^{*}+k_{3}p\beta S^{*}{L}^{*}+k_2(1-p)\beta S^{*}{B}^{*}\hfill \\ & & +k_{3}p\beta S^{*}{L}^{*}]+(k_2\epsilon B^{*}+k_3\gamma L^{*})\hfill \\ & =& 2k_2(1-p)\beta S^{*}{L}^{*}+2k_{3}p\beta S^{*}{B}^{*}+4S^{*}\beta {\left[k_2{k}_{3}p(1-p)L^{*}{B}^{*}\right]}^{{\displaystyle \frac{1}{2}}}\hfill \\ & & +2{\left(k_2{k}_3\epsilon \gamma L^{*}{B}^{*}\right)}^{{\displaystyle \frac{1}{2}}}.\hfill \end{array}$$

(9)

$$\begin{array}{ccc}\hfill W_3& =& \left[-k_2(1-p)\beta {\displaystyle \frac{S^{{*}^2}{L}^{*}}{S}}-k_2(1-p)\beta S{L}^{*}\right]+\left[-k_{3}p\beta {\displaystyle \frac{S^{{*}^2}{B}^{*}}{S}}-k_{3}p\beta S{B}^{*}\right]+[-k_2\epsilon {\displaystyle \frac{B{L}^{*}}{L}}\hfill \\ & & -k_3\gamma {\displaystyle \frac{L{B}^{*}}{B}}]+\left[-k_2(1-p)\beta {\displaystyle \frac{S^{{*}^2}{B}^{*}}{S}}-k_{3}p\beta {\displaystyle \frac{S^{{*}^2}{L}^{*}}{S}}-k_2(1-p)\beta {\displaystyle \frac{S{L}^{*}B}{L}}-k_{3}p\beta {\displaystyle \frac{SL{B}^{*}}{B}}\right]\hfill \\ & ⩽& 2k_2(p-1)\beta S^{*}{L}^{*}-2k_{3}p\beta S^{*}{B}^{*}-2{\left(k_2{k}_3\epsilon \gamma L^{*}{B}^{*}\right)}^{{\displaystyle \frac{1}{2}}}-4S^{*}\beta {\left[k_2{k}_{3}p(1-p)L^{*}{B}^{*}\right]}^{{\displaystyle \frac{1}{2}}}.\hfill \end{array}$$

(10)

Collecting (6) and (10), we arrive at $W_2+W_3\le 0.$ Thanks to ${\mathcal{R}}_0^{\alpha }>1$ and $(1-p)\gamma =p\epsilon ,$ we can obtain that ${}_0^{c}D{\left(W\right)}_t^{\alpha }⩽0,$ for all $t\ge 0$. The equality holds only at the virus-present equilibrium point $(S^{*},L^{*},B^{*}).$ Moreover, the largest invariant set of ${\mathsf{\Lambda }}_2=\left\{(S,L,B)\in {\mathbb{R}}_{+}^3:{}_0^c{D}_t^{\alpha }W\left(t\right)=0\right\}$ is the singleton $\left\{E_2^{*}\right\}.$ By LaSalle’s invariance principle, the virus-present equilibrium point $E_2^{*}$ is globally asymptotical stability if ${\mathcal{R}}_0^{\alpha }>1$ and $(1-p)\gamma =p\epsilon .$ The proof is complete. □

Remark 1.

We would like to mention that the construction method of Lyapunov function in Theorem 5.4 is different from that Theorem 5 in Ref. [19]. The sufficient conditions for the global stability of the virus-present equilibrium point are unlike those in Ref. [19]. For model (3), we guess that the virus-present equilibrium point $E_2^{*}$ is also globally asymptotically stable as ${\mathcal{R}}_0^{\alpha }>1.$

6. FOCP Formulation

In this section, we formulate the FOCP to reduce the spread of infection, simultaneously reducing the related cost. To achieve our goal, we introduce two control measures: (i) improvement measures to raise public safety awareness $u_1\left(t\right)$; (ii) treatment measures that optimize the anti-virus software to kill the virus $u_2\left(t\right).$ Based on the control variables stated above the new FOCP of the Caputo fractional order model (3) is reformulated as

$$\left\{\begin{array}{c}{}_0^c{D}_t^{\alpha }S\left(t\right)=b-\left(1-u_1\left(t\right)\right)\beta S\left(t\right)(L\left(t\right)+B\left(t\right))-\mu S\left(t\right),\hfill \\ {}_0^c{D}_t^{\alpha }L\left(t\right)=(1-p)\left(1-u_1\left(t\right)\right)\beta (L\left(t\right)+B\left(t\right))S\left(t\right)-\gamma L\left(t\right)+\epsilon B\left(t\right)-\mu L\left(t\right),\hfill \\ {}_0^c{D}_t^{\alpha }B\left(t\right)=p(1-u_1\left(t\right))\beta (L\left(t\right)+B\left(t\right))S\left(t\right)+\gamma L\left(t\right)-\epsilon B\left(t\right)-u_2\left(t\right)B\left(t\right)-\mu B\left(t\right),\hfill \\ {}_0^c{D}_t^{\alpha }R\left(t\right)=u_2\left(t\right)B\left(t\right)-\mu R\left(t\right).\hfill \end{array}\right.$$

(11)

with initial data $S\left(0\right)=S^0,L\left(0\right)=L^0,B\left(0\right)=B^0,R\left(0\right)=R^0$ and the Lebesgue measurable control set is $\mathcal{U}=\left\{(u_1\left(t\right),u_2\left(t\right)):0\le u_i\left(t\right)\le u_{imax}\le 1,i=1,2,t\in [0,T]\right\},$ where T is the final time of implementing control measures. The objective of the control problem is to minimize the number of infected nodes under the cost of incorporating control strategies. In mathematical perspective, for a fixed terminal time T, the problem is to minimize the objective functional

$$\begin{array}{c}\hfill J(u_1,u_2)={\int }_0^T\left[AB\left(t\right)+{\displaystyle \frac{B_1}{2}}{u}_1^2\left(t\right)+{\displaystyle \frac{B_2}{2}}{u}_2^2\left(t\right)\right]dt,\end{array}$$

(12)

where the coefficients $A,B_1,B_2$ are weight constants of infected nodes and control measures. The following theorem follows from a direct application of Theorem 4.1 in Ref.[35].

Theorem 9.

Let the control function $u=(u_1,u_2)\in \mathcal{U}$ be measurable on $[0,T]$ with value of each of $u_1,u_2$ lies in $[0,1].$ Then there exist adjoint variables ${\lambda }_i,i=1,2,3,4$ satisfy

$\begin{array}{cc}\hfill {}_0^c{D}_t^{\alpha }{\lambda }_1(T-t)=& \left[-(1-u_1(T-t))\beta (L(T-t)+B(T-t))-\mu \right]{\lambda }_1(T-t)\hfill \\ & +\left[(1-p)(1-u_1(T-t))\beta (L(T-t)+B(T-t))\right]{\lambda }_2(T-t)\hfill \\ & +\left[p(1-u_1(T-t))\beta (L(T-t)+B(T-t))\right]{\lambda }_3(T-t),\hfill \\ \hfill {}_0^c{D}_t^{\alpha }{\lambda }_2(T-t)=& \left[-(1-u_1(T-t))\beta S(T-t)\right]{\lambda }_1(T-t)+[(1-p)(1-u_1(T-t))\beta S(T-t)\hfill \\ & -(\gamma +\mu )]{\lambda }_2(T-t)+\left[p(1-u_1(T-t))\beta S(T-t)+\gamma \right]{\lambda }_3(T-t),\hfill \\ \hfill {}_0^c{D}_t^{\alpha }{\lambda }_3(T-t)=& A-\left[(1-u_1(T-t))\beta S(T-t)\right]{\lambda }_1(T-t)+[(1-p)(1-u_1(T-t))\beta S(T-t)\hfill \\ & +\epsilon ]{\lambda }_2(T-t)+\left[p(1-u_1(T-t))\beta S(T-t)-\epsilon -u_2(T-t)-\mu \right]{\lambda }_3(T-t)\hfill \\ & +u_2(T-t){\lambda }_4(T-t),\hfill \\ \hfill {}_0^c{D}_t^{\alpha }{\lambda }_4(T-t)=& -\mu {\lambda }_4(T-t),\hfill \end{array}$ with the terminal conditions:

$S\left(0\right)=S^0,L\left(0\right)=L^0,B\left(0\right)=B^0,R\left(0\right)=R^0,$${\lambda }_1\left(T\right)={\lambda }_2\left(T\right)={\lambda }_3\left(T\right)={\lambda }_4\left(T\right)=0.$ Moreover, the optimal controls $u_1^{*}$ and $u_2^{*}$ are given as

$\begin{array}{cc}& u_1{}^{*}=min\left\{u_{1max},max\left(0,{\displaystyle \frac{\left[-{\lambda }_1\left(t\right)+{\lambda }_2\left(t\right)(1-p)+{\lambda }_3\left(t\right)p\right]\beta (L\left(t\right)+B\left(t\right))S\left(t\right)}{B_1}}\right)\right\},\hfill \\ & u_2{}^{*}=min\left\{u_{2max},max\left(0,{\displaystyle \frac{({\lambda }_3\left(t\right)-{\lambda }_4\left(t\right))B\left(t\right)}{B_2}}\right)\right\}.\hfill \end{array}$

7. Numerical Results and Discussion

In this section, the system (3) and system (11)-(12) are numerically simulated, respectively. All the simulation results are helped by MATLAB software. In the upcoming subsections, simulations for various parameters and different control strategies will be given comparatively.

7.1. Numerical Results of Fractional SLBR Model

In this subsection, in order to illustrate our analytical results, we adopt the numerical solver called the Adams-Bashforth-Moulton predictor corrector algorithm for getting approximate solutions of fractional ordinary differential equations [36,37].

7.1.1. Stability Analysis of Virus-Free Equilibrium

In what follows, we demonstrate numerically the analytical findings obtained in Section 5 by some graphical representation. In Figure 1-4, we plot the 2D system’s solutions using $a=0.5,b=2,\mu =0.3,\beta =0.003,p=0.5,\gamma =0.002,\alpha =0.9,\epsilon =0.5$ and using three sets of different initial values, which are $\left(S\right(0),L(0),B(0),R(0\left)\right)=(55,80,79,69),\left(S\right(0),L(0),B(0),R(0\left)\right)=(32,65,40,45)$ and $\left(S\right(0),L(0),B(0),R(0\left)\right)=$$(18,47,25,18),$ respectively. For this set of parameter values, ${\mathcal{R}}_0^{\alpha }\approx 0.5473<1$. We demonstrate the comparison of four state functions $S\left(t\right),L\left(t\right),B\left(t\right),R\left(t\right)$ under various initial conditions. It can be seen from Figure 3 that the value of all the curves eventually converges to 0, which verifies the stability of the virus-free equilibrium under different initial value conditions.

7.1.2. Stability Analysis of Virus-Present Equilibrium

Consider system (3) with $b=0.6,a=0.9,\mu =0.3,\beta =0.3,p=0.5,\gamma =0.002,\alpha =0.9,\epsilon =0.002.$ For this set of parameter values, $(1-p)\gamma =p\epsilon ,$ and ${\mathcal{R}}_0^{\alpha }\approx 1.2483>1$. Here we take three sets of different initial values, which are $\left(S\right(0),L(0),B(0),R(0\left)\right)=(50,300,200,150),\left(S\right(0),L(0),B(0),R(0\left)\right)=(100,425,300,195)$ and $\left(S\right(0),L(0),B(0),R(0\left)\right)$$=(200,500,310,210),$ respectively. As is shown in Figure 6-7, the number of L nodes and B nodes are gradually decreases and eventually keeps at a certain level. This observation is consistent with Theorem 8.

Figure 5. Variation of S over time.

Figure 5. Variation of S over time.

Figure 6. Variation of L over time.

Figure 6. Variation of L over time.

Figure 7. Variation of B over time.

Figure 7. Variation of B over time.

Figure 8. Variation of R over time.

Figure 8. Variation of R over time.

7.1.3. Model comparison

In this subsection, we present numerical comparison of several parameters in system (3). First, we show the effect of fractional order $\alpha$ on the system’s dynamics. Here, we use parameter values given as follows: $b=0.2,a=0.9,\mu =0.3,\beta =0.0002,p=0.5,\gamma =0.002,\epsilon =0.5.$ Figure 9-12 show the solution curves of system (3) with initial value $\left(S\right(0),L(0),B(0),R(0\left)\right)=(50,350,550,150),$ when $\alpha =0.90,\alpha =0.93$ and $\alpha =0.5.$ We can see from Figure 9-12 that fractional order $\alpha$ is a significant factor which affects the convergence speed of the solutions of system (3).

Figure 9. Variation of S over time.

Figure 9. Variation of S over time.

Figure 10. Variation of L over time.

Figure 10. Variation of L over time.

Figure 11. Variation of B over time.

Figure 11. Variation of B over time.

Figure 12. Variation of R over time.

Figure 12. Variation of R over time.

Next, we investigate the impact of parameter $\mu$ on the system’s dynamics. Here, we choose $b=2,a=0.9,\beta =0.0003,p=0.5,\gamma =0.002,\alpha =0.9,\epsilon =0.5.$ Figure 13-16 depict the solution curves of system (3) with initial value $\left(S\right(0),L(0),B(0),R(0\left)\right)=(75,350,550,150),$ when $\mu =0.3,\mu =0.25$ and $\mu =0.2.$ We observe that the higher the probability of the nodes leaving the SCADA system, the faster the number of L and B nodes decrease, which indicates that the nodes leaving the system is conducive to inhibiting the spread of the virus during the virus outbreak.

Figure 13. Variation of S over time.

Figure 13. Variation of S over time.

Figure 14. Variation of L over time.

Figure 14. Variation of L over time.

Figure 15. Variation of B over time.

Figure 15. Variation of B over time.

Figure 16. Variation of R over time.

Figure 16. Variation of R over time.

Finally, we explore the impact of parameter $\beta$ on the system’s dynamics. Here, we assume $b=0.2,a=0.9,\mu =0.3,p=0.5,\gamma =0.002,\alpha =0.9,\epsilon =0.5.$ Figure 17-21 describe the solution curves of system (3) with initial value $\left(S\right(0),L(0),B(0),R(0\left)\right)=(50,350,550,150),$ when $\beta =0.002,\beta =0.3$ and $\beta =0.13$. The results show that the smaller the value of $\beta$ is, the easier the virus propagation is inhibited. As is shown in Figure 17-21, when $\beta =0.002$, the number of S nodes increase significantly, while the number of L, B and R nodes decrease to varying degrees.

Figure 17. Variation of S over time.

Figure 17. Variation of S over time.

Figure 18. Variation of L over time.

Figure 18. Variation of L over time.

Figure 19. Variation of B over time.

Figure 19. Variation of B over time.

Figure 20. Variation of R over time.

Figure 20. Variation of R over time.

7.2. FOCP of SLBR transmission

In this subsection, we apply the fractional Euler method by combining it with the forward-backward predict-evaluate-correct-evaluate method [38]. To depict the effect of optimal control on the response of the system under study, we set $b=0.5,a=0.5,\mu =0.3,\beta =0.3,p=0.7,\gamma =0.002,\epsilon =0.5,A=100,B_1=0.6,B_2=0.4.$ We assume a final time $T=5$ for optimal control. Figure 21-24 show the time series plots of control signals that are applied in the model with initial value $\left(S\right(0),L(0),B(0),R(0\left)\right)=(53,350,325,480),$ when $\alpha =0.94,\alpha =0.92,\alpha =0.88,\alpha =0.98$ and $\alpha =0.90.$ The number of latent, breaking and recovered compartments is reduced in the case where the optimal control scheme is utilized. Figure 25-26 show the time series plots of two optimal control strategies $u_1$ and $u_2$ with different fractional order $\alpha =0.94,\alpha =0.92,\alpha =0.88,\alpha =0.98$ and $\alpha =0.90.$ We observe that the two control measures $u_1$ and $u_2$ should maintain maximum efforts for almost 0.5 year and 1 year respectively before decreasing to zero. The order of derivative can differ from range to range. If we vary the order of derivatives while keeping other parametric values fixed, the results are very close. This demonstrates that the spread of industrial virus is effectively controlled after raising public safety awareness and optimizing anti-virus software no matter in what range.

Figure 21. Variation of S over time.

Figure 21. Variation of S over time.

Figure 22. Variation of L over time.

Figure 22. Variation of L over time.

Figure 23. Variation of B over time.

Figure 23. Variation of B over time.

Figure 24. Variation of R over time.

Figure 24. Variation of R over time.

Figure 25. Variation of $u_1$ over time.

Figure 25. Variation of $u_1$ over time.

Figure 26. Variation of $u_2$ over time.

Figure 26. Variation of $u_2$ over time.

8. Conclusions

In this paper, we have studied a generalized mathematical model for the transmission dynamics of industrial virus based on SCADA system in the form of fractional-order model with a well-known Caputo fractional derivative. Theoretically, we have discussed the well-posedness of the given fractional-order model such as existence and uniqueness of solutions, non-negativity and uniform boundedness as well as stability of the equilibrium points. In addition, we have considered the control parameters of the given fractional-order model as time-dependent controls in order to formulate a FOCP constrained by the proposed fractional-order model, where improvement and treatment measures have been considered as two control efforts in the given FOCP. Also, conditions for fractional optimal control of industrial virus have been derived and analyzed. Numerically, graphical representations have been used to demonstrate the analytical outcomes, providing confidence in the results obtained. It follows from the research results that we should raise public safety awareness and optimize the anti-virus software for the aim of stopping industrial virus infection. The method employed in this paper can be easily adapted to solve other mathematical models related to the spread of infectious diseases such as plant diseases, human infectious diseases and so on.

Undoubtedly, dynamics of industrial virus is far more complicated and varied than one captured by the current mathematical model. In practical applications, the computer virus may spread once a user opens an infected file. Like human infectious disease, the computer virus usually takes a so-called incubation period to exhibit symptoms or have an impact on other individuals. As a future direction, we will incorporate time delay into the present model so as to obtain more dynamic behaviors, especially the conditions of Hopf bifurcation occurring in the system we propose.