1. Introduction

Due to their continued global prevalence, infectious diseases have been receiving great attention. While some developed nations have been affected by the adverse consequences of infectious diseases, their prevalence and impact are more profound in developing nations where prevention and treatment are not readily available. For instance, the spread of Sexually Transmitted Diseases (STDs) such as Human Immunodeficiency Syndrome (HIV) have had devastating effects on the socio-economic structure of many developing nations. In the past, it was believed that, with the introduction of effective antibiotics and vaccination programs as well as improved sanitation, infectious diseases will soon be eradicated. However, the world has witnessed the propagation and continued presence of infectious diseases at a global scale in spite of humanity's relentless effort to get rid of them. As a result of this alarming reality, the study of infectious diseases has become an area of significant scientific research. In particular, there has been immense scientific investigation focused in developing critical comprehension of the conditions or factors that contribute to the epidemic of diseases and the controlling measures that can be employed to curb this epidemic. These investigations attempt to address a number of vital scientific questions; how effective will an introduction of vaccination help decrease the impact of the epidemic? How will the behaviours of individuals affect the spread of the disease when vaccines are introduced?

As a common tool in mathematical modelling, system of equations has always been effective in comprehending disease dynamics among population. See for instance [19,20,28,32,33]. However, most of these models do not explicitly take into account the impacts of behavioural responses of individuals on disease dynamics including whether there will be epidemic or not. Along with the rapid development of Internet and social network applications, social behaviour has become a new challenge in public health. As a result, infectious disease models incorporating social influence, risk perception and decision-making have attracted more and more interests of researchers from multiple areas. For example, in [16] (Science 2013), it was shown that when a social contagion is coupled to a biological contagion, the disease-behaviour system exhibits complex dynamics and social impact can be either positive or negative. In [34], injunctive social norms were added to an existing behaviour-incidence model to study the dynamics of vaccinating behaviour. To study the pandemic potential of influenza such as H1N1 or H5N1, an epidemiological game-theoretic model of an influenza pandemic was developed in [41]. The model compared the perspectives for antiviral coverage at the individual level (individual behaviour) and the population level to determine the optimal [41]. As examples, risk prediction in decision-making from psychological point of view was studied using the fuzzy trace theory [39] and later, risk perception for HIV transmission response through multiple pathways is discussed in [42]. Some recent progress on population dynamics involving decision-making under the consideration of both human-environment and biological conditions can be found in [1,2,17,34,42,43]. In particular, a review on coupled disease-behavior dynamics including social and disease perspectives from the approach of complex networks was given in [43].

At the meantime, as an emerging class of models, Economic Epidemiology (EE) is interdisciplinary and utilizes economic concepts to explicitly incorporate behavioural related responses. For some previous efforts in this direction see [8,10,18,35,37]. Also, a detailed account of some results in this field maybe found in [37] and recently in [5,36]. Different from the traditional mathematical epidemiology, EE models apply incentives for healthy behaviour and associated behavioural responses to offer unique insights to transmissibility of infectious diseases and thereby recommending optimal control strategies to contain their spread.

As classical models for epidemiology, SIR (susceptible-infected-removed), SI, SIS and SIRS models and their various extensions have been extensively studied [11]. For instance, [23] discusses continuous time mathematical epidemiological models integrated with utility functions and decision making and a discrete time EE models were discussed in [5], where equilibrium dynamics of EE under rational expectations were investigated.

In this paper, we introduce a discrete time EE model of vaccination incorporating social consideration and decision-making. Individuals are assumed to have control over the contacts they make. That is, it is assumed that contacts are made in order to maximize utility (satisfaction) subject to disease dynamics, health stock and probability of infection. It is also assumed that newborns and older susceptible individuals are vaccinated. The approach used in this study is similar to the one used in [4] but more detailed investigation on properties of system bifurcation and sensitivity analysis for the infection parameter and contact rates provides in-depth understanding to system dynamics. We also present results obtained from sensitivity analysis with specific disease parameters to show the behaviour of some particular diseases. Our results show that the system exhibits saddle-point and period-doubling bifurcation when older susceptible individuals are vaccinated. The converse also holds. That is the model does not show these bifurcation properties if there is no available vaccination for older susceptible individuals.

The organization of the paper is as follows: Section 2 describes the model, the various parameters and the underlining assumptions governing the formulation of the system. Section 3 describes the optimal behaviour of susceptible individuals. Parameter analysis is carried out in Section 4, where the endemic equilibria (when diseases is present in the population) of the prevalence rate of disease and sensitivity analysis on the equilibria of the number of contacts is carried out. Also, local stability analysis around the disease endemic equilibria as well as bifurcation analysis are discussed in this section. In Section 5, simulation of the system is presented and discussed. Conclusion on the paper is presented in Section 6.

---

2. Model development

The model is set in a discrete time interval such that individuals make decisions in discrete time given their disease status. The model considers a total population, say $N$, divided into three mutually exclusive disease groups: Susceptible ($S$), Infected ($I$) and Vaccinated ($V$) category. The disease category is written as the proportion of the total population. Let $s_{t},i_{t}$ and $v_{t}$ denote the proportion of the disease categories (Susceptible, Infected and Vaccinated) to the total population at a given time, $t$ respectively, hence $s_{t}+i_{t}+v_{t}=1$. We derive the epidemiological model by making the following assumptions:

● The population size is assumed constant with equal birth and death rates per year. This is denoted as $\mu$. The death rate is assumed the same for all disease categories. The new births enter the susceptible group (we assumed no infection among the new births). The constant $m>0$ is the proportion of children who are vaccinated in their first year. Adult susceptible individuals are assumed to receive vaccination. Let $n$ denotes the proportion of these vaccinated individuals. We assume that vaccination does not confer permanent immunity. This assumption is consistent with epidemiological data since for many infectious diseases immunity wanes following natural infection. While the duration of immunity provided by vaccines varies, live vaccines generally induce longer lived immunity than sub unit vaccines [50]. For example for Tetanus, the estimated duration of protection from vaccine after receiving all the recommended doses is $25$ years for $72 %$ of the vaccinated individuals while for measles permanent immunity is obtained for $96 %$ of the vaccinated individuals [38]. Therefore, a fraction of the vaccinated individuals are assumed to lose their immunity. Let $\nu$ denote this fraction.

● Following the work in [26] and [44], the treatment function of infected individuals is modelled as

$T(i_{t})= \begin{cases} \sigma i_{t} & \text{if} 0\leq i_{t} \leq i_{c},\\ \nonumber \sigma i_{c} & \text{if} i_{t} > i_{c}. \end{cases}$

Where $\sigma$ and $i_{c}$ are the respective treatment parameter and the disease prevalence level beyond which the health care system can not operate at a particular time. $\sigma$ is the percentage of individuals who receive treatment per year. The treatment function indicates that the treatment increases linearly with infection prevalence and then reaches a peak at which it becomes constant.

● There is no disease related death.

From the above, we have the following system of equations explaining the epidemiological model:

${s_{t + 1}} - {s_t} = \mu - m\mu - ({p_t} + \mu + n){s_t} + T({i_t}) + \nu {v_t},$

(1)

${i_{t + 1}} - {i_t} = {p_t}{s_t} - \mu {i_t} - T({i_t}),$

(2)

${v_{t + 1}} - {v_t} = m\mu + n{s_t} - (\mu + \nu ){v_t},$

(3)

where $p_{t}$ is the probability of infection. By intuition, individuals can either be in the susceptible or infected or vaccinated disease category at a time. Also, a contact made by a susceptible individual can either result in infection or not. The probability of a susceptible individual transitioning to the infected group is modelled by following the work in [3]-[7]. We have the probability of infection as

$p_{t}=1-(1-\lambda i_{t})^{c_{t}},\label{probinf}\nonumber$

where $c_{t}$ and $\lambda$ are the number of contacts made by an individual at time $t$ and the probability of a susceptible individual contracting the disease from a single infected contact per unit time respectively. $\lambda$ can also be viewed as the rate (force) of infection per unit time. Since $s_{t}+i_{t}+v_{t}=1$, Eq. (1) -(3) can be reduced to the following system of equations:

$\begin{cases} s_{t+1}-s_{t} =\mu-m\mu-(p_{t}+\mu+n)s_{t}+T(i_{t})+\nu(1-s_{t}-i_{t}),\\ i_{t+1}-i_{t} = p_{t}s_{t}-\mu i_{t}-T(i_{t}).\label{sni2} \end{cases}$

(4)

To introduce behavioural influence of individuals into the dynamics of the model, suppose an individual k¹ independently makes a decision by choosing a number of contacts $c$. Associated with each contact comes a risk of infection. Individuals are faced with the decision to choose $c$ in any period $t$ such that his/her utility is maximized. We have the following as the individual's objective function:

$\sum\limits_{j = 0}^\infty {{\beta ^j}} U({c_{k,t + j}},{h_{k,t + j}}),$

¹k = 1; 2; 3:: can be interpreted as individual 1, 2 and 3 etc..

where $0<\beta<1$ and $h_{k,t}$ are the discount factor and individual $k$'s health stock in period $t$. $\beta$ measures the estimated worth individuals place on their future utility. If they place less value on their future utility, the discount factor takes up high value. On the other hand, the discount factor will take low value when more value is placed on future utility. In the analysis, we specify the utility function $U$ as follows:

$U(c_{k,t},h_{i,t})=c_{k,t}-\delta c_{k,t}^{2}+\phi h_{k,t}\label{utilfun},$

(5)

where $\phi > 0$ and $\delta$ measures the relative importance of the health stock [3] and the level of contacts that yields the possible maximum utility respectively and $h$ denotes the health stock. Thus, we have $c_{k,t}=1/2\delta$ as the level of contact that yields maximum utility. Fig. 1 illustrates this. Notice that the graph for $\delta=0.05$ has a maximum utility less than the graph for $\delta=0.025$. The graph is plotted by dropping the indexes attached to $c \text{ and } h$. Furthermore, low values of $\delta$ yield high level of utility and vice-versa (that is there is a case of disutility when $\delta$ takes high values). Note that $\phi=0$ implies that individuals do not place importance on $h$ and thus an increase in the value of $\phi$ implies that much importance is placed on $h$ [3]. We define the dynamics of the health stock by following the work by Grossman [25]. We have the equation for the health stock as

$h_{k,t+1}=h^{'}+(1-\epsilon)h_{k,t}-i_{t+1},\label{healthstock}$

(6)

where $\epsilon$ captures the depreciation rate of health $h$. $h'$ is the autonomous health stock. It implies from Eq. (6) that the steady state for the health stock of uninfected individuals is $h^{*}=h'/\epsilon$. Hence, $h^*$ serves as the ceiling for $h$.

---

3. Optimal behaviour of individuals

Suppose that individuals are not altruistic (that is, they are not concerned about the welfare of the general public), then each individual will be willing to go for a contact level that will yield maximum satisfaction. Since individuals in the infected and vaccinated category are already infected and vaccinated respectively, we further assumed that unless they are altruistic, they will opt for the maximum number of contacts (that is $c_{k,t}=1/2\delta=\bar{c}$) in order to gain the possible optimal level of utility. Beyond this level, the utility of individuals will decrease (that is individuals will be faced with disutility for any additional contact they make. This phenomenon is due to the quadratic nature of the utility function. Therefore, to maintain a strictly increasing utility function the contact level should not exceed $\bar{c}$. This is shown in Fig. 2. Note that as $c$ moves beyond the vertical line ($c=1/2 \delta$), the utility $U(c)$ begins to decrease. The figure also exhibits how the variation in the health stock affects the level of utility for individuals. In our illustration, we assumed an exogenous increase in the health stock from 2 to 4. This causes an upward shift of the utility curve resulting in the maximum level of utility for individuals moving from 7 to 9². To carry out our analysis, we assume that individuals are identical in behaviour with the exception that they have different health stock. Hence, we can drop the subscripts and analyze in terms of a single individual. We further assume that the infected and the vaccinated individuals are not altruistic and that they go for a level of contact that yields maximum utility. The decision faced by the susceptible individual is to maximize the objective function subject to Eq. (6) and the dynamics of the disease. The following Euler equation was used to explain the optimal behaviour of susceptible individuals [3]:

$\frac{\partial U}{\partial c_{t}}=-\beta\frac{\partial U}{\partial h_{t+1}}\frac{\partial h_{t+1}}{\partial i_{t+1}}\frac{\partial p_t}{\partial c_{t}}.\label{Euler}$

(7)

²The maximum utility was calculated by substituting the value for $c=1/2\delta$ into the utility function.

Eq. (7) measures the trade-off in the model in that an additional contact made by a susceptible individual may or may not result in an infection. Therefore, in order for the individual to be in equilibrium, the individual should make contacts that will satisfy the above equation.

The right hand side of Eq. (7) measures the marginal benefit associated with an additional contact made by a susceptible individual while the expression on the left hand side measures the damage the individual incurs as a result of an additional contact made. See [3] for detailed explanation on the above. By employing Eq. (4), (5) and (6) we have the Euler equation reducing to

$1-2\delta c_{t}=\beta \phi p_{c,t},\label{euler2}$

(8)

where

$p_{c,t}=\frac{\partial p_t}{\partial c_{t}}=-(1-\lambda i_{t})^{c_{t}}\ln(1-\lambda i_{t})=-\frac{(1-p_{t})}{c_{t}}\ln(1-p_t).$

The value of $p_{c,t}$ shows that unlike the classical mathematical epidemiology, the probability of infection is dependent on how susceptible individuals responds to disease risk. Therefore, we can interpret $p_{c,t}$ as the magnitude in which the rate of probability of infection changes with respect to optimal number of contacts.

---

4. Parameter analysis

---

4.1. Equilibria analysis

This section discusses the existence of the disease steady state. Let $E(s^*,i^*,c^*,h^*)$ be the equilibrium(equilibria) for the disease categories. The endemic steady state is found by solving the time invariant version of System (4) and Eq. (6) and (8) simultaneously. The following equations represent the endemic steady state (equilibrium at which the disease is presence in the population) for $0<i^*\leq i_{c}$:

$\begin{cases} s^*=\frac{\mu(1-m)+(\sigma-\nu)i^* +\nu}{p+\mu+n+\nu},\\ i^*=\frac{ps^*}{\mu+\sigma},\\ 1-2 \delta c^* =\beta \phi p_{c},\\ h^* =h'+(1-\epsilon)h^*-i^*. \label{stdyeqnforsi1} \end{cases}$

(9)

And for $i^*>i_{c}$, $s^* \text{ and } i^*$ satisfies

$\begin{cases} s^*=\frac{\mu(1-m)+\nu+\sigma i_{c}-\nu i^*}{p+\mu+n+\nu},\\ \nonumber i^*=\frac{ps^*}{\mu+\sigma i_{c}}\label{stdyeqnforsi2}, \end{cases}$

where

$\begin{cases} p=1-(1-\lambda i^*)^{c^*},\\ p_{c}=-(1-\lambda i^*)^{c^*}\ln(1-\lambda i^*).\label{probsteady} \end{cases}$

(10)

It can be verified that there is no explicit solution for the endemic steady state. Numerical method is employed in examining the existence of positive endemic steady state(s). For example, Fig. 3 shows the relationship between the prevalence rate $i$ and the number of contacts $c$. Note that, for all values of $c$, there is a positive relationship between $i$ and $c$. The graph shown is for the condition $0<i^{*}<i_{c}$. The same result holds for $i^*>i_{c}$. It further shows that as individuals increase their number of contacts, it reaches a level beyond which the disease prevalence remains constant.

The disease free equilibrium, $E_{0}\left(\frac{\mu(1-m)+\nu}{\mu+n+\nu},0,\frac{1}{2\delta}, \frac{h'}{\epsilon}\right)$ is obtained from Eq. (9). The expression for $E_0$ confirms the expectation that, in the absence of infectious disease(s) susceptible individuals will go for maximum number of contacts and yet maintain their highest health stock. It further shows that if the proportion of vaccinated newborns and adult susceptible individuals are zero (that is both newborns and adult susceptible individuals are not vaccinated), $s$ will equal 1.

---

4.2. Sensitivity analysis of number of contacts

This section discusses the sensitivity analysis of the number of contacts to the parameter $\lambda$ and $\phi$, denoted by $S(c^*,\lambda)$ and $S(c^*,\phi)$ respectively. By employing the method used in [29], we have

$S(c^*,\lambda)=\frac{\partial c^*}{\partial \lambda}\frac{\lambda}{c^*}=\frac{\zeta_{1}}{\zeta_{2}},\nonumber$

where

$\begin{cases} \zeta_{1}=\beta \phi (1-\lambda i^*)^{c^*} (c^*\ln (1-\lambda i^*)+1)\lambda i^*,\\ \nonumber \zeta_{2}=(1-\lambda i^*)(\Delta \ln(1-\lambda i^*)-2\delta)c^*, \end{cases}$

and $\Delta= \beta \phi (1-\lambda i^*)^{c^*} \ln (1-\lambda i^*)$.

From the above expression, we can obtain results on the relationship between number of contacts and the infection parameter.

Theorem 4.1. (Positive relationship) If $\sqrt{M_{1}}<\lambda i^{*}< M_{2} \text{ or } c^*>\text{max}\left(M_{3},\right. \left.\frac{2}{\sqrt{e M_{1}}}\right)$ then $S(c^*,\lambda)>0$. Where

$\begin{cases} M_{1}=\frac{2\delta e}{\beta \phi},\\ \nonumber M_{2}=1-e^{-\frac{1}{c^*}},\\ M_{3}=-\frac{1}{\ln (1-\lambda i^*)}. \label{Sclambda1>0} \end{cases}$

Proof.

Case 1. Assume that $\sqrt{M_1}<\lambda i^*<M_2$. Clearly $\zeta_1>0$ since $\lambda i^*<M_2$. This implies that

$c^*\ln (1-\lambda i^{*})+1>0. \label{z11>0} \nonumber$

Similarly,

$\begin{array}{l} {\zeta _2} = (1 - \lambda {i^*})(\beta \phi {(1 - \lambda {i^*})^{{c^*}}}{(\ln (1 - \lambda {i^*}))^2} - 2\delta ){c^*}\\ > \left( {\frac{{\beta \phi }}{e}{{(\ln (1 - \lambda {i^*}))}^2} - 2\delta } \right)(1 - \lambda {i^*}){c^*} \end{array}$

(11)

since $\lambda i^*<M_2$ implies that $(1-\lambda i^*)^{c^*}>\left(e^{-\frac{1}{c^*}}\right)^{c*}=\frac{1}{e}$.

However, $(\ln(1-\lambda i^*))^{2}>(\lambda i^*)^{2}$ since $0\leq \lambda i^*\leq 1$. Hence Eq. (11) implies that

$\zeta_{2}>\left(\frac{\beta \phi}{e}(\lambda i)^{2}-2\delta\right)(1-\lambda i^{*})c^{*}>0 \label{z21>0}$

(12)

since $\lambda i^{*}>\sqrt{M_1}$.

Therefore, if $\sqrt{M_{1}}<\lambda i^*<M_{2}$, then $\zeta_{1} >0\text{ and } \zeta_{2}>0$ and thus $S(c^*,\lambda)>0$.

Case 2. Suppose that $c^*>\text{max}\left(M_3,\frac{2}{\sqrt{e M_1}}\right)$.

Since $c^{*}>M_3$, we have $c^{*}\ln (1-\lambda i^*)+1<0$. Thus $\zeta_{1}<0$. On the other hand,

$\begin{array}{l} {\zeta _2} = (1 - \lambda {i^*})(\beta \phi {(1 - \lambda {i^*})^{{c^*}}}{(\ln (1 - \lambda {i^*}))^2} - 2\delta ){c^*}\\ \quad \quad < (1 - \lambda {i^*})\left( {\frac{{4\beta \phi }}{{{e^2}{c^*}^2}} - 2\delta } \right){c^*} \end{array}$

due to $f(\lambda i^*)=(1-\lambda i^*)^{c^*}(\ln(1-\lambda i^*))^{2}$ has a maximum at $1-\lambda i^*= e^{-\frac{2}{c^*}}$, therefore $\text{max}(f(\lambda i^*))=f(1-e^{-\frac{2}{c^*}})=\frac{4}{e^{2} {c^*}^{2}}$. Since $c^*>\frac{2}{\sqrt{e M_{1}}}$,

$\zeta_{2}<(1-\lambda i^*)\left(\frac{\beta \phi}{e}M_{1}-2 \delta \right)=0 \label{z31>0}.$

(13)

Hence, $S(c^*,\lambda)>0$.

Theorem 1 implies that an increase in the probability of infection with each infected contact $\lambda$ will result in an increase in the number of contacts $c$ on the part of the individual. This may arise as a result of fatalistic behaviour on the part of susceptible individuals. Also, a decrease in $\lambda$ will result in a decrease in $c$. This situation can be attributed to public education attempted at encouraging individuals to reduce contact even though $\lambda$ maybe low.

Theorem 4.2. (Negative relationship) If $\lambda i^*>\;max\;\left(\sqrt{M_{1}},M_{2}\right)\; or\; \frac{2}{\sqrt{e M_{1}}} <c^{*}< M_{3} \;then \; S(c^{*},\lambda)<0$. Where

$\begin{cases} M_{1}=\frac{2\delta e}{\beta \phi},\\ \nonumber M_{2}=1-e^{-\frac{1}{c^*}},\\ M_{3}=-\frac{1}{\ln (1-\lambda i^*)}. \end{cases}$

Proof.

Case 3. Assume $\lambda i^*>\text{max}(\sqrt{M_1},M_2)$. This implies that $\lambda i^*>\sqrt{M_1} \text{ and } \lambda i^* >M_2$.

Since $\lambda i^*>M_2$, we have

$c^*\ln(1-\lambda i^*)+1<0, \nonumber$

which implies that $\zeta_{1}<0$. On the other hand,

$\zeta_{2}=(1-\lambda i^*)(\beta\phi(1-\lambda i^*)^{c*}(\ln(1-\lambda i^*))^{2}-2\delta)c^*>0$

(14)

since $\lambda i^*>\sqrt{M_{1}}$ as shown in Eq. (12) above.

Hence, if $\lambda i^*>\text{max}(\sqrt{M_1},M_2)$, then $\zeta_1<0 \text{ and } \zeta_2>0$. Thus, $S(c^*,\lambda)<0$.

Case 4. Suppose $\frac{2}{\sqrt{e M_1}}<c^*<M_3$. Clearly, $\zeta_1>0$ since $c^*<M_3$ gives us $c^*\ln(1-\lambda i^*)+1>0$, which implies that $\zeta_1>0$.

Similarly,

$\zeta_{2}=(1-\lambda i^*)(\beta\phi(1-\lambda i^*)^{c*}(\ln(1-\lambda i^*))^{2}-2\delta)c^*<0$

(15)

as shown in Eq. (13) above.

Therefore, if $\frac{2}{\sqrt{e M_1}}<c^*<M_3$, then $\zeta_1>0$ and $\zeta_2<0$ and thus $S(c^*,\lambda)<0$.

Theorem 4.2 implies that a reduction in the value of $\lambda$ (probability of infection with each infected contact) will result in an increase in the number of contacts and the opposite holds. That is an increase in the value of $\lambda$ will result in a decrease in the number of contacts by a susceptible individual. The resulting behaviour implies that individuals are risk-averse and thus are very much concerned about their health status. Therefore, public health policy directed at reducing $\lambda$ may in the long-run result in high prevalence rate of the disease, in that the reduction in $\lambda$ may cause individuals to increase their number of contacts.

Also, we have

$S(c^*,\phi)=\frac{\partial c^*}{\partial \phi }\frac{\phi}{c^*}=-\frac{\bar{c}\beta \phi p_{c}}{c^*}.\label{scphi}$

(16)

From Eq. (10), we have $p_{c}$ as positive. Also, since $\phi$, $\beta$ and $c^*$ are positive, we have Eq. (16) as negative. This means that the more importance individuals place on their health status, the less contacts they will be willing to make. The converse holds; the less importance they place on their health status, the more contacts they are willing to make. Of course, we can have a counter intuition to the above, where there exist a positive relationship between $c^*$ and $\phi$. For instance, the Ebola outbreak has shown that even though individuals may place much importance on their health status, people may still will to expose themselves to the disease. This maybe due to people putting less priority on their health status as compared to that of the health of the infected.

---

4.3. Analysis around the endemic equilibrium

Stability analysis of the system is carried out under two cases: When older susceptible individuals are not vaccinated (the proportion of vaccinated susceptible adults, $n=0$) and when they are vaccinated (the proportion of the vaccinated susceptible adults, $n\neq0$). Given the relation $s+i+v=1$, we have established in Section 2 that the model can be reduced to two discrete equations involving susceptible ($S$) and infected ($I$) categories. Therefore, we employed System(4) in our analysis. We investigated the stability properties of the model by first linearizing System (4), Eq. (6) and (8) around the endemic steady states (that is equilibrium at which disease is present in the population) by employing first order Taylor series approximation. Below is the linearized system and equations:

$\hat{s}_{t+1} = (1-p-\mu-\nu-n)\hat{s}_{t}+T(\hat{i_{t}})-(sp_{i}+\nu)\hat{ i_{t}}-sp_{c} \hat{c}_{t},\\ \hat{i}_{t+1} = p\hat{s}_{t}+(1-\mu+sp_{i})\hat{i}_{t}-T(\hat{i}_{t})+sp_{c} \hat{c}_{t},\\ \hat{h}_{t+1} = (1-\epsilon)\hat{h}_{t}-\hat{i}_{t+1},\\ \hat{c}_{t} = (\kappa \frac{c^{*}}{i^{*}})\hat{i}_{t},$

where

$T(\hat{i}_{t})= \begin{cases} \sigma \hat{i}_{t} & \text{if} 0\leq i_{t} \leq i_{c},\\ \nonumber 0 & \text{if} i_{t} > i_{c}, \end{cases}$

and

$\kappa = \frac{\partial c^*}{\partial i^*}\frac{i^*}{c^*}=\frac{\beta \phi(1-\lambda i^*)^{c^*}(c^*\ln (1-\lambda i^*)+1)\lambda i^*}{(1-\lambda i^*)(\Delta(\ln(1-\lambda i^*))-2\delta)c^* }. \nonumber$

$\kappa$ in the linearized system measures the elasticity or the sensitivity of the number of contacts to the infection or prevalence rate of disease. This parameter has its expression equal to that of $S(c^*,\lambda)$. Thus, the conditions for which $\kappa$ is positive or negative are the same as $S(c^*,\lambda)$. $\kappa>0$ implies individuals are not mindful of their health stock or disease status, in that they increase their number of contacts even in the presence of increasing prevalence. This behaviour according to Aadland et al. [4] can be due to fatalistic behaviour on the part of individuals. And for the individual to make an optimal choice, he or she has to increase his or her number of contacts. If $\kappa=0$, the model reduces to the standard Mathematical epidemiological model (see for instance [9,31,45]). By substituting out the control variable, the number of contacts $c$, and making use of the condition $0\leq i_{t}\leq i_{c}$, we have the linearized system in matrix form as follows:

$\begin{bmatrix} \hat{s}_{t+1}\\ \hat{i}_{t+1} \end{bmatrix} = \underbrace{\begin{bmatrix} 1-p-\mu-\nu-n & \sigma-\nu-\theta\\ p & 1-\mu-\sigma+\theta \end{bmatrix}}_\text{A} \begin{bmatrix} \hat{s}_{t}\\ \hat{i}_{t} \end{bmatrix},\label{Matrixwithvacc}$

(17)

where $\theta=s\left(p_{i}+p_{c}\kappa c^*/i^*\right)$ and $\frac{\partial p}{\partial i^*}= p_{i}=c^* \lambda(1- \lambda i^*)^{c^*-1}$. $\theta$ captures the summation of the effect of a change in disease prevalence on the probability of infection where susceptible individuals do not have control over contacts levels (that is the number of contacts $c$ is assumed fixed) and the effect of a change in disease prevalence on the probability of infection as a result of a change in the optimal number of contacts made by susceptible individuals [4]. From System (17), we have the eigenvalues for matrix $A$ as follows:

$\lambda_{1,2}=\frac{X_{1}\pm \sqrt{\psi}}{2},\label{eigennot0}$

(18)

where

$\begin{cases} X_{1}=2(1-\mu)+\theta-(p+\nu+n+\sigma),\\ \psi=(p+\sigma-\nu-\theta-n)^{2}+4pn.\label{psi} \end{cases}$

(19)

From expression (18) the system is locally stable provided that $|\lambda_{1,2}|$ are less than one. The system will exhibit a stable cycle if these eigenvalues have an imaginary part. That is, if $\psi<0$ and the norm of the eigenvalues are less than one. This implies that the system exhibits a dampened cycle (that is the oscillation in the system decay or dies out after a disturbance) for some parameter values, which can be amplified when $\kappa>0$ as shown for instance in [4]. If the proportion of vaccinated susceptible adults ($n=0$) is zero (older susceptible individuals are not receiving vaccination), we have the eigenvalues as

$\begin{cases} \lambda_{1}=1-\sigma-p-\mu+\theta,\\ \lambda_{2}=1-\mu-\nu. \end{cases}$

(20)

Proposition 1. The system does not exhibit stable or dampened cycle for both $n=0 \text{ or } n\neq0$.

The proof is straight forward by observing that the expression for $\psi$ in equation (19) can not be less than zero, making the eigenvalues for both cases to be real numbers. The significance of this is that the introduction of vaccination into the model does not cause the system to exhibit dampened cycle(s).

Proposition 2. If $n=0$, System (17) is locally stable if

$\sigma+p+\mu-2<\theta<\sigma+p+\mu. \nonumber$

Proof. Suppose $|\lambda_{1,2}|<1$, then from the expression for $\lambda_{1}$ we have

$-1<1-\sigma-p-\mu+\theta<1.\nonumber$

This implies that

$\sigma+p+\mu-2<\theta<\sigma+p+\mu.\label{eigproof1} \nonumber$

Also for $\lambda_{2}$, we have

$-\nu<\mu<2-\nu.\label{eigproof2}$

(21)

Since $0<\mu<1$ and $0<\nu<1$, Eq. (21) holds always, therefore we have $\sigma+p+\mu-2<\theta<\sigma+p+\mu$ as the stability condition.

Proposition 2 implies that the system is stable if the sum effect of the change in disease prevalence on the probability of infection (number of contacts is held fixed) and the effect of a change in prevalence on the probability of infection due to a change in optimal number of contacts is less than the sum of the treatment rate ($\sigma$) of disease, the probability of infection ($p$) and $\mu$.

It follows from case 2 that for $n=0$ and $i_{t}>i_{c}$, the system is locally stable if

$p+\mu-2<\theta<p+\mu.$

(22)

Furthermore, if the endemic steady state for which disease prevalence is less than treatment capacity is stable it suffices to conclude that the endemic steady-state equilibrium for the system will be stable if the disease prevalence is greater than the treatment capacity.

Proposition 3. If $n \neq 0$, (that is older susceptible individuals are vaccinated) and $0\leq i_{t} \leq i_{c}$, System (17) is locally stable if

$\begin{cases} \theta_1<\theta<\frac{H}{L},\\ \nonumber L<2, \end{cases}$

where

$\begin{cases} \theta_{1}=\frac{2(F+L)-H-4}{2-L},\\ \nonumber H=LF-np,\\ L=\nu+\mu+n,\\ F=p+\mu+\sigma. \end{cases}$

Proof. Stability requires that

$|Tr(A)|<\det(A)+1<2.\label{stacon}$

(23)

Eq. (23) implies that

$\begin{cases} -1-\det(A)<Tr(A)<\det(A)+1,\\ \det(A)<1, \end{cases}$

(24)

where

$\begin{cases} Tr(A)=2+\theta-F-L,\\ \det(A)=1-F-L+H-\theta L+\theta.\nonumber \end{cases}$

Case 5. $1<L<2$: We have

$Tr(A) < \det(A)+1,\nonumber$

implies

$\theta < \frac{H}{L}. \label{Con1}$

(25)

And

$-1-\det(A)< Tr(A), \nonumber$

implies

$\theta > \theta_1 = \frac{2(F+L)-H-4}{2-L}. \label{Con2}$

(26)

Also,

$\det (A) < 1, \nonumber$

implies

$\theta > \theta_2 = \frac{H-F-L}{L-1}. \label{Con3}$

(27)

For $\theta_{2}<\theta_{1}$, we have

$\frac{H}{L}<F+\frac{(L-2)^{2}}{L}. \label{Con4}$

(28)

Eq. (28) always holds, since

$FL+L^{2}-4L+4-H=pn+(n+\mu+\det(A)-2)^{2}>0. \label{con5}$

Furthermore,

$\theta_{2}-\theta_{1}=\frac{pn+(L-2)^{2}}{(L-2)(L-1)}<0$

Hence, from Eq. (25), (26) and (27) we have

$\theta_{1}<\theta<\frac{H}{L}.$

Case 6. $L<1$: We have

$Tr(A) < \det(A)+1,\nonumber$

implies

$\theta < \frac{H}{L}. \nonumber$

And

$\begin{cases} -1-\det(A)< Tr(A),\\ \nonumber \det (A) < 1, \end{cases}$

implies that $\theta_{1} < \theta < \theta_{2}$. Also, for

$\theta_{2}-\frac{H}{L}=\frac{pn+L^{2}}{L(1-L)}>0.\nonumber$

Hence, $\theta_{2}>\frac{H}{L}$. Therefore, we have

$\theta_{1}<\theta<\frac{H}{L}. \nonumber$

Case 7. $L>2$: From Eq. (26) and (27), $\theta_{2}<\theta<\theta_{1}$. However,

$\theta_{2}-\theta_{1}=\frac{pn+(L-2)^{2}}{(L-2)(L-1)}>0.\nonumber$

Hence, $\theta_{2}<\theta<\theta_{1}$ can not be satisfied. Therefore, for $L>2$ the system will always be unstable. Thus, we conclude that Proposition 3 holds.

Proposition 3 implies that in case of full vaccination, we expect $\mu+\nu<1$. That is birth/death rate and those of the vaccinated who lose immunity can not exceed 100 %. Furthermore, from Proposition 3, we can establish that $pn<L(2-L)$: implying that in case of full vaccination, the condition reduces to $p<1-(\mu+\nu)^{2}$. Which implies that the probability of infection should have a bound that reduces (increases) with birth/death rate and those with lost immunity. The gray area in Fig. 4 shows the stability region of System (17). Notice that the stability region indicates that $pn<1$.

It follows from Proposition 3 that System (17) is locally stable for $i_{t}>i_{c}$ and $n \neq 0$ if

$\begin{cases} \theta_{1}<\theta<\frac{H}{L},\\ L<2, \end{cases}$

(29)

where

$\begin{cases} \theta_{1}=\frac{(\mu+\nu-2)(p+\mu)+n\mu-2(L-2)}{L-2},\\ H=(\nu+\mu)(p+\mu)+n\mu,\\ L=\nu+\mu+n. \end{cases}$

(30)

---

4.4. Bifurcation analysis

This section discusses conditions under which the model will exhibits saddle-point and period-doubling bifurcation. We have a saddle-point bifurcation if $\lambda_{1} \text{ and } \lambda_{2}$ equal positive one and a period-doubling bifurcation if $\lambda_{1} \text{ and } \lambda_{2}$ equal negative one [22]. Thus, for $0\leq i_{t} \leq i_{c}$, the system exhibits saddle-point bifurcation if

$\theta=\frac{(\nu+\mu)(p+\sigma+\mu)+n(\sigma+\mu)}{\nu+\mu+n},\nonumber$

and period-doubling bifurcation if

$\theta=\frac{(\mu+\nu-2)(\sigma+p+\mu)+n(\mu+\sigma)-2(\nu+\mu+n-2)}{\nu+\mu+n-2}.\nonumber$

Proposition 4. The model does not exhibit both saddle-point and period-doubling bifurcation if $n = 0$.

Proof. It follows from the corresponding eigenvalues for case $n=0$ that the condition for saddle-point bifurcation is

$\begin{cases} \theta-(\sigma+p+\mu)=0,\\ \mu+\nu=0, \label{saddlebif} \end{cases}$

(31)

and that of period-doubling bifurcation is

$\begin{cases} \theta-(\sigma+p+\mu)=-2,\\ \mu+\nu=2. \label{periodbif} \end{cases}$

(32)

Since $0<\mu<1 \text{ and } 0<\nu<1$, it implies that conditions (31) and (32) can not be satisfied.

---

5. Simulation and discussion

This section discusses some simulation results on the theorems and propositions of Section 4.

---

5.1. Sensitivity analysis of number of contacts

Fig. 6 confirms Theorems 1 and 2. Fig. 6(a) -6b satisfies theorem 1. Fig. 6(a) confirms the condition for which both $\zeta_{1} \text{ and } \zeta_{2}$ are positive; which corresponds to the condition $\sqrt{M_{1}}<\lambda i^{*}< M_{2}$. We chose $\phi= 1, \beta = 0.5, \delta = 0.05, \text{ and } c^{*}= 5$. Fig. 6(b) confirms the case when both $\zeta_{1} \text{ and } \zeta_{2}$ are negative; which gives the condition $c^*>\text{max}\left(M_{3},\frac{2}{\sqrt{e M_{1}}}\right)$. $\phi, \beta , \delta \text{ and } \lambda i^{*} = 0.4$ are chosen as 1, 0.5, 0.05 and 0.4 respectively.

Fig. 6(c)-6(d) confirms theorem 2. Fig. 6(c) confirms the condition for which $\lambda i^{*}>\text{max}(\sqrt{M_{1}},M_{2})$; this corresponds to the case where $\zeta_{1}<0 \text{ and } \zeta_{2}>0$. We chose $\phi=1,\beta=0.5,c^{*}=2 \text{ and } \delta=0.0001$. Fig. 6(d) satisfies the condition $\frac{2}{\sqrt{e M_{1}}}<c^{*}<M_{3}$, which $\zeta_{1}>0 \text{ and } \zeta_{2}<0$. We chose the values for $\phi, \beta, \lambda i^{*}$ and $\delta$ as 1, 0.5, 0.3 and 0.1 respectively.

---

5.2. Sensitivity analysis with specific disease parameters

In this section we present a general numerical sensitivity analysis in which model parameters are motivated by the outbreak of measles virus among 0-12 month old babies. Measles is a highly contagious and a serious respiratory disease caused by a virus. In spite of the availability of a safe and effective vaccine, the disease has remained one of the leading causes of death among young children globally [46,48,49]. Measles is prevalent in developing countries where per capita incomes are low and where their health care system is weak [49]. In this analysis, we consider parameter values related to children who received measles vaccination by the time they celebrate their first birthdays in the Kissii county in Kenya. Table 1 shows the relevant values for the model parameters employed:

Note that we assume $n=0$ since we are strictly considering infectivities among babies. Even though there is no specific antiviral treatment for the disease, it is shown that the number of deaths associated with the disease has been reduced through supportive care and Vitamin A supplements [49]. Thus, we assume $\sigma=40 %$. The relative importance of health stock, $\phi$ is set at 1. Furthermore, the immunity confer by a single dose of measles vaccine given to 12 or 15 months old babies is estimated to be between 85 % and 95 % with a second dose conferring 100 % immunity [47]. Therefore the percentage of babies who lose immunity, $\nu$ is set at 10 %.

The resulting numerical solution for the above parameters yields eigenvalues $\lambda_1 = 0.7356$ and $\lambda_2 =0.8725$. This is indicative of stability of the system around the disease endemic steady state. Table 2 summarizes the results. The results show that at the endemic steady state, the disease prevalence among the babies, $i^*$ is about 26.6 %, proportion of babies susceptible to the disease, $s^*$ is about 53.3 % and the proportion of babies who maintained their immunity against the disease is about 20 %. The probability of infection, $p$ is about 21.4 %.

Figures 5(a)-5(d) show the simulation output for the various disease categories and the number of contacts by babies. $(i,s,v)= (0.1, 0.9,0.5)$ is the initial point for the numerical solution. Fig. 5(a)shows that as time increases, the proportion of susceptible babies declines and eventually stabilizes. The decrease may partly be attributed to the high rate of mass vaccination of new born babies in the county. Also, since measles is very contagious, the interplay of the disease dynamics and the number of contacts may also contribute to this phenomenon. Also, Fig. 5(b) indicates the disease prevalence increases sharply, reaches a threshold and then decreases to the endemic steady state value. The initial increase of the disease may be due to the associated high level of contacts as well as the highly infectious nature of the disease. Conversely, the decline of the disease prevalence may be due to early diagnoses of the disease and control mechanisms put in place to curb the outbreak of the disease. This means that a timely intervention to curb an outbreak or resurgence of infectious disease can go a long way to mitigate the damage that could be triggered by the disease. Furthermore, fig. 5(c) and 5(d) are the graphs of the proportion of babies vaccinated and the associated number of contacts respectively. Both graphs exhibit an initial decline, upward movements and then stability. These dynamics may be explained by the effectiveness of measles vaccines. Therefore, it is vital to increase the awareness of parents about the availability of a preventive vaccine that will protect the babies as well as the possible damage the disease can cause to new born babies [21,40].

---

5.3. Analysis around the endemic steady state

For confirmation of some of the results on the stability of the system we chose the following parameters:

The parameters in Table 3 indicate that 80 % of the newborns in the population are vaccinated. The treatment rate is 60 % of the infected population, probability of a susceptible individual becoming infected by a single infected contact is 0.6, the value for $\beta$ implies that the annual discount rate is 4 %, the values for $\delta$ and $\phi$ indicates that the maximum level of contacts is $10$ and the relative importance individuals attach to their health status is 3 respectively.

Proposition 2. As a form of demonstration, we chose the values for death and birth rates as 5 %. The proportion of vaccinated individuals whose vaccination wear out per annum is chosen as 10 %. The proportion of older susceptible individuals who receive vaccination is set at 0 (that is $n=0$). Fig. 7 shows the simulation result confirming this proportion. The initial points given in the order of $(i,s,v,c)$ are given as $(0.01,0.01,0.01,1)$, $(0.1,0.2,0.4,3)$, $(0.08,0.5,0.6,10)$ for the red, green and blue paths respectively.

Proposition 3. Table 4 contains the respective parameters that satisfy the conditions under which the system is stable or unstable. As a form of demonstration, Fig. 8 shows the numerical simulation for System (4) and Eq. (8) for condition $L<1$. The initial points are given as those given for Proposition 2.

Furthermore, for $L<1$, Table 5 indicates that at equilibrium, approximately 19 %, 19.4 % and 61.6 % of the population will remain susceptible, infected and vaccinated respectively. The probability of infection is approximately 66.3 % and individuals are willing to go for a number of contacts that is approximately 9.

Based on the values for condition $L<1$, we calculate the corresponding value of $\kappa$ as approximately 0.0151. This indicates a positive relationship between disease prevalence and contacts made by susceptible individuals. This behaviour can be attributed to fatalistic behaviour on the part of susceptible individuals as pointed out in [4,13]. This can also be attributed to the fact that, the vaccinated group is large as compared to the infected group and thus rational individuals have incentive to increase contact levels. For instance, notice that the steady state value for the proportion of the vaccinated is about 62 %, which is much higher than that of the disease prevalence (19.4 %).

To illustrate the condition where there is negative relationship between the contact levels and disease prevalence rate, we set $\nu=0.1$. This gives the steady state values in Table 6.

The values in Table 6 give the value for $\kappa$ as approximately -0.0558. This is indicates that susceptible individuals will opt to reduce their number of contacts $c$ in the presence of an increasing disease prevalence. Notice that the value for $\kappa$ can be said to be small. This can be attributed to the fact that as the fraction of those who lose immunity from the disease decreases, it decreases the probability of infection in the long run and thus causing disease prevalence to decrease. Hence, rational susceptible individuals have the incentive to increase their number of contacts. This maybe why the steady state of the prevalence rate is about 5.9 % and that of the number of contact is approximately 9. Fig. 9 shows the simulation results for this case.

Also, Fig. 10 shows the simulation result for $1<L<2$. The initial points are given in the order of $(i,s,v,c)$ as $(0.01,0.01,0.01,1)$, $(0.1,0.2,0.4,3)$, $(0.05,0.4,0.4,8)$ for the red, green and blue paths respectively.

---

5.4. Bifurcation analysis

The parameter values in Table 7 are chosen to confirm the conditions for which the system exhibits bifurcation properties:

Fig. 11(a) and 11(b) show period-doubling bifurcation diagram for the system. The infection parameter $\lambda$ is the bifurcation parameter. Also, the number of contacts is varying in this case. It shows that as the values of $\lambda$ increases, the system makes a transition from a single equilibrium path to multiple equilibria paths. The system changes qualitative behaviour around 0.56, 0.64 and 0.75. This behaviour may be due to the unpredictability nature of behavioural responses of individuals. This behaviour is clearly shown by 11(a) (the bifurcation diagram for disease prevalence).

Fig. 12(a) and 12(b) show the period-doubling bifurcation diagram for the case where the number of contacts is fixed. In carrying out the numerical simulation, we fixed the number of contacts at 8. It is clear from the figures that the system makes a smooth transition of the system from an equilibrium path to a double equilibria paths is observed.

You will notice from Fig. 11 and 12 that there is a difference in the transitioning process. The implication of the above phenomenon is that if rational individuals are allowed to make choices on contact levels, the system can behave in a chaotic manner. This may be due to the unpredictable nature of the behaviours of these individuals. From policy perspective, policy makers should take into consideration the behaviour of individuals, in order to have a clear picture of the dynamics of the epidemiology of diseases, the infection parameter ($\lambda$) should be kept within the range that gives a clear picture of the dynamic of the system.

Fig. 13(a) and 13(b) show the period-doubling bifurcation diagram for $n$ being the bifurcation parameter and individuals in control of the number of contacts. That is $c$ is varying. This diagram shows no chaotic behaviour. This may be due to the fact that as the proportion of vaccinated adults increases, individuals are faced with low risk of infection and thus showing some level of predictability in their behaviours. The case for which $c$ taking up a fixed value shows similar diagram but with a bifurcation value less than its counterpart ($c$ varying).

The epidemiological implication of the above dynamic is that the system has the tendency to switch from equilibria path to chaotic path. In other words, there are regions of multistability in which the disease can have a stabilizing effect as well as chaotic effect. This dynamic can become more complex as disease parameters are varied across some range of values. In [15] for instance, it is shown that as disease induced death rate with higher transmissibility, the system exhibits more complexities that give rise to period-doubling cascades couple with other dynamics. For other studies on possible secondary infection, period-doubling bifurcations and chaotic behaviour in epidemic models, we refer [12,14,30].

---

6. Conclusion

Decisions made by individuals in the presence of infectious disease(s) are most often done in a selfish manner. This is due to the risk and benefits associated with such decisions. Individuals are faced with either forgoing a number of contacts just to maintain their health stock or risk their health stock by making some number of contacts. Associated with this choices are benefits(utilities). Also, associated with this is the cost of becoming infected. These phenomena affect disease dynamics. That is individual decisions determine whether a population will be faced with an epidemic or not.

It is therefore imperative that decision makers take into consideration the effects of private choices on epidemiological processes in order to better understand the disease dynamics for better policy on "disease-epidemic-control". In this paper, behavioural responses by individuals is incorporated into a "disease vaccination model". This gave an explicit fashion in which we can analyze how individuals respond to diseases with available vaccination. In our analyses, we were able to establish that, with the introduction of vaccination the system did not exhibit dampened cycles. And this also holds when adults are not vaccinated. We have established that the system exhibits period-doubling and saddle-point bifurcation where there is transitioning from unstable to stable endemic steady-state paths when adult individuals are vaccinated. With regards to the case where adult susceptible individuals are not vaccinated, the system exhibits either stability or instability; there is no bifurcation.

In carrying out the numerical simulations, the parameter values are chosen to reflect the case where susceptible individuals increase their number of contacts for an initial decrease in disease prevalence. Also, we chose the values of the parameters to reflect the opposite: susceptible individuals respond positively in terms of choices on contact levels to disease prevalence, a situation Aadland et. al. refers to as dynamic resonance [7]. In our case, by choosing high values for the proportion of the vaccinated individuals who lose vaccination (all other parameters are constant), we established a positive relationship between contacts levels and disease prevalence; A situation that causes disease prevalence to increase in the long run. This behaviour can be attributed to social factors such as those that contributed to the Ebola outbreak; where people for instance, consider health condition of their infected relatives is more important than their own health condition and thus make contacts with these relatives. In [4] for instance where "syphilis cycles" is studied, Aadland et. al. pointed out that this dynamics produces cycles in syphilis infection such that an increase in disease prevalence can cause susceptible individuals to increase their partners and vice-versa, and a decrease in the disease prevalence will opt to reduce partners. This phenomenon amplifies the disease cycles in the population.

In a nutshell, this paper is an attempt to further emphasize the importance of considering private choices in formulating policies on tackling infectious diseases. And that instead of public policy makers imposing policies on individuals in the advent of disease spread, they should rather understand how their behaviours can affect the epidemiological process of the disease.

---

Acknowledgments

The authors thank the referee for valuable comments. The project was partly supported by the Natural Sciences and Engineering Research Council of Canada (NSERC).

---