# Mathematical Analysis of the Transmission Dynamics of Skin Cancer Caused by UV Radiation

Authors: Tahera Parvin; Md. Haider Ali Biswas; Bimal Kumar Datta
Journal: Journal of Applied Mathematics
 (2022)

---

## Abstract

Nowadays, skin cancer is a worldwide panic. It is related to ultraviolet radiation. In this paper, we have formulated a SIRS type mathematical model to show the effects of ultraviolet radiation on skin cancer. At first, we have showed the boundedness and positivity of the model solutions to verify the model’s existence and uniqueness. The boundedness and positivity gave the solutions of our model bounded and positive, which was very important for real-world situation because in real world, population cannot be negative. Then, we have popped out all the equilibrium points of our model and verified the stability of the equilibrium points. This stability test expressed some physical situation of our model. The disease-free equilibrium point is locally asymptotically stable ifR0<1 and if R0>1, then it is unstable. Again, the endemic equilibrium point is stable, if R0>1 and unstable if R0<1. In order to understand the dynamical behavior of the model’s equilibrium points, we examined the phase portrait. We also have observed the sensitivity of the model parameters. After this, we have investigated the different scenarios of bifurcations of the model’s parameters. At the set of Hopf bifurcation parameters when infection rate due to UV rays is less than α1=0.01, proper control may eradicate the existence of disease. From transcritical bifurcation, we can say when recovery rate greater than 1.9, then the disease of skin cancer can be eliminated and when recovery rate less than 1.9 then the disease of skin cancer cannot be eradicated. Finally, numerical analysis is done to justify our analytical findings.


---

## Body

## 1. Introduction

Cancer is currently the leading cause of death worldwide. Error mutations in DNA are the most common cause of cancer. UV light, pollution, and other environmental factors primarily damage DNA, which can lead to uncontrolled cell development. There is no one who has not heard about cancer’s horrors. With a fatality rate of one in every six persons, cancer is one of the most common causes of death in contemporary times. Skin cancer is the most dangerous of them all. Skin cancer is a common occurrence in the United States. By the age of 70, one out of every five Americans will have developed skin cancer. Over two Americans die from skin cancer every hour. The two most common kinds of skin cancer are melanoma and nonmelanoma. Every month, more than 5,400 people die from nonmelanoma skin cancer throughout the world. Reference [1] shows, there were 19 292 789 new cancer diagnoses in 2020 and 10 million cancer-related deaths. Melanoma is the deadlier among all kinds of skin cancer [2]. Day by day the incidence rate of skin cancer is increasing (see Figures 1–5).Figure 1
Different cancer cases in 2020 worldwide (Source: [1]).Figure 2
New cases of melanoma from 2017 to 2020 worldwide (Source: [1]).Figure 3
New cases of nonmelanoma from 2017 to 2020 worldwide (Source: [1]).Figure 4
Percentage of new cases for melanoma skin cancer w.r.t. age structure groups in United State from 2014 to 2018 (Source: [3]).Figure 5
Percentage of deaths for melanoma skin cancer w.r.t. age structure groups in United State from 2015 to 2019 (Source: [3]).Numerous mathematical models have been created during the past ten years to explain the real-world situation among other topics. Different phenomena have been explained using these theories. The models that have been suggested are mostly ordinary differential equation models, both with and without delay components, that are linear and nonlinear [5]. Research of Newman et al. [4] discovered the UV index is increasing for the depletion of the ozone layers that is shown in Figure 6. According to the most recent WHO statistics published in 2018, the number of melanoma skin cancers death in Bangladesh is 320 out of 472 new cases of melanoma.Figure 6
Behavior of UV index with ozone layer (Source: [4]).Main cause of skin cancer is over exposure of UV rays. UVR is the part of electromagnetic spectrum which wavelengths is 100-400 nm. It is ejected by the sun and other manmade sources. Reduced stratospheric ozone layer will permit more UVB to reach the atmosphere. As a result, increased UV radiation from the sun and sunbeds may cause DNA damage in skin cells. When enough DNA damage accumulates over time, it can develop to skin tumors, which can then progress to skin cancer. Many scholars have studied theory-based and statistics-based studies on skin cancer. Fears et al. [6] formulated a mathematical model of skin cancer. In their study they had discussed the effects of age and UV rays on skin cancer. They had considered only the population who are fair skinned in the United States. De Gruul and Leun [7] developed a model considering the dose of UV rays for skin cancer incidence on human taking the results of previously discussed animal research. In their theory-based investigation, Moan et al. [8] and Shore [9] observed the relation between skin cancer and ultraviolet radiation. They also talked about skin cancer treatment alternatives. Bharath and Turner [10] showed the effect of climate change on skin cancer in their research. Besides these, Newton-Bishop et al. [11], Kim and He [12], Greinert et al. [13], and Berwick et al. [14] also discussed skin cancer model where UV radiation was the main risk factor. Biswas et al. has used mathematical modeling to observe the most depletory infectious disease [15]. We refer readers to [16–19] for more details of simple mathematical model.In our study, we proposed a four compartmental model based on skin cancer transmission characteristics. However, this is the skin cancer mathematical models in terms of system of nonlinear differential equations based on certain fundamental assumptions. Then, we have analyzed different types of analytical analysis of our propose model. Finally, we have observed the numerical simulations to validate our model and analytical findings.


## 2. Formulation of Mathematical Model

Although skin cancer is noncommunicable, but in very rare cases, cancer is transmitted by organ transplant [20]. For this reason, we have considered skin cancer as infectious. Assume the total number of populations is fixed entire the whole process which is defined by Nt. We consider four compartments with some fundamental assumptions. Skin cancer is very slow process which is caused due to ultraviolet radiation. People who are work at outside and remain in contact with sunlight are define as susceptible individuals and denoted by S1t. The progression of illness transmission is crucial to the disease’s dynamics. There are usually varied ranges of incubation duration for most noncommunicable diseases. Skin cancer is a noncommunicable disease mainly caused by long-term ultraviolet radiation’s exposure. Besides this there are several factors which are also responsible for skin cancer. So, considering the real phenomena another category we examine the infected individuals which are denoted by I1t. Those who have survived from skin cancer and are immune to it are denoted by Rt. This compartmental model’s flow chart is provided in Figure 7.Figure 7
Flow diagram of skin cancer model.According to the flowchart of the model in Figure7, the mathematical model of skin cancer can be written in the form of following nonlinear system of ordinary differential equations:
(1)dS1dt=k1−α1U+γ1I1S1−μ0S1+ηR,(2)dI1dt=α1U+γ1I1S1−ε1+μ0+ψ1I1,(3)dRdt=ψ1I1−μ0+ηR,(4)dUdt=r−μU.With initial conditionsS10=S10>0, I10=I10≥0, R0=R0≥0, and U0=U0>0.In Table1, we have described the parameters of our model.Table 1
Description of model parameters and state variables.
SymbolDescriptionk1Constant source rateα1Infection rateγ1Transmission rateμoNatural death rateηRelapse rateε1Cancer induced death rateψ1Recovery raterThe index of ultraviolet radiationμThe decay rate of ultraviolet radiationS1Susceptible individualsI1Infected individualsRRecovered individualsUUV radiation


## 3. Verification of the Properties of the Model’s Solution

Boundedness and positivity of the solution is very important properties of the solution of the autonomous system. Mainly it is used to define the well-posed system. So, at first, we examined the boundedness and positivity of the model solutions.


### 3.1. Boundedness of the Solutions

We have examined the boundedness of the system (1)–(4)'s solutions in this part. The system (1)–(4) describes the dynamical behavior of population during skin cancer period. Thus, the population size should be finite and nonnegative.Lemma 1.
The solutions of the system (1)–(4) are uniformly bounded within the regionΔ=S1,I1,R,U.
(5)∈ℝ+4:Pt=S1t+I1t+Rt+Ut,0<Pt≤μk1+μ0rμμ0.Proof.
A uniformly bounded family of functions is a family of bounded functions that can all be bounded by the same constant. This constant is larger than or equal to the absolute value of any value of any of the functions in the family.
Since the total population size isNt. So, we can write Nt=S1t+I1t+Rt.
Then,dNt/dt=dS1t/dt+dI1t/dt+dRt/dt.
From eqn. (1)-(3) of the system (1)–(4), we obtain,(6)Ntdt=k1−μ0N−ε1I1,⇒dNtdt+μ0N=k1−ε1I1,⇒dNtdt+μ0N≤k1.
Here, integrating factor,I.F=e∫μ0dt=eμ0t.
Multiplying both sides byeμ0t, we get
(7)eμ0tdNtdt+μ0Nteμ0t≤k1eμ0t,⇒dNteμ0t≤k1eμ0tdt.
Integrating both sides, we get(8)Nteμ0t≤k1eμ0tμ0+c11.
Using initial condition, we get(9)Nteμ0t≤k1eμ0tμ0+N0−k1μ0,⇒Nt≤k1μ0+N0e−μ0t−k1e−μ0tμ0,i.e.Nt≤k1μ01−e−μ0t+N0e−μ0t.
Att⟶∞, 0<Nt≤k1/μ0.
Now, from eqn. (4) of the system (1)–(4), we obtain(10)Utdt=r−μUt,⇒dUtdt+μUt=r.
Here, integrating factor,I.F=e∫μdt=eμt.
Multiplying both sides byeμt, we get
(11)eμtdUtdt+μUteμt=reμt,⇒dUteμt=reμtdt.
Integrating both sides, we get(12)Uteμt=reμtμ+c21.
Using initial condition, we get,(13)Uteμt=reμtμ+U0−rμ,⇒Ut=rμ+U0e−μt−re−μtμ,i.e.Ut=rμ1−e−μt+U0e−μt.
Att⟶∞, Ut=r/μ.
AssumePt=Nt+Ut.
So,(14)0+0<Nt+Ut≤k1μ0+rμ,(15)i.e.0<Pt≤μk1+μ0rμμ0.
Hence, the solutions of the system (1)–(4) are bounded in the regionΔ.


### 3.2. Positivity of the Solutions

Here, we check the positivity of the compartmentsS1t, I1t, Rt, and Ut. To investigate the positivity of these compartments we use lemma 2.Lemma 2.
IfS1t>0,I1t≥0,Rt≥0,Ut>0 and St,It,Mt,Ut∈ℝ+4, then the solutions S1t,I1t,Rt,Ut of the system (1)-(4) are positively invariant.Proof.
From equation (1), we get,(16)dS1dt=k1−α1U+γ1I1S1−μ0S1+ηR,dS1dt+μ0S1≥k1.
Here, integrating factor,I.F=e∫μ0dt=eμ0t.Multiplying both sides byeμ0t, we get
(17)eμ0tdS1dt+μ0S1eμ0t≥k1eμ0t,⇒dS1eμ0t≥k1eμ0t.Integrating both sides, we get(18)S1eμ0t≥k1eμ0tμ0+c31.Using initial condition, we get,(19)S1eμ0t≤k1eμ0tμ0+S10−k1μ0,⇒S1≤k1μ0+S10e−μ0t−k1e−μ0tμ0,i.e.S1≤k1μ01−e−μ0t+S10e−μ0t,Att⟶0,S1t>0.And also att⟶∞,S1t>0.Similarly, we can verify the positivity ofI1t, Rt, and Ut under the initial conditions. Therefore, the solutions S1t,I1t,Rt,Ut of the system (1)–(4) are positively invariant. Hence, the Lemma 2 is proved.


## 3.1. Boundedness of the Solutions

We have examined the boundedness of the system (1)–(4)'s solutions in this part. The system (1)–(4) describes the dynamical behavior of population during skin cancer period. Thus, the population size should be finite and nonnegative.Lemma 1.
The solutions of the system (1)–(4) are uniformly bounded within the regionΔ=S1,I1,R,U.
(5)∈ℝ+4:Pt=S1t+I1t+Rt+Ut,0<Pt≤μk1+μ0rμμ0.Proof.
A uniformly bounded family of functions is a family of bounded functions that can all be bounded by the same constant. This constant is larger than or equal to the absolute value of any value of any of the functions in the family.
Since the total population size isNt. So, we can write Nt=S1t+I1t+Rt.
Then,dNt/dt=dS1t/dt+dI1t/dt+dRt/dt.
From eqn. (1)-(3) of the system (1)–(4), we obtain,(6)Ntdt=k1−μ0N−ε1I1,⇒dNtdt+μ0N=k1−ε1I1,⇒dNtdt+μ0N≤k1.
Here, integrating factor,I.F=e∫μ0dt=eμ0t.
Multiplying both sides byeμ0t, we get
(7)eμ0tdNtdt+μ0Nteμ0t≤k1eμ0t,⇒dNteμ0t≤k1eμ0tdt.
Integrating both sides, we get(8)Nteμ0t≤k1eμ0tμ0+c11.
Using initial condition, we get(9)Nteμ0t≤k1eμ0tμ0+N0−k1μ0,⇒Nt≤k1μ0+N0e−μ0t−k1e−μ0tμ0,i.e.Nt≤k1μ01−e−μ0t+N0e−μ0t.
Att⟶∞, 0<Nt≤k1/μ0.
Now, from eqn. (4) of the system (1)–(4), we obtain(10)Utdt=r−μUt,⇒dUtdt+μUt=r.
Here, integrating factor,I.F=e∫μdt=eμt.
Multiplying both sides byeμt, we get
(11)eμtdUtdt+μUteμt=reμt,⇒dUteμt=reμtdt.
Integrating both sides, we get(12)Uteμt=reμtμ+c21.
Using initial condition, we get,(13)Uteμt=reμtμ+U0−rμ,⇒Ut=rμ+U0e−μt−re−μtμ,i.e.Ut=rμ1−e−μt+U0e−μt.
Att⟶∞, Ut=r/μ.
AssumePt=Nt+Ut.
So,(14)0+0<Nt+Ut≤k1μ0+rμ,(15)i.e.0<Pt≤μk1+μ0rμμ0.
Hence, the solutions of the system (1)–(4) are bounded in the regionΔ.


## 3.2. Positivity of the Solutions

Here, we check the positivity of the compartmentsS1t, I1t, Rt, and Ut. To investigate the positivity of these compartments we use lemma 2.Lemma 2.
IfS1t>0,I1t≥0,Rt≥0,Ut>0 and St,It,Mt,Ut∈ℝ+4, then the solutions S1t,I1t,Rt,Ut of the system (1)-(4) are positively invariant.Proof.
From equation (1), we get,(16)dS1dt=k1−α1U+γ1I1S1−μ0S1+ηR,dS1dt+μ0S1≥k1.
Here, integrating factor,I.F=e∫μ0dt=eμ0t.Multiplying both sides byeμ0t, we get
(17)eμ0tdS1dt+μ0S1eμ0t≥k1eμ0t,⇒dS1eμ0t≥k1eμ0t.Integrating both sides, we get(18)S1eμ0t≥k1eμ0tμ0+c31.Using initial condition, we get,(19)S1eμ0t≤k1eμ0tμ0+S10−k1μ0,⇒S1≤k1μ0+S10e−μ0t−k1e−μ0tμ0,i.e.S1≤k1μ01−e−μ0t+S10e−μ0t,Att⟶0,S1t>0.And also att⟶∞,S1t>0.Similarly, we can verify the positivity ofI1t, Rt, and Ut under the initial conditions. Therefore, the solutions S1t,I1t,Rt,Ut of the system (1)–(4) are positively invariant. Hence, the Lemma 2 is proved.


## 4.Model Analysis

Since it is impossible to find the exact solution of the nonlinear autonomous system (1)–(4), we have to analyze the qualitative behavior of the solutions in the neighborhood of the equilibrium points. So, in this section the nonlinear system of equation (1)–(4) has qualitatively analyzed to find the local and global stability of the different equilibrium points.


### 4.1. Equilibrium Points

The equilibrium points of the system (1)–(4) are obtained by equating(20)dS1dt=dI1dt=dRdt=dUdt=0Thus, we have(21)k1−α1U+γ1I1S1−μ0S1+ηR=0,(22)α1U+γ1I1S1−ε1+μ0+ψ1I1=0,(23)ψ1I1−μ0+ηR=0,(24)r−μU=0.


#### 4.1.1. Ultraviolet Free Equilibrium Points

For this caseU=0, then from (21)–(24) we get two equilibrium points, one is disease free and the another one is endemic equilibrium point in absent of ultraviolet radiation.E11=k1/μ0,0,0,0 and E12=S¯1,I¯1,R,0.Here, Complete disease and risk of disease-free equilibrium isE11=k1/μ0,0,0,0.And endemic equilibrium point isE12=S¯1,I¯1,R,0.Here,(25)S¯1=ϵ1+μ0+ψ1γ1,I¯1=−η+μ0ϵ1μ0−γ1k1+μ0ψ1+μ02γ1ϵ1η+ϵ1μ0+ημ0+μ0ψ1+μ02,R=−ψ1ϵ1μ0−γ1k1+μ0ψ1+μ02γ1ϵ1η+ϵ1μ0+ημ0+μ0ψ1+μ02.


#### 4.1.2. Disease Free Equilibrium Point WhenU≠0

At disease free equilibrium point (DFE),I1=R=0. Thus, the system (21)–(24) reduces to
(26)k1−α1rμS1−μ0S1=0,U=rμ,which implies
(27)S1=k1μα1r+μμ0.So, disease-free equilibrium point (DFE) is as follows:(28)E13:k1μα1r+μμ0,0,0,rμ.


#### 4.1.3. Endemic Equilibrium Points WhenU≠0

If all populations exist, the system (1)–(4) present endemic equilibrium (EE) point given by E1∗=S1∗,I1∗,R∗,U∗.Equation (24) gives U=r/μ.From equation (23), we get
(29)I1=μ0+ηψ1R.Putting the values ofUandI1 from equation (22), we obtain
(30)S1=με1+μ0+ψ1μ0+ηRα1ψ1r+γ1μμ0+ηR.Then equation (21) reduces to
(31)a0R2+a1R+a2=0.Here,(32)a0=ηψ1γ1μμ0+η−γ1με1+μ0+ψ1μ0+η2,a1=k1ψ1γ1μμ0+η−α1rψ1ε1+μ0+ψ1μ0+η−μμ0ψ1ε1+μ0+ψ1μ0+η+α1rηψ12,a2=k1α1rηψ12.Only real positive solutions of the quadratic equation (31) provide biological relevant steady state. Based on parameters values of system (1)-(4), we can have between zero and two endemic equilibria. Among them at least one will be positive using Descurt’s rule of sign if(1)
a0,a1>0anda2<0(2)
a0,a1<0anda2>0(3)
a0>0anda1,a2<0(4)
a0<0anda1,a2>0Using MATLAB, we get two positive endemic equilibrium points, whereS1∗,I1∗,R∗,U∗ represent the number of susceptible, infected, and recovered individuals and the last one is the index of ultraviolet radiation.
(33)S1∗=με1+μ0+ψ1μ0+ηR∗α1ψ1r+γ1μμ0+ηR∗,I1∗=μ0+ηψ1R∗,R∗=R∗,U∗=rμ.So, the endemic equilibrium point is(34)E1∗:με1+μ0+ψ1μ0+ηR∗α1ψ1r+γ1μμ0+ηR∗,μ0+ηψ1R∗,R∗,rμ.


### 4.2. Basic Reproduction Number

By focusing on the critical components of a disease, determining threshold values for illness survival, and evaluating the effect of various control techniques, mathematical modeling can play an important role in helping to quantify feasible disease control strategies. The basic reproduction number, also known as the basic reproductive number or basic reproductive ratio [21], is a critical threshold variable. It is generally denoted by R0. The epidemiological definition of R0 is the average number of secondary cases produced by one infected individual introduced into a population of susceptible individuals, where an infected individual has acquired the disease, and susceptible individuals are healthy but can acquire the disease. It is a key epidemiological quantity, because it determines the size and duration of epidemics. If R0>1, the occurrence of the disease will increase. If R0<1, the occurrence of the disease will decrease and the disease will ultimately be eliminated. When R0=1, the disease will be constant. Using the Van Den Driesseche and Watmough next generation approach and Blower et al. [22] concepts, we calculated the basic infection reproduction number of the systems (1)-(4), for more details see also [23, 24]. The vectors Fi and Vi are filled with appropriate terms from the infected class equations using this method. Terms that describe the appearance of new illnesses belong in the Fi category, while terms that describe the transmission of existing infections are in the Vi category and should be avoided. The matrices F and V are constructed and evaluated at a nontrivial disease-free equilibrium using the Jacobian matrices generated by differentiating Fi and Vi with respect to the relevant subset of variables.Here according to our model of skin cancer, we consider fast skin cancer which is caused by the contact of ultraviolet radiation and SLOW skin cancer refers to that skin cancer which is caused by other reasons. So, forR0Fast,F=γ1S10, and V=ε1+μ0+ψ1.So, we obtain(35)V−1=1ε1+μ0+ψ1.Thus,(36)FV−1=γ1S10ε1+μ0+ψ1(37)∴ρFV−1=γ1S10ε1+μ0+ψ1.So, forR0Slow, F=α1S10 and V=ε1+μ0+ψ1.So, we obtain(38)V−1=1ε1+μ0+ψ1.Thus,(39)FV−1=α1S10ε1+μ0+ψ1,(40)∴ρFV−1=α1S10ε1+μ0+ψ1.So, the basic reproduction number for our total model isR0=R0Fast+R0Slow.
(41)i.e.,R0=γ1S10ε1+μ0+ψ1+α1S10ε1+μ0+ψ1=α1+γ1S10ε1+μ0+ψ1.So, at disease free equilibrium point(42)R0=α1+γ1S10ε1+μ0+ψ1=α1+γ1k1μ0ε1+μ0+ψ1.For the base line parameters of our system the basic reproduction number is1.005..


## 4.1. Equilibrium Points

The equilibrium points of the system (1)–(4) are obtained by equating(20)dS1dt=dI1dt=dRdt=dUdt=0Thus, we have(21)k1−α1U+γ1I1S1−μ0S1+ηR=0,(22)α1U+γ1I1S1−ε1+μ0+ψ1I1=0,(23)ψ1I1−μ0+ηR=0,(24)r−μU=0.


### 4.1.1. Ultraviolet Free Equilibrium Points

For this caseU=0, then from (21)–(24) we get two equilibrium points, one is disease free and the another one is endemic equilibrium point in absent of ultraviolet radiation.E11=k1/μ0,0,0,0 and E12=S¯1,I¯1,R,0.Here, Complete disease and risk of disease-free equilibrium isE11=k1/μ0,0,0,0.And endemic equilibrium point isE12=S¯1,I¯1,R,0.Here,(25)S¯1=ϵ1+μ0+ψ1γ1,I¯1=−η+μ0ϵ1μ0−γ1k1+μ0ψ1+μ02γ1ϵ1η+ϵ1μ0+ημ0+μ0ψ1+μ02,R=−ψ1ϵ1μ0−γ1k1+μ0ψ1+μ02γ1ϵ1η+ϵ1μ0+ημ0+μ0ψ1+μ02.


### 4.1.2. Disease Free Equilibrium Point WhenU≠0

At disease free equilibrium point (DFE),I1=R=0. Thus, the system (21)–(24) reduces to
(26)k1−α1rμS1−μ0S1=0,U=rμ,which implies
(27)S1=k1μα1r+μμ0.So, disease-free equilibrium point (DFE) is as follows:(28)E13:k1μα1r+μμ0,0,0,rμ.


### 4.1.3. Endemic Equilibrium Points WhenU≠0

If all populations exist, the system (1)–(4) present endemic equilibrium (EE) point given by E1∗=S1∗,I1∗,R∗,U∗.Equation (24) gives U=r/μ.From equation (23), we get
(29)I1=μ0+ηψ1R.Putting the values ofUandI1 from equation (22), we obtain
(30)S1=με1+μ0+ψ1μ0+ηRα1ψ1r+γ1μμ0+ηR.Then equation (21) reduces to
(31)a0R2+a1R+a2=0.Here,(32)a0=ηψ1γ1μμ0+η−γ1με1+μ0+ψ1μ0+η2,a1=k1ψ1γ1μμ0+η−α1rψ1ε1+μ0+ψ1μ0+η−μμ0ψ1ε1+μ0+ψ1μ0+η+α1rηψ12,a2=k1α1rηψ12.Only real positive solutions of the quadratic equation (31) provide biological relevant steady state. Based on parameters values of system (1)-(4), we can have between zero and two endemic equilibria. Among them at least one will be positive using Descurt’s rule of sign if(1)
a0,a1>0anda2<0(2)
a0,a1<0anda2>0(3)
a0>0anda1,a2<0(4)
a0<0anda1,a2>0Using MATLAB, we get two positive endemic equilibrium points, whereS1∗,I1∗,R∗,U∗ represent the number of susceptible, infected, and recovered individuals and the last one is the index of ultraviolet radiation.
(33)S1∗=με1+μ0+ψ1μ0+ηR∗α1ψ1r+γ1μμ0+ηR∗,I1∗=μ0+ηψ1R∗,R∗=R∗,U∗=rμ.So, the endemic equilibrium point is(34)E1∗:με1+μ0+ψ1μ0+ηR∗α1ψ1r+γ1μμ0+ηR∗,μ0+ηψ1R∗,R∗,rμ.


## 4.1.1. Ultraviolet Free Equilibrium Points

For this caseU=0, then from (21)–(24) we get two equilibrium points, one is disease free and the another one is endemic equilibrium point in absent of ultraviolet radiation.E11=k1/μ0,0,0,0 and E12=S¯1,I¯1,R,0.Here, Complete disease and risk of disease-free equilibrium isE11=k1/μ0,0,0,0.And endemic equilibrium point isE12=S¯1,I¯1,R,0.Here,(25)S¯1=ϵ1+μ0+ψ1γ1,I¯1=−η+μ0ϵ1μ0−γ1k1+μ0ψ1+μ02γ1ϵ1η+ϵ1μ0+ημ0+μ0ψ1+μ02,R=−ψ1ϵ1μ0−γ1k1+μ0ψ1+μ02γ1ϵ1η+ϵ1μ0+ημ0+μ0ψ1+μ02.


## 4.1.2. Disease Free Equilibrium Point WhenU≠0

At disease free equilibrium point (DFE),I1=R=0. Thus, the system (21)–(24) reduces to
(26)k1−α1rμS1−μ0S1=0,U=rμ,which implies
(27)S1=k1μα1r+μμ0.So, disease-free equilibrium point (DFE) is as follows:(28)E13:k1μα1r+μμ0,0,0,rμ.


## 4.1.3. Endemic Equilibrium Points WhenU≠0

If all populations exist, the system (1)–(4) present endemic equilibrium (EE) point given by E1∗=S1∗,I1∗,R∗,U∗.Equation (24) gives U=r/μ.From equation (23), we get
(29)I1=μ0+ηψ1R.Putting the values ofUandI1 from equation (22), we obtain
(30)S1=με1+μ0+ψ1μ0+ηRα1ψ1r+γ1μμ0+ηR.Then equation (21) reduces to
(31)a0R2+a1R+a2=0.Here,(32)a0=ηψ1γ1μμ0+η−γ1με1+μ0+ψ1μ0+η2,a1=k1ψ1γ1μμ0+η−α1rψ1ε1+μ0+ψ1μ0+η−μμ0ψ1ε1+μ0+ψ1μ0+η+α1rηψ12,a2=k1α1rηψ12.Only real positive solutions of the quadratic equation (31) provide biological relevant steady state. Based on parameters values of system (1)-(4), we can have between zero and two endemic equilibria. Among them at least one will be positive using Descurt’s rule of sign if(1)
a0,a1>0anda2<0(2)
a0,a1<0anda2>0(3)
a0>0anda1,a2<0(4)
a0<0anda1,a2>0Using MATLAB, we get two positive endemic equilibrium points, whereS1∗,I1∗,R∗,U∗ represent the number of susceptible, infected, and recovered individuals and the last one is the index of ultraviolet radiation.
(33)S1∗=με1+μ0+ψ1μ0+ηR∗α1ψ1r+γ1μμ0+ηR∗,I1∗=μ0+ηψ1R∗,R∗=R∗,U∗=rμ.So, the endemic equilibrium point is(34)E1∗:με1+μ0+ψ1μ0+ηR∗α1ψ1r+γ1μμ0+ηR∗,μ0+ηψ1R∗,R∗,rμ.


## 4.2. Basic Reproduction Number

By focusing on the critical components of a disease, determining threshold values for illness survival, and evaluating the effect of various control techniques, mathematical modeling can play an important role in helping to quantify feasible disease control strategies. The basic reproduction number, also known as the basic reproductive number or basic reproductive ratio [21], is a critical threshold variable. It is generally denoted by R0. The epidemiological definition of R0 is the average number of secondary cases produced by one infected individual introduced into a population of susceptible individuals, where an infected individual has acquired the disease, and susceptible individuals are healthy but can acquire the disease. It is a key epidemiological quantity, because it determines the size and duration of epidemics. If R0>1, the occurrence of the disease will increase. If R0<1, the occurrence of the disease will decrease and the disease will ultimately be eliminated. When R0=1, the disease will be constant. Using the Van Den Driesseche and Watmough next generation approach and Blower et al. [22] concepts, we calculated the basic infection reproduction number of the systems (1)-(4), for more details see also [23, 24]. The vectors Fi and Vi are filled with appropriate terms from the infected class equations using this method. Terms that describe the appearance of new illnesses belong in the Fi category, while terms that describe the transmission of existing infections are in the Vi category and should be avoided. The matrices F and V are constructed and evaluated at a nontrivial disease-free equilibrium using the Jacobian matrices generated by differentiating Fi and Vi with respect to the relevant subset of variables.Here according to our model of skin cancer, we consider fast skin cancer which is caused by the contact of ultraviolet radiation and SLOW skin cancer refers to that skin cancer which is caused by other reasons. So, forR0Fast,F=γ1S10, and V=ε1+μ0+ψ1.So, we obtain(35)V−1=1ε1+μ0+ψ1.Thus,(36)FV−1=γ1S10ε1+μ0+ψ1(37)∴ρFV−1=γ1S10ε1+μ0+ψ1.So, forR0Slow, F=α1S10 and V=ε1+μ0+ψ1.So, we obtain(38)V−1=1ε1+μ0+ψ1.Thus,(39)FV−1=α1S10ε1+μ0+ψ1,(40)∴ρFV−1=α1S10ε1+μ0+ψ1.So, the basic reproduction number for our total model isR0=R0Fast+R0Slow.
(41)i.e.,R0=γ1S10ε1+μ0+ψ1+α1S10ε1+μ0+ψ1=α1+γ1S10ε1+μ0+ψ1.So, at disease free equilibrium point(42)R0=α1+γ1S10ε1+μ0+ψ1=α1+γ1k1μ0ε1+μ0+ψ1.For the base line parameters of our system the basic reproduction number is1.005..


## 5. Stability Analysis

The behavior of solutions that start near the equilibrium solution is addressed by the physical stability of an equilibrium solution to a system of differential equations. There are two types of physical stability—local and global stability. The local stability of an equilibrium point means that if the system is placed near the point, it will eventually migrate to the equilibrium point. The term “global stability” refers to the system’s ability to reach equilibrium from any possible starting point.


### 5.1. Local Stability of Equilibrium Points

In this section, we observed the stability of all of the equilibrium points of the system.Theorem 1.
The equilibrium pointsE11=k1/μ0,0,0,0 of the system (1)-(4) are locally asymptotically stable if R0<1 and if R0>1, then it is unstable.Proof.
Let(43)dS1dt=P1,dI1dt=Q1,dRdt=R1,dUdt=T1.
Thus, we have(44)P1=k1−α1U+γ1I1S1−μ0S1+ηR,(45)Q1=α1U+γ1I1S1−ε1+μ0+ψ1I1,(46)R1=ψ1I1−μ0+ηR,(47)T1=r−μU.The Jacobian matrix of the system (44)–(46) can be written as
(48)J=∂P1,Q1,R1,T1∂S1,I1,R,U=∂P1∂S1∂P1∂I1∂P1∂R∂P1∂U∂Q1∂S1∂Q1∂I1∂Q1∂R∂Q1∂U∂R1∂S1∂R1∂I1∂R1∂R∂R1∂U∂T1∂S1∂T1∂I1∂T1∂R∂T1∂U.So that(49)J=−μ0−I1γ1−Uα1−S1γ1η−S1α1I1γ1+Uα1S1γ1−μ0−ψ1−ϵ10S1α10ψ1−η−μ00000−μ,At pointE11=k1/μ0,0,0,0 the Jacobian matrix becomes
(50)JE11=−μ0−γ1k1μ0η−α1k1μ00γ1k1μ0−μ0−ψ1−ϵ10α1k1μ00ψ1−η−μ00000−μ.The characteristic equation of matrixJE11 is detJE11−ρI=0.So, we get(51)−μ0−ρ−γ1k1μ0η−α1k1μ00γ1k1μ0−μ0−ψ1−ρ−ϵ10α1k1μ00ψ1−η−μ0−ρ0000−μ−ρ=0.So, the eigen values are(52)ρ1=−μ0,ρ2=γ1k1μ0+α1k1μ0−μ0−ψ1−ϵ1−α1k1μ0,=α1+γ1k1μ0−μ0−ψ1−ϵ1−α1k1μ0,=ϵ1+μ0+ψ1α1+γ1k1μ0ϵ1+μ0+ψ1−1−α1k1μ0,=ϵ1+μ0+ψ1R0−1−−α1k1μ0,ρ3=−η−μ0,ρ4=−μ.The eigen values ρ1=−μ0,ρ3=−η−μ0andρ4=−μ are negative, and ρ2 will be negative if R0<1.So,E11 is stable, ifR0<1, otherwise unstable. Hence, Theorem 1 is proved.Theorem 2.
The equilibrium pointsE12=S¯1,I¯1,R¯,0 of the system (1)-(4) are stable, if ϵ1μ0<γ1k1+μ0ψ1+μ02 otherwise unstable.Proof.
From equation (43) the Jacobian matrix at pointE12=S¯1,I¯1,R¯,0 becomes
(53)JE12=η+μ0ϵ1μ0−γ1k1+μ0ψ1+μ02ϵ1η+ϵ1μ0+ημ0+μ0ψ1+μ02−μ0−ϵ1−μ0−ψ1η−α1ϵ1+μ0+ψ1γ1−η+μ0ϵ1μ0−γ1k1+μ0ψ1+μ02ϵ1η+ϵ1μ0+ημ0+μ0ψ1+μ0200α1ϵ1+μ0+ψ1γ10ψ1−η−μ00000−μ.
The characteristic equation of matrixJE12 is detJE12−ρI=0.So, we get(54)η+μ0ϵ1μ0−γ1k1+μ0ψ1+μ02ϵ1η+ϵ1μ0+ημ0+μ0ψ1+μ02−ρ−μ0−ϵ1−μ0−ψ1η−α1ϵ1+μ0+ψ1γ1−η+μ0ϵ1μ0−γ1k1+μ0ψ1+μ02ϵ1η+ϵ1μ0+ημ0+μ0ψ1+μ02−ρ0α1ϵ1+μ0+ψ1γ10ψ1−η−μ0−ρ0000−μ−ρ=0,⇒a11−ϵ1−μ0−ψ1η−α1ϵ1+μ0+ψ1γ10a22a23a240ψ1−η−μ0−ρ0000−μ−ρ=0,a11=η+μ0ϵ1μ0−γ1k1+μ0ψ1+μ02ϵ1η+ϵ1μ0+ημ0+μ0ψ1+μ02−ρ−μ0,a22=η+μ0ϵ1+μ0+ψ1ϵ1μ0−γ1k1+μ0ψ1+μ02μ0+ρ−η+μ0ϵ1μ0−γ1k1+μ0ψ1+μ02/ϵ1η+ϵ1μ0+ημ0+μ0ψ1+μ02ϵ1η+ϵ1μ0+ημ0+μ0ψ1+μ02−ρ,a23=−ηη+μ0ϵ1μ0−γ1k1+μ0ψ1+μ02μ0+ρ−η+μ0ϵ1μ0−γ1k1+μ0ψ1+μ02/ϵ1η+ϵ1μ0+ημ0+μ0ψ1+μ02ϵ1η+ϵ1μ0+ημ0+μ0ψ1+μ02,a24=α1ϵ1+μ0+ψ1γ1+α1η+μ0ϵ1+μ0+ψ1ϵ1μ0−γ1k1+μ0ψ1+μ02γ1μ0+ρ−η+μ0ϵ1μ0−γ1k1+μ0ψ1+μ02/ϵ1η+ϵ1μ0+ημ0+μ0ψ1+μ02ϵ1η+ϵ1μ0+ημ0+μ0ψ1+μ02,⇒a11−ϵ1−μ0−ψ1η−α1ϵ1+μ0+ψ1γ10a22a23a2400a33a34000−μ−ρ=0,a33=ψ1μ0+ρ−η+μ0ϵ1μ0−γ1k1+μ0ψ1+μ02/ϵ1η+ϵ1μ0+ημ0+μ0ψ1+μ02ϵ1η+ϵ1μ0+ημ0+μ0ψ1+μ02ηη+μ0ϵ1μ0−γ1k1+μ0ψ1+μ02−μ−ρ−η,a34=−ψ1a24.So, the eigen values are(55)ρ1=η+μ0ϵ1μ0−γ1k1+μ0ψ1+μ02ϵ1η+ϵ1μ0+ημ0+μ0ψ1+μ02−μ0,ρ2=η+μ0ϵ1+μ0+ψ1ϵ1μ0−γ1k1+μ0ψ1+μ02μ0+ρ−η+μ0ϵ1μ0−γ1k1+μ0ψ1+μ02/ϵ1η+ϵ1μ0+ημ0+μ0ψ1+μ02ϵ1η+ϵ1μ0+ημ0+μ0ψ1+μ02,ρ3=η+μ0ϵ1μ0−γ1k1+μ0ψ1+μ02ϵ1η+ϵ1μ0+ημ0+μ0ψ1+μ02−μ0,ρ4=−μHere,ρ4 is negative.The eigen valuesρ1,ρ2,andρ3 will be negative, if ϵ1μ0<γ1k1+μ0ψ1+μ02.So,E12 is stable, if ϵ1μ0<γ1k1+μ0ψ1+μ02 otherwise unstable. Hence, Theorem 2 is proved.Theorem 3.
The equilibrium pointsE13:k1μ/α1r+μμ0,0,0,r/μ of the system (1)-(4) are locally asymptotically stable if R0<1 and if R0>1, then it is unstable.Proof.
From equation (49) at point E13:k1μ/α1r+μμ0,0,0,r/μ the Jacobian matrix becomes
(56)JE13=−μ0−α1rμ−γ1k1μα1r+μμ0η−α1k1μα1r+μμ0α1rμγ1k1μα1r+μμ0−μ0−ψ1−ϵ10α1k1μα1r+μμ00ψ1−η−μ00000−μ.
The characteristic equation of matrixJE13 is detJE13−ρI=0.So, we get(57)−μ0−ρ−α1rμ−γ1k1μα1r+μμ0η−α1k1μα1r+μμ0μγ1k1μα1r+μμ0−μ0−ψ1−ρ−ϵ10α1k1μα1r+μμ00ψ1−η−μ0−ρ0000−μ−ρ=0,⇒−μ0−ρ−α1rμ−γ1k1μα1r+μμ0η−α1k1μα1r+μμ00b22b23b2400b33b34000−μ−ρ=0,b22=γ1k1μα1r+μμ0−μ0−ψ1−ρ−ϵ1−α1γ1k1rα1r+μμ0μ0+ρ+α1r/μ,b23=α1ηrμμ0+ρ+α1r/μ,b24=α1k1μα1r+μμ0−α12k1rα1r+μμ0μ0+ρ+α1r/μ,b33=α1ηψ1rμμ0+ρ+α1r/μϵ1+μ0+ψ1+ρ−γ1k1μ/α1r+μμ0+α1γ1k1r/α1r+μμ0μ0+ρ+α1r/μ−μ0−η,b34=ψ1α1k1μ/α1r+μμ0−α12k1r/α1r+μμ0μ0+ρ+α1r/μϵ1+μ0+ψ1+ρ−γ1k1μ/α1r+μμ0+α1γ1k1r/α1r+μμ0μ0+ρ+α1r/μ.So, the eigen values are(58)ρ1=−μ0−α1rμ,ρ2=α1+γ1k1μ0−μ0−ψ1−ϵ1−α1γ1k1rα1r+μμ0μ0+α1r/μ−α1k1μ0=μ0+ψ1+ϵ1γ1k1μ0μ0+ψ1+ϵ1−1−α1γ1k1rα1r+μμ0μ0+α1r/μ−α1k1μ0=μ0+ψ1+ϵ1R0−1−α1γ1k1rα1r+μμ0μ0+α1r/μ−α1k1μ0,ρ3=α1γ1k1rα1r+μμ0μ0+α1r/μα1+γ1k1μ0ϵ1+μ0+ψ1−1−α1ηψ1rμ0+η=α1γ1k1rα1r+μμ0μ0+α1r/μR0−1−α1ηψ1rμ0+η,ρ4=−μ.The eigen valuesρ1 and ρ4 are negative, and ρ2 and ρ3 will be negative if R0<1.So,E13 is locally asymptotically stable if R0<1, and if R0>1, then it is unstable. Hence, Theorem 3 is proved.Figure8 represents the stability of the equilibrium point E13. Here, the parameters values are k1=0.01,α1=0.003,γ1=0.0018,μ0=0.00561,η=0.0092,ε1=0.001,ψ1=0.93,r=10,μ=0.09.In this case, neglecting infected and recovered class the equilibrium point becomes (0.03,111). For this parameter values the basic reproduction number R0is less than unity.Figure 8
Phase portrait for equilibrium point (0.03,111) whenR0<1.Figure9 shows unstable phase portrait for equilibrium point E13. Here, the parameters values are k1=10.1, α1=0.3, γ1=0.0018, μ0=0.561, η=0.0092, ε1=0.001, ψ1=0.93, r=10, and μ=0.9. In this case, neglecting infected and recovered class the equilibrium point becomes (2.6,11). For this parameter values, the basic reproduction number R0 is greater than unity. So, we can say Figures 8 and 9 justify our analytical findings for the equilibrium point E13.Theorem 4.
The equilibrium pointsE1∗=S1∗,I1∗,R∗,U∗ of the system (1)-(4) is stable, if R0>1 and unstable if R0<1.Proof.
From equation (49) the Jacobian matrix at point E1∗=S1∗,I1∗,R∗,U∗ becomes,
(59)JE1∗=−μ0−I1∗γ1−U∗α1−S1∗γ1η−S1∗α1I1∗γ1+U∗α1S1∗γ1−μ0−ψ1−ϵ10S1∗α10ψ1−η−μ00000−μ.
The characteristic equation of matrixJE1∗ is detJE1∗−ρI=0.Figure 9
Unstable phase portrait for equilibrium point (2.6,11) whenR0>1.So, we get(60)−μ0−I1∗γ1−U∗α1−ρ−S1∗γ1η−S1∗α1I1∗γ1+U∗α1S1∗γ1−μ0−ψ1−ϵ1−ρ0S1∗α10ψ1−η−μ0−ρ0000−μ−ρ=0,−μ0−I1∗γ1−U∗α1−ρ−S1∗γ1η−S1∗α10c22c23c2400c33c34000−μ−ρ=0,c22=S1∗γ1−μ0−ψ1−ρ−ϵ1−S1∗γ1I1∗γ1+U∗α1μ0+ρ+I1∗γ1+U∗α1,c23=S1∗γ1−μ0−ψ1−ρ−ϵ1+ηI1∗γ1+U∗α1μ0+ρ+I1∗γ1+U∗α1,c24=S1∗α1−S1∗α1I1∗γ1+U∗α1μ0+ρ+I1∗γ1+U∗α1,c33=−η−μ0−ρ−ψ1ϵ1+μ0+ψ1+ρ−S1∗γ1−ηI1∗γ1+U∗α1/μ0+ρ+I1∗γ1+U∗α1ϵ1+μ0+ψ1+ρ−S1∗γ1+S1∗γ1I1∗γ1+U∗α1/μ0+ρ+I1∗γ1+U∗α1,c34=ψ1S1∗α1−S1∗α1I1∗γ1+U∗α1/μ0+ρ+I1∗γ1+U∗α1ϵ1+μ0+ψ1+ρ−S1∗γ1+S1∗γ1I1∗γ1+U∗α1/μ0+ρ+I1∗γ1+U∗α1.Here, one of the eigen value isρ=−μ, and according to Routh-Hurwitz criteria the remaining roots will be negative if R0>1.So,E1∗ is stable, if R0>1 otherwise unstable. Hence, Theorem 4 is proved.


## 5.1. Local Stability of Equilibrium Points

In this section, we observed the stability of all of the equilibrium points of the system.Theorem 1.
The equilibrium pointsE11=k1/μ0,0,0,0 of the system (1)-(4) are locally asymptotically stable if R0<1 and if R0>1, then it is unstable.Proof.
Let(43)dS1dt=P1,dI1dt=Q1,dRdt=R1,dUdt=T1.
Thus, we have(44)P1=k1−α1U+γ1I1S1−μ0S1+ηR,(45)Q1=α1U+γ1I1S1−ε1+μ0+ψ1I1,(46)R1=ψ1I1−μ0+ηR,(47)T1=r−μU.The Jacobian matrix of the system (44)–(46) can be written as
(48)J=∂P1,Q1,R1,T1∂S1,I1,R,U=∂P1∂S1∂P1∂I1∂P1∂R∂P1∂U∂Q1∂S1∂Q1∂I1∂Q1∂R∂Q1∂U∂R1∂S1∂R1∂I1∂R1∂R∂R1∂U∂T1∂S1∂T1∂I1∂T1∂R∂T1∂U.So that(49)J=−μ0−I1γ1−Uα1−S1γ1η−S1α1I1γ1+Uα1S1γ1−μ0−ψ1−ϵ10S1α10ψ1−η−μ00000−μ,At pointE11=k1/μ0,0,0,0 the Jacobian matrix becomes
(50)JE11=−μ0−γ1k1μ0η−α1k1μ00γ1k1μ0−μ0−ψ1−ϵ10α1k1μ00ψ1−η−μ00000−μ.The characteristic equation of matrixJE11 is detJE11−ρI=0.So, we get(51)−μ0−ρ−γ1k1μ0η−α1k1μ00γ1k1μ0−μ0−ψ1−ρ−ϵ10α1k1μ00ψ1−η−μ0−ρ0000−μ−ρ=0.So, the eigen values are(52)ρ1=−μ0,ρ2=γ1k1μ0+α1k1μ0−μ0−ψ1−ϵ1−α1k1μ0,=α1+γ1k1μ0−μ0−ψ1−ϵ1−α1k1μ0,=ϵ1+μ0+ψ1α1+γ1k1μ0ϵ1+μ0+ψ1−1−α1k1μ0,=ϵ1+μ0+ψ1R0−1−−α1k1μ0,ρ3=−η−μ0,ρ4=−μ.The eigen values ρ1=−μ0,ρ3=−η−μ0andρ4=−μ are negative, and ρ2 will be negative if R0<1.So,E11 is stable, ifR0<1, otherwise unstable. Hence, Theorem 1 is proved.Theorem 2.
The equilibrium pointsE12=S¯1,I¯1,R¯,0 of the system (1)-(4) are stable, if ϵ1μ0<γ1k1+μ0ψ1+μ02 otherwise unstable.Proof.
From equation (43) the Jacobian matrix at pointE12=S¯1,I¯1,R¯,0 becomes
(53)JE12=η+μ0ϵ1μ0−γ1k1+μ0ψ1+μ02ϵ1η+ϵ1μ0+ημ0+μ0ψ1+μ02−μ0−ϵ1−μ0−ψ1η−α1ϵ1+μ0+ψ1γ1−η+μ0ϵ1μ0−γ1k1+μ0ψ1+μ02ϵ1η+ϵ1μ0+ημ0+μ0ψ1+μ0200α1ϵ1+μ0+ψ1γ10ψ1−η−μ00000−μ.
The characteristic equation of matrixJE12 is detJE12−ρI=0.So, we get(54)η+μ0ϵ1μ0−γ1k1+μ0ψ1+μ02ϵ1η+ϵ1μ0+ημ0+μ0ψ1+μ02−ρ−μ0−ϵ1−μ0−ψ1η−α1ϵ1+μ0+ψ1γ1−η+μ0ϵ1μ0−γ1k1+μ0ψ1+μ02ϵ1η+ϵ1μ0+ημ0+μ0ψ1+μ02−ρ0α1ϵ1+μ0+ψ1γ10ψ1−η−μ0−ρ0000−μ−ρ=0,⇒a11−ϵ1−μ0−ψ1η−α1ϵ1+μ0+ψ1γ10a22a23a240ψ1−η−μ0−ρ0000−μ−ρ=0,a11=η+μ0ϵ1μ0−γ1k1+μ0ψ1+μ02ϵ1η+ϵ1μ0+ημ0+μ0ψ1+μ02−ρ−μ0,a22=η+μ0ϵ1+μ0+ψ1ϵ1μ0−γ1k1+μ0ψ1+μ02μ0+ρ−η+μ0ϵ1μ0−γ1k1+μ0ψ1+μ02/ϵ1η+ϵ1μ0+ημ0+μ0ψ1+μ02ϵ1η+ϵ1μ0+ημ0+μ0ψ1+μ02−ρ,a23=−ηη+μ0ϵ1μ0−γ1k1+μ0ψ1+μ02μ0+ρ−η+μ0ϵ1μ0−γ1k1+μ0ψ1+μ02/ϵ1η+ϵ1μ0+ημ0+μ0ψ1+μ02ϵ1η+ϵ1μ0+ημ0+μ0ψ1+μ02,a24=α1ϵ1+μ0+ψ1γ1+α1η+μ0ϵ1+μ0+ψ1ϵ1μ0−γ1k1+μ0ψ1+μ02γ1μ0+ρ−η+μ0ϵ1μ0−γ1k1+μ0ψ1+μ02/ϵ1η+ϵ1μ0+ημ0+μ0ψ1+μ02ϵ1η+ϵ1μ0+ημ0+μ0ψ1+μ02,⇒a11−ϵ1−μ0−ψ1η−α1ϵ1+μ0+ψ1γ10a22a23a2400a33a34000−μ−ρ=0,a33=ψ1μ0+ρ−η+μ0ϵ1μ0−γ1k1+μ0ψ1+μ02/ϵ1η+ϵ1μ0+ημ0+μ0ψ1+μ02ϵ1η+ϵ1μ0+ημ0+μ0ψ1+μ02ηη+μ0ϵ1μ0−γ1k1+μ0ψ1+μ02−μ−ρ−η,a34=−ψ1a24.So, the eigen values are(55)ρ1=η+μ0ϵ1μ0−γ1k1+μ0ψ1+μ02ϵ1η+ϵ1μ0+ημ0+μ0ψ1+μ02−μ0,ρ2=η+μ0ϵ1+μ0+ψ1ϵ1μ0−γ1k1+μ0ψ1+μ02μ0+ρ−η+μ0ϵ1μ0−γ1k1+μ0ψ1+μ02/ϵ1η+ϵ1μ0+ημ0+μ0ψ1+μ02ϵ1η+ϵ1μ0+ημ0+μ0ψ1+μ02,ρ3=η+μ0ϵ1μ0−γ1k1+μ0ψ1+μ02ϵ1η+ϵ1μ0+ημ0+μ0ψ1+μ02−μ0,ρ4=−μHere,ρ4 is negative.The eigen valuesρ1,ρ2,andρ3 will be negative, if ϵ1μ0<γ1k1+μ0ψ1+μ02.So,E12 is stable, if ϵ1μ0<γ1k1+μ0ψ1+μ02 otherwise unstable. Hence, Theorem 2 is proved.Theorem 3.
The equilibrium pointsE13:k1μ/α1r+μμ0,0,0,r/μ of the system (1)-(4) are locally asymptotically stable if R0<1 and if R0>1, then it is unstable.Proof.
From equation (49) at point E13:k1μ/α1r+μμ0,0,0,r/μ the Jacobian matrix becomes
(56)JE13=−μ0−α1rμ−γ1k1μα1r+μμ0η−α1k1μα1r+μμ0α1rμγ1k1μα1r+μμ0−μ0−ψ1−ϵ10α1k1μα1r+μμ00ψ1−η−μ00000−μ.
The characteristic equation of matrixJE13 is detJE13−ρI=0.So, we get(57)−μ0−ρ−α1rμ−γ1k1μα1r+μμ0η−α1k1μα1r+μμ0μγ1k1μα1r+μμ0−μ0−ψ1−ρ−ϵ10α1k1μα1r+μμ00ψ1−η−μ0−ρ0000−μ−ρ=0,⇒−μ0−ρ−α1rμ−γ1k1μα1r+μμ0η−α1k1μα1r+μμ00b22b23b2400b33b34000−μ−ρ=0,b22=γ1k1μα1r+μμ0−μ0−ψ1−ρ−ϵ1−α1γ1k1rα1r+μμ0μ0+ρ+α1r/μ,b23=α1ηrμμ0+ρ+α1r/μ,b24=α1k1μα1r+μμ0−α12k1rα1r+μμ0μ0+ρ+α1r/μ,b33=α1ηψ1rμμ0+ρ+α1r/μϵ1+μ0+ψ1+ρ−γ1k1μ/α1r+μμ0+α1γ1k1r/α1r+μμ0μ0+ρ+α1r/μ−μ0−η,b34=ψ1α1k1μ/α1r+μμ0−α12k1r/α1r+μμ0μ0+ρ+α1r/μϵ1+μ0+ψ1+ρ−γ1k1μ/α1r+μμ0+α1γ1k1r/α1r+μμ0μ0+ρ+α1r/μ.So, the eigen values are(58)ρ1=−μ0−α1rμ,ρ2=α1+γ1k1μ0−μ0−ψ1−ϵ1−α1γ1k1rα1r+μμ0μ0+α1r/μ−α1k1μ0=μ0+ψ1+ϵ1γ1k1μ0μ0+ψ1+ϵ1−1−α1γ1k1rα1r+μμ0μ0+α1r/μ−α1k1μ0=μ0+ψ1+ϵ1R0−1−α1γ1k1rα1r+μμ0μ0+α1r/μ−α1k1μ0,ρ3=α1γ1k1rα1r+μμ0μ0+α1r/μα1+γ1k1μ0ϵ1+μ0+ψ1−1−α1ηψ1rμ0+η=α1γ1k1rα1r+μμ0μ0+α1r/μR0−1−α1ηψ1rμ0+η,ρ4=−μ.The eigen valuesρ1 and ρ4 are negative, and ρ2 and ρ3 will be negative if R0<1.So,E13 is locally asymptotically stable if R0<1, and if R0>1, then it is unstable. Hence, Theorem 3 is proved.Figure8 represents the stability of the equilibrium point E13. Here, the parameters values are k1=0.01,α1=0.003,γ1=0.0018,μ0=0.00561,η=0.0092,ε1=0.001,ψ1=0.93,r=10,μ=0.09.In this case, neglecting infected and recovered class the equilibrium point becomes (0.03,111). For this parameter values the basic reproduction number R0is less than unity.Figure 8
Phase portrait for equilibrium point (0.03,111) whenR0<1.Figure9 shows unstable phase portrait for equilibrium point E13. Here, the parameters values are k1=10.1, α1=0.3, γ1=0.0018, μ0=0.561, η=0.0092, ε1=0.001, ψ1=0.93, r=10, and μ=0.9. In this case, neglecting infected and recovered class the equilibrium point becomes (2.6,11). For this parameter values, the basic reproduction number R0 is greater than unity. So, we can say Figures 8 and 9 justify our analytical findings for the equilibrium point E13.Theorem 4.
The equilibrium pointsE1∗=S1∗,I1∗,R∗,U∗ of the system (1)-(4) is stable, if R0>1 and unstable if R0<1.Proof.
From equation (49) the Jacobian matrix at point E1∗=S1∗,I1∗,R∗,U∗ becomes,
(59)JE1∗=−μ0−I1∗γ1−U∗α1−S1∗γ1η−S1∗α1I1∗γ1+U∗α1S1∗γ1−μ0−ψ1−ϵ10S1∗α10ψ1−η−μ00000−μ.
The characteristic equation of matrixJE1∗ is detJE1∗−ρI=0.Figure 9
Unstable phase portrait for equilibrium point (2.6,11) whenR0>1.So, we get(60)−μ0−I1∗γ1−U∗α1−ρ−S1∗γ1η−S1∗α1I1∗γ1+U∗α1S1∗γ1−μ0−ψ1−ϵ1−ρ0S1∗α10ψ1−η−μ0−ρ0000−μ−ρ=0,−μ0−I1∗γ1−U∗α1−ρ−S1∗γ1η−S1∗α10c22c23c2400c33c34000−μ−ρ=0,c22=S1∗γ1−μ0−ψ1−ρ−ϵ1−S1∗γ1I1∗γ1+U∗α1μ0+ρ+I1∗γ1+U∗α1,c23=S1∗γ1−μ0−ψ1−ρ−ϵ1+ηI1∗γ1+U∗α1μ0+ρ+I1∗γ1+U∗α1,c24=S1∗α1−S1∗α1I1∗γ1+U∗α1μ0+ρ+I1∗γ1+U∗α1,c33=−η−μ0−ρ−ψ1ϵ1+μ0+ψ1+ρ−S1∗γ1−ηI1∗γ1+U∗α1/μ0+ρ+I1∗γ1+U∗α1ϵ1+μ0+ψ1+ρ−S1∗γ1+S1∗γ1I1∗γ1+U∗α1/μ0+ρ+I1∗γ1+U∗α1,c34=ψ1S1∗α1−S1∗α1I1∗γ1+U∗α1/μ0+ρ+I1∗γ1+U∗α1ϵ1+μ0+ψ1+ρ−S1∗γ1+S1∗γ1I1∗γ1+U∗α1/μ0+ρ+I1∗γ1+U∗α1.Here, one of the eigen value isρ=−μ, and according to Routh-Hurwitz criteria the remaining roots will be negative if R0>1.So,E1∗ is stable, if R0>1 otherwise unstable. Hence, Theorem 4 is proved.


## 6. Method of Parameters Estimation

We use the least-square method to carry out the parameter estimation, which is implemented by the command fmincon, a part of the optimization toolbox in MATLAB. The least-square estimation is to find the parameter values to minimize the following objective function(61)fΘ,n=∑j=1nI1t−I∧1t2,where Θ is a parameter vector that is estimated by this method, n is the number of data points, I1t is the actual skin cancer infected person and I^1t is the number of skin cancer patient from the simulation.To estimate the parameters, we fit our model to the yearly new cases of skin cancer patient. In Figure10, yearly global cases of skin cancer incidence are represented by pink colored dash line and the cases from mathematical simulation are represented by green colored solid line.Figure 10
Yearly global cases of skin cancer incidence are represented by pink colored dash line and the cases from mathematical simulation are represented by green colored solid line.The estimated parameters are given in the result and discussion section.


## 7. Sensitivity Analysis of Model Parameters

To determine the model’s robustness to parameter values, we ran a sensitivity analysis. This will aid us in determining the parameters that have a significant impact on cancer invasion, such as the number of infected reproductions (R0). We used the normalized forward sensitivity index of a variable to a parameter approach described by Omoloye et al. [25] to do the sensitivity study. The ratio of relative change in the variable to relative change in the parameter is defined as this. When the variable is a differentiable function of the parameter, the sensitivity index can also be computed using partial derivatives.


### 7.1. Local Sensitivity Indices forR0

Definition.
The normalized forward sensitivity index of a variable,Q, that depends differentiably on a parameter, w, is defined as
(62)℘wQ=∂Q∂w×wQ.In particular, sensitivity indices of the basic reproduction numberR0 with respect to the model parameters are computed as follows:
(63)∂R0∂α1×α1R0=k1μ0ϵ1+μ0+ψ1×α1R0=α1α1+γ1,∂R0∂γ1×γ1R0=k1μ0ϵ1+μ0+ψ1×γ1R0=γ1α1+γ1,∂R0∂k1×k1R0=α1+γ1μ0ϵ1+μ0+ψ1×k1R0=+1,∂R0∂ϵ1×ϵ1R0=−k1α1+γ1μ0ϵ1+μ0+ψ12×ϵ1R0=−ϵ1ϵ1+μ0+ψ1,∂R0∂μ0×μ0R0=−k1α1+γ1ϵ1+2μ0+ψ1μ02ϵ1+μ0+ψ12×μ0R0=−ϵ1+2μ0+ψ1ϵ1+μ0+ψ1,∂R0∂ψ1×ψ1R0=−−k1α1+γ1μ0ϵ1+μ0+ψ12×ψ1R0=−ψ1ϵ1+μ0+ψ1.The sensitivity index of Table2 gives the idea about how basic reproduction number R0 changes with the changes of the model parameters. According to Table 2, 10% increase or reduction of k1 causes 10% increase or reduction the value of R0, 10% increase or reduction of α1 causes 6.2% increase or reduction the value of R0, 10% increase or reduction of γ1 causes 3.7% increase or reduction the value of R0, 10% increase or reduction of ε1 causes 0.01% reduction or increase the value of R0, 10% increase or reduction of μ0 causes 10% reduction or increase the value of R0, and 10% increase or reduction of ψ1 causes 9.9% reduction or increase the value of R0.Table 2
Sensitivity indices ofR0 on parameters for system (1)-(4).
ParametersDescriptionSensitivity indexk1Constant source rate+1α1Infection rate due to ultraviolet radiation+0.625γ1Infection rate due to other risk factor+0.375ε1Skin cancer induced death rate−0.001μ0Natural death rate−1.006ψ1Recovery rate−0.99


### 7.2. Global Sensitivity Analysis

Local sensitivity only gives some parameters which are differentially dependent to the output. But global sensitivity gives idea about the significant of all parameters on the baseline output. We observed the monotonicity of the parameters of our model to calculate the global sensitivity according to Mckay et al. [26]. Figure 11 represents the monotonicity plot.Figure 11
Figure of monotonicity test to the parameters of system (1)–(4) vs. infected individuals.
(a)(b)(c)(d)(e)(f)(g)(h)(i)When infected individual reaches its maximum density then using MATLAB simulation, we have found effects of model parameters on infected individual. For the details calculation techniques of sensitivity see the works of O’Hara et al. [27]. The sensitivity indices of infected individual with respect to model parameters are k1=0.005,α1=0.4,γ1=1.2,μ0=−0.027,η=0.0001,ε1=−0.001,ψ1=−1.2,r=0.4,μ=−0.04. Figure 12 shows the effects of parameters variation on R0. Figure 13 represents the Tornado plot of the sensitivity. From Figure 13, we can conclude that k1,α1,γ1,η,andr are positively related to the density of infected individuals whereas ε1,ψ1,μ,andμ0 are negatively related. Among all parameters γ1,ψ1,α1,andr are most sensitive.Figure 12
Effects of parameters variation onR0.Figure 13
Effects of model parameters on infected individuals when the density of infected individual is maximum after 3 years.


## 7.1. Local Sensitivity Indices forR0

Definition.
The normalized forward sensitivity index of a variable,Q, that depends differentiably on a parameter, w, is defined as
(62)℘wQ=∂Q∂w×wQ.In particular, sensitivity indices of the basic reproduction numberR0 with respect to the model parameters are computed as follows:
(63)∂R0∂α1×α1R0=k1μ0ϵ1+μ0+ψ1×α1R0=α1α1+γ1,∂R0∂γ1×γ1R0=k1μ0ϵ1+μ0+ψ1×γ1R0=γ1α1+γ1,∂R0∂k1×k1R0=α1+γ1μ0ϵ1+μ0+ψ1×k1R0=+1,∂R0∂ϵ1×ϵ1R0=−k1α1+γ1μ0ϵ1+μ0+ψ12×ϵ1R0=−ϵ1ϵ1+μ0+ψ1,∂R0∂μ0×μ0R0=−k1α1+γ1ϵ1+2μ0+ψ1μ02ϵ1+μ0+ψ12×μ0R0=−ϵ1+2μ0+ψ1ϵ1+μ0+ψ1,∂R0∂ψ1×ψ1R0=−−k1α1+γ1μ0ϵ1+μ0+ψ12×ψ1R0=−ψ1ϵ1+μ0+ψ1.The sensitivity index of Table2 gives the idea about how basic reproduction number R0 changes with the changes of the model parameters. According to Table 2, 10% increase or reduction of k1 causes 10% increase or reduction the value of R0, 10% increase or reduction of α1 causes 6.2% increase or reduction the value of R0, 10% increase or reduction of γ1 causes 3.7% increase or reduction the value of R0, 10% increase or reduction of ε1 causes 0.01% reduction or increase the value of R0, 10% increase or reduction of μ0 causes 10% reduction or increase the value of R0, and 10% increase or reduction of ψ1 causes 9.9% reduction or increase the value of R0.Table 2
Sensitivity indices ofR0 on parameters for system (1)-(4).
ParametersDescriptionSensitivity indexk1Constant source rate+1α1Infection rate due to ultraviolet radiation+0.625γ1Infection rate due to other risk factor+0.375ε1Skin cancer induced death rate−0.001μ0Natural death rate−1.006ψ1Recovery rate−0.99


## 7.2. Global Sensitivity Analysis

Local sensitivity only gives some parameters which are differentially dependent to the output. But global sensitivity gives idea about the significant of all parameters on the baseline output. We observed the monotonicity of the parameters of our model to calculate the global sensitivity according to Mckay et al. [26]. Figure 11 represents the monotonicity plot.Figure 11
Figure of monotonicity test to the parameters of system (1)–(4) vs. infected individuals.
(a)(b)(c)(d)(e)(f)(g)(h)(i)When infected individual reaches its maximum density then using MATLAB simulation, we have found effects of model parameters on infected individual. For the details calculation techniques of sensitivity see the works of O’Hara et al. [27]. The sensitivity indices of infected individual with respect to model parameters are k1=0.005,α1=0.4,γ1=1.2,μ0=−0.027,η=0.0001,ε1=−0.001,ψ1=−1.2,r=0.4,μ=−0.04. Figure 12 shows the effects of parameters variation on R0. Figure 13 represents the Tornado plot of the sensitivity. From Figure 13, we can conclude that k1,α1,γ1,η,andr are positively related to the density of infected individuals whereas ε1,ψ1,μ,andμ0 are negatively related. Among all parameters γ1,ψ1,α1,andr are most sensitive.Figure 12
Effects of parameters variation onR0.Figure 13
Effects of model parameters on infected individuals when the density of infected individual is maximum after 3 years.


## 8. Numerical Results and Discussions

In this section, we have executed the numerical simulations of system (1)-(4) using ODE45-solver in MATLAB programming to verify our analytical findings. To solve the epidemic model (1)–(4), we consider the initial values asS10=1000,I10=2,R0=1,U0=2, and all the values of the parameters estimated from statistical data that are given in Table 3. The values of parameters are given in Table 3. We perform simulations for the fixed final time 10 years.Table 3
Parameters and values used in model (1)–(4).
ParameterDescriptionValuek1Source rate of susceptible individuals1.1α1Infection rate for the ultraviolet radiation0.003γ1Transmission rate0.0018μ0Natural death rate0.00561ηRelapse rate0.0092ε1Skin cancer induced death rate0.001ψ1Recovery rate0.93rConstant infusing rate of index of ultraviolet radiation10μUV ray’s decay rate0.09From local sensitivity, we getα1andγ1which are the most sensitive except constant source rate and from global sensitivity analysis we obtained, α1,γ1,r,andψ1 are most sensitive parameters. For this reason, we have observed the effect of these parameters on different compartments in Figures 14–25. Figures 26 and 27 represent the solution trajectories of the system (1)–(4).Figure 14
Transformation of susceptible individuals for varied infection rateα1.Figure 15
Transformation of infected individuals for varied infection rateα1.Figure 16
Transformation of recovered individuals for varied infection rateα1.Figure 17
Transformation of susceptible individuals for constant injecting rate of UV indexr with various values.Figure 18
Transformation of infected individuals for constant injecting rate of UV indexr with various values.Figure 19
Transformation of recovered individuals for constant injecting rate of UV indexr with various values.Figure 20
Transformation of susceptible individuals for recovery rateψ1 with various values.Figure 21
Transformation of infected individuals for recovery rateψ1 with various values.Figure 22
Transformation of recovered individuals for recovery rateψ1 with various values.Figure 23
Transformation of susceptible individuals for transmission rateγ1 with various values.Figure 24
Transformation of infected individuals for transmission rateγ1 with various values.Figure 25
Transformation of recovered individuals for transmission rateγ1 with various values.Figure 26
Numerical solution of system (1)-(4).
(a)(b)(c)(d)Figure 27
Solution trajectory of susceptible, infected and recovered individual of system (1)-(4).In Figures28–31, we have asserted the effects of different parameters values on basic reproduction number R0.Figure 28
Variation of different individuals whenR0>1.Figure 29
Variation of different individuals whenR0<1.Figure 30
The impact ofα1onR0.Figure 31
The impact ofμ0onR0.Figure26 shows the solution trajectories of the system (1)-(4). The number of susceptible individuals is diminishing at a significant rate with time which is shown in Figure 26(a). The number of infected individuals is rising at a significant rate from starting time to first 4 years and reaches its peak point after 3.5 years. About 3.5 to 4 years later, the number of infected individuals is becoming reducing gradually which is shown in Figure 26(b). The number of recovered individuals is consecutively rising at a significant rate with time which is shown in Figure 26(c). The amount of UV ray’s index is successively increasing at an effective rate with time which is shown in Figure 26(d).Figure27 represents the combined disease behavior of susceptible, infected, and recovered individual with time in a one window, where green line represents the number of susceptible, red line represents the number of infected, and blue line represents the number of recovered individuals.Figure14 displays when the values of α1 fall and other parameters remain constant, the number of persons who are susceptible rises. The saffron-colored line represents the number of susceptible when α1=0.003, the red-colored line represents the number of susceptible when α1=0.0025, and the blue-colored line represents the number of susceptible when α1=0.002..Figure15 reveals when the values of α1 fall and other parameters remain constant, the number of persons who are infected reduces from starting time to around 4.5 or 5 years, after 5 years the number of infected individual rises when α1 fall. The saffron-colored line represents the number of infected when α1=0.003, the red-colored line represents the number of infected when α1=0.0025, and the blue-colored line represents the number of infected when α1=0.002.Figure16 represents with the values of α1 increases and other parameters remains constant, the number of persons who are recovered gradually decreases. This is caused due to lack of treatment options, limitation of money, and many reasons. The saffron-colored line represents the number of recovered when α1=0.003, the red-colored line represents the number of recovered when α1=0.0025, and the blue-colored line represents the number of recovered when α1=0.002.Figure17 llustrates that when r increases, the number of susceptible individuals decreases gradually.Figure18 demonstrates that as r grows, so does the number of infected individuals unto 4.8 years, and then gradually decrease.The number of recovered individuals steadily reduces asr increases, as seen in Figure 19.Figure20 represents the effects of recovery rate on susceptible individual. Due to skin cancer people loses their immunity. As a result, recovered individuals again become susceptible. This figure gives idea about the reinfection of skin cancer.Figure21 demonstrates that when the recovery rate ψ1 improves, the number of infected people decreases.Figure22 represents with the values of ψ1 increases and other parameters remains constant, the number of persons who are recovered gradually increases. The saffron-colored line represents the number of recovered when ψ1=0.93, the red-colored line represents the number of recovered when ψ1=0.73, and the blue-colored line represents the number of recovered when ψ1=0.53..Figure23 shows that when the value of γ1 decreases while the other parameters stay constant, the number of people who are susceptible increases. The number of susceptible when γ1=0.0018 is shown by the saffron-colored line, the number of susceptible when γ1=0.0016 is represented by the red-colored line, and the number of susceptible when γ1=0.0014 is represented by the blue-colored line.Figure24 shows that when the value of γ1 decreases while the other parameters stay constant, the number of infected individuals decreases from the beginning to roughly 4.5 or 5 years, and then rises after 5 years when γ1 decreases. In Figure 24, the number of infected when γ1=0.0018 is shown by the saffron-colored line, the number of infected when γ1=0.0016 is represented by the red-colored line, and the number of infected when γ1=0.0014 is represented by the blue-colored line.In Figure25, we see that the number of recovered is decreasing with the increases of transmission rate.Figures28 and 29 represent the stable and unstable conditions of skin cancer whenR0>1andR0<1, respectively.Figure30 gives the relation between α1andR0. This figure shows R0 is increasing proportional to α1..Figure31 gives the relation between μ0andR0. From this figure, we can say, R0 is increasing with μ0 inversely. Figures 28–31 give the idea about how the basic reproduction number is changes for different values of parameters which are associated to this threshold number. So, these figures will help us to take decision about control strategy.


## 9. Bifurcation Analysis

Bifurcation analysis is very important to understand the dynamical behavior of a mathematical model. There are many types of bifurcation to understand the model behavior. We have studied three types of bifurcation. These are (i) Transcritical bifurcation, (ii) Hopf bifurcation, and (iii) Saddle-node bifurcation. These types of bifurcation have significant impact on epidemiology. We have observed Transcritical bifurcation considering the effect of recovery rateψ1 for a set of parameters which is considered as Transcritical bifurcating parameters. From Transcritical bifurcation, we obtained a critical value of recovery rate ψ1∗. When recovery rate will be above this critical point, the disease will eradicate from the system (1)–(4) and when recovery rate will be below this critical point the disease will continue to exist. We have studied the Hopf bifurcation considering the infection rate due to UV radiationα1. This gives the idea about stability of endemic equilibrium at the critical values of α1. We also studied the behavior of saddle-node bifurcation for a set of parameters.


### 9.1. Transcritical Bifurcation

In this section, we have observed the behavior of the Transcritical bifurcation using Theorem5.Theorem 5.
The system (1)–(4) experience a Transcritical bifurcation at DFE pointE11 when
(64)ψ1=μ0+ϵ1−γ1k1μ0.Proof.
Suppose our system is represented byHS1,I1,R,U represent the system (1)–(4). From equation (16) at point E11=k1/μ0,0,0,0, the Jacobian matrix becomes
(65)JE11=−μ0−γ1k1μ0η−α1k1μ00γ1k1μ0−μ0−ψ1−ϵ10α1k1μ00ψ1−η−μ00000−μ.Puttingψ1=μ0+ϵ1−γ1k1/μ0 from equation (65) we obtained
(66)JE11=−μ0−γ1k1μ0η−α1k1μ0000α1k1μ00ψ1−η−μ00000−μ.One of the eigen value of matrixJE11 is zero. Consider Corresponding zero eigen value, the eigen vectors of JE11andJE11T are XandY, respectively. After numerical YTHψ1E11;ψ1∗=0,YTDHψ1E11;ψ1∗X≠0, and YTD2Hψ1E11;ψ1∗X,X≠0.So, according to Sotomayor’s theorem [28, 29], at DFE point E11, the system (1)–(4) experience Transcritical bifurcation at critical value ψ1∗=μ0+ϵ1−γ1k1/μ0.After numerical simulation we getψ1∗=1.9.


### 9.2. Hopf Bifurcation

Hopf bifurcation is very important to understand the dynamic behavior of endemic equilibrium. To understand the Hopf bifurcation, we used Theorem6.Theorem 6.
A necessary and sufficient condition for Hopf bifurcation of the system (1)–(4) atE1∗=S1∗,I1∗,R∗,U∗ is that there must be a parameter α1=α1∗ such that
(67)iReρα1∗=0anddReρα1dα1α1=α1∗≠0.where ρ is the complex roots of Jacobian matrix of system (1)–(4) at E1∗=S1∗,I1∗,R∗,U∗.Proof.
For observing Hopf bifurcation we considerα1 as bifurcation parameters. From equation (49), the Jacobian matrix at point E1∗=S1∗,I1∗,R∗,U∗ becomes
(68)JE1∗=−μ0−I1∗γ1−U∗α1−S1∗γ1η−S1∗α1I1∗γ1+U∗α1S1∗γ1−μ0−ψ1−ϵ10S1∗α10ψ1−η−μ00000−μ.
The characteristic equation of matrixJE1∗ is detJE1∗−ρI=0.After MATLAB simulation, we obtain forα1=0.01, one pair of eigen values of (68) is imaginary. According to Saha et al. [29] which can be written as
(69)iReρα1∗=0anddReρα1dα1α1=α1∗≠0.So, Hopf bifurcation occurs at pointE1∗=S1∗,I1∗,R∗,U∗ of system (1)-(4) for α1=0.01.Hence Theorem6 is proved.Figure32 shows the bifurcation diagram for Transcritical bifurcation. Transcritical bifurcation is very significant for verify disease behavior.Figure 32
For Transcritical bifurcating parametersψ1, bifurcation diagram with respect to infected individuals.Figure33 represents the Transcritical bifurcation for the parameters set k1=1,α1=0.0001,γ1=0.01,μ0=0.005,η=0.003,ε1=0.08,r=0,andμ=0.09. From Figure 33, we can say, when recovery rate greater than 1.9 then the disease of skin cancer can be eliminated and when recovery rate less than 1.9, then the disease of skin cancer cannot be eradicated.Figure 33
Phase portrait for Transcritical bifurcation (a) and (b) forψ1=0.9 (c) for ψ1=1.9 and (d) for ψ1=2.
(a)(b)(c)(d)Figure34 reveals that at the set of Hopf bifurcation parameters when infection rate due to UV rays is less than α1=0.01, proper control may eradicate the existence of disease.Figure 34
Diagram of Hopf bifurcation.Figure35 represents the Hopf bifurcation for the parameters set k1=1,γ1=0.01,μ0=0.05,η=0.3,ε1=0.08,ψ1=0.9,r=5, and μ=0.9. Figure 35 gives stable limit cycles for α1=0.01andα1=0.1 as well as unstable phase portrait for α1=0.001.Figure 35
Phase portrait for Hopf bifurcation with different values of infection rateα1.
(a)
α1=0.001(b)
α1=0.01(c)
α1=0.1


## 9.1. Transcritical Bifurcation

In this section, we have observed the behavior of the Transcritical bifurcation using Theorem5.Theorem 5.
The system (1)–(4) experience a Transcritical bifurcation at DFE pointE11 when
(64)ψ1=μ0+ϵ1−γ1k1μ0.Proof.
Suppose our system is represented byHS1,I1,R,U represent the system (1)–(4). From equation (16) at point E11=k1/μ0,0,0,0, the Jacobian matrix becomes
(65)JE11=−μ0−γ1k1μ0η−α1k1μ00γ1k1μ0−μ0−ψ1−ϵ10α1k1μ00ψ1−η−μ00000−μ.Puttingψ1=μ0+ϵ1−γ1k1/μ0 from equation (65) we obtained
(66)JE11=−μ0−γ1k1μ0η−α1k1μ0000α1k1μ00ψ1−η−μ00000−μ.One of the eigen value of matrixJE11 is zero. Consider Corresponding zero eigen value, the eigen vectors of JE11andJE11T are XandY, respectively. After numerical YTHψ1E11;ψ1∗=0,YTDHψ1E11;ψ1∗X≠0, and YTD2Hψ1E11;ψ1∗X,X≠0.So, according to Sotomayor’s theorem [28, 29], at DFE point E11, the system (1)–(4) experience Transcritical bifurcation at critical value ψ1∗=μ0+ϵ1−γ1k1/μ0.After numerical simulation we getψ1∗=1.9.


## 9.2. Hopf Bifurcation

Hopf bifurcation is very important to understand the dynamic behavior of endemic equilibrium. To understand the Hopf bifurcation, we used Theorem6.Theorem 6.
A necessary and sufficient condition for Hopf bifurcation of the system (1)–(4) atE1∗=S1∗,I1∗,R∗,U∗ is that there must be a parameter α1=α1∗ such that
(67)iReρα1∗=0anddReρα1dα1α1=α1∗≠0.where ρ is the complex roots of Jacobian matrix of system (1)–(4) at E1∗=S1∗,I1∗,R∗,U∗.Proof.
For observing Hopf bifurcation we considerα1 as bifurcation parameters. From equation (49), the Jacobian matrix at point E1∗=S1∗,I1∗,R∗,U∗ becomes
(68)JE1∗=−μ0−I1∗γ1−U∗α1−S1∗γ1η−S1∗α1I1∗γ1+U∗α1S1∗γ1−μ0−ψ1−ϵ10S1∗α10ψ1−η−μ00000−μ.
The characteristic equation of matrixJE1∗ is detJE1∗−ρI=0.After MATLAB simulation, we obtain forα1=0.01, one pair of eigen values of (68) is imaginary. According to Saha et al. [29] which can be written as
(69)iReρα1∗=0anddReρα1dα1α1=α1∗≠0.So, Hopf bifurcation occurs at pointE1∗=S1∗,I1∗,R∗,U∗ of system (1)-(4) for α1=0.01.Hence Theorem6 is proved.Figure32 shows the bifurcation diagram for Transcritical bifurcation. Transcritical bifurcation is very significant for verify disease behavior.Figure 32
For Transcritical bifurcating parametersψ1, bifurcation diagram with respect to infected individuals.Figure33 represents the Transcritical bifurcation for the parameters set k1=1,α1=0.0001,γ1=0.01,μ0=0.005,η=0.003,ε1=0.08,r=0,andμ=0.09. From Figure 33, we can say, when recovery rate greater than 1.9 then the disease of skin cancer can be eliminated and when recovery rate less than 1.9, then the disease of skin cancer cannot be eradicated.Figure 33
Phase portrait for Transcritical bifurcation (a) and (b) forψ1=0.9 (c) for ψ1=1.9 and (d) for ψ1=2.
(a)(b)(c)(d)Figure34 reveals that at the set of Hopf bifurcation parameters when infection rate due to UV rays is less than α1=0.01, proper control may eradicate the existence of disease.Figure 34
Diagram of Hopf bifurcation.Figure35 represents the Hopf bifurcation for the parameters set k1=1,γ1=0.01,μ0=0.05,η=0.3,ε1=0.08,ψ1=0.9,r=5, and μ=0.9. Figure 35 gives stable limit cycles for α1=0.01andα1=0.1 as well as unstable phase portrait for α1=0.001.Figure 35
Phase portrait for Hopf bifurcation with different values of infection rateα1.
(a)
α1=0.001(b)
α1=0.01(c)
α1=0.1


## 10. Conclusion

In this study we worked on mathematical model of skin cancer. Mainly by our study, we have tried to answer the four key questions: (1) how UV radiation is related to skin cancer? (2) Are our models valid? (3) Which parameters are most significant to causes skin cancer? (4) What values of parameters give idea about how the disease can be minimized?We used next generation matrix method to find out basic reproduction number which is very important for a disease model. Basic reproduction number of the model for the baseline parameters isR0=1.005>1. We have checked the stability of that equilibrium points analytically and numerically to prove the existence of the model. The global and local sensitivity analysis of parameters gives that infection rate α1, transmission rate γ1, constant infusing rate of UV rays r, and the recovery rate ψ1 are the most sensitive parameters. We have showed these parameters effects to our model outputs in Figures 14–31.The bifurcation analysis gives some critical values of parameters in which our system (1)–(4) experience Transcritical and Hopf bifurcation. For Transcritical bifurcation the critical value isψ1=1.9, and for the Hopf bifurcation the critical value is α1=0.01 which we obtained from the close look of Figures 32–35.


## 11. Future Work

We would like to apply optimal control measures (awareness and treatment) to reduce the transmission dynamics of skin cancer. Therefore, we want to extend our work in future to observe the above cases by introducing a specific objective function.



---

*Source: 5445281-2022-12-07.xml*