1 Introduction

Many compartmental models like SIS, SIR, SIRS, SEIRS, MSEIR, MSEIRS have been introduced since 1927 in order to analyze the infectious diseases in detail by means of mathematical tools. When analyzed such models, it is supposed that the population is in a homogeneous structure, people react the same to infectious diseases and transform from one situation to another collectively. We should emphasize that this significant assumption is essential for formulating a complex system. Infectious diseases that cannot be fully defined have taken a crucial role in shaping world history while reducing the world population. For instance, in the second century AD, Antonine Plague gave rise to a sharp decline in the population and economic trouble in the Roman Empire and so this disease that caused the disruption of order led to the collapse of the Roman Empire. On the other hand, the Han empire also collapsed owing to similar reasons in the third century AD. The plague epidemic, which started in the Eastern Roman Empire in 541, affected the Mediterranean region adversely until 750 and caused an estimated 20-30 percent decrease in the population. With infectious diseases such as smallpox, measles, and diphtheria, 90 percent of the indigenous people died, and their population, which was 30 million in Mexico, decreased to 3 million between 1519–1530. Moreover, the plague called Black Death that ended up with great destruction between 1347 and 1351 destroyed about a third of the European population, which was 70 million in the 1340s and it repeated regularly more than 300 years in many parts of Europe. Also in the nineteenth century, cholera killed millions of people, and in 1918-1919 influenza killed more than 20 million people in the United States [1]. Between June 1918 and December 1920, Spanish Flu, the deadliest pandemic, caused the death of at least 50 million people, which constituted about 3 percent of the world population, and 500 million people were affected by this virus in almost every region of the world [2].

The increase in the effectiveness of sanitation programs, antibiotics, and vaccination programs in the 1960s led to the idea that infectious diseases will be overcome. Furthermore, infectious diseases adapt and develop in new environments and also new ecosystems, global warming, rise in international travel, changes in the economic structure allow novel infectious diseases to emerge easier. Resistance to antibiotics utilized in gonorrhea, tuberculosis, and pneumonia has developed and some diseases like dengue fever, yellow fever, malaria have become widespread especially in regions where climate change is experienced [3, 4]. In the last 40 years, various new infectious diseases such as Lyme (1975), legionary (1976), toxic shock syndrome (1978), hepatitis C (1989), hepatitis E (1990), hantavirus (1993), hepatitis G (1995), Nipah (1998), SARS (2003) have been diagnosed [5]. The Human Immunodeficiency Virus (HIV), the etiological agent (infectious substances) for Acquired Immune Deficiency Syndrome (AIDS), was first diagnosed in 1981 and surrounded the world. According to 2009 data, AIDS was a crucial infectious disease that caused about 25 million deaths and 33 million people had to live together with this dangerous disease and fight to death [6].

Numerous deaths due to plague, smallpox, cholera, and influenza outbreaks have led to the creation of various models for the monitoring and prevention of the spread of infectious diseases since the 17th century. Initially, it was aimed at monitoring only the simple data collection and the increase and decrease in the number of deaths, but today, the models are set using all the possibilities of information and communication technologies. In this context, system dynamics and agent-based simulation models, geographic information systems, spatial data mining, complex networks, heuristic methods, and dynamic system approaches play an important role. The basics of all these methods are based on mathematical models that seek solutions with differential equations laid in the 1920s. The first mathematical model related to infectious diseases belongs to Daniel Bernoulli published in 1766. In Bernoulli’s study, variolation discussions, a vaccination method developed for smallpox, are tackled by means of a mathematical model [7]. In 1906, Hamer designed and analyzed a discrete-time model to examine the recurrence of measles outbreaks. The model is important because it is the first model that assumes that the number of new cases that occur per unit time depends on the number of susceptible and disease-carrying individuals. After these studies, A. G. McKendrick and W. O. Kermack made the biggest contribution, which was shown as the first example in the prediction of the spread of infectious diseases by means of mathematical models and published in 1927 [8].

As an example of the SIS model, let us make mention of the spread of staph infection in hospitals. Methicillin-resistant Staphylococcus aureus (MRSA) is a bacterium that can cause serious infections in humans and is often referred to as staph. In fact, this bacterium is resistant to the typically used antibiotic methicillin. In hospitals, patients with compromised immune systems or elderly patients can easily catch the MRSA bacterium and get bloodstream infectious. MRSA, traditionally a crucial problem in hospitals, accounts for the significant part of hospital fatalities and brings about more deaths than AIDS each year [9]. Recently, a different genetic type of MRSA has been identified and it has been found that this novel strain CA-MRSA can also infect healthy and young individuals while traditional strain HA-MRSA has occasionally. In some studies, it is asserted that CA-MRSA may surpass HA-MRSA and so deaths may increase yearly owing to the severity of the problem handled. Therefore, a compartmental mathematical model has been developed to determine whether the CA-MRSA overtakes HA-MRSA. For more details on the above-mentioned model, we refer the reader to [9].

Fractional calculus, a venerable branch of mathematical analysis, examines the possibilities of the order of the derivative and integral operator, which is generally denoted by D and I, to be real or complex number. Although the idea of fractional derivative was first introduced when L’Hopital asked a question about fractional derivative to Leibniz in 1695, the definition of fractional derivative was not presented by Leibniz at that time. Hence, the fractional derivative definition has been put forward by subsequent scientists and today has many applications in biology, physics, chemistry, engineering, and so on. The topic of fractional calculus has attracted great attention from several authors due to applicability to various real-world problems and has the potential to achieve better results than the classical derivative [10, 11]. The most common way of explaining fractional operators is to present the Riemann-Liouville (RL) integral and derivative of an arbitrary complex analytic function. These substantial definitions are given as follows

$$\begin{aligned} {}_{RL}{\mathscr {I}}^\alpha f(t)= & {} \frac{1}{\Gamma (\alpha )}\int _{a}^{t} \frac{f(\tau )}{(t-\tau )^{1-\alpha }}d\tau ,\quad Re(\alpha )>0; \end{aligned}$$
(1)
$$\begin{aligned} {}_{RL}{\mathscr {D}}^\alpha f(t)= & {} \frac{d^n}{dx^n}_{RL}{\mathscr {I}}^{n-\alpha }f(t),\quad n=\lfloor Re(\alpha )\rfloor +1,\quad Re(\alpha )\ge 0. \end{aligned}$$
(2)

respectively. However, the definition of RL is not sufficient to fully explain fractional calculus because it has some shortcomings in practice related to initial conditions. So, various fractional derivatives with a wide variety of properties have been defined to compensate for these shortcomings. It should be emphasized that in fractional calculus, unlike classical analysis, there are dozens of differintegral definitions instead of a specific definition. This active and popular field of research has great importance for the modeling of certain biological, physical, etc processes. For more information on fractional operators, see [12,13,14,15,16,17,18,19,20,21,22,23] and the references given there. On the other hand, Caputo fractional derivative obtained by modifying the Riemann-Liouville definition is of paramount significance because of the applicability to various problems in daily life [24]. So, we present the definition of Caputo fractional operator as below:

$$\begin{aligned} _C{\mathscr {D}}^\alpha f(t)=\frac{1}{\Gamma (n-\alpha )}\int _{a}^{t}\frac{f^{(n)}(\tau )}{(t-\tau )^{\alpha -n+1}}d\tau , \end{aligned}$$
(3)

where \(Re(\alpha )\ge 0\) and \(n=\lfloor Re(\alpha )\rfloor +1\). From the applications point of view, this useful fractional derivative has been widely interpreted in connection with different areas such as finance, diffusion, control engineering, optics, non-Markovian physical processes, etc. One of the best ways to comprehend the fractional calculus is to consider the tautochrone problem in which non-integer order is needed to get the desired results. It can be predicted that in the future, a definition covering all definitions of fractional operators available in the literature may be introduced. To learn more about the content of this study, we refer the readers to [25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43].

In hospitals, the infection of MRSA, a gram-positive bacterium, has traditionally been a crucial problem causing many deaths each year. Hospital-acquired MRSA strain (HA-MRSA) is usually seen in elderly people or patients with a weak immune system while community-acquired MRSA strain can also infect young and healthy individuals. Thus in order to critically enunciate the actual characteristics of the situation, an expansive study needs to be carried out. It can be noted that, in the recent past, multifarious scientists have shown that the fractional calculus can more precisely explain natural phenomena than the integer calculus. In light of these facts, fractional calculus has become increasingly important and common in modeling realistic cases, particularly those with memory effects. We get inspired to investigate and evaluate the fractional variant of the governing model with an effective fractional operator called Caputo fractional derivative. To the best of our knowledge, this is the first time a fractional operator has been employed for the analysis and investigation of the underlying model. Consequently, several properties of the model such as existence and uniqueness of the solution, positivity, stability analysis, sensitivity analysis, basic reproduction number have been analyzed. Furthermore, by means of an effective method called Laplace-Adomian decomposition method (LADM), simulation and comparative analysis are carried out on the basis of the model being considered and the values of parameters to determine the efficiency and effectiveness of the current fractional variant of the model in the Caputo sense.

The remainder of this study is as follows: In Sect 2, model description and necessary prerequisites are given in both classical and fractional manner. We show the existence and uniqueness of solutions related to the Caputo fractional derivative and some other crucial theoretical information such as stability analysis, sensitivity analysis, and so on are furnished in Sect. 3. Moreover, we present the solution for the MRSA model using the method of LADM related to Caputo’s definition and plot various graphs to get comparative results in Sect. 4. Finally, we give some important concluding remarks in Sect. 5.

2 Model description and prerequisites

In this segment, we present the methicillin-resistant Staphylococcus aereus (MRSA) model with two strains, CA-MRSA and HA-MRSA, by making use of standard works on the current model in [9, 23, 25]. Hence, the population exposed to an infectious disease in the model is divided into three groups:

  • \(H(t)=\)patients colonized with HA-MRSA,

  • \(C(t)=\)patients colonized with CA-MRSA,

  • \(S(t)=\)susceptible patients who are not colonized with HA-MRSA or CA-MRSA.

The parameters of the model we studied are given as follows:

  • \({\mathscr {N}}=\) the total number of the individuals who are the patients in the hospital,

  • \(\Lambda =\) the daily entrance rate of patients to the hospital,

  • \(\delta _C=\) the daily exit rate of patients who are colonized with CA-MRSA through death or discharge,

  • \(\delta _H=\) the daily exit rate of patients who are colonized with HA-MRSA through death or discharge,

  • \(\delta _S=\) the daily exit rate of susceptible patients through death or discharge,

  • \(\beta _C=\) the daily transfer rate of CA-MRSA among patients,

  • \(\beta _H=\) the daily transfer rate of HA-MRSA among patients,

  • \(\gamma _C=\) the daily rate of undergoing decolonization measures of patients who are colonized with CA-MRSA,

  • \(\gamma _H=\) the daily rate of undergoing decolonization measures of patients who are colonized with HA-MRSA,

Susceptible patients denoted by S(t) can become contaminated when visited by the colonized healthcare team or as a result of failure to pay attention to hand hygiene, the healthcare team can become colonized by coming into contact with infected patients. The susceptible patient is the person who does not have enough immunity to prevent the infection that will occur with this factor when it is exposed to a certain pathogen. Also, C(t) and S(t) are infectious individuals who carry the disease and can transmit it to susceptible individuals. Although some people have a natural immunity to certain infectious diseases, individuals in this situation can be neglected in modeling as they will constitute a very small proportion of the population. In the model, the concept of population refers to the number of people living in a particular geographic area. The population N and the number of individuals susceptible to the disease S are considered as two different concepts in order to avoid confusion in other models in which the entire population is not sensitive to a particular infectious disease. In the model, it is assumed that the individuals in both compartments will be in contact with each other and the level of contagiousness will occur during each contact so that individuals susceptible to the disease can be infected by infectious individuals. The SIS model is one of the compartmental models in which individuals are not immune and are susceptible to illness with recovery. As patients become colonized or decolonized, they move between compartments S(t), C(t) and H(t). It should be noted that HA-MRSA occurs generally in elderly patients while CA-MRSA can also occur in healthy young people.

The transition between states is modeled by

$$\begin{aligned} \begin{aligned} \frac{dS(t)}{dt}&=\Lambda -\frac{\beta _H S(t)H(t)}{{\mathscr {N}}} -\frac{\beta _C S(t)C(t)}{{\mathscr {N}}}+{\gamma _H} H(t)+{\gamma _C} C(t)-{\delta _S} S(t), \\ \frac{dH(t)}{dt}&=\frac{{\beta _H}S(t)H(t)}{{\mathscr {N}}}-{\gamma _H}H(t)-{\delta _H}H(t), \\ \frac{dC(t)}{dt}&=\frac{\beta _C S(t)C(t)}{{\mathscr {N}}}-{\gamma _C}C(t)-{\delta _C}C(t). \end{aligned} \end{aligned}$$
(4)

Assuming the hospital is always full, \(\Lambda =\delta _S S(t)+\delta _H H(t)+\delta _C C(t)\) and \({\mathscr {N}}=C(t)+H(t)+S(t)\), then the system (4) can be also shown by the following system:

$$\begin{aligned} \begin{aligned} \frac{dH(t)}{dt}&=\frac{\beta _H}{{\mathscr {N}}}({\mathscr {N}}-C(t)-H(t)) -(\delta _H+\gamma _H)H(t), \\ \frac{dC(t)}{dt}&=\frac{\beta _C}{{\mathscr {N}}}({\mathscr {N}}-C(t)-H(t))C(t) -(\delta _C+\gamma _C)C(t). \end{aligned} \end{aligned}$$
(5)

Here, S(t) is determined by the equation \(S(t)={\mathscr {N}}-H(t)-C(t)\).

The essential objective is to investigate transmission dynamics in the hospital of MRSA strains: HA-MRSA and CA-MRSA. When it is assumed that patients can only be colonized with HA-MRSA or CA-MRSA, according to global results, competitive exclusion occurs between the two strains of MRSA [23, 25]. In this case, the strain having a larger reproduction number become endemic and the other one is extinguished because of competition. As a consequence, some novel studies present that patients can be colonized with both HA-MRSA and CA-MRSA, that is, it is allowed patients to be co-colonized with two strain. The extended version of the single-colonization model (4) to the co-colonization model has been presented in [25].

2.1 Mathematical formulation of the fractional-order model

In current years, fractional calculus has attracted great attention, which enables us to take derivative or integral of arbitrary order. The main advantage of the fractional differential equations systems is to ensure the memory-effect in the model. Taking into account all these benefits, we present a novel system of fractional differential equations to model the CA-MRSA and HA-MRSA . This model is given by Caputo fractional derivative as follows:

$$\begin{aligned} \begin{aligned} _C{\mathscr {D}}^\alpha S(t)&=\Lambda ^\alpha -\frac{\beta _H^\alpha S(t)H(t)}{{\mathscr {N}}} -\frac{\beta _C^\alpha S(t)C(t)}{{\mathscr {N}}}+{\gamma _H^\alpha } H(t) +{\gamma _C^\alpha } C(t)-{\delta _S^\alpha } S(t), \\ _C{\mathscr {D}}^\alpha H(t)&=\frac{{\beta _H^\alpha }S(t)H(t)}{{\mathscr {N}}} -{\gamma _H^\alpha }H(t)-{\delta _H^\alpha }H(t), \\ _C{\mathscr {D}}^\alpha C(t)&=\frac{\beta _C^\alpha S(t)C(t)}{{\mathscr {N}}}-{\gamma _C^\alpha }C(t)-{\delta _C^\alpha }C(t). \end{aligned} \end{aligned}$$
(6)

If we presume \(\Lambda =\delta _S S(t)+\delta _H H(t)+\delta _C C(t)\) and \({\mathscr {N}}=C(t)+H(t)+S(t)\), we can rewrite the (6) as below:

$$\begin{aligned} \begin{aligned} _C{\mathscr {D}}^\alpha H(t)&=\frac{\beta _H^\alpha }{{\mathscr {N}}}({\mathscr {N}} -C(t)-H(t))-(\delta _H^\alpha +\gamma _H^\alpha )H(t), \\ _C{\mathscr {D}}^\alpha C(t)&=\frac{\beta _C^\alpha }{{\mathscr {N}}}({\mathscr {N}}-C(t) -H(t))C(t)-(\delta _C^\alpha +\gamma _C^\alpha )C(t). \end{aligned} \end{aligned}$$
(7)

S(t) can be calculated by the \(S(t)={\mathscr {N}}-H(t)-C(t)\).

3 Mathematical analyses of the MRSA model

3.1 Existence and uniqueness of solutions via Caputo fractional derivative

Here, by utilizing the theory of fixed-point, we furnish the existence and uniqueness of the solution for the non-linear system under the account by means of Caputo fractional derivative. Let us assume that \({\mathscr {B}}(J)\) is a Banach space for the continuous real-valued functions defined on the interval \(J=[0,a]\) with sub norm and \({\mathscr {Q}}={\mathscr {B}}(J)\times {\mathscr {B}}(J)\times {\mathscr {B}}(J)\) with the norm \(\parallel (S,H,C)\parallel =\parallel S \parallel +\parallel H \parallel +\parallel C\parallel \), \(\parallel S \parallel =\sup _{t\in J}|S|\), \(\parallel H \parallel =\sup _{t\in J}|H|\), \(\parallel C \parallel =\sup _{t\in J}|C|\). If we apply Caputo fractional integral on equation (6), we have

$$\begin{aligned} \begin{aligned} S(t)-S(0)&={_C{\mathscr {D}}^\alpha }\left[ \Lambda ^\alpha -\frac{\beta _H^\alpha S(t)H(t)}{{\mathscr {N}}}-\frac{\beta _C^\alpha S(t)C(t)}{{\mathscr {N}}}+{\gamma _H^\alpha } H(t)+{\gamma _C^\alpha } C(t)-{\delta _S^\alpha } S(t)\right] , \\ H(t)-H(0)&={_C{\mathscr {D}}^\alpha }\left[ \frac{{\beta _H^\alpha }S(t)H(t)}{{\mathscr {N}}}-{\gamma _H^\alpha }H(t)-{\delta _H^\alpha }H(t)\right] , \\ C(t)-C(0)&={_C{\mathscr {D}}^\alpha }\left[ \frac{\beta _C^\alpha S(t)C(t)}{{\mathscr {N}}}-{\gamma _C^\alpha }C(t)-{\delta _C^\alpha }C(t)\right] . \end{aligned} \end{aligned}$$
(8)

Assuming that

$$\begin{aligned} \begin{aligned} {\mathbf {K}}_1&=\Lambda ^\alpha -\frac{\beta _H^\alpha S(t)H(t)}{{\mathscr {N}}}-\frac{\beta _C^\alpha S(t)C(t)}{{\mathscr {N}}}+{\gamma _H^\alpha } H(t)+{\gamma _C^\alpha } C(t)-{\delta _S^\alpha } S(t), \\ {\mathbf {K}}_2&=\frac{{\beta _H^\alpha }S(t)H(t)}{{\mathscr {N}}}-{\gamma _H^\alpha }H(t)-{\delta _H^\alpha }H(t), \\ {\mathbf {K}}_3&=\frac{\beta _C^\alpha S(t)C(t)}{{\mathscr {N}}}-{\gamma _C^\alpha }C(t)-{\delta _C^\alpha }C(t), \end{aligned} \end{aligned}$$
(9)

then the system (8) can be written by means of Caputo operator as follows

$$\begin{aligned} \begin{aligned} S(t)-S(0)&={\mathscr {M}}(\alpha )\int _0^t \frac{{\mathbf {K}}_1(\alpha ,\Theta ,S(\Theta ))}{(t-\Theta )^\alpha }d\Theta , \\ H(t)-H(0)&={\mathscr {M}}(\alpha )\int _0^t \frac{{\mathbf {K}}_2(\alpha ,\Theta ,H(\Theta ))}{(t-\Theta )^\alpha }d\Theta , \\ C(t)-C(0)&={\mathscr {M}}(\alpha )\int _0^t \frac{{\mathbf {K}}_3(\alpha ,\Theta ,C(\Theta ))}{(t-\Theta )^\alpha }d\Theta . \end{aligned} \end{aligned}$$
(10)

It should be emphasized that \({\mathbf {K}}_1(S,\Theta )\), \({\mathbf {K}}_2(H,\Theta )\) and \({\mathbf {K}}_3(C,\Theta )\) satisfy the Lipschitz condition if and only if S(t), H(t) and C(t) have an upper bound. Let S(t) and \(S^*(t)\) be couple functions, then we have

$$\begin{aligned} \parallel {\mathbf {K}}_1(\alpha ,t,S(t))-{\mathbf {K}}_1(\alpha ,t,S^*(t)) \parallel =\left| \left| \left( -\frac{\beta _H^\alpha H(t)}{{\mathscr {N}}} -\frac{\beta _C^\alpha C(t)}{{\mathscr {N}}}\right) (S(t)-S^*(t))\right| \right| .\nonumber \\ \end{aligned}$$
(11)

Supposing \(\nu _1 := \left| \left| -\frac{\beta _H^\alpha H(t)}{{\mathscr {N}}}-\frac{\beta _C^\alpha C(t)}{{\mathscr {N}}} \right| \right| \), we obtain

$$\begin{aligned} ||{\mathbf {K}}_1(\alpha ,t,S(t))-{\mathbf {K}}_1(\alpha ,t,S^*(t))||\le \nu _1||S(t)-S^*(t)||, \end{aligned}$$
(12)

and in a similar way, it can be obtained

$$\begin{aligned} \begin{aligned} ||{\mathbf {K}}_2(\alpha ,t,H(t))-{\mathbf {K}}_2(\alpha ,t,H^*(t))||&\le \nu _2||H(t)-H^*(t)|| , \\ ||{\mathbf {K}}_3(\alpha ,t,C(t))-{\mathbf {K}}_3(\alpha ,t,C^*(t))||&\le \nu _3||C(t)-C^*(t)||,\\ \end{aligned} \end{aligned}$$
(13)

where \(\nu _2=-\frac{\beta _H^\alpha S(t)}{{\mathscr {N}}}-\gamma _H^\alpha -\delta _H^\alpha \) and \(\nu _3=-\frac{\beta _C^\alpha S(t)}{{\mathscr {N}}}-\gamma _C^\alpha -\delta _C^\alpha \). So, this indicates that the Lipschitz condition is satisfied for \({\mathbf {K}}_1\), \({\mathbf {K}}_2\) and \({\mathbf {K}}_3\).

Recursively, (10) can be expressed as

$$\begin{aligned} \begin{aligned} S_n(t)&={\mathscr {M}}(\alpha )\int _0^t \frac{{\mathbf {K}}_1(\alpha ,\Theta ,S_{n-1} (\Theta ))}{(t-\Theta )^\alpha }d\Theta , \\ H_n(t)&={\mathscr {M}}(\alpha )\int _0^t \frac{{\mathbf {K}}_2(\alpha ,\Theta ,H_{n-1} (\Theta ))}{(t-\Theta )^\alpha }d\Theta , \\ C_n(t)&={\mathscr {M}}(\alpha )\int _0^t \frac{{\mathbf {K}}_3(\alpha ,\Theta ,C_{n-1}(\Theta ))}{(t-\Theta )^\alpha }d\Theta , \end{aligned} \end{aligned}$$
(14)

associated with the initial conditions \(S_0(t)=S(0)\), \(H_0(t)=H(0)\), \(C_0(t)=C(0)\). After subtracting the successive terms, we have

$$\begin{aligned} \begin{aligned} \Psi _{S,n}(t)=S_n(t)-S_{n-1}(t)&={\mathscr {M}}(\alpha )\int _0^t \frac{{\mathbf {K}}_1(\alpha ,\Theta ,S_{n-1}(\Theta ))-{\mathbf {K}}_1(\alpha ,\Theta ,S_{n-2}(\Theta ))}{(t-\Theta )^\alpha }d\Theta , \\ \Psi _{H,n}(t)=H_n(t)-H_{n-1}(t)&={\mathscr {M}}(\alpha )\int _0^t \frac{{\mathbf {K}}_2(\alpha ,\Theta ,H_{n-1}(\Theta ))-{\mathbf {K}}_2(\alpha ,\Theta ,H_{n-2}(\Theta ))}{(t-\Theta )^\alpha }d\Theta , \\ \Psi _{C,n}(t)=C_n(t)-C_{n-1}(t)&={\mathscr {M}}(\alpha )\int _0^t \frac{{\mathbf {K}}_3(\alpha ,\Theta ,C_{n-1}(\Theta ))-{\mathbf {K}}_3(\alpha ,\Theta ,C_{n-2}(\Theta ))}{(t-\Theta )^\alpha }d\Theta . \end{aligned}\nonumber \\ \end{aligned}$$
(15)

Considering as follows

$$\begin{aligned} \begin{aligned} S_n(t)&=\sum _{j=0}^{n}\Psi _{S,j}(t) , \\ H_n(t)&= \sum _{j=0}^{n}\Psi _{H,j}(t), \\ C_n(t)&=\sum _{j=0}^{n}\Psi _{C,j}(t), \end{aligned} \end{aligned}$$
(16)

and in addition, by employing the equations (12), (13) and taking into account \(\Psi _{S,n-1}(t)=S_{n-1}(t)-S_{n-2}(t)\), \(\Psi _{H,n-1}(t)=H_{n-1}(t)-H_{n-2}(t)\), \(\Psi _{C,n-1}(t)=C_{n-1}(t)-C_{n-2}(t)\), we can reach

$$\begin{aligned} \begin{aligned} ||\Psi _{S,n}(t)||&={\mathscr {M}}(\alpha )\nu _1\int _0^t\frac{||\Psi _{S,n-1}(\Theta )||}{(t-\Theta )^\alpha }d\Theta , \\ ||\Psi _{H,n}(t)||&={\mathscr {M}}(\alpha )\nu _2\int _0^t\frac{||\Psi _{H,n-1}(\Theta )||}{(t-\Theta )^\alpha }d\Theta , \\ ||\Psi _{C,n}(t)||&={\mathscr {M}}(\alpha )\nu _3\int _0^t\frac{||\Psi _{C,n-1}(\Theta )||}{(t-\Theta )^\alpha }d\Theta . \end{aligned} \end{aligned}$$
(17)

Now, let us prove the theorem below:

Theorem 1

The fractional community-acquired and hospital-acquired MRSA model (6) has a unique solution under the condition that

$$\begin{aligned} \frac{{\mathscr {M}}(\alpha )}{\alpha }r^\alpha \nu _i <1,\,\,\,\,\,i=1,2,3 \end{aligned}$$
(18)

when \(t\in [0,r]\).

Proof

As shown above, the functions S(t), H(t) and C(t) are bounded and \({\mathbf {K}}_1\), \({\mathbf {K}}_2\), \({\mathbf {K}}_3\) satisfy the Lipschitz condition. Hence, by means of the recursive principle and (17), we get

$$\begin{aligned} \begin{aligned} ||\Psi _{S,n}(t)||&\le ||S_0(t)|| \left( \frac{{\mathscr {M}}(\alpha )}{\alpha }r^\alpha \nu _1\right) ^n, \\ ||\Psi _{H,n}(t)||&\le ||H_0(t)|| \left( \frac{{\mathscr {M}}(\alpha )}{\alpha }r^\alpha \nu _2\right) ^n, \\ ||\Psi _{C,n}(t)||&\le ||C_0(t)|| \left( \frac{{\mathscr {M}}(\alpha )}{\alpha }r^\alpha \nu _3\right) ^n. \end{aligned} \end{aligned}$$
(19)

So, it can be clearly considered that \(||\Psi _{S,n}(t)||\rightarrow 0\), \(||\Psi _{H,n}(t)||\rightarrow 0\) and \(||\Psi _{C,n}(t)||\rightarrow 0\) when \(n\rightarrow \infty \). Additionally, by using the triangle inequality and the system (19) for any p, we attain

$$\begin{aligned} \begin{aligned} ||S_{n+p}(t)-S_n(t)||&\le \sum _{j=n+1}^{n+p}k_1^j=\frac{k_1^{n+1}-k_1^{n+p+1}}{1-k_1},\\ ||H_{n+p}(t)-H_n(t)||&\le \sum _{j=n+1}^{n+p}k_2^j=\frac{k_2^{n+1}-k_2^{n+p+1}}{1-k_2},\\ ||C_{n+p}(t)-C_n(t)||&\le \sum _{j=n+1}^{n+p}k_3^j=\frac{k_3^{n+1}-k_3^{n+p+1}}{1-k_3}, \end{aligned} \end{aligned}$$
(20)

where \(k_i=\frac{{\mathscr {M}}(\alpha )}{\alpha }r^\alpha \nu _i<1\). Accordingly, \(S_n\), \(H_n\), \(C_n\) are Cauchy sequences in \({\mathscr {B}}(J)\). For this reason, they are uniformly convergent. Through the limit theorem, it can be said that the limit of the sequences (14) is the unique solution of the system (6). This completes the proof. \(\square \)

3.2 Stability analysis and iterative solutions via Caputo fractional derivative

In this portion, iterative solutions are presented by making use of the Laplace transform of Caputo fractional derivative. On the other hand, stability criteria for the fractional MRSA model (6) is given with the help of fixed point theorem. Let \(({\mathscr {B}},||.||)\) be a Banach space and \({\mathscr {Q}}^*\) be a self-map of \({\mathscr {B}}\). Also, let us consider the recursive procedure in the form of the \(y_{n+1}=h({\mathscr {Q}}^*,y_n)\) and \({\mathscr {G}}({\mathscr {Q}}^*)\) be a fixed point set of non-empty \({\mathscr {Q}}^*\). Here, the sequence \(y_n\) converges to the point of \(q^* \in {\mathscr {G}}({\mathscr {Q}}^*)\). Moreover, we define \(||z_{n+1}^*-h({\mathscr {Q}}^*,z_n^*)||\) such that \(\{z_n^*\subseteq {\mathscr {B}}\}\). The iterative approach, \(y_{n+1}=h({\mathscr {Q}}^*,y_n)\) is \({\mathscr {Q}}^*\)-stable if \(\lim _{n\rightarrow \infty }c_n=0\), that is, \(\lim _{n\rightarrow \infty }c_n^*=p^*\). For the sequence \(z_n\) to be convergent, it must have an upper limit. If all the conditions mentioned above are met for \(y_{n+1}={\mathscr {Q}}^*\) where n is considered as Picard’s iteration as in [26], then the iteration is \({\mathscr {Q}}^*\)-stable. So, we can express the theorem below:

Theorem 2

Let \(({\mathscr {B}},||.||)\) be a Banach space and \({\mathscr {Q}}^*\) be a self-map on \({\mathscr {B}}\), then for all \(x,y\in {\mathscr {B}}\) we have

$$\begin{aligned} \Vert {\mathscr {Q}}^{*}_{x}- {\mathscr {Q}}^{*}_{y}\Vert \le K\Vert x-{\mathscr {Q}}^{*}_{x}\Vert +k\Vert x-y \Vert \end{aligned}$$
(21)

where \(K \ge 0\), \(0\le k <1\). Assuming \({\mathscr {Q}}^*\) is Picard \({\mathscr {Q}}^*\)-stable, the recursive formula can be presented as follows

$$\begin{aligned} \begin{aligned} S_{n+1}(t)=&S_{n}(t)+ {\mathscr {L}}^{-1}\left\{ \frac{1}{s^\alpha }{\mathscr {L}}\left\{ \Lambda ^\alpha - \frac{\beta _H^\alpha S_{n}(t) H_{n}(t)}{{\mathscr {N}}} - \frac{\beta _C^\alpha S_{n}(t) C_{n}(t)}{{\mathscr {N}}}\right. \right. \\&\left. \left. +\gamma _H^\alpha H_{n}(t)+\gamma _C^\alpha C_{n}(t)-\delta _S^\alpha S_{n}(t)\right\} \right\} , \\ H_{n+1}(t)=&H_{n}(t)+{\mathscr {L}}^{-1}\left\{ \frac{1}{s^\alpha }{\mathscr {L}}\left\{ \frac{{\beta _H^\alpha }S_n(t)H_n(t)}{{\mathscr {N}}}-{\gamma _H^\alpha }H_n(t)-{\delta _H^\alpha }H_n(t)\right\} \right\} , \\ C_{n+1}(t)=&C_{n}(t)+ {\mathscr {L}}^{-1}\left\{ \frac{1}{s^\alpha }{\mathscr {L}}\left\{ \frac{\beta _C^\alpha S_n(t)C_n(t)}{{\mathscr {N}}}-{\gamma _C^\alpha }C_n(t)-{\delta _C^\alpha }C_n(t)\right\} \right\} . \end{aligned}\nonumber \\ \end{aligned}$$
(22)

Theorem 3

Let \({\mathscr {F}}\) be a self-map, \(P>0\), \(Q>0\) and \(R>0\) are three different constants, then

$$\begin{aligned} \begin{aligned} {\mathscr {F}}[S_n(t)]=&S_{n+1}(t)=S_{n}(t)+ {\mathscr {L}}^{-1}\left\{ \frac{1}{s^\alpha }{\mathscr {L}}\left\{ \Lambda ^\alpha - \frac{\beta _H^\alpha S_{n}(t) H_{n}(t)}{{\mathscr {N}}}- \frac{\beta _C^\alpha S_{n}(t) C_{n}(t)}{{\mathscr {N}}}\right. \right. \\&\left. \left. +\gamma _H^\alpha H_{n}(t)+\gamma _C^\alpha C_{n}(t)-\delta _S^\alpha S_{n}(t)\right\} \right\} , \\ {\mathscr {F}}[H_n(t)]=&H_{n+1}(t)=H_{n}(t)+{\mathscr {L}}^{-1}\left\{ \frac{1}{s^\alpha }{\mathscr {L}}\left\{ \frac{{\beta _H^\alpha }S_n(t)H_n(t)}{{\mathscr {N}}}-{\gamma _H^\alpha }H_n(t)-{\delta _H^\alpha }H_n(t)\right\} \right\} , \\ {\mathscr {F}}[C_n(t)]=&C_{n+1}(t)=C_{n}(t)+ {\mathscr {L}}^{-1}\left\{ \frac{1}{s^\alpha }{\mathscr {L}}\left\{ \frac{\beta _C^\alpha S_n(t)C_n(t)}{{\mathscr {N}}}-{\gamma _C^\alpha }C_n(t)-{\delta _C^\alpha }C_n(t)\right\} \right\} , \end{aligned}\nonumber \\ \end{aligned}$$
(23)

which is \({\mathscr {F}}\)-stable in \(L^1(a,b)\) if the following conditions are satisfied

$$\begin{aligned} \left\{ \begin{aligned}&\left\{ 1-\beta _H^\alpha (P+Q)f_1(\rho )-\beta _C^\alpha (P+R)g_1(\rho )+(\gamma _H^\alpha +\gamma _C^\alpha +\delta _S^\alpha )k_1(\rho )\right\}<1,\\ {}&\left\{ 1-\beta _H^\alpha (P+Q)f_2(\rho )-(\gamma _H^\alpha +\delta _H^\alpha )k_2(\rho )\right\}<1, \\ {}&\left\{ 1-\beta _C^\alpha (P+R)g_3(\rho )-(\gamma _C^\alpha +\delta _C^\alpha )k_3(\rho )\right\} <1. \end{aligned} \right. \end{aligned}$$
(24)

Proof

It is clear that \({\mathscr {F}}\) is a fixed point. So, we determine the following iterations for all \((m,n)\in {\mathbb {N}}\times {\mathbb {N}}\).

$$\begin{aligned} \begin{aligned} {\mathscr {F}}(S_{n}(t))-{\mathscr {F}}(S_{m}(t)) =&S_{n}(t)-S_{m}(t)\\&+ {\mathscr {L}}^{-1}\left\{ \frac{1}{s^\alpha }{\mathscr {L}}\left\{ \Lambda ^\alpha - \frac{\beta _H^\alpha S_{n}(t) H_{n}(t)}{{\mathscr {N}}}- \frac{\beta _C^\alpha S_{n}(t) C_{n}(t)}{{\mathscr {N}}}\right. \right. \\&\left. \left. +\gamma _H^\alpha H_{n}(t)+\gamma _C^\alpha C_{n}(t)-\delta _S^\alpha S_{n}(t)\right\} \right\} \\&-{\mathscr {L}}^{-1}\left\{ \frac{1}{s^\alpha }{\mathscr {L}} \left\{ \Lambda ^\alpha - \frac{\beta _H^\alpha S_{m}(t) H_{m}(t)}{{\mathscr {N}}}- \frac{\beta _C^\alpha S_{m}(t) C_{m}(t)}{{\mathscr {N}}}\right. \right. \\&\left. \left. +\gamma _H^\alpha H_{m}(t)+\gamma _C^\alpha C_{m}(t)-\delta _S^\alpha S_{m}(t)\right\} \right\} ,\\ {\mathscr {F}}(H_{n}(t))-{\mathscr {F}}(H_{m}(t)) =&H_{n}(t)-H_{m}(t)\\&+{\mathscr {L}}^{-1}\left\{ \frac{1}{s^\alpha }{\mathscr {L}}\left\{ \frac{\beta _H^\alpha S_{n}(t) H_{n}(t)}{N}-\gamma _H^\alpha H_{n}(t)-\delta _H^\alpha H_{n}(t)\right\} \right\} \\&-{\mathscr {L}}^{-1}\left\{ \frac{1}{s^\alpha }{\mathscr {L}}\left\{ \frac{\beta _H^\alpha S_{m}(t) H_{m}(t)}{{\mathscr {N}}}-\gamma _H^\alpha H_{m}(t)-\delta _H^\alpha H_{m}(t)\right\} \right\} ,\\ {\mathscr {F}}(C_{(n)}(t))-{\mathscr {F}}(C_{(m)}(t)) =&C_{(n)}(t)-C_{(m)}(t)\\&+ {\mathscr {L}}^{-1}\left\{ \frac{1}{s^\alpha }{\mathscr {L}}\left\{ \frac{\beta _C^\alpha S_{n}(t) C_{n}(t)}{{\mathscr {N}}}-\gamma _C^\alpha C_{n}(t)-\delta _C^\alpha C_{n}(t)\right\} \right\} \\&-{\mathscr {L}}^{-1}\left\{ \frac{1}{s^\alpha }{\mathscr {L}}\left\{ \frac{\beta _C^\alpha S_{m}(t) C_{m}(t)}{{\mathscr {N}}}-\gamma _C C_{m}(t)-\delta _C C_{m}(t)\right\} \right\} . \end{aligned}\nonumber \\ \end{aligned}$$
(25)

Taking the norm of both sides of the first equation in (25), we have

$$\begin{aligned} \begin{aligned} \Vert {\mathscr {F}}(S_{n}(t))-{\mathscr {F}}(S_{m}(t))\Vert =&\Vert S_{n}(t)-S_{m}(t)\\&+{\mathscr {L}}^{-1}\left\{ \frac{1}{s^\alpha }{\mathscr {L}}\left\{ \Lambda ^\alpha - \frac{\beta _H^\alpha S_{n}(t) H_{n}(t)}{{\mathscr {N}}}- \frac{\beta _C^\alpha S_{n}(t) C_{n}(t)}{{\mathscr {N}}}\right. \right. \\&\left. \left. +\gamma _H^\alpha H_{n}(t)+\gamma _C^\alpha C_{n}(t)-\delta _S^\alpha S_{n}(t)\right\} \right\} \\&-{\mathscr {L}}^{-1}\left\{ \frac{1}{s^\alpha } {\mathscr {L}}\left\{ \Lambda ^\alpha - \frac{\beta _H^\alpha S_{m}(t) H_{m}(t)}{{\mathscr {N}}}- \frac{\beta _C^\alpha S_{m}(t) C_{m}(t)}{{\mathscr {N}}}\right. \right. \\&\left. \left. +\gamma _H^\alpha H_{m}(t)+\gamma _C^\alpha C_{m}(t)-\delta _S^\alpha S_{m}(t)\right\} \right\} \bigg \Vert , \end{aligned}\nonumber \\ \end{aligned}$$
(26)

and if we use the triangular inequality, we can write

$$\begin{aligned} \begin{aligned} \Vert {\mathscr {F}}(S_{n}(t))-{\mathscr {F}}(S_{m}(t))\Vert \le&\Vert S_{n}(t)-S_{m}(t)\Vert \\&+\bigg \Vert {\mathscr {L}}^{-1}\left\{ \frac{1}{s^\alpha }{\mathscr {L}}\left\{ \Lambda ^\alpha - \frac{\beta _H^\alpha S_{n}(t) H_{n}(t)}{{\mathscr {N}}}- \frac{\beta _C^\alpha S_{n}(t) C_{n}(t)}{{\mathscr {N}}}\right. \right. \\&\left. \left. +\gamma _H^\alpha H_{n}(t)+\gamma _C^\alpha C_{n}(t)-\delta _S^\alpha S_{n}(t)\right\} \right\} \\&-{\mathscr {L}}^{-1}\left\{ \frac{1}{s^\alpha }{\mathscr {L}}\left\{ \Lambda ^\alpha - \frac{\beta _H^\alpha S_{m}(t) H_{m}(t)}{{\mathscr {N}}}- \frac{\beta _C^\alpha S_{m}(t) C_{m}(t)}{{\mathscr {N}}}\right. \right. \\&\left. \left. +\gamma _H^\alpha H_{m}(t)+\gamma _C^\alpha C_{m}(t)-\delta _S^\alpha S_{m}(t)\right\} \right\} \bigg \Vert . \end{aligned}\nonumber \\ \end{aligned}$$
(27)

By some necessary simplifications, (27) takes the form of

$$\begin{aligned} \begin{aligned} \Vert {\mathscr {F}}(S_{n}(t))-{\mathscr {F}}(S_{m}(t))\Vert \le&\Vert S_{n}(t)-S_{m}(t)\Vert \\&+ {\mathscr {L}}^{-1}\bigg \{\frac{1}{s^\alpha }{\mathscr {L}}\bigg \{\Vert -\beta _H^\alpha \frac{S_{n}(t)}{{\mathscr {N}}}(H_{n}(t)-H_{m}(t))\Vert +\Vert \\&-\beta _H^\alpha \frac{H_{m}}{{\mathscr {N}}}(S_{n}(t)-S_{m}(t))\Vert + \Vert -\beta _C^\alpha \frac{S_{n}}{{\mathscr {N}}}(C_{n}(t)-C_{m}(t))\Vert \\&+\Vert -\beta _C^\alpha \frac{C_{m}}{{\mathscr {N}}}(S_{n}(t)-S_{m}(t)) \Vert +\Vert \gamma _H^\alpha (H_{n}(t)-H_{m}(t)) \Vert +\Vert \gamma _C^\alpha (C_{n}(t)\\&-C_{m}(t)) \Vert +\Vert -\delta _S^\alpha (S_{n}(t)-S_{m}(t)) \Vert \bigg \}\bigg \}. \end{aligned}\nonumber \\ \end{aligned}$$
(28)

Owing to the same behavior of functions inside the system handled, it can be assumed that

$$\begin{aligned} \begin{aligned}&\Vert H_{n}(t)-H_{m}(t) \Vert \cong \Vert S_{n}(t)-S_{m}(t) \Vert \\&\Vert C_{n}(t)-C_{m}(t) \Vert \cong \Vert S_{n}(t)-S_{m}(t) \Vert . \end{aligned} \end{aligned}$$
(29)

Substituting (29) into the relation (28), we reach

$$\begin{aligned} \begin{aligned}&\Vert {\mathscr {F}}(S_{n}(t))-{\mathscr {F}}(S_{m}(t))\Vert \le \Vert S_{n}(t)-S_{m}(t)\Vert \\&\qquad + {\mathscr {L}}^{-1}\bigg \{\frac{1}{s^\alpha }{\mathscr {L}}\bigg \{\Vert -\beta _H^\alpha \frac{S_{n}(t)}{{\mathscr {N}}}(S_{n}(t)-S_{m}(t))\Vert +\Vert -\beta _H^\alpha \frac{H_{m}(t)}{{\mathscr {N}}}(S_{n}(t)-S_{m}(t))\Vert + \\&\qquad \Vert -\beta _C^\alpha \frac{S_{n}(t)}{{\mathscr {N}}}(S_{n}(t)-S_{m}(t))\Vert +\Vert -\beta _C^\alpha \frac{C_{m}(t)}{{\mathscr {N}}}(S_{n}(t)-S_{m}(t)) \Vert \\&\qquad +\Vert \gamma _H^\alpha (H_{n}(t)-H_{m}(t)) \Vert +\Vert \gamma _C^\alpha (C_{n}(t)-C_{m}(t)) \Vert +\Vert -\delta _S^\alpha (S_{n}(t)-S_{m})(t) \Vert \bigg \}\bigg \}. \end{aligned}\nonumber \\ \end{aligned}$$
(30)

Because the sequences \(S_n(t)\), \(H_m(t)\) and \(C_m(t)\) are convergent and bounded, there exist three different constants \(P>0\), \(Q>0\) and \(R>0\) for all t. Hence, we have

$$\begin{aligned} \begin{aligned}&\Vert S_{n}(t) \Vert< P, \Vert H_{m}(t)\Vert< Q, \Vert C_{m}(t)\Vert < R, \ (m,n)\in {\mathbb {N}}\times {\mathbb {N}}. \end{aligned} \end{aligned}$$
(31)

By the relations (30) and (31), one can attain

$$\begin{aligned} \begin{aligned}&\Vert {\mathscr {F}}(S_{n}(t))-{\mathscr {F}}(S_{m}(t))\Vert \\&\quad \le \bigg [1-\beta _H^\alpha (P+Q) f_1(\rho )-\beta _C^\alpha (P+R)g_1(\rho ) +(\gamma _H^\alpha +\gamma _C^\alpha +\delta _S^\alpha )k_1(\rho )\bigg ] \Vert (S_{n}(t)-S_{m}(t))\Vert , \end{aligned}\nonumber \\ \end{aligned}$$
(32)

where \(f_1\), \(g_1\) and \(k_1\) are the functions obtained by the inverse Laplace transform in (30). In a similar manner, we reach

$$\begin{aligned} \Vert {\mathscr {F}}(H_n(t))-{\mathscr {F}}(H_m(t))\Vert \le [1-\beta _H^\alpha (P+Q)f_2(\rho )-(\gamma _H^\alpha +\delta _H^\alpha )k_2(\rho )]\Vert H_n(t)-H_m(t)\Vert ,\nonumber \\ \end{aligned}$$
(33)

and

$$\begin{aligned} \Vert {\mathscr {F}}(C_n(t))-{\mathscr {F}}(C_m(t))\Vert \le [1-\beta _C^\alpha (P+R)g_3(\rho )-(\gamma _C^\alpha +\delta _C^\alpha )k_3(\rho )]\Vert C_n(t)-C_m(t)\Vert ,\nonumber \\ \end{aligned}$$
(34)

where the condition (24) is satisfied. Therefore, it can be said that \({\mathscr {F}}\) has a fixed-point. In order to prove that \({\mathscr {F}}\) satisfies the conditions of Theorem 3.3, we presume that (33) and (34) hold and also the following expressions are valid

$$\begin{aligned} k=(0,0,0),\,\,\,\,\,K=\left\{ \begin{aligned}&\left\{ 1-\beta _H^\alpha (P+Q)f_1(\rho )-\beta _C^\alpha (P+R)g_1(\rho )+(\gamma _H^\alpha +\gamma _C^\alpha +\delta _S^\alpha )k_1(\rho )\right\}<1,\\&\left\{ 1-\beta _H^\alpha (P+Q)f_2(\rho )-(\gamma _H^\alpha +\delta _H^\alpha )k_2(\rho )\right\}<1,\\&\left\{ 1-\beta _C^\alpha (P+R)g_3(\rho )-(\gamma _C^\alpha +\delta _C^\alpha )k_3(\rho )\right\} <1. \end{aligned} \right. \nonumber \\ \end{aligned}$$
(35)

So, the desired result is obtained. \(\square \)

3.3 The positiveness of solutions for MRSA model

We wish to determine the invariant region and showing that all solutions of the fractional differential equations system (6) are positive for all \(t\ge 0\). The main purpose is to furnish the feasibility of the solutions for the model investigated looking at whether they enter the invariant region \(\Upsilon \). Benefiting from the Caputo fractional operator’s advantages, we assume that

$$\begin{aligned} \Upsilon =(S,H,C)\in {\mathbf {R}}_+^3, {\mathbf {R}}_+^3=\{n\in {\mathbf {R}}_+^3 : n\ge 0\} \end{aligned}$$
(36)

be any solution of the model (6) having non-negative initial conditions. Additionally, we have \(n=(S(t),H(t),C(t))^T.\)

We must also demonstrate that the vector field points to \({\mathbf {R}}_+^3\) upon each hyperplane which is bounded by the non-negative hyperoctant. We can write

$$\begin{aligned} \begin{aligned}&_C{\mathscr {D}}^\alpha S(t)=\Lambda ^\alpha +\gamma _H^\alpha H(t)+\gamma _C^\alpha C(t)\ge 0, \\&_C{\mathscr {D}}^\alpha H(t)=\frac{\beta _H^\alpha S(t)H(t)}{{\mathscr {N}}}\ge 0, \\&_C{\mathscr {D}}^\alpha C(t)=\frac{\beta _C^\alpha S(t)C(t)}{{\mathscr {N}}}\ge 0. \end{aligned} \end{aligned}$$
(37)

Hence, the convenient region can be stated as

$$\begin{aligned} \Upsilon =\{(S,H,C): S\ge 0, H\ge 0, C\ge 0, {\mathscr {N}}\le 1\}. \end{aligned}$$
(38)

Therefore, the underlying model is biologically appropriate and mathematically well-posed in the region \(\Upsilon \) when \(t>0\). Besides, this region is positively invariant, that is, solutions of the system (6) are positive for all t.

3.4 Disease-free equilibrium point (DFE)

To verify the existence of the equilibrium points, let us assume that \({\mathbf {E}}(S^*,H^*,C^*)\) be the equilibrium points of the model (6). The equilibrium points can be obtained by setting the right-hand side of the differential equations in (6) equal to zero. That is to say

$$\begin{aligned} \begin{aligned}&\Lambda ^\alpha -\frac{\beta _H^\alpha S(t)H(t)}{{\mathscr {N}}}-\frac{\beta _C^\alpha S(t)C(t)}{{\mathscr {N}}}+{\gamma _H^\alpha } H(t)+{\gamma _C^\alpha } C(t)-{\delta _S^\alpha } S(t)=0,\\&\frac{{\beta _H^\alpha }S(t)H(t)}{{\mathscr {N}}}-{\gamma _H^\alpha }H(t)-{\delta _H^\alpha }H(t)=0, \\&\frac{\beta _C^\alpha S(t)C(t)}{{\mathscr {N}}}-{\gamma _C^\alpha }C(t)-{\delta _C^\alpha }C(t)=0. \end{aligned} \end{aligned}$$
(39)

If the disease does not exist (\(H^*=0\), \(C^*=0\)), the system (39) reduces to \(\Lambda ^\alpha -\delta _S^\alpha S=0\). After solving, we have \(S^*=\frac{\Lambda ^\alpha }{\delta _S^\alpha }\). Thus, one can attain

$$\begin{aligned} {\mathbf {E}}^0=(S^*,H^*,C^*)=\left( \frac{\Lambda ^\alpha }{\delta _S^\alpha },0,0\right) . \end{aligned}$$
(40)

This represents the disease-free equilibrium point which means there is no infection [25].

3.5 Reproduction number

The basic reproduction number denoted by \({\mathscr {R}}_0=\rho (FV^{-1})\) where \(\rho (.)\) stands for the spectral radius of the matrix \(FV^{-1}\) can be obtained by the next-generation matrix approach [31]. The matrix F of transmission and matrix V of transformation for our fractional model (6) are presented by

$$\begin{aligned} F=\begin{bmatrix} \frac{\beta _H^\alpha S}{{\mathscr {N}}} &{} 0 \\ 0 &{} \frac{\beta _C^\alpha S}{{\mathscr {N}}} \ \end{bmatrix} \end{aligned}$$
(41)

and

$$\begin{aligned} V=\begin{bmatrix} \gamma _H^\alpha +\delta _H^\alpha &{} 0 \\ 0 &{} \gamma _C^\alpha +\delta _C^\alpha \ \end{bmatrix}. \end{aligned}$$
(42)

So, we can attain the reproduction ratio for both HA-MRSA and CA-MRSA respectively, as follows

$$\begin{aligned} {\mathscr {R}}_0^H=\frac{\Lambda ^\alpha \beta _H^\alpha }{\delta _S^\alpha {\mathscr {N}}(\gamma _H^\alpha +\delta _H^\alpha )},\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, {\mathscr {R}}_0^C=\frac{\Lambda ^\alpha \beta _C^\alpha }{\delta _S^\alpha {\mathscr {N}}(\gamma _C^\alpha +\delta _C^\alpha )}. \end{aligned}$$
(43)

It can be expressed that when both \({\mathscr {R}}_0^H<1\) and \({\mathscr {R}}_0^C<1\), the disease-free equilibrium \({\mathbf {E}}^0_{SC}=(C,H)=(0,0)\) for the single-colonization fractional model (6) is globally asymptotically stable. Moreover, if \(1<{\mathscr {R}}_0^C<{\mathscr {R}}_0^H\), then the boundary equilibrium \({\mathbf {E}}_H:(C,H)=\left( 0,{\mathscr {N}}\left( 1-\frac{1}{{\mathscr {R}}_0^H}\right) \right) \) is stable and \({\mathbf {E}}_C:(C,H)=\left( {\mathscr {N}}\left( 1-\frac{1}{{\mathscr {R}}_0^C}\right) ,0\right) \) is unstable. Similarly, if \(1<{\mathscr {R}}_0^H<{\mathscr {R}}_0^C\), \({\mathbf {E}}_C\) is stable and \({\mathbf {E}}_H\) is unstable [25].

Table 1 Sensitivity indices of the reproduction number \(R_0\) against mentioned parameters
Full size table

3.6 Sensitivity analysis

The sensitivity analysis of \(R_0\) has drawn a lot of attention in various scientific areas. As the parameters of a dynamical model are estimated, it is possible to have some uncertainty about their values employed to draw conclusions about the system studied. To decrease the spread of the infectious disease, it can be performed the sensitivity analysis by determining the parameters. Sensitivity analysis is an important part of the disease model analysis although computation of it can become exhaustive for complex dynamical systems. For this reason, it is very important to evaluate the effects of each parameter on the spread of the disease and therefore finds the parameters that have the most important effect on the reduction and spread of the outbreak. To this end, we carry out the sensitivity analysis with the help of the sensitivity index for the parameters of the underlying model. This technique helps to measure the most sensitive parameters inside the system for the reproduction number \(\textit{R}_0\). The following formula is used to calculate the sensitivity index of the reproduction number \(\textit{R}_0\) of the model (6) presented by the Caputo derivative.

$$\begin{aligned} S_\omega =\frac{\omega }{|R_0|}\times \frac{\partial R_0}{\partial \omega }. \end{aligned}$$
(44)

Three methods are normally used to calculate the sensitivity indices, (i) by direct differentiation, (ii) by a Latin hypercube sampling method (iii) by the linearizing system under consideration, and then solving the obtain set of linear algebraic equations. We will apply the direct differentiation method as it gives analytical expressions for the indices. The indices not only show us the influence of various aspects associated with the spreading of infectious disease but also give us important information regarding the comparative change between \(R_0\) and different parameter. Consequently, it helps in developing control strategies. Table 1 shows that the parameters \(\Lambda \) and \(\beta _H\) have a positive influence on the reproduction number \(R_0\), which describes that the growth or decay of these parameters say by 10 percent will increase or decrease the reproduction number by 9.99 percent and 9.99 percent, respectively. But on the other hand, the index for parameters \(\delta _S\), N, \(\gamma _H\), and \(\delta _H\) illustrates that increasing their values by 10 percent will decrease the values of reproduction number \(R_0\) by 10.0 percent, 9.99 percent, 4.16 percent, and 5.83 percent respectively. On the other hand, sensitivity analysis of different parameters is carried out in Figs. 1, 2, 3, 4, 5 and 6.

Fig. 1
figure 1

\({\mathscr {R}}_0\) versus sensitive parameters \(\beta _H\) and \({\mathscr {N}}\)

Full size image
Fig. 2
figure 2

\({\mathscr {R}}_0\) versus sensitive parameters \(\delta _H\) and \(\beta _H\)

Full size image
Fig. 3
figure 3

\({\mathscr {R}}_0\) versus sensitive parameters \(\delta _S\) and \(\beta _H\)

Full size image
Fig. 4
figure 4

\({\mathscr {R}}_0\) versus sensitive parameters \(\gamma _H\) and \(\delta _H\)

Full size image
Fig. 5
figure 5

\({\mathscr {R}}_0\) versus sensitive parameters \(\beta _H\) and \(\Lambda \)

Full size image
Fig. 6
figure 6

\({\mathscr {R}}_0\) versus sensitive parameters \(\delta _S\) and \({\mathscr {N}}\)

Full size image

4 Discussions

In this section, we present the solution for the system (6) employing the Laplace-Adomian decomposition method (LADM). The fractional-order shows the realistic decline behavior of infection of disease. Therefore, here, the advantage of the fractional model is shown by the method of LADM on the various graphs. Let us take into account the following fractional system

$$\begin{aligned} \begin{aligned} _C{\mathscr {D}}^\alpha S(t)&=\Lambda ^\alpha -\frac{\beta _H^\alpha S(t)H(t)}{{\mathscr {N}}} -\frac{\beta _C^\alpha S(t)C(t)}{{\mathscr {N}}}+{\gamma _H^\alpha } H(t) +{\gamma _C^\alpha } C(t)-{\delta _S^\alpha } S(t), \\ _C{\mathscr {D}}^\alpha H(t)&=\frac{{\beta _H^\alpha }S(t)H(t)}{{\mathscr {N}}} -{\gamma _H^\alpha }H(t)-{\delta _H^\alpha }H(t), \\ _C{\mathscr {D}}^\alpha C(t)&=\frac{\beta _C^\alpha S(t)C(t)}{{\mathscr {N}}}-{\gamma _C^\alpha }C(t)-{\delta _C^\alpha }C(t). \end{aligned} \end{aligned}$$
(45)

with the initial conditions \(S(0) = n_1 ,~ H(0) = n_2 ,~ C(0) = n_3\). Applying Laplace transform to both sides of the system (45), we have

$$\begin{aligned} \begin{array}{rcl} {\mathscr {L}}\{_C{\mathscr {D}}^\alpha S(t)\} &{}=&{} {\mathscr {L}}\{\Lambda ^\alpha -\frac{\beta _H^\alpha S(t)H(t)}{{\mathscr {N}}}-\frac{\beta _C^\alpha S(t)C(t)}{{\mathscr {N}}}+{\gamma _H^\alpha } H(t)+{\gamma _C^\alpha } C(t)-{\delta _S^\alpha } S(t)\}, \\ \\ {\mathscr {L}}\{_C{\mathscr {D}}^\alpha H(t)\} &{}=&{} {\mathscr {L}}\{\frac{{\beta _H^\alpha }S(t)H(t)}{{\mathscr {N}}}-{\gamma _H^\alpha }H(t)-{\delta _H^\alpha }H(t)\}, \\ \\ {\mathscr {L}}\{_C{\mathscr {D}}^\alpha C(t)\} &{}=&{} {\mathscr {L}}\{\frac{\beta _C^\alpha S(t)C(t)}{{\mathscr {N}}}-{\gamma _C^\alpha }C(t)-{\delta _C^\alpha }C(t)\}. \end{array} \end{aligned}$$
(46)

By the definition of Caputo Laplace transform, it can be written

$$\begin{aligned} \begin{array}{rcl} s^\alpha {\mathscr {L}}\{S(t)\}-s^{\alpha -1} S(0) &{}=&{} {\mathscr {L}}\{\Lambda ^\alpha -\frac{\beta _H^\alpha S(t)H(t)}{{\mathscr {N}}}-\frac{\beta _C^\alpha S(t)C(t)}{{\mathscr {N}}}+{\gamma _H^\alpha } H(t)+{\gamma _C^\alpha } C(t)-{\delta _S^\alpha } S(t)\}, \\ \\ s^\alpha {\mathscr {L}}\{H(t)\}-s^{\alpha -1} H(0) &{}=&{} {\mathscr {L}}\{\frac{{\beta _H^\alpha }S(t)H(t)}{{\mathscr {N}}}-{\gamma _H^\alpha }H(t)-{\delta _H^\alpha }H(t)\}, \\ \\ s^\alpha {\mathscr {L}}\{C(t)\}-s^{\alpha -1} C(0) &{}=&{} {\mathscr {L}}\{\frac{\beta _C^\alpha S(t)C(t)}{{\mathscr {N}}}-{\gamma _C^\alpha }C(t)-{\delta _C^\alpha }C(t)\}, \end{array} \end{aligned}$$
(47)

and so we reach

$$\begin{aligned} \begin{array}{rcl} {\mathscr {L}} \{S(t)\} &{}=&{} \frac{n_1}{s} + \frac{1}{s^\alpha }{\mathscr {L}}\{\Lambda ^\alpha -\frac{\beta _H^\alpha S(t)H(t)}{{\mathscr {N}}}-\frac{\beta _C^\alpha S(t)C(t)}{{\mathscr {N}}}+{\gamma _H^\alpha } H(t)+{\gamma _C^\alpha } C(t)-{\delta _S^\alpha } S(t)\}, \\ \\ {\mathscr {L}}\{H(t)\} &{}=&{} \frac{n_2}{s} +\frac{1}{s^\alpha }{\mathscr {L}}\{\frac{{\beta _H^\alpha }S(t)H(t)}{{\mathscr {N}}}-{\gamma _H^\alpha }H(t)-{\delta _H^\alpha }H(t)\}, \\ \\ {\mathscr {L}}\{C(t)\} &{}=&{} \frac{n_3}{s} +\frac{1}{s^\alpha }{\mathscr {L}}\{\frac{\beta _C^\alpha S(t)C(t)}{{\mathscr {N}}}-{\gamma _C^\alpha }C(t)-{\delta _C^\alpha }C(t)\}. \end{array} \end{aligned}$$
(48)

Let us assume that the solution S(t), H(t) and C(t) is in the form of infinite series as follows

$$\begin{aligned} \begin{aligned}&S(t)=\sum _{n=0}^{\infty } S_{n}, ~~~ H(t) = \sum _{n=0}^{\infty } H_{n}, ~~~ C(t) = \sum _{n=0}^{\infty } C_{n}. \end{aligned} \end{aligned}$$
(49)

The non-linear terms SH and SC involved in the system (45) are decomposed by Adomian polynomial as below

$$\begin{aligned} S H = \sum _{n=0}^{\infty } A_n,~~~ S C = \sum _{n=0}^{\infty } B_n, \end{aligned}$$
(50)

where \(A_n\) and \(B_n\) are Adomian polynomials. They are given by

$$\begin{aligned} A_n = \frac{1}{\Gamma (n+1)} \frac{d^n}{d\lambda ^n} \bigg [\sum _{k=0}^{n} \lambda ^k S{}_k \sum _{k=0}^{n} \lambda ^k H{}_k \bigg ]\bigg |_{\lambda =0}, \nonumber \\ B_n = \frac{1}{\Gamma (n+1)} \frac{d^n}{d\lambda ^n} \bigg [\sum _{k=0}^{n} \lambda ^k S{}_k \sum _{k=0}^{n} \lambda ^k C{}_k \bigg ]\bigg |_{\lambda =0}. \end{aligned}$$
(51)

The first six polynomials are as follows

$$\begin{aligned} A_0= & {} S_0H_0, \nonumber \\ A_1= & {} S_0H_1+H_0S_1, \nonumber \\ A_2= & {} S_0H_2+S_1H_1+S_2H_0, \nonumber \\ A_3= & {} S_0H_3+S_1H_2+S_2H_1+S_3H_0, \nonumber \\ A_4= & {} S_0H_4+S_1H_3+S_2H_2+S_3H_1+S_4H_0, \nonumber \\ A_5= & {} S_0H_5+S_1H_4+S_2H_3+S_3H_2+S_4H_1+S_5H_0. \end{aligned}$$
(52)

and

$$\begin{aligned} B_0= & {} S_0C_0, \nonumber \\ B_1= & {} S_0C_1+C_0S_1, \nonumber \\ B_2= & {} S_0C_2+S_1C_1+S_2C_0, \nonumber \\ B_3= & {} S_0C_3+S_1C_2+S_2C_1+S_3C_0, \nonumber \\ B_4= & {} S_0C_4+S_1C_3+S_2C_2+S_3C_1+S_4C_0, \nonumber \\ B_5= & {} S_0C_5+S_1C_4+S_2C_3+S_3C_2+S_4C_1+S_5C_0. \end{aligned}$$
(53)

Utilizing the (52) and (53), one can attain

$$\begin{aligned}&\begin{array}{rcl} {\mathscr {L}} \{S_0\} = \frac{n_1}{s}, {\mathscr {L}}\{H_0\} = \frac{n_2}{s}, {\mathscr {L}}\{C_0\} = \frac{n_3}{s}, \end{array} \end{aligned}$$
(54)
$$\begin{aligned}&\begin{array}{rcl} {\mathscr {L}} \{S_1\} &{}=&{} (\Lambda ^\alpha - \frac{\beta _{H}^\alpha A_0}{{\mathscr {N}}} - \frac{\beta _{C}^\alpha B_0}{{\mathscr {N}}} + \gamma _{H}^\alpha H_0 + \gamma _{C}^\alpha C_0 - \delta _{S}^\alpha S_0)\frac{1}{s^{\alpha +1}}, \\ \\ {\mathscr {L}}\{H_1\} &{}=&{} (\frac{\beta _{H}^\alpha A_0}{{\mathscr {N}}} -\gamma _{H}^\alpha H_0 - \delta _{H}^\alpha H_0)\frac{1}{s^{\alpha +1}}, \\ \\ {\mathscr {L}}\{C_1\} &{}=&{} (\frac{\beta _{C}^\alpha B_0}{{\mathscr {N}}} -\gamma _C^\alpha C_0 - \delta _C^\alpha C_0)\frac{1}{s^{\alpha +1}}, \end{array} \end{aligned}$$
(55)

and

$$\begin{aligned} \begin{array}{rcl} {\mathscr {L}} \{S_2\} &{}=&{} (- \frac{\beta _{H} A_1}{N} - \frac{\beta _{C} B_1}{N} + \gamma _{H}H_1 + \gamma _{C} C_1 - \delta _{S} S_1)\frac{1}{s^{\alpha +1}}, \\ \\ {\mathscr {L}}\{H_2\} &{}=&{} (\frac{\beta _{H} A_1}{N} -\gamma _{H}H_1 - \delta _{H} H_1)\frac{1}{s^{\alpha +1}}, \\ \\ {\mathscr {L}}\{C_2\} &{}=&{} (\frac{\beta _{C} B_1}{N} -\gamma _C C_1 - \delta _C C_1)\frac{1}{s^{\alpha +1}}, \end{array} \end{aligned}$$
(56)

by continuing in a similar manner, we get

$$\begin{aligned} \begin{array}{rcl} {\mathscr {L}} \{S_{n+1}\} &{}=&{} (- \frac{\beta _{H}^\alpha A_n}{{\mathscr {N}}} - \frac{\beta _{C}^\alpha B_n}{{\mathscr {N}}} + \gamma _{H}^\alpha H_n + \gamma _{C}^\alpha C_n - \delta _{S}^\alpha S_n)\frac{1}{s^{\alpha +1}}, \\ \\ {\mathscr {L}}\{H_{n+1}\} &{}=&{} (\frac{\beta _{H}^\alpha A_n}{{\mathscr {N}}} -\gamma _{H}^\alpha H_n - \delta _{H}^\alpha H_n)\frac{1}{s^{\alpha +1}}, \\ \\ {\mathscr {L}}\{C_{n+1}\} &{}=&{} (\frac{\beta _{C}^\alpha B_n}{{\mathscr {N}}} -\gamma _C^\alpha C_n - \delta _C^\alpha C_n)\frac{1}{s^{\alpha +1}}. \end{array} \end{aligned}$$
(57)

Now applying inverse Laplace transformation on system (54), (55) and (56) we get

$$\begin{aligned}&\begin{array}{rcl} S_0 = n_1, H_0 = n_2, C_0 = n_3, \end{array} \end{aligned}$$
(58)
$$\begin{aligned}&\begin{array}{rcl} S_1 &{}=&{} (\Lambda ^\alpha - \frac{\beta _{H}^\alpha n_1n_2}{{\mathscr {N}}} - \frac{\beta _{C}^\alpha n_1n_3}{{\mathscr {N}}} + \gamma _{H}^\alpha n_2 + \gamma _{C}^\alpha n_3 - \delta _{S}^\alpha n_1)\frac{t^{\alpha }}{\Gamma ({\alpha +1})}, \\ \\ H_1 &{}=&{} (\frac{\beta _{H}^\alpha n_1n_2}{{\mathscr {N}}} -\gamma _{H}^\alpha n_2 - \delta _{H}^\alpha n_2)\frac{t^{\alpha }}{\Gamma ({\alpha +1})}, \\ \\ C_1 &{}=&{} (\frac{\beta _{C}^\alpha n_1n_3}{{\mathscr {N}}} -\gamma _C^\alpha n_3 - \delta _C^\alpha n_3)\frac{t^{\alpha }}{\Gamma ({\alpha +1})}. \end{array} \end{aligned}$$
(59)

and

$$\begin{aligned} S_2= & {} \bigg [-(\frac{\beta _{H}^\alpha n_2}{{\mathscr {N}}} + \frac{\beta _{C}^\alpha n_3}{{\mathscr {N}}} + \delta _S^\alpha )(\Lambda ^\alpha - \frac{\beta _{H}^\alpha n_1n_2}{{\mathscr {N}}} - \frac{\beta _{C}^\alpha n_1n_3}{{\mathscr {N}}} + \gamma _{H}^\alpha n_2 + \gamma _{C}^\alpha n_3 - \delta _{S}^\alpha n_1) \nonumber \\&+ (\gamma _{H}^\alpha -\frac{\beta _{H}^\alpha n_1}{{\mathscr {N}}} )(\frac{\beta _{H}^\alpha n_1n_2}{{\mathscr {N}}} -\gamma _{H}^\alpha n_2 - \delta _{H}^\alpha n_2)\nonumber \\&+ (\gamma _{C}^\alpha -\frac{\beta _{C}^\alpha n_1}{{\mathscr {N}}} )(\frac{\beta _{C}^\alpha n_1n_3}{{\mathscr {N}}} -\gamma _C^\alpha n_3 - \delta _C^\alpha n_3)\bigg ]\frac{t^{2\alpha }}{\Gamma ({2\alpha +1})}, \nonumber \\ H_2= & {} \bigg [(\frac{\beta _{H}^\alpha n_1}{{\mathscr {N}}}-\gamma _{H}^\alpha -\delta _H^\alpha )(\frac{\beta _{H}^\alpha n_1n_2}{{\mathscr {N}}} -\gamma _{H}^\alpha n_2 - \delta _{H}^\alpha n_2)\nonumber \\&+ \frac{\beta _{H}^\alpha n_2}{{\mathscr {N}}}(\Lambda ^\alpha - \frac{\beta _{H}^\alpha n_1n_2}{{\mathscr {N}}} - \frac{\beta _{C}^\alpha n_1n_3}{{\mathscr {N}}} + \gamma _{H}^\alpha n_2 + \gamma _{C}^\alpha n_3 - \delta _{S}^\alpha n_1)\bigg ]\frac{t^{2\alpha }}{\Gamma ({2\alpha +1})}, \nonumber \\ C_2= & {} \bigg [(\frac{\beta _{C}^\alpha n_1}{{\mathscr {N}}}-\gamma _{C}^\alpha -\delta _C^\alpha )(\frac{\beta _{C}^\alpha n_1n_3}{{\mathscr {N}}} -\gamma _{C}^\alpha n_3 - \delta _{C}^\alpha n_3) \nonumber \\&+ \frac{\beta _{C}^\alpha n_3}{{\mathscr {N}}}(\Lambda ^\alpha - \frac{\beta _{H}^\alpha n_1n_2}{{\mathscr {N}}} - \frac{\beta _{C}^\alpha n_1n_3}{{\mathscr {N}}} + \gamma _{H}^\alpha n_2 + \gamma _{C}^\alpha n_3 - \delta _{S}^\alpha n_1)\bigg ]\frac{t^{2\alpha }}{\Gamma ({2\alpha +1})}. \end{aligned}$$
(60)

By using the parameter values in Table 2 when \(n_1=90\), \(n_2=80\) and \(n_3=80\), the solution for three terms can be given by

$$\begin{aligned} S(t)= & {} 90-17.299\frac{t^{\alpha }}{\Gamma ({\alpha +1})}+6.2467\frac{t^{2\alpha }}{\Gamma ({2\alpha +1})} ,\nonumber \\ H(t)= & {} 80-16.8\frac{t^{\alpha }}{\Gamma ({\alpha +1})}+2.144\frac{t^{2\alpha }}{\Gamma ({2\alpha +1})}, \nonumber \\ C(t)= & {} 80-11.26\frac{t^{\alpha }}{\Gamma ({\alpha +1})}+0.497\frac{t^{2\alpha }}{\Gamma ({2\alpha +1})}. \end{aligned}$$
(61)

In particular case let \(\alpha =1\), so the solution after three terms is as below

$$\begin{aligned} S(t)= & {} 90-17.299t+6.2467t^2 , \nonumber \\ H(t)= & {} 80-16.8t+2.144t^2, \nonumber \\ C(t)= & {} 80-11.26t++0.497t^2. \end{aligned}$$
(62)
Table 2 Parameter values for the transmission dynamics of CA-MRSA and HA-MRSA obtained from the Beth Israel Deaconess Medical Center
Full size table
Fig. 7
figure 7

Comparison of Caputo and classical derivative for susceptible class S(t)

Full size image
Fig. 8
figure 8

Comparison of Caputo and classical derivative for susceptible class S(t)

Full size image
Fig. 9
figure 9

Comparison of Caputo and classical derivative for HA-MRSA class H(t)

Full size image
Fig. 10
figure 10

Comparison of Caputo and classical derivative for HA-MRSA class H(t)

Full size image
Fig. 11
figure 11

Comparison of Caputo and classical derivative for CA-MRSA class C(t)

Full size image
Fig. 12
figure 12

Comparison of Caputo and classical derivative for CA-MRSA class C(t)

Full size image
Fig. 13
figure 13

Comparison with Caputo derivative for S(t), H(t) and C(t) when \(\alpha =0.8\)

Full size image

For analysis purposes in the present work, the values of the parameters are obtained from the Beth Israel Deaconess Medical Center and values of the fractional-order parameter \(\alpha \) are chosen differently as can be seen in the figures. Initial conditions are given for obtaining the solutions of fractional-order staph infection model (6) as [90, 80, 80]. The behavior of susceptible patients who are not colonized with HA-MRSA or CA-MRSA S(t) can be seen in Figs. 7, 8. Patients colonized with HA-MRSA H(t) and patients colonized with CA-MRSA C(t) can be seen in Figs. 9, 10, 11, 12, respectively. Also, the behavior of the solution curves of all three classes can be observed together in Fig. 13 when arbitrary order \(\alpha =0.8\). Figures 7, 8, 9, 10, 11, 12 are depicted with the varying values of the fractional-order \(\alpha \). The fractional and classical solution curves are shown separately for \(\alpha =1,0.9,0.8,0.7,0.6\) and \(\alpha =1,0.98,0.96,0.94\). In this way, it can be understood that as the \(\alpha \) values approach 1, the solution curves approach the classical solution. It should be noted that the classical solution is represented by \(\alpha =1\). For susceptible patients who are not colonized with HA-MRSA or CA-MRSA S(t), from Figs. 7 and 8, one can see that such choices of values \(\alpha \) indicate the decreasing rate for decreasing values of the fractional-order \(\alpha \). When susceptible subclass S(t) decrease as in Figs. 7 and 8, the danger in the hospital will also get decreased due to the potential decrease in the infectious classes H(t) and C(t). As clearly seen, the choices of arbitrary order \(\alpha \) can bring about a strongly decreasing rate in S(t), which can be superb and better for the individuals in the hospital. For the patients colonized with HA-MRSA H(t), Figs. 9 and 10 affirm a decreasing rate with a slight form of the increasing rate as the \(\alpha \) values keep varying lower. On the other hand, Figs. 11 and 12 for the patients colonized with CA-MRSA behaves similarly. This is a well-motivated impact as the number of H(t) and C(t) get decreased vehemently with the decrease of the susceptible class S(t). Also, this is a realistic situation which suggests that if the S(t) decreases then its density would decrease and consequently the infected classes H(t) and C(t) would reduce at the speedier rate for smaller values of the fractional-order \(\alpha \). Further observations affirm that for smaller values of \(\alpha \), the stability for S(t), H(t) and C(t) does not take long duration as it does for \(\alpha \rightarrow 1\).

5 Concluding remarks

Some important conclusions that can be drawn from this study can be listed as follows:

  • The governing model has been proposed for the very first time under Caputo fractional operator. This fractional model (6) analyzes community-acquired and hospital-acquired methicillin-resistant Staphylococcus aures which when assuming different values of the fractional-order \(\alpha \) allowed us to observe the behavior of the solution curves in detail.

  • The key explanation for advantages of the non-local fractional derivatives is owing to the presence of many orders of freedom for the operators which is not possible in the classical case.

  • The current research is crucial since the model under consideration is of a non-linear type.

  • By employing the concept of fixed point theory, the existence and uniqueness of solutions have been achieved by Caputo fractional derivative.

  • Furthermore, the iterative solutions have been shown by utilizing the Laplace transform of Caputo fractional derivative. Also, the stability criteria for the fractional model (6) has also been presented by using the fixed point theorem.

  • On the other hand, the sensitivity analysis of the basic reproduction number of the aforementioned model has been carried out in order to determine the parameters that have the most crucial effect on the spread of the disease.

  • With the help of an efficient method of LADM, the numerical characteristics of the fractional model (6) have been depicted by the parameter values obtained from the Beth Israel Deaconess Medical Center.

  • Multifarious graphical diagrams in this study suggest that the introduced fractional model is sufficiently capable to the entire system.

  • Generally, it can be expressed that fractional models are thus not only an important tool for a more advantageous representation of the behavior of physical organisms and epidemiology but also capable of better prediction for ecosystems.

  • It should be noted that the main goal of the present work is to employ the LADM and see its effectiveness on the model under consideration. However, as a future direction, the proposed model can be analyzed by utilizing some other useful numerical methods and compared with each other.