A model of bi-mode transmission dynamics of hepatitis C with optimal control.

In this paper, we present a rigorous mathematical analysis of a deterministic model for the transmission dynamics of hepatitis C. The model is suitable for populations where two frequent modes of transmission of hepatitis C virus, namely unsafe blood transfusions and intravenous drug use, are dominant. The susceptible population is divided into two distinct compartments, the intravenous drug users and individuals undergoing unsafe blood transfusions. Individuals belonging to each compartment may develop acute and then possibly chronic infections. Chronically infected individuals may be quarantined. The analysis indicates that the eradication and persistence of the disease is completely determined by the magnitude of basic reproduction number R(c). It is shown that for the basic reproduction number R(c) < 1, the disease-free equilibrium is locally and globally asymptotically stable. For R(c) > 1, an endemic equilibrium exists and the disease is uniformly persistent. In addition, we present the uncertainty and sensitivity analyses to investigate the influence of different important model parameters on the disease prevalence. When the infected population persists, we have designed a time-dependent optimal quarantine strategy to minimize it. The Pontryagin's Maximum Principle is used to characterize the optimal control in terms of an optimality system which is solved numerically. Numerical results for the optimal control are compared against the constant controls and their efficiency is discussed.

Introduction

Hepatitis C virus (HCV) is the common cause of liver diseases worldwide and a major public health problem (Bisceglie 1998). The disease hepatitis C was first recognized in 1975 and the causative agent HCV was identified in 1989. Hepatitis C is characterized by an acute, often asymptomatic stage, followed in most cases by chronic infection that can lead to cirrhosis and liver cancer. HCV, which causes hepatitis C in humans, is a small, enveloped, positive-sense single-stranded RNA virus of the family Flaviviridae. Replication of the RNA-based virus involves the use of the enzyme RNA-dependent RNA polymerase (RdRP) that has a high error rate. Consequently, the virus mutates very rapidly and has no single genotype. Currently, there are seven known genotypes of HCV (Jawaid and Khuwaja 2008; Torresi et al. 2011). According to World Health Organization, nearly 3 % of the world population has been infected with HCV. An estimated 170 million people globally are infected with chronic HCV and are at risk of developing liver cancer/cirrhosis, while nearly 350,000 people die worldwide as a result of hepatitis C-related liver diseases each year (Waheed et al. 2009; Jiwani and Gul 2011).

Hepatitis C is characterized by an acute and a chronic stage. Initial infection by HCV results in acute hepatitis C, which is largely asymptomatic. Only about 15 % of cases display mild symptoms such as decreased appetite, fatigue, nausea, muscle or joint pains, and weight loss. The infection resolves spontaneously in 20% of the cases. About 80% of the people exposed to HCV, however, eventually develop a chronic infection, thus progressing to the chronic stage of the disease. This stage can last for decades. Most people experience minimal or no symptoms during the initial few years of the infection. After several years of living with the disease hepatitis C becomes the primary cause of cirrhosis and liver cancer. Nearly 5–20 % of chronic hepatitis C patients develop cirrhosis over 30 years and 1–5 % die from cirrhosis or liver cancer. Those who develop cirrhosis have a 20-fold greater risk of hepatocellular carcinoma, a rate of 13 % per year. It is estimated that hepatitis C is the cause of 27 % of cirrhosis cases and 25 % of hepatocellular carcinoma worldwide (Bisceglie et al. 1991; Fattovich 1997; Kiyosawa 1990; Hutin 2004).

The standard treatment of infected individuals comprises a combination of pegylated interferon alpha (Peg IFN-α) and the potent antiviral drug Ribavirin for a period of 24 or 48 weeks, depending on the HCV genotype (Jawaid and Khuwaja 2008). The response to treatment also varies by genotype and ranges from 70–80 % for genotypes 2 and 3 to almost non-existent for genotype 6. Also, recent treatment advances of genotype 1 infection using directly acting antiviral agents are encouraging; there is still a need to develop an effective vaccination strategy capable of preventing infection. The vaccines used in clinical and pre-clinical trials include recombinant proteins, DNA-based proteins, synthetic peptide vaccines, etc. The implementation and success of the above-mentioned vaccines along with future designs of vaccine strategies have been discussed in detail by Toressi et al. (2011). Furthermore, the high cost of the therapy results in low clinical usage. In USA, the median HCV health care expenditure of Interferon therapy exceeds (Armstrong and Charland 2004). Recovery from infection with HCV does not result in long-term immunity. Therefore, any model for hepatitis C must reflect this lack of acquired immunity by allowing recovered individuals to become susceptible again.

HCV is primarily spread by blood-to-blood contact through blood transfusions, the use of poorly sterilized medical equipment and intravenous drug use (IDU). Screening at blood transfusion has resulted in a significant decline in HCV transmission, especially in the developed world. However there are countries lacking resources where screening costs have posed a hinderance to blood screening for HCV (Alter 2007). The primary route of transmission in the developed world is IDU, while in the developing world HCV is predominantly spread through unsafe blood transfusions and therapeutic procedures (Maheshwari and Thuluvath 2009). The HCV prevalence among IDUs is greater than HIV [67, 72.5 and 73.4 % among the IDUs in China, Russia and the USA, respectively (Nelson et al. 2011)]. The World Health Organization declares IDUs as an important target group for HCV treatment and prevention (Nelson et al. 2011). Among IDUs in USA and Europe, transmissions occur via sharing injection equipments such as syringes, needles and other paraphernalia (Mathei et al. 2006). Both modes of HCV transmission (unsafe blood transfusion and IDUs) are significant for certain communities [studies indicate prevalence percentages of 48.67 ± 1.75 % and 57 ± 17.7 %, respectively, in Pakistan (Waheed et al. 2009) see also Ali et al. (2009)].

Developing countries of South Asia (especially Pakistan) have high HCV prevalence rates. Pakistan is a highly populated, developing country with low health and education standards and therefore the population is particularly vulnerable to infection with the hepatitis C virus. In Pakistan 10 million people are presumed to be infected with HCV. Percentage prevalence of HCV in the adult population of Pakistan is estimated to be about 4.95 ± 0.53 % (Waheed et al. 2009). It was estimated that there were about 5 million drug users in Pakistan, out of which 15 % were regular IDUs. In contrast to the US, where the majority of HCV infections are caused by genotype 1 of HCV, an overwhelming majority of hepatitis C patients in Pakistan are infected with genotype 3a of the virus. In 2004, a panel of gastroenterologists reported that 75–90 % of HCV patients in Pakistan had genotype 3a (Waheed et al. 2009; Jiwani and Gul 2011). The therapy of the genotype 3a HCV is more expensive than other genotypes of the virus which poses a great threat of affordability in a developing country such as Pakistan.

Other HCV transmission modes, though less efficient (and not included in the current model), include occupational exposures where HCV prevalence among healthcare workers is roughly 1–2 % (Alter 2007) [about 5.2 ± 0.63 % in Pakistan (Waheed et al. 2009)] with an average transmission rate of 0.5 % to patients (Alter 2007); perinatal exposure with a transmission rate of 4–7 % per pregnancy (Alter 2007), which increases with HIV co-infection; sexual exposure, where the transmission rate is controversial with inconsistent results reported by various studies (Alter 2007). However the risk increases with unsafe, high risk sex practices, multiple sex partners, etc. as this increases the probability of a sexual relationship with an HCV positive partner (Alter 2007). Other human activities, least rarely associated with HCV transmission, include percutaneous exposure to blood such as during tattooing, body-piercing, intranasal drug use, and acupuncture (Alter 2007).

Several studies on modeling epidemics in populations already exist (Hutin 2004; Dontwi et al. 2010) and they are pertinent to our work. Most of these papers classify individuals in the population into different states and then formulate a system of ordinary differential equations (ODE) to analyze the time evolution of each of these population states. Reade et al. (1998) has discussed an ODE model of infections with acute and chronic stages. Similarly, Luo and Xiang (2012) analyzed a four state system with exposed, acute and chronic states. Suna et al. (2012) has analyzed an SEIRS model where it was assumed that recovered individuals lose their infection-acquired immunity. Martcheva and Castillo-Chavez (2003) have formulated a model for hepatitis C lacking an exposed class and have discussed the stability of the equilibrium states. Dontwi et al. (2010) considered the transmission of hepatitis C through IDUs in an acute, chronic and recovered model. Other deterministic models pertinent to our work have also been proposed and analysed (Busenberg and Hadeler (1990);  Esteva and Vargas (1999); Hadeler and Castillo-Chavez (1995); Hethcote and Thieme (1985); and Hethcote et al. (2002)).

We present a deterministic model with two distinct susceptible population groups comprising IDUs and those getting blood transfusion (unsafe needles in healthcare setting etc.), respectively. So we are considering the transmission of hepatitis C among these two population groups (representing two distinctive modes of transmission). Individuals in each group are classified as susceptible, acute, chronic, quarantined and recovered. The acute and chronic states represent the individuals in the acute and chronic stage of the infections, respectively. The quarantine state represents the chronically infected individuals getting quarantined. The isolation of those with disease symptoms is among the first infection control measures ever recorded (Hethcote 2000). Over the course of time, quarantine strategy has been used to combat the spread of many emerging and re-emerging human diseases such as leprosy, plague, cholera, typhus, yellow fever, smallpox, diphtheria, tuberculosis, measles, Ebola, pandemic influenza and, more recently, SARS (Chowell et al. 2004; Yan et al. 2007; Lloyd-Smith et al. 2003). Luo and Xiang (2012) has performed a global analysis of a four-state ‘susceptible’, ‘exposed’, ‘acute’ and ‘chronic’ model and proven certain theorems regarding the stability of the steady states of the model. However, almost no analysis of the effects of a quarantine class on a disease with a chronic stage has been done and therefore, our paper will be one of the first attempts to study the effect of quarantine (education campaigns, self quarantine etc.) on the spread of a disease with a chronic stage. It should be pointed out that quarantine or isolation here is not being considered in the traditional sense. Quarantine involves only the isolation of individuals who have been exposed to the disease. Since hepatitis C is largely asymptomatic, and chronically infected individuals can live for years without being identified, a quarantine class can be used to consider individuals who have been clearly identified as being infected with HCV. Such individuals can then undergo treatment for HCV and take extra precautions while interacting with the people around them (self quarantine). The notion of quarantine introduced also includes the educational campaigns to spread awareness about HCV, its transmission modes and how to proceed if infection is identified, free supply of needle to intravenous drug users. Such quarantine strategy can potentially play a significant role in hepatitis C control by decreasing the effective contact rate.

Section 2 of the paper presents the mathematical formulation of the model, its basic important properties and the analysis of a reduced model (of one mode, as the dynamics of the two modes are identical). Section 3 includes the rigorous mathematical analysis of the full model. The effect of using quarantine on population is discussed using a threshold quantity along with sensitivity and uncertainty analyses presented with simulations to investigate the dependence of the reproduction number on some crucial parameters. In Sect. 4, we design an optimal control strategy and numerically solve the optimality system to illustrate the effects of an optimal control strategy. Finally Sect. 5 presents the conclusions drawn from the analysis.

Model formulation

The primary route of transmission of HCV in the developed world is IDU, while in the developing world the main source is through blood transfusions. Our model will consider transmission among both of these susceptible groups as they are statistically significant. S 1 represents the susceptible population who are intravenous drug users and similarly S 2 corresponds to the susceptible population who undergo blood transfusion (frequently subjected to unsafe needles in healthcare setting etc.). The susceptible individuals from either group can get infected and move to the respective compartment.

The asymptomatic nature of the hepatitis C and its slow progression make it difficult to characterize the natural history of disease (Bisceglie et al. 1991). The following assumptions are made in the construction of the model.

All infected individuals develop the acute form of hepatitis C first.

Individuals with either the acute or chronic form of hepatitis C are capable of transmitting the disease.

Individuals with the acute form of the disease either progress to the chronic form or recover naturally. Since the acute form of the disease is largely asymptomatic, there is little chance of treatment at this stage.

There is no life long immunity against HCV after recovering, thus the recovered individuals move back to the susceptible class.

The model assumes that the susceptible drug user population S 1, has a constant recruitment rate and natural death rate μ. Susceptible drug user individuals who get infected suffer from the acute form of hepatitis C and move to the compartment A 1 with the force of infection given by λ 1. Individuals in A 1, in addition to the natural death rate μ, die at a disease-induced death rate δ a. They also have a natural recovery rate of κ 1. Individuals with the acute form of the infection progress to the chronic form of the disease at a rate ξ 1, in which case the individual is shifted to compartment C 1. Individuals in C 1, in addition to the natural death rate μ, also die at a disease-induced death rate δ c. Furthermore, these individuals recover at a rate ψ1 and thus move to the recovered compartment R 1. Also, the individuals in compartment C 1 are quarantined and moved to compartment Q 1 at a rate α 1. Individuals in Q 1, in addition to the natural death rate μ, also die at a disease-induced death rate δ . Quarantined individuals can either become susceptible once more at a rate γ 1 f 1 or regress to become acutely infected with HCV at a rate γ 1 (1 − f 1). Recovery from HCV does not result in long-term immunity. Therefore, recovered patients in R become susceptible at a rate ω1. Also the infected individual from one group can infect the individual from the same group as well as the individuals from other group and vice versa. This aspect has been captured in the force of infections λ 1 and λ 2. However the relative ability to infect the other group might be different and this will be discussed as the analysis is carried out later (Table 1).

Table 1

State variables

Variable	Description
N(t)	Total population
S 1, 2(t)	Population of susceptible individuals
A 1, 2(t)	Population with acute hepatitis C
C 1, 2(t)	Population with chronic hepatitis C
Q 1, 2(t)	Population of quarantined individuals
R 1, 2(t)	Population of Recovered individuals

Subscript 1 denotes intravenous drug users; subscript 2 represents individuals who undergo blood transfusions

Mathematically, the model is as follows: where

The description and values of the model are presented in Table 2. The order of magnitudes of the assumed values is discussed in the “Appendix”.

Table 2

Description and values of the model parameters

	Description	Values	References
Π_1	Recruitment rate of drug users (IDUs)	10	Conservative estimate
Π_2	Recruitment rate for those undergoing blood transfusions etc	10	Conservative estimate
μ	Natural death rate	1/(12 × 60)	Conservative estimate
β 1	Effective contact rate	0.3	Zhang and Zhou (2012)
β 2	Effective contact rate	0.2	Zhang and Zhou (2012)
γ 1	Recovery rate of quarantined	1/(3 × 12)	Conservative estimate
γ 2	Recovery rate of quarantined	1/(3 × 12)	Conservative estimate
f 1,2	Fraction of quarantined that becomes susceptible	0.9	Conservative estimate
ξ 1,2	Progression rate from acute to chronic	2/12	Corson et al. (2013)
α 1	Proportion of chronically infected being quarantined	1/7	Conservative estimate
α 2	Proportion of chronically infected being quarantined	1/10	Conservative estimate
κ 1,2	Proportion of acute infection recovering spontaneously	0.26	Corson et al. (2013)
ψ 1,2	Proportion of chronic infection recovering spontaneously	0.05/12	Zhang and Zhou (2012)
ω 1,2	Proportion of recovered who lost immunity (both)	0.75	Corson et al. (2013)
η 1,2	Modification parameter for infectiousness of acute infection	1.25	Zhang and Zhou (2012)
ζ_1,2	Modification parameter for infectiousness of quarantined	0.2	Conservative estimate
θ 12, θ 21	Modification parameter for cross infectiousness	0.01	Conservative estimate
δ a	Disease-induced death rate for individuals with acute infection	0.001	Conservative estimate
δ c	Disease-induced death rate for chronically infected individuals	0.001	Conservative estimate
δ q	Disease-induced death rate for quarantined individuals	0.0005	Conservative estimate

Subscripts 1 and 2 denote drug users and individuals who undergo blood transfusions, respectively

Basic properties

Since the model (1) monitors human populations, all its associated parameters are non-negative. Further, the following non-negativity result about population holds.

Lemma 1

The variables of the model (1) are non-negative for all time t > 0. In other words, solutions of the model (1) with positive initial data will remain positive for all t > 0.

Proof is presented in the “Appendix”.

Lemma 2

The closed set is positively invariant and attracting.

Proof

If we add all the equations of the above model (1) we will have, Since N(t) ≥ 0, a standard comparison theorem can be used to show that Particularly, if Thus, the region D is positively invariant.

Further if then either the solution enters D in finite time, or N(t) approaches asymptotically. Hence, the region D attracts all solutions in R +10 and the solutions remain bounded.

Since the region D is positively invariant and attracting, it is sufficient to consider the dynamics of the flow generated by the model (1) in D where the model is epidemiologically and mathematically well-posed.

Analysis of reduced model

Before continuing with the analysis of the model (1), we will first discuss a reduced version (2), of main model (1), comprising of just one group of population corresponding to one mode of transmission (either IDUs or those getting blood transfusion). The model is as follows (with subscripts 1 or 2 dropped for simplicity)

The basic properties (positivity of the states and positive invariance of a corresponding region in R +5) of the reduced model are identical to those of (1). This model has a DFE where, N = S + A + C + Q + R. The local stability property of ℵ0 will be determined using the next generation operator method described in van den Driessche and Watmough (2002). The non-negative matrix F, of the new infection terms, and the M-matrix, V, of the transition terms associated with the model are given by and,

The eigenvalues of matrix FV −1 are

It follows that the basic reproduction number R 0 = ρ(FV −1), is given by, where,

Now we discuss the existence of an endemic equilibrium of the reduced model. We define endemic equilibria to be those fixed points of the system in which at least one of the infected compartments of the model is non-zero.

Let ℵ1 = (S **, A **, C **, Q **, R **) denote an arbitrary endemic equilibrium of the reduced model so that

Solving the equations of the reduced model at steady state gives

Note that if C ** = 0, then A ** = Q ** = R ** = 0 and we obtain the disease-free equilibrium solution. Thus, we may assume that Therefore, where,

Finally the endemic steady states are given by:

Thus, we have established the following result:

Lemma 3

The reduced model has endemic equilibria, given by ℵ1, whenever R 0 > 1.

In next section we will perform the analysis of our full model (1)

Equilibrium states and sensitivity analysis

Disease-free equilibrium (DFE)

The model (1) has a DFE, obtained by setting the right hand sides of the equations in (1) to zero, given by

Local stability

The local stability property of ℵ0 will be determined using the next generation operator method (van den Driessche and Watmough 2002). The drug users are unlikely to infect the individuals having blood transfusions, mainly because the drug users need to follow the screening procedure (previous health record, smoking/drug habits etc.) used at blood centers. Therefore, we will take θ 21 = 0. However, we will discuss the case where θ 21 > 0. The non-negative matrix F, of the new infection terms, and the M-matrix, V, of the transition terms associated with the model (1) are given in the “Appendix”.

The eigenvalues of (FV −1) are where It is easy to verify (by expanding the terms) that the denominators of R 0 and R 0′ are greater than zero. Clearly F − V is reducible and the equations of the infected compartment decouple near the disease-free equilibrium (DFE). The two non-zero eigenvalues correspond to the basic reproduction numbers for each mode of transfer. Since the reproduction number is the spectral radius of (FV −1) (van den Driessche and Watmough 2002), we have

The basic reproduction number is interpreted as the average number of new infections that one infectious individual can produce if introduced into a population composed of susceptibles. Since R c is the maximum of the two basic reproduction numbers for each mode, it is sufficient to discuss one of them as their expressions are completely identical, except that they represent different modes of transfer of HCV. Assume, without loosing generality, that R c = R 0 (because the expressions are identical to the case R c = R 0′, and there is no biological bias to choose one over the other). Susceptible individuals acquire infection following contact with either an acute (A 1), chronic (C 1) or quarantined (Q 1) individual. The number of infections produced by an acutely infected individual (near the DFE) is given by the product of the infection rate of an acute individual (β 1 η 1) and the average duration in the acute class (). Furthermore, the number of infections produced by a chronically infected individual (near the DFE) is given by the product of the infection rate of a chronic individual (β 1), the average duration in the chronic class () and the probability that an acute individual survives and progresses to the chronic stage (). Similarly, the number of infections produced by a quarantined individual (near the DFE) is given by the product of the infection rate of a quarantined individual (), the average duration in the quarantined class () and the probability that an acute individual survives and progresses to the quarantined stage ( Finally, we observe that a fraction of newly infected individuals will re-enter the acute class (A 1). Thus, the average number of new infections generated by a single infectious individual is given by

The local stability of the DFE holds due to Theorem 2 of van den Driessche and Watmough (2002).

Theorem 1

The disease-free equilibrium DFE, ℵ0, of the model (1) is locally asymptotically stable if R c < 1 and unstable if R c > 1

Theorem 1 implies that, with R c < 1, a small influx of infectious individuals will not lead to large outbreaks of the disease. To ensure that the disease elimination is independent of the initial sizes of sub-populations, it is necessary to show that the DFE is globally asymptotically stable if R c < 1. This is explored below. Since R c is calculated with θ 21 = 0, we continue with this assumption.

Global stability

Theorem 2

The disease-free equilibrium DFE of the model (1), given by (4) is globally asymptotically stable whenever R c < 1.

The epidemiological implication of the above result is that the disease can be eliminated from the population if the basic reproduction number R c can be brought down to (and maintained at) a value less than unity (that is, the condition R c < 1 is sufficient and necessary for disease elimination). Figure 2 depicts numerical results by simulating the model (1) using various initial conditions with R c < 1. It is evident from the simulation that all solutions converge to DFE, ℵ0, in line with Theorem 2. Proof is presented in the “Appendix”.

Fig. 2

Disease-free equilibrium: simulations showing the total chronically infected population eventually dying out for different sets of initial conditions

Flow diagram of the model

Lemma 4

The disease is uniformly persistent in D if and only if R  > 1: there exists a δ > 0 such that where X represents the infected states of (1).

Proof

The necessity for R  > 1 follows from Theorem 1 and the fact the global stability of the (DFE) precludes any kind of persistence for R  < 1. The theorem can be proved using the approach used to prove Proposition 3.3 by Li et al. (1999), by applying a uniform persistence result in Freedman et al. 1994) and noting that the DFE of the model (1) is unstable whenever R  > 1 (Theorem 1). □

Endemic equilibrium

In this section we discuss the endemic equilibrium, the steady state for which at least one of the infected components of the model is non-zero. We assume, as before, θ 21 = 0.

Theorem 3

The model (1) has an endemic equilibrium whenever R  > 1.

A proof is given in the “Appendix”.

After solving the system for the endemic states explicitly, we came across some interesting results. In case of R c = R 0, the infected population of blood transfusion will vanish and infected population of drug users will prevail. This implies there is a boundary endemic equilibrium where only drug users infected population exists. This can be verified easily from the expressions of λ 1** and λ 2** presented in the “Appendix”. On the other hand R c = R 0′ will result in endemic equilibrium in both blood transfusion and drug users infected population. This is an example of co-existence E c = (S 1**, A 1**, C 1**, Q 1**, R 1**, S 2**, A 2**, C 2**, Q 2**, R 2**) where each infected population of both modes prevails. The difference in notation ( for E b against S 1** for E c) of the equilibria only conveys the difference in magnitude. A discussion is presented in the “Appendix”.

These results can be related to the nature of the infection functions (λ 1, λ 2). Since the blood transfusion infected can infect the drug user population, so when R c = R 0′ this results not only in an endemic among blood transfusions population but also among the drug users. On the other hand, the assumption θ 21 = 0 restricts the drug user to infect blood transfusion population. Therefore, in case of R c = R 0 the endemic equilibrium exists only among drug users. However, if we take θ 21 > 0 we will have three endemic equilibria. Two of these are similar to E b and E c, while the third one is of the form But we could not find the exact conditions in which either of them exists, primarily because with θ 21 > 0, explicit expressions for R cannot be found.

Sensitivity analysis

The asymptotic dynamics of the model are completely determined by the threshold quantity R c, which determines the prevalence of the disease. Since we have a deterministic model, the only uncertainty is generated by the input variation (initial conditions and model parameters). Model parameters are the most integral part of the input data. Therefore, in this section we present parameter-related global uncertainty and sensitivity analyses on R c. Measurement errors or imperfect measurement techniques and natural variations are among factors of uncertainty in parameter estimates. To qualitatively decide which parameters influence the model output (R c) most, uncertainty analysis is carried out, and the degree of confidence on the available parameter estimates is quantified. Critical model parameters for the model are identified, and their impact on the model output, in combination with the other model parameters, is quantified using the sensitivity analysis.

We quantify the uncertainty and sensitivity of R c (using the assumption R c = R 0) as a function of 12 model parameters ( and δ q). For the sensitivity analysis, Partial Rank Correlation Coefficient (PRCC) measures the impact of the parameters on the output variable. To reduce the nonlinearity effects, PRCC method uses the rank transformation by rearranging the data in ascending order and then replacing the values with their array ranks. The Rank Correlation Coefficient (RCC) is used to measure the amount of monotonicity between the input and output variables.

The assumed distributions and mean values of the model parameters are mentioned in the “Appendix”. Based on the uncertainty analysis, our estimates of R 0 for hepatitis C is 2.47 with 95 % CI (1.60, 3.71) (Fig. 3).

Fig. 3

Uncertainty analysis: the probability that R 0 > 1 is 99 % with 95 % confidence interval (1.60, 3.71). This suggests that hepatitis C will get endemic under the present conditions. However, the time taken to reach that state could be large. 10,000 values were generated for each parameter according to their assumed distributions and mean values. These values (presented in “Appendix”) were used to calculate R 0 and its central tendency

The most significant (PRCC values above 0.5 or below −0.5 in Fig. 4) sensitivity parameters to R 0 are This implies that even a small error in the estimation of these parameters can greatly affect the value of R 0 and hence, the analysis of our model. Therefore, these parameters need to be estimated with utmost precision and accuracy to capture the transmission dynamics of the hepatitis C. Out of these four parameters, α 1, β 1 and are the ones which can be controlled and we can try to keep theses values within a range so that the R 0 value does not exceed 1. This is explored below using simulations of R 0 plotted against these parameters. The analyses further suggest the quarantine strategy aimed to reduce the infected population yields the desired result as evident from the fact that α and R 0 have a negative correlation.

Fig. 4

Sensitivity analysis: the proportion of chronically infected being quarantined α 1, proportion of acute infections recovering spontaneously κ 1, effective contact rate β 1 and modification parameter for infectiousness of quarantined are the most significant parameters. This means that even a small error in the estimation of these parameters can greatly affect the value of R 0 and hence, the analysis of our model. Partial Rank Correlation Coefficients (PRCC) are calculated with respect to R 0. Parameters with modulus of PRCC values in excess of 0.5 are declared sensitive to R 0

In order to qualitatively measure the effect of quarantine on the transmission dynamics of hepatitis C, a threshold analysis of the parameter α associated with the quarantine of chronically infected individuals is discussed. We compute the partial derivative of R 0 with respect to α 1: where, This quadratic polynomial (in k 3) has a negative and a positive root (k 3+), given by with It is easy to see that

Thus, the quarantine of chronically infected individuals will reduce R 0 and, therefore, reduce disease burden (new infections, mortality etc.) if k 3 exceeds the threshold k 3+ (i.e., if the per capita rate of individuals leaving quarantine k 3 is greater than the threshold value k 3+). This case is presented in Fig. 5. The parameter values in Table 2 ensure that k 3 exceeds k 3+. Negative correlation of quarantine rate α 1 with R 0 also reinforces the claim that quarantining has a positive effect on reducing disease spread.

Fig. 5

Positive effect of quarantine measures on the infected population with k 3+= 0.0112 < k 3 = 0.0361

Lemma 5

The use of quarantine of the chronically infected individuals will have positive (negative) population-level impact if k 3 > k 3+ (< k 3+).

Figure 6 presents the dependence of the basic reproduction number on the parameters α 1 and β 1, where α 1 denotes the quarantine rate of chronic and β 1 denotes the effective contact rate. These parameters were chosen because of high correlation with the R 0 as shown by the sensitivity analysis. From the contour plot, we see that if β 1 is larger, then R 0 is always greater than one, which implies that it is important to control the effective contact rate. Figure 6b shows that the basic reproductive number may be less than one if α 1 and β 1 can be restricted to a range, leading to the potential extinction of the disease. Also as β 1 increases the rate of change R 0 is high.

Fig. 6

Plots of the basic reproduction number R 0 with respect to infected being quarantined rate α 1 and effective contact rate β 1; a a contour plot of the surface R 0, showing higher quarantine rate α 1 will reduce R 0. b Surface plot of R 0, higher quarantine rate α 1 and low contact rate β 1 will keep R 0 < 1

Optimal control

Pontryagin and Boltyanskii (1986) formulated the optimal control theory for the models with underlying dynamics defined by a system of ordinary differential equations. The theory, application areas and the numerical methods have progressed considerably. The Pontryagin’s Maximum Principle allows us to adjust the control in a model to achieve the desired results. The control parameters are mostly functions of time, mainly appearing as the coefficients in the model.

Optimal control theory has been employed to make decisions involving epidemic and biological models. The desired results and performance of the control functions depend on different situations. Fister et al. (1998), Kirschner et al. (1997) in their study of HIV models (1998); (1997) used optimal control to design the treatment strategies. Agusto (2009) used optimal control strategies of a tuberculosis transmission model. Jung et al. (2002) provide a very good example of deciding how to divide the efforts between two treatment strategies (case holding and case finding) of the two strain TB model. Joshi (2002) formulated two control functions as coefficients of the ODE system representing treatment effects in a two drug regime in an HIV immunology model. The goal was to maximize the concentration of T cells while minimizing the toxic effects of the drug. The analytic and numerical results illustrated the level of two drugs to be used over the chosen time interval. The required balancing effect between two competing goals was well predicted by optimal control theory. Behncke (2000) studied SIR models including vaccination, quarantine and health promotion campaign and obtained analytical results for optimal control. The optimal control intervention policies for stochastic epidemic models were treated by Clancy (1999).

Pontryagin’s Maximum Principle appends an adjoint system of differential equations with terminal boundary conditions, to the original model (state system) of differential equations, in the attempt to characterize an optimal control. The optimality system, which characterizes the optimal controls, consists of the differential equations of the original model (state system) along with the adjoint differential equations (adjoint system). The number of equations in the adjoint system is same as that of the state system. The adjoint functions behave very similar to the Lagrange multipliers (appending constraints to the function of several variables to be maximized or minimized). The adjoint variables maximize or minimize the state variables with respect to the desired objective functional. The details of the necessary conditions for the adjoint and optimal controls are presented in Pontryagin and Boltyanskii (1986), Fleming and Rishel (1975). For the application of these results see the work of Fister et al. (1998).

Since an effective vaccine is not available against all the genotypes of HCV, we have to look for alternate strategies to control the spread of HCV. The isolation of those with disease symptoms, constitute what is probably the first infection control measure since the beginning of recorded human history (Hethcote 2000). Over the decades, quarantine has been used to reduce the transmission of numerous emerging and re-emerging human diseases such as leprosy, plague, cholera, typhus, yellow fever, smallpox, diphtheria, tuberculosis, measles, ebola, pandemic influenza and, more recently, SARS. In our model of HCV, the quarantine compartment was introduced to investigate the effect on the infected population size and results were discussed in the last section. Now we attempt to control the quarantine rate of the chronically infected individuals to control the HCV transmission. This section will explore the effects of quarantine control rate, of chronically infected individuals, on the total size of the infected population. In addition to the dynamics of the original model (1), the quarantine rates labeled as α 1 and α 2 will now be considered as time dependent control parameters. Let this model, with time-dependent control parameters, be labeled as (4.1).

Yan et al. (2007) discuss the application of optimal and sub-optimal control for SARS outbreak, a pair of control parameters were introduced representing quarantine and isolation strategies. The use of quarantine and isolation control has also been studied by Yan and Zou (2008, 2009). The cost of the quarantine facility will strongly influence the policy and we have to consider it in our analysis.

Now we design an optimal control strategy to minimize an objective functional that takes into account both the cost and the number of infectious individuals. The control set U is where b are positive constant which are fixed.

The goal is to minimize the cost function defined as

The coefficients X and W are balancing cost factors due to scales and importance of the all parts of the objective function. This specification involves the numbers of individuals with chronic infection as well as the cost for maintaining quarantine control facilities for drug users as well as blood transfusion individuals. The total cost also includes the cost related of organizational, management, and cooperation, etc. Hence, the cost function should be nonlinear. In this paper, a quadratic function is implemented for measuring the control cost with reference to literature in epidemics control (Fister et al. 1998; Kirschner et al. 1997; Jung et al. 2002).

We need to find an optimal control pair (α 1*(t),α 2*(t)) such that The existence of a solution to the optimal control problem can be obtained by verifying sufficient conditions. We refer to the conditions in Theorem III.4.1 and its corresponding Corollary in Fleming and Rishel (1975). The boundedness of solutions to the system (4.1) for the finite time interval is needed to establish these conditions. Pontryagin’s Maximum Principle (Pontryagin and Boltyanskii 1986) provides the necessary conditions to be satisfied by the optimal control pair. This principle reduces (4.1), (5) and (6) into a problem of minimizing pointwise a Hamiltonian, H, with respect to α 1(t) and α 2(t) where k represents the right hand side of the ith equation of the model (4.1). Using Pontryagin’s Maximum Principle (Pontryagin and Boltyanskii 1986) and the optimal control existence result from Fleming and Rishel (1975), we have the following result:

Theorem 4

There exists a unique optimal control pair (α 1*(t),α 1*(t)) which minimizes J over U. Also, there exists an adjoint system of ϕ 's (see Eq. 9) such that the optimal treatment control pair is characterized as

The proof is presented in the “Appendix”.

The following optimality system, consisting of 20 equations, characterizes the optimal vaccination control as defined in (8)

Next, we discuss the numerical solutions of the optimality system and the corresponding optimal control pairs and the parameters. Solving the optimality system (9), using an iterative method, will result in the optimal quarantine strategy. First solve the state equations with a guess for the control pair (α 1(t), α 2(t)) over the simulated time using a forward fourth order Runge–Kutta scheme. The adjoint functions have final time conditions. Because of this transversality conditions on the adjoint functions (9), the adjoint equations are then solved by a backward fourth order Runge–Kutta scheme using the current iteration solution of the state equations. Then, the controls are updated using a convex combination of the previous control and the value from the characterizations (8). This process is repeated and iteration is stopped if the values are converging.

Numerical solutions to the optimal system (9) are carried out using MATLAB and are presented here. The parameter values are used from Table 2 and the initial conditions are given in the “Appendix”. The parameter values used have R c > 1 when the model without time dependent control is considered. Thus, the disease is not expected to die out without intervention strategies.

Figure 7 represents the optimal quarantine strategy to be employed to minimize the cost and the infected population. Considering the practical constraints, an upper bound of 0.5 was chosen for the optimal quarantine control pair (α 1(t), α 2(t)).The optimal control α 1 is at its upper bound value at the start of time before reaching the minimum value at a slower rate. On the other hand, α 2 is at the maximum value initially and then it sharply decreases before getting to the minimum value of 0. In fact, at the beginning of simulated time, the optimal controls are staying at the upper bound to quarantine as many chronically infected individuals as possible to prevent the infected population from increasing. The steadily decreasing of the control pairs is determined by the balance between the cost of the infected individuals and the cost of the control facilities. Also the assumption of θ 21 = 0 appears to be reasonable as a positive value of θ 21 does not change the optimal control strategy.

Fig. 7

Optimal quarantine control: simulation presents the quarantine strategy to be followed to prevent the epidemic and disease spread. Also the positive value of the cross infectiousness parameter θ 21 is not making any significant difference in the outcome as assumed earlier

Figure 8a shows the total infected population for the optimal control and constant control l. It is clear that with the use of an optimal control strategy disease epidemic can be prevented and disease remains under control at all times. Figure 9 shows the cost associated with the optimal and constant control strategy. It is clear the costs of optimal strategy are much less than the cost of relatively low (but practically feasible) constant control and in fact differ by order of magnitude of tens. It is important to note that high constant quarantine rate (α 1,α 2 = 0.4) incurs almost the same cost as of optimal control. However, practically it is highly unlikely to implement these high constant controls primarily due to the lack of required resources and facilities.

Fig. 8

Chronically infected population: simulation presents comparison of the total chronically infected individuals under optimal and constant control. Clearly optimal strategy prevents the epidemic and retains the infected population to a minimum

Fig. 9

Accumulated cost: simulation presents comparison of the cost incurred to implement optimal and different constant control strategies to control hepatitis C. Optimal strategy is considerably cheaper than different feasible constant control strategies

Figure 10 captures the effect of change in effective contact rate over the optimal control strategy. It is clear from the simulation that an increase in the contact rate will lead to higher rates of quarantine. This result is in line with the sensitivity analysis where it was shown that contact rate β 1 have a positive correlation with R c (which determines the disease prevalence).

Fig. 10

Optimal control: simulation presents the change in the optimal quarantine strategy as the effective contact rate is changed for intravenous drug users and for blood transfusion group

Conclusions

In this paper we have discussed an ODE-based deterministic model for the transmission dynamics of the Hep C virus. The main contributions of this work are the inclusion of two distinct susceptible population groups; individuals associated with the unsafe blood transfusions and intravenous drug users (IDUs) along with looking at the effects of quarantine on the dynamics. The two different susceptible groups were considered as these are the major transmission pathways for Hep C, in developing countries such as Pakistan. Further in the absence of an effective vaccine against all known genotypes of hepatitis C formulating optimal quarantine strategies would provide insight to public health officials.

Dynamical system techniques were used to analyze the model. We observed the following points

The disease-free equilibrium of the model is globally asymptotically stable whenever the basic reproduction number R c is less than unity;

The disease persists uniformly if and only if R c > 1, in which case the model has an endemic equilibrium;

The infected population corresponding to the mode (blood transfusions or IDUs) with the higher basic reproduction number surely persists;

To ascertain the relative importance of various parameters, sensitivity analysis was performed, this would help epidemiologists and public health officials to focus on the more important parameters in formulating a disease control policy. Our analysis led to the following observations

The model is most sensitive to the control variables α (proportion of infected population being quarantined), β (effective contact rate) and (relative infectiousness of quarantined individuals);

Given the nature of the disease, controlling quarantine parameter (i.e. devising an effective quarantine strategy) seems to the most workable solution;

Quarantining of the chronically infected individuals has a positive population-level impact if a certain condition is met;

Finally we considered control strategies to prevent the spread of the disease. In our analysis we assumed a quadratic cost function due to obvious non linearity of the cost as briefly discussed in the optimal control section, and the fact that convexity of the function allows one to apply established results from optimal control theory, as has been done in similar work in the literature. Using techniques from optimal control theory we concluded that

An optimal control (rather than a high constant control) is preferable, where quarantining rate of chronically infected individuals is a function of time. That is, the proportion that is quarantined optimally with respect to time has a higher favorable impact (as compared to implementing a high but constant quarantine rate) in keeping the cost of disease control low.

However, it should be pointed out that the ideal time-varying optimal strategy might not be applied easily. Still, it does provide a basis on which to design practical quasi-optimal control strategies

There are several directions in which to extend the present work, we outline some of these here. In our analysis we took the modification parameter for cross infectiousness, θ 21 = 0, this was justified as the numerics showed no qualitative difference between the zero and non-zero θ 21 case. It would be interesting to analytically verify this numerical observation. We considered two modes of transmission of the disease, one can also look at a third important mode of transmission, the spread of virus via facilities such as the dentists or barber shops. Finally an interesting extension of this work would be looking at a Markov Chain based stochastic model, involving the two transmission pathways.

While early diagnosis and treatment of HCV might be able to reduce the progression of disease, the majority of infected individuals are asymptomatic and most infected persons are unaware of their exposure to the virus. Increase in the public awareness of HCV can play a significant role in increasing the recovery from acute and chronic stages of the disease. Therapy along with an effective quarantine strategy can greatly reduce the prevalence of hepatitis C.

Parameter	Distribution	Mean	SD
μ	(N)	1.4E−03	2E−04
δ a	(G)	3E−03	7E−04
δ c	(G)	9.3E−04	9E−04
δ q	(G)	2E−04	4E−04
γ	(N)	2.8E−02	9.9E−05
ξ	(N)	1.67E−01	1.5E−02
α	(U)	1.5E−01	2.9E−02
κ	(N)	2.6E−01	4.9E−02
ψ	(N)	4E−03	1E−04
β	(G)	3E−01	7E−04
η	(G)	9.3E−04	9E−04
ζ	(G)	2E−04	4E−04