Modelling the spread of atherosclerosis considering relapse and linear treatment

Atherosclerosis is a narrowing of the arteries due to a build-up of plaque in the artery walls. One of the reasons for the spread of atherosclerosis is the widespread of improper lifestyle in a population because of the tendency to follow the improper lifestyle of people in the surrounding environment. Atherosclerosis can cause complications such as heart disease, stroke, etc. This study aims to construct a model for the spread of atherosclerosis among human population by considering relapse and linear treatment rates. The model was then analyzed analytically and numerically. Analytical study reveals the existence and local stability criteria of the equilibrium points, determine the basic reproduction numbers and investigate the existence of bifurcations of the constructed model using the center-manifold theorem. Some numerical simulations were conducted for some possible prevention scenarios in the field.


Introduction
The significant trigger for vascular disease worldwide and a global threat is atherosclerosis. Atherosclerosis is a narrowing of the arteries [1] due to a buildup of plaque on the artery walls so that it will block blood flow to the body's organs. These deposits are consist of cholesterol, fatty substances, cellular waste products, calcium, fibrin, and other materials found in the blood. Atherosclerosis increased by several risk factors, including high blood pressure, high blood cholesterol, smoking, obesity, diabetes, physical inactivity [2], alcohol consumption [3], unhealthy eating patterns, insulin resistance, increasing age. Epidemiology of atherosclerosis is quite challenging to determine because, at first, atherosclerosis disease often causes no symptoms. Symptoms appear when blood flow to the organ is obstructed and takes years.
Atherosclerosis can occur in the brain, heart, kidneys, arms, and legs. If atherosclerosis is left and not treated immediately, it will cause dangerous complications, such as ischemic heart disease and stroke. One of the leading causes of the increase in cases of patients with atherosclerosis is an unhealthy lifestyle. Based on data from the World Health Organization, complications of atherosclerosis, namely heart disease and ischemic stroke, rank first, and the second cause of global death in 2016 [4]. WHO also mentions cardiovascular disease (CVD) is estimated to be the cause of death of 17.9 million people every year [5]. CVD is a group of heart and blood vessel disorders and includes coronary heart disease, cerebrovascular disease, rheumatic heart disease, and other conditions. Heart attacks and strokes cause 85% of deaths due to CVD, and one-third of these deaths occur prematurely in people under 70 years [5]. The most appropriate treatment for atherosclerosis involves three things: a healthy lifestyle, medication, and medical procedures. However, this cannot cure a person of atherosclerosis, but Several mathematical models have been introduced by many authors to understand some types of self-harm diseases, such as obesity [6,7], smoking [8], diabetes [9,10], and many more. These articles model the spread of bad habit, which caused several diseases based on a wellknown SIR compartmental model. This SIR model (and it modification) is already used by many authors to model transmissible diseases [11,12,13,14]. On the mentioned references, they model the spread of disease through social contact between susceptible and infected individuals, where infected individuals can make a persuasive approach to make susceptible individuals follow a bad habit, leading to the related disease.
Using a similar approach for the transmission mechanism, this article introduced a novel mathematical model for the transmission of atherosclerosis among the human population. The infection occurs during social contact between susceptible and infected individuals. We include three important factors in our model, namely the treatment, reinfection, and relapse. Mathematical analysis regarding the stability of the equilibrium and the sensitivity analysis on the basic reproduction number is our analysis's focus. We obtain that our model may exhibit a backward bifurcation phenomenon. Furthermore, we also find the sensitivity analysis that the treatment rate is promising to reduce the spread of atherosclerosis significantly.
This manuscript presented as follows. We introduce our model in section 2, followed by mathematical analysis regarding the equilibrium points and the basic reproduction number in section 3. Bifurcation analysis is given in section 4 and followed with some numerical experiments in section 5. Finally, some concluding remarks and conclusions are given in section 6.

Model construction
In this section, we model the spread of atherosclerosis disease that focused on social interaction habits with individuals who have unhealthy lifestyles. The population is divided into three compartments, namely susceptible (S), untreated infected (I) and treated infected individuals (I t ) . We consider the following assumptions to construct our model in system (1). A1. Population is closed. Migrations are neglected. A2. Number of newborn is the same with the number of natural death. A3. All newborns are susceptible to atherosclerosis. A4. Individuals get infected with atherosclerosis due to intensive social interaction with infected individuals, which leads to unhealthy lifestyles and finally to the disease. A5. People with atherosclerosis cannot recover. A6. Reinfection occurs when atherosclerosis individuals undergoing treatment conduct an unhealthy lifestyle caused by the habit of social interactions with people with atherosclerosis. A7. Relapse occurs when treated infected individuals return to an unhealthy lifestyle. A8. Only infected individuals without treatment can exhibit death due to disease.
Based on the above assumptions and transmission diagram in Figure 1, a mathematical model of atherosclerosis is given by  Figure 1: The transmission diagram of atherosclerosis considering relapse and treatment with a non-negative initial conditions: Due to the assumption A2, A = µ N then dN dt = −γ I < 0, which indicate that the total of human population is decreasing. The description of each parameters in system (1) is given in Table 1. Table 1: Description of parameters in system (1).

Parameters Description Condition Unit
A Number of newborn per unit time.
The rate of social interaction success to infect susceptible individual The rate of social interaction success to re-infect treated infected individual Death rate induced by atherosclerosis Natural death rate. µ > 0 1 time Before we analyzed system (1), it is essential to guarantee the positivity and boundedness of our model. The following theorem state that our system always has a positive solution for all t > 0 as long as the initial condition is non-negative. Theorem 1. Solution of system (1) is always non-negative, i.e S(t) ≥ 0, I(t) ≥ 0, I t ≥ 0, for all t > 0 with a non-negative initial condition. Furthermore, our system is attracted to a positive invariant set : is non-negative over the boundary field of R 3 + . Therefore, the direction of the vector field tends to be inward from its boundary. As a result, by having a non-negative initial condition, all solutions of system (1) remain positive for all time t > 0.
Furthermore, adding all equation of (1) gives Hence, we have that dN dt ≤ A−µN . This result tells us that even though the number of newborns and the number of natural death is the same, due to death by disease make the number of the human population becoming not constant. Solving the last equation respect to N (t), we have that lim t→∞ N ≤ A µ .

Model Analysis
Using the next generation matrix (NGM) approach [15], model (1) has a basic reproduction number given by In many mathematical models [16,17,18,19,20,21,22,23,24], the basic reproduction number become an important threshold which determine whether the disease may exist or disappeared from the population. In our case, we state the relation of our basic reproduction number in (4) respect to the local stability of DF E in the following theorem.
Theorem 2. The disease-free equilibrium of system (1) is locally asymptotically stable if and only if R 0 < 1, and unstable if R 0 > 1.
Proof. To analyze the local stability of the DF E, we linearize system (1) at DF E. The characteristic polynomial of this linearized matrix is given by  (5), we get an eigenvalue to appear explicitly namely, λ 1 = −µ which is negative. The following quadratic equation contains the other two eigenvalues Hence, based on the Routh-Hurwitz criteria, the condition that the system (1) is said to be stable if a 1 > 0 and a 2 > 0 which implies u+2µ+α+γ and R 0 given in equation (4). Since R 1 − R 0 < 0, it can be written that R 1 < R 0 < 1. Therefore, the disease-free equilibrium point will be asymptotically stable if and only if the value of R 0 < 1. Otherwise, if R 0 > 1, then DFE is unstable.

Existence of the endemic equilibrium point
The endemic equilibrium point of system (1) can not be calculated and shown explicitly. Hence, we will show the existence of the endemic equilibrium, where S and I t present as a function of I, while I is taken from the positive roots of the polynomial which coming from backward substitution on system (1) respect to each variables.
It is easy to show that equation (7) always has a unique positive root is R 0 > 1 by using the standard analysis of a quadratic equation. For two positive roots, it is more challenging. We use a similar method which used by authors in [25,26,27,28]. We analyze the existence of the  (7) using the sign of partial derivative of I respect R 0 at R 0 = 1, since R 0 = 1 is the threshold number of our epidemic model. Taking the partial derivative ∂I ∂R 0 from equation (7) yields Substitute I = 0 on the last equation gives us .
Hence, we can conclude that ∂I ∂R 0 < 0 if and only ifr(R 0 ) < 0 which equivalent with : From above analysis, we can conclude that it is possible to have indicate the existence of other endemic equilibrium point even though R 0 < 1. We write this results in the following theorem.

Bifurcation Analysis
To perform the bifurcation analysis of system (1) in this article conducted using the centermanifold approach by implementing the Castillo-Song theorem [29]. Hence, let With this assumption, system (1) can be written as follows. For the analysis, we choose β 1 as the bifurcation parameter. Solving R 0 = 1 respect to β 1 give us Substitute β 1 = β # 1 and DF E to the Jacobian matrix of system (1), the eigenvalues of this matrix are Since λ 1 = 0 is one of the eigenvalues, and the others are negative, hence we can use the centermanifold theory by Castillo and Song [29] to analyze the bifurcation analysis of system (1). The right vector eigen of the Jacobian Matrix respect to λ = 0 is given by w = (w 1 , w 2 , w 3 ), where On the other hand, the left eigen vector v = (v 1 , v 2 , v 3 ) is given by : Since we have that v 1 = 0, we do not have to calculate the partial derivative of g 1 . Hence, we only have to calculate the partial derivative of g 2 and g 3 to calculate A and B. First, we calculate the following partial derivatives.
It is clear to see that B always positive, while A could be positive or negative. Using β 2 as the indicator to determine the existence of forward or backward bifurcation, solving B = 0 yield Hence, we have that A < 0 ⇐⇒ β 2 < β critical 2 , which indicates system (1) exhibit a forward bifurcation at R 0 = 1. On the other hand, if A > 0 ⇐⇒ β 2 > β critical 2 , then system (1) exhibit a backward bifurcation at R 0 = 1. This results summarized in the following theorem. The explanation of Figure 2 is as follows.
a. When x < x 1 , There is only one equilibrium point, i.e the disease free equilibrium, and it is locally stable. b. When x ∈ (x 1 , x 0 ), there exist a stable disease free equilibrium, and two endemic equilibrium. The small endemic equilibrium is unstable, while the bigger one is locally stable. c. When x > x 0 , there exist only one stable endemic equilibrium which locally stable, while the disease-free equilibrium become unstable. (ii) Forward bifurcation By choosing β 2 = 0.17 < β critical 2 , according to Theorem 4, forward bifurcation will arises at R 0 = 1. The results illustrated in Figure 3. The explanation of Figure 3 is as follows.
a. When x < x 0 , there is only exist one equilibrium, i.e the disease-free equilibrium. This equilibrium is locally stable. b. When x > x 0 , the the disease free equilibrium become unstable, and the endemic equilibrium arises, and it is locally stable.

Numerical experiments
In this section, we conduct two kind of numerical experiments, namely the sensitivity analysis to determine the most influential parameter to the basic reproduction number, and the autonomous simulation to predict the long-time behavior of our model respect to the change of parameters.

Sensitivity analysis
As mentioned before, the sensitivity analysis is conducted to determine the most influential parameter to the change of the basic reproductive number, locally. To conduct this, we use the local-sensitivity analysis used by authors in [30] with the following formula where ϑ is a parameter in R 0 . If the sign of E ϑ R 0 is positive, then we have that R 0 increases whenever ϑ increases. On the other hand, if E ϑ R 0 < 0, then R 0 decreases whenever ϑ increases. We use the same parameters value as in previous section, while β 1 chosen to be 0.15 to have R 0 < 1, and β 1 = 0.163 such that R 0 > 1. The results for the local sensitivity is given in Table 2.   Table 2, it can be seen that A and β 2 do not affect R 0 . On the other hand, R 0 decreases whenever µ, u or γ increases. Furthermore, R 0 increases whenever β 1 and α increases. The value of this sensitivity index represent the change of R 0 respect to the change of parameters. For an example, whenever u increases 10%, then R 0 decreases 9.9%. Hence, we can see that β 1 is the most influential parameter on R 0 , followed with u, µ, α and γ, respectively.
Next, since we have that β 1 and u is the most influential parameter to R 0 , we continue our analysis to see the contour plot of R 0 respect to β 1 and u. The results is given in Figure 4. It is easy to see how R 0 increases respect to β 1 , and decreases respect to u.

Autonomous simulation on the variation of treatment rate
The treatment rate related to the effort that given by the government, which in this case is the effort by the hospital to accommodate number of new patient of atherosclerosis. To do simulation in this section, we use the following parameters value