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 if and if , then it is unstable. Again, the endemic equilibrium point is stable, if and unstable if . 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 , 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.
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).
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.
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 . 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 . 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 . Those who have survived from skin cancer and are immune to it are denoted by . This compartmental model’s flow chart is provided in Figure 7.
According to the flowchart of the model in Figure 7, the mathematical model of skin cancer can be written in the form of following nonlinear system of ordinary differential equations:
With initial conditions , , , and .
In Table 1, we have described the parameters of our model.
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 .
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 is . So, we can write .
Then, .
From eqn. (1)-(3) of the system (1)–(4), we obtain,
Here, integrating factor, .
Multiplying both sides by , we get
Integrating both sides, we get
Using initial condition, we get
At , .
Now, from eqn. (4) of the system (1)–(4), we obtain
Here, integrating factor, .
Multiplying both sides by , we get
Integrating both sides, we get
Using initial condition, we get,
At , .
Assume
So,
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 compartments , , , and . To investigate the positivity of these compartments we use lemma 2.
Lemma 2. If and , then the solutions of the system (1)-(4) are positively invariant.
Proof. From equation (1), we get, Here, integrating factor, .
Multiplying both sides by , we get
Integrating both sides, we get
Using initial condition, we get,
And also at .
Similarly, we can verify the positivity of , , and under the initial conditions. Therefore, the solutions 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
Thus, we have
4.1.1. Ultraviolet Free Equilibrium Points
For this case , 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.
and .
Here, Complete disease and risk of disease-free equilibrium is .
And endemic equilibrium point is .
Here,
4.1.2. Disease Free Equilibrium Point When
At disease free equilibrium point (DFE), . Thus, the system (21)–(24) reduces to which implies
So, disease-free equilibrium point (DFE) is as follows:
4.1.3. Endemic Equilibrium Points When
If all populations exist, the system (1)–(4) present endemic equilibrium (EE) point given by .
Equation (24) gives .
From equation (23), we get
Putting the values of from equation (22), we obtain
Then equation (21) reduces to
Here,
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)(2)(3)(4)
Using MATLAB, we get two positive endemic equilibrium points, where represent the number of susceptible, infected, and recovered individuals and the last one is the index of ultraviolet radiation.
So, the endemic equilibrium point is
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 . The epidemiological definition of 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 , the occurrence of the disease will increase. If , the occurrence of the disease will decrease and the disease will ultimately be eliminated. When , 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 and are filled with appropriate terms from the infected class equations using this method. Terms that describe the appearance of new illnesses belong in the category, while terms that describe the transmission of existing infections are in the category and should be avoided. The matrices and are constructed and evaluated at a nontrivial disease-free equilibrium using the Jacobian matrices generated by differentiating and 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, for ,, and .
So, we obtain
Thus,
So, for , and .
So, we obtain
Thus,
So, the basic reproduction number for our total model is .
So, at disease free equilibrium point
For the base line parameters of our system the basic reproduction number is .
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 points of the system (1)-(4) are locally asymptotically stable if and if , then it is unstable.
Proof. Let Thus, we have
The Jacobian matrix of the system (44)–(46) can be written as
So that
At point the Jacobian matrix becomes
The characteristic equation of matrix is .
So, we get
So, the eigen values are The eigen values are negative, and will be negative if .
So, is stable, otherwise unstable. Hence, Theorem 1 is proved.
Theorem 2. The equilibrium points of the system (1)-(4) are stable, if otherwise unstable.
Proof. From equation (43) the Jacobian matrix at point becomes The characteristic equation of matrix is .
So, we get
So, the eigen values are
Here, is negative.
The eigen values will be negative, if .
So, is stable, if otherwise unstable. Hence, Theorem 2 is proved.
Theorem 3. The equilibrium points of the system (1)-(4) are locally asymptotically stable if and if , then it is unstable.
Proof. From equation (49) at point the Jacobian matrix becomes The characteristic equation of matrix is .
So, we get
So, the eigen values are
The eigen values and are negative, and and will be negative if .
So, is locally asymptotically stable if , and if , then it is unstable. Hence, Theorem 3 is proved.
Figure 8 represents the stability of the equilibrium point . Here, the parameters values are In this case, neglecting infected and recovered class the equilibrium point becomes (0.03,111). For this parameter values the basic reproduction number is less than unity.
Figure 9 shows unstable phase portrait for equilibrium point . Here, the parameters values are , , , , , , , , and . In this case, neglecting infected and recovered class the equilibrium point becomes (2.6,11). For this parameter values, the basic reproduction number is greater than unity. So, we can say Figures 8 and 9 justify our analytical findings for the equilibrium point .
Theorem 4. The equilibrium points of the system (1)-(4) is stable, if and unstable if .
Proof. From equation (49) the Jacobian matrix at point becomes, The characteristic equation of matrix is .
So, we get
Here, one of the eigen value is , and according to Routh-Hurwitz criteria the remaining roots will be negative if .
So, is stable, if 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 where is a parameter vector that is estimated by this method, is the number of data points, is the actual skin cancer infected person and 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 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 (). 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 for
Definition. The normalized forward sensitivity index of a variable, , that depends differentiably on a parameter, , is defined as
In particular, sensitivity indices of the basic reproduction number with respect to the model parameters are computed as follows:
The sensitivity index of Table 2 gives the idea about how basic reproduction number changes with the changes of the model parameters. According to Table 2, 10% increase or reduction of causes 10% increase or reduction the value of , 10% increase or reduction of causes 6.2% increase or reduction the value of , 10% increase or reduction of causes 3.7% increase or reduction the value of , 10% increase or reduction of causes 0.01% reduction or increase the value of , 10% increase or reduction of causes 10% reduction or increase the value of , and 10% increase or reduction of causes 9.9% reduction or increase the value of .
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.
(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 Figure 12 shows the effects of parameters variation on . Figure 13 represents the Tornado plot of the sensitivity. From Figure 13, we can conclude that are positively related to the density of infected individuals whereas are negatively related. Among all parameters are most sensitive.
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 as 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.
From local sensitivity, we get which are the most sensitive except constant source rate and from global sensitivity analysis we obtained, 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).
(a)
(b)
(c)
(d)
In Figures 28–31, we have asserted the effects of different parameters values on basic reproduction number .
Figure 26 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).
Figure 27 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.
Figure 14 displays when the values of fall and other parameters remain constant, the number of persons who are susceptible rises. The saffron-colored line represents the number of susceptible when , the red-colored line represents the number of susceptible when , and the blue-colored line represents the number of susceptible when .
Figure 15 reveals when the values of 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 fall. The saffron-colored line represents the number of infected when , the red-colored line represents the number of infected when , and the blue-colored line represents the number of infected when
Figure 16 represents with the values of 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 , the red-colored line represents the number of recovered when , and the blue-colored line represents the number of recovered when
Figure 17 llustrates that when increases, the number of susceptible individuals decreases gradually.
Figure 18 demonstrates that as grows, so does the number of infected individuals unto 4.8 years, and then gradually decrease.
The number of recovered individuals steadily reduces as increases, as seen in Figure 19.
Figure 20 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.
Figure 21 demonstrates that when the recovery rate improves, the number of infected people decreases.
Figure 22 represents with the values of 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 , the red-colored line represents the number of recovered when , and the blue-colored line represents the number of recovered when .
Figure 23 shows that when the value of decreases while the other parameters stay constant, the number of people who are susceptible increases. The number of susceptible when is shown by the saffron-colored line, the number of susceptible when is represented by the red-colored line, and the number of susceptible when is represented by the blue-colored line.
Figure 24 shows that when the value of 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 decreases. In Figure 24, the number of infected when is shown by the saffron-colored line, the number of infected when is represented by the red-colored line, and the number of infected when is represented by the blue-colored line.
In Figure 25, we see that the number of recovered is decreasing with the increases of transmission rate.
Figures 28 and 29 represent the stable and unstable conditions of skin cancer when, respectively.
Figure 30 gives the relation between This figure shows is increasing proportional to .
Figure 31 gives the relation between From this figure, we can say, is increasing with 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 for a set of parameters which is considered as Transcritical bifurcating parameters. From Transcritical bifurcation, we obtained a critical value of recovery rate . 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. This gives the idea about stability of endemic equilibrium at the critical values of . 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 Theorem 5.
Theorem 5. The system (1)–(4) experience a Transcritical bifurcation at DFE point when
Proof. Suppose our system is represented by represent the system (1)–(4). From equation (16) at point , the Jacobian matrix becomes
Putting from equation (65) we obtained
One of the eigen value of matrix is zero. Consider Corresponding zero eigen value, the eigen vectors of are , respectively. After numerical and
So, according to Sotomayor’s theorem [28, 29], at DFE point , the system (1)–(4) experience Transcritical bifurcation at critical value .
After numerical simulation we get .
9.2. Hopf Bifurcation
Hopf bifurcation is very important to understand the dynamic behavior of endemic equilibrium. To understand the Hopf bifurcation, we used Theorem 6.
Theorem 6. A necessary and sufficient condition for Hopf bifurcation of the system (1)–(4) at is that there must be a parameter such that where is the complex roots of Jacobian matrix of system (1)–(4) at .
Proof. For observing Hopf bifurcation we consider as bifurcation parameters. From equation (49), the Jacobian matrix at point becomes The characteristic equation of matrix is .
After MATLAB simulation, we obtain for , one pair of eigen values of (68) is imaginary. According to Saha et al. [29] which can be written as
So, Hopf bifurcation occurs at point of system (1)-(4) for .
Hence Theorem 6 is proved.
Figure 32 shows the bifurcation diagram for Transcritical bifurcation. Transcritical bifurcation is very significant for verify disease behavior.
Figure 33 represents the Transcritical bifurcation for the parameters set 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.
(a)
(b)
(c)
(d)
Figure 34 reveals that at the set of Hopf bifurcation parameters when infection rate due to UV rays is less than , proper control may eradicate the existence of disease.
Figure 35 represents the Hopf bifurcation for the parameters set , and . Figure 35 gives stable limit cycles for as well as unstable phase portrait for .
(a)
(b)
(c)
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 is . 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 , transmission rate , constant infusing rate of UV rays , and the recovery rate 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 and for the Hopf bifurcation the critical value is 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.
Data Availability
All data are available within the paper cited in the text.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The first author is supported by the M.Sc. (NST) fellowship, bearing ID: MSc-201203, Serial: 291, Merit: 77, No: 39.00.0000.012.002.06.21. Session: 2021-2022 which is provided by the Ministry of Science and Technology, Government of the People’s Republic of Bangladesh. The NST fellowship is gratefully acknowledged for financial support.