An SIR epidemic model for COVID-19 spread with fuzzy parameter: the case of Indonesia

The aim of this research is to construct an SIR model for COVID-19 with fuzzy parameters. The SIR model is constructed by considering the factors of vaccination, treatment, obedience in implementing health protocols, and the corona virus-load. Parameters of the infection rate, recovery rate, and death rate due to COVID-19 are constructed as a fuzzy number, and their membership functions are used in the model as fuzzy parameters. The model analysis uses the generation matrix method to obtain the basic reproduction number and the stability of the model’s equilibrium points. Simulation results show that differences in corona virus-loads will also cause differences in the transmission of COVID-19. Likewise, the factors of vaccination and obedience in implementing health protocols have the same effect in slowing or stopping the transmission of COVID-19 in Indonesia.


Introduction
A new virus that can cause an increase in pneumonia first appeared in Wuhan, China, at the beginning of December 2019. The virus was called SARS-COV-2 and the disease associated with the virus was called COVID-19 [1]. The rapid spread of the disease throughout the world has made the World Health Organization declare the COVID-19 outbreak a global pandemic on March 12, 2020 [2]. COVID-19 is transmitted by a person infected with COVID-19 through physical contact with another person, or through small droplets from the mouth of a person with COVID-19 that is touched by another person. According to the data collected by John Hopkins University, as of the beginning of November 2020, COVID-19 has spread to 219 countries in the world with the total number of cases infected with COVID-19 reaching 53,699,160, with 1,308,261 deaths, and 37,469,072 people have been declared cured.
The spread of COVID- 19 to Indonesia has been evident since the first confirmed case on March 2, 2020. COVID-19 is continuing to spread to all provinces, and as of the beginning of November 2020 it has reached a total of 457,735 infected cases, 15,037 deaths, and 385,094 people have been declared cured [3].
A number of mathematicians have performed various studies to make a model and predict the spread of COVID-19 since it is becoming a pandemic. The first study [4] to pre-dict COVID-19 in China was conducted using GLM method and Richard's model. Then [5] predicted COVID-19 in Indonesia based on early endemic data using Richard's curve. Other models and predictions using statistical approaches [6][7][8][9][10] or those using SIR, SEIR, and their extensions [11][12][13][14][15][16] have been widely performed. Moreover, some researchers used fractional order in epidemic models to model the spread of COVID-19. Fractional order derivative with Mittag-Leffler function as a nonsingular kernel type [17] and Caputo derivative [18][19][20] are used in modeling the transmission of COVID-19. Then [21,22] considered the fractal-fractional derivative in the Atangana-Baleanu sense to obtain the stability of the model, and [23] presented the existence and uniqueness solution of the model via fractal-fractional operators. However, the parameters used in any existing SIR epidemic model and its extensions employ crisp numbers, whereas uncertainty in parameters and heterogeneity in the population are very possible to occur. Therefore the use of uncertain parameters or fuzzy parameters is very important, because the model will reflect the real world problems. Supporting studies could be used as references, such as the fuzzy epidemic models for human infectious diseases [24][25][26], other models considering the uncertainty of parameter space and heterogeneity of a population [27,28], fuzzy dynamic systems [29], and dynamic behavior of an epidemic model with fuzzy transmission [30].
In this study, we consider a mathematical model of SIR in a normalized form with three control parameters, namely vaccination control, treatment control and implementation of health protocols. The parameters of infection rate, recovery rate, and death rate due to COVID-19 are treated as fuzzy numbers that depend on individual virus-load.

Method
The method used to construct the model is the SIR model by considering vaccination, treatment, and the implementation of health protocols as control parameters. The parameters of the infection rate, recovery rate, and death rate due to COVID-19 are constructed as a membership function of a fuzzy number [28]. These parameters depend on the corona virus-load in an individual and control parameters. The model analysis uses the generation matrix method to obtain the basic reproduction number and stability of the SIR model for the spread of COVID-19. The numerical simulation of the model uses data on the number of COVID-19 cases in Indonesia [3]. The parameters of vaccination effectiveness, treatment effectiveness, the level of obedience in implementing health protocols, and corona virus-load are assumed in this simulation.

SIR model of COVID-19
Consider an SIR model for COVID-19 that describes the dynamics of direct transmission of COVID-19 with interactions between suspected and infected, change from being infected to recovering, pure birth/death rates, vaccine effectiveness, treatment effectiveness, obedience in implementing health protocols, and deaths due to the COVID-19 infection. The schematic diagram of the COVID-19 transmission flow is given in Fig. 1, and the model is mathematically described as follows: where S is the proportion of susceptible individuals, I is the proportion of infected individuals, R is the proportion of recovered individuals, β is the infection rate parameter; γ is the recovery rate parameter; μ is the natural birth/death rate parameter, τ is the vaccine effectiveness parameter, θ is the treatment effectiveness parameter, π is the effectiveness of obedience in implementing health protocols, μ C is the death rate parameter due to COVID-19. Now, we can extend the SIR model (1)

The SIR fuzzy model of COVID-19 spread
Consider the SIR model for COVID-19 in (1) Let β = β( ) be the chance of transmission between a suspected and an infected individual with the amount of the corona virus-load . Some values of β are more reasonable compared to some others, and it turns β into a membership function of fuzzy numbers. To construct the membership function, we assume that if the number of corona virus-loads in an individual is relatively low, the chance of transmission is negligible, and there is a minimum corona virus-load min needed to be able to transmit to other individuals. Furthermore, there is a certain amount of corona virus-load 0 , where the transmission rate is maximum and equal to one. However, we assume that the total amount of corona virusload on a person is limited by max [28]. We can also consider that the vaccination and the discipline to follow health protocols will affect the infection rate of COVID-19. Let τ and π be the parameters representing the vaccine effectiveness and the level of discipline in implementing health protocols, respectively. Then fuzzy membership function of the infectivity contact rate is given as follows: The graphic of β( ) is given in Fig. 2. The death rate due to COVID-19 infection can also be assumed as a membership function of a fuzzy number. The function is an increasing function of corona virus-load . However, due to some reasons, such as a person infected with COVID-19 suffering from other diseases, immunity power, availability of medicine, etc., the function might not reach its maximum value equal to one. Likewise, treatment for COVID-19 will affect the death rate due to COVID-19 infection. Therefore, we assume that the maximum value of the . Thus, we can define the function μ C ( ) as follows (depicted in Fig. 3): where μ C 0 ; (0 < μ C 0 < 1) is the lowest death rate due to COVID-19 infection and θ is the treatment effectiveness.
The recovery rate of the COVID-19 infection group γ = γ ( ) is also a function of the corona virus-load . The higher the corona virus-load , the longer the recovery process will take from infection. So, γ ( ) is a decreasing function. Moreover, we can also consider where γ 0 is the lowest recovery rate.
In the fuzzy SIR model, the membership function of the infection rate β( ), the recovery rate γ ( ), and the death rate μ C ( ) due to COVID-19 infection are treated as fuzzy parameters of the model.

The equilibrium points and the basic reproduction number
There are two equilibrium points in model (4)-(6), namely the disease-free equilibrium point and the endemic equilibrium point. To determine these two equilibrium points, each of the equations in the equations must be equal to zero, that is, dS dt = 0, dI dt = 0, dan dR dt = 0 so that: θ + γ ( ) I + (π + τ )S -μR = 0, then the equilibrium points for S, I, and R are as follows.

The disease-free equilibrium point for the SIR fuzzy model
The points of equilibrium for disease free are conditions where there is no spread of COVID-19, namely I = I 0 = 0. Thus, from Eq. (4), we obtain from Eq. (6) and Eq. (13), we obtain Thus, the disease-free equilibrium point for the SIR fuzzy model (4)- (6) is
It is easy to see that λ 1 = -μ is one of the eigenvalues P 1 (λ). The other eigenvalues are the solutions of P 2 (λ) = 0. Based on the Routh-Hurwitz condition, P 2 (λ) has two roots which have a negative real part if j 1j 4 > 0 and j 2 j 3j 1 j 4 > 0.
Since the disease-free equilibrium is stable if 0 ( ) < 1 and unstable for 0 ( ) > 1, then system (4)-(6) is at a bifurcation point when 0 ( ) = 1. Let * be the bifurcation value of the system, then * is the solution of the equation , where * ≤ 0 .
In this way, we can think of * as a parameter related to the corona virus control in the sense that if a corona virus is transmitted in some number of individuals, it should be noted that is not higher than * .
Corollary The disease-free equilibrium and the endemic equilibrium of system (4)-(6) are locally asymptotically stable for < * and > * , respectively.

Numerical simulation for transmission of COVID-19 in Indonesia
Numerical simulations are carried out using the initial values for N , S, I, and R as in Table 1. The parameters β, γ , and μ C are calculated based on Eqs. (7)-(9) and the parameters       Table 2 and Table 3. The values of 0 , min , and should be determined by a virologist, while the value of parameter ϑ should be determined by a physician. For each of the controlled parameters, namely the parameters of vaccine effectiveness, level of obedience in implementing health protocols, treatment effectiveness, and corona virus-load, simulations are carried out three times for each parameter with different values, as in Table 3. Based on Table 3, it can be explained that if a treatment, vaccination, and health protocols have not been implemented, the basic reproduction number is 9.32, which means that a person with COVID-19 can infect nine other people. Meanwhile, if the effectiveness of can be controlled. Furthermore, if the effectiveness of vaccination or the effectiveness of implementing health protocols is more than 10%, then the basic reproductive number is less than 0.44, meaning that the COVID-19 outbreak will disappear in the population. The simulation results for the corona virus-load are presented in Fig. 5, Fig. 6, and Fig. 7.
Based on these figures, by giving 100 the corona virus-load, the COVID-19 outbreak will never vanish in the population, and the recovery rate tends to remain, but it is smaller than those infected with COVID-19. Meanwhile, if the virus corona-load is 50, the epidemic tends to decrease, and if = 20, the COVID-19 outbreak will vanish in the population.
The simulation results of treatment effectiveness are presented in Fig. 8, Fig. 9, and Fig. 10. From the figures, it can be seen that if treatment is only 2% effective, COVID-19 will become endemic, while if the treatment effectiveness is 47% and 85%, the COVID-19 outbreak will disappear in the population.
The simulation results of the effectiveness of vaccination are presented in Fig. 11, Fig. 12, and Fig. 13. Based on these figures, if the effectiveness of vaccination is only 2%, then COVID-19 will become endemic, while if the effectiveness is more than 4.5%, the COVID-19 outbreak tends to vanish in the population. The simulation results of obedience in implementing health protocols are presented in Fig. 14, Fig. 15, and Fig. 16. Based on these figures, if obedience in implementing health protocols is only 2%, then COVID-19 will become endemic in the population. However, if obedience to follow health protocols is more than 4.5%, then the COVID-19 outbreak will vanish in the population. From the simulation results of vaccine effectiveness and the level of adherence to follow health protocols, it can be seen that, for the same π and τ values, the results will be the same. This shows that vaccination and the implementation of health protocols have the same effect on the spread of COVID-19 in Indonesia.

Conclusion
A study using an SIR model for the spread of COVID-19 in Indonesia taking into account the factors of vaccination, treatment, implementation of health protocols and the corona virus-load has been performed. In this study, the parameters β, γ , and μ C are treated as membership functions of fuzzy numbers and are represented as fuzzy parameters. Those parameters depend on the corona virus load . The points of disease-free equilibrium and endemic equilibrium are locally asymptotically stable for 0 ( ) < 1 and 0 ( ) > 1, respectively. Based on the simulation results, it is found that vaccination and the implementation of health protocols have a significant effect in slowing or stopping the spread of COVID-19 in Indonesia. Likewise, treatment has an effect in slowing or stopping the rate of infection of COVID-19 but not as much as the effect of vaccination and the implementation of health protocols.