Parameter Identification in a Tuberculosis Model for Cameroon

A deterministic model of tuberculosis in Cameroon is designed and analyzed with respect to its transmission dynamics. The model includes lack of access to treatment and weak diagnosis capacity as well as both frequency- and density-dependent transmissions. It is shown that the model is mathematically well-posed and epidemiologically reasonable. Solutions are non-negative and bounded whenever the initial values are non-negative. A sensitivity analysis of model parameters is performed and the most sensitive ones are identified by means of a state-of-the-art Gauss-Newton method. In particular, parameters representing the proportion of individuals having access to medical facilities are seen to have a large impact on the dynamics of the disease. The model predicts that a gradual increase of these parameters could significantly reduce the disease burden on the population within the next 15 years.


Introduction
Tuberculosis (TB) is a preventable and curable disease caused by Mycobacterium tuberculosis (Mtb) that most often affects the lungs. To date, TB claims the second largest number of fatalities due to a single infectious agent right after Human Immunodeficiency Virus and Acquired Immune Deficiency Syndrome (HIV/AIDS) [1]. Young adults in their most productive years are especially affected. Worldwide, 8.8 million people were infected in 2010. Sub-Saharan Africa carries the greatest burden with over 260 new cases per 100,000 population in 2011 [1]. Structurally speaking, Mtb's infection can remain latent, become active, or progress from latent TB to active TB either by endogenous re-activation and/or exogenous re-infection. Active TB is often acquired through co-infection of Mtb with other diseases (such as diabetes and HIV/ AIDS) or abuse of alcohol and tobacco.
Mathematical analysis of disease transmission models can significantly help to understand the underlying mechanisms of these processes and to design potential therapies [2,3]. The earliest mathematical models describing TB dynamics were constructed in the 1960's by the statistician H. T. Waaler [4]. His models focus on prediction and control strategies using simulation approaches. However, he considered exponential population dynamics in the absence of TB, which is certainly not realistic. Further simulation models have been suggested by Waaler [5], Revelle and coworkers [6], and Ferebee [7]. Revelle [6] first introduced a TB dynamics model consisting of systems of nonlinear ordinary differential equations (ODEs), which he used to study a minimal cost strategy against TB. Later, Blower and colleagues [8][9][10] discussed the persistence condition of TB inside the population and determined the characteristic reproduction ratio R 0 , the average number of new infectious cases caused by a single infectious case in a fully susceptible population over the course of the entire infectious period. A sensitivity analysis of R 0 has been performed by several authors [8]. However, the sensitivity analysis of parameters with respect to R 0 alone does not fully exhibit the impact of these parameters on the global time dependent behavior of the system, especially in the presence of backward bifurcation. From their model, Blower and colleagues derived that 1 < R 0 < 9. Moreover, the infection rate, the probability of fast progression, the re-activation rate, and the TB related death rate appeared to be the most important parameters. Chavez and colleagues [11] developed a mathematical analysis of a TB model without fast progression. Since then, most publications included sophisticated mathematical theories to study the dynamics of tuberculosis, such as center manifold theory and Lyapunov functions [12][13][14][15][16].
In the reality of developing countries the true challenge of TB control is the high level of undiagnosed infectious cases as well as lost sight cases compared to diagnosed infectious cases. In the following, the term "undiagnosed infectious population" means those people that have not yet been to a hospital for diagnosis or have not been detected, but have a pulmonary TB [17,18]; the term "lost sight population" includes those people that have been diagnosed as having active TB and begun their treatment, but quit before the end. Compared to existing results in [11-14, 16, 19-23] (and references therein), the model presented here considers, in addition to undiagnosed infectious and lost sight population, the aspects of exogenous re-infections, disease relapse as well as primary active TB cases, natural recovery and traditional medicine or self-medication (practiced in Sub-Saharan Africa). Moreover, our model takes undiagnosed population, lost sight population, and exogenous re-infections as important factors of TB transmission in Sub-Saharan Africa into account. As a consequence, we divide the infective class into the following three subgroups: i) diagnosed infectious population, ii) undiagnosed infectious population and iii) lost sight population. According to the National Committee for the Fight against TB in Cameroon (NCFT) [24], about 8% of diagnosed infectious people who begin a therapy, never return to the hospital for the rest of sputum examinations and treatment, and then become part of the lost sight population. The design of the model is accompanied by a qualitative analysis to gain insight into the transmission dynamics of TB.
The main difficulty in model construction is the large number of unknown model parameters. Usually, the goal is to determine parameter values that minimize the difference between experimental measurement values and model predictions in a least-squares sense. Often, however, there are not enough data available to determine all parameter values. In this case, the least-squares problems are rank-deficient, which means that there exists a continuum of possible solutions all of which give nearly the same "nice fits". Standard optimization methods typically ignore this detail, failing to take the structure of the underlying inverse problem into account. These algorithms tend to simply return a "solution", but do not consider the question of the identifiability of the parameters or the uniqueness of the solutions. This can yield misleading results, especially when parameters are co-regulated. In view of this situation, we solve such nonlinear least-squares problems by an error-oriented Gauss-Newton method [25], which monitors the numerical rank of the Jacobian matrix and converges locally, for a well-defined class of statistically "adequate" problems, to a solution that is unique within the subspace of identifiable parameters. As usual, special attention must be paid to initial guesses for the parameter values to be estimated. Since the algorithm converges to a local optimum, these values should reflect good approximations to the "true" values. In fact, this methodology has not been applied to the modelling of TB dynamics before and represents an important motivation for writing this paper. Alternatively, one might think of developing simpler models with less parameters. Such models, however, are prone to lose their predictive ability, if important mechanisms of disease transmission are missing. On the contrary, larger models can help in determining the importance of assumed mechanisms for the observed effects.
The paper is organized as follows. At the beginning, we derive the model equations and briefly review the mathematical algorithm used for sensitivity analysis and parameter identification. Then we show that the model is mathematically well-posed and epidemiologically reasonable. We present the values of model parameters and discuss the results of the sensitivity analysis. Finally, by means of a controlled dynamical system we simulate an increasing access to TB treatment and discuss the effect of this strategy. Even though the used data originate from Cameroon, the class of TB epidemiological models presented here can be extended to many classes of infective individuals and data for further African countries.

Methods
This section presents our newly suggested mathematical model and briefly reviews the mathematical algorithms for sensitivity analysis and parameter identification.

Suggested epidemiological model
Based on the available data, a finite total population at time t denoted by N(t) is considered and sub-divided into the following mutually exclusive sub-populations ("compartments"): S susceptible: healthy people not yet exposed to TB E latently infected: exposed to TB but not infectious I diagnosed infectious: have active TB confirmed after a sputum examination in a hospital J undiagnosed infectious: have not yet been to a hospital for diagnosis but are active for confirmation by a sputum examination L lost sight: people who have been diagnosed as having active TB, begun their treatment and quitted before the end R recovered: people cured after treatment in the hospital With these definitions, we are now ready to derive the details of our compartment model. In Africa, reliable TB tests [26] are often missing or too expensive. Hence, TB diagnosis based on a single sputum examination can often only be classified as "probable" or "presumed", and cannot detect cases of less infectious forms of TB [27].
Transmission of Mtb occurs due to sufficient contacts among susceptible people and infectious TB cases. Diagnosed infectious people are in most cases hospitalized for at least two months or are advised to lessen their infectiousness in their neighborhood. Their distribution in the population is not necessarily homogeneous. People from class L (lost sight) are first diagnosed before quitting the treatment. Even though they are more likely to transmit the disease than diagnosed cases who continue their treatment, since they are aware of their infectiousness they are less likely to transmit TB than undiagnosed cases. We consider a density dependent force of infection for both diagnosed infectious people as well as lost sight infectious people, but with different contact rates for the infection [28]. Since undiagnosed infectious people remain inside the population, there is an unlimited possibility of contacts with the susceptible population [17]. Thus, susceptible individuals acquire Mtb infection from individuals with active TB (classes I and J) and lost sight people (class L) at a rate ν(I, J, L) given by Here, β i , i = 1, 2, 3, are the effective contact rates with diagnosed, lost sight and undiagnosed infectious population sufficient to transmit infection to susceptible people. The effective contact rates β 1, 2 in a given population are measured in effective contacts per unit time. This may be expressed as the product of the total contact rate per unit time (η i ) by the risk of infection (ϕ i ) given contact between an infectious and a susceptible individual, β i = η i ϕ i . All recruitment is into the susceptible class and occurs at an average scale Λ. The fixed survey for non-disease related death is μ, thus 1/μ is the average lifetime. Diagnosed infectious, undiagnosed infectious and lost sight population have additional constant death rates due to the disease, defined by d 1 , d 2 and d 3 , respectively.
A proportion p of the latently-infected individuals develop fast active TB and the remainder (1−p) develop latent TB and enter the latent class E. Among latently-infected individuals developing active TB, a fraction f is assumed to undergo a fast progression directly to the diagnosed infectious class I, while the remainder (1−f) enters the undiagnosed infectious class J. We set Once latently infected with Mtb, an individual will remain so for life, unless reactivation occurs. Latently infected individuals are assumed to acquire some immunity as a result of infection, which reduces the risk of subsequent infection but does not fully prevent it.
Due to endogenous reactivation, a fraction 1−r 1 of latently infected individuals who did not receive effective chemoprophylaxis become infectious with a constant rate k. They reinfect after effective contact with individuals in the active TB classes or lost sight class at a rate l e ¼ s 1 nðI; J; LÞ; where σ 1 is the factor reducing the risk of infection as a result of acquiring immunity for latently infected individuals. Among latently infected individuals who become infectious, the fraction h is diagnosed and treated under the "Stop TB" program [29], while the remaining fraction 1−h is not diagnosed and becomes undiagnosed infectious. We assume that, after some time suffering from TB, some undiagnosed infectious people decide to go to hospital with a rate θ. Also, we assume that among diagnosed infectious people who begin their treatment therapy, a fraction r 2 of I take the full dose and make all the sputum examinations and will be declared cured from the disease. Some diagnosed infectious people who do not finish their dose of drugs and sputum examinations or whose treatment was unsuccessful, will not return to the hospital for the rest of sputum examinations and check-up. They will enter the class L of lost sight population at a constant rate α. Lost sight population can return to the hospital at a constant rate δ.
As suggested by Murray et al. [30], recovered individuals can only have partial immunity. Hence, they can undergo a TB reactivation or relapse with a constant rate γ. The remaining people can be reinfected (exogenously) after an effective contact with individuals in the active TB classes and lost sight class at a rate l r ¼ s 2 nðI; J; LÞ; where σ 2 is the factor reducing the risk of infection as a result of acquiring partial immunity for recovered individuals. Due to their own immunity, traditional medicine, natural recovery and drugs bought in the street (practiced in sub-Saharan Africa), a fraction of lost sight and undiagnosed infectious population can spontaneously recover at constant rates ρ and ω and enter the latent class E and recovery class R, respectively. The whole compartment model is schematically presented in a flow diagram in Fig. 1. This diagram gives rise to the following mathematical model of ordinary differential equations: where

Sensitivity analysis and parameter identification
In order to be able to assess targeted public health education strategies and chemoprophylaxis against TB spread in a population, we want to test the suitability of the model by fitting it to data from Cameroon [27]. For this purpose, we briefly describe the mathematical techniques that we use for parameter identification. In mathematical short-hand notation, the system of differential Equations (2) can be formally written as where p is the vector of parameters, and the right-hand side f depends on both the states, y 2 R n , and the parameter vector, p 2 R q . The initial condition vector, y 0 , has the same dimension as the state vector y. The TB model (2) can be written in the form of Equation (3) Assume there are given m experimental measurement time-point τ 1 , Á Á Á, τ m , and corresponding data values z j 2 R n , j = 1, Á Á Á, m associated with corresponding measurement tolerances δz j 2 R n . For ease of presentation, these tolerances are assumed to be positive, but the algorithm to be described is also able to tackle strictly zero tolerances (that indicate equality constraints); details are omitted here, but can be found, e.g., in the book [25]. Parameter identification consists of solving the least-squares minimization problem with diagonal weighting That means we want to minimize the relative deviation of model and data at the measurement time points τ j . Again in short-hand notation, the minimization problem (4) can be written as where F(p) = (F 1 (p), . . .,F m (p)) is a vector of length N = mÁn with entries defined by F:R q ! R N is a non-linear mapping and structured as a stacked vector. If not all components of a measurement z j are given, the number N is accordingly smaller, N < nm. The above problem (6), which is highly nonlinear in p, can be solved by affine covariant Gauss-Newton iteration, see [25], where each iteration step k requires the solution of a linear least-squares problem, where J(p k ) = F 0 (p k ) 2 R N×q denotes the Jacobian matrix or, equivalently, the sensitivity matrix. Its elements are the sensitivities computed at all data time-points. These entries represent the sensitivity of the solution y with respect to the parameters p at the time points of measurements. In model (2), values of model parameters and population classes vary over orders of magnitude. To achieve comparability, the sensitivity values have to be normalised by the absolute values of species and parameters to obtain scaled sensitivities whereby p scal and y scal are scaling values specified by the user. An analysis of the scaled matrix J(p) gives some hints whether the current combination of model and data will allow an identification of a given parameter. Parameters with very small sensitivity have nearly no influence on the solution at the measurement time-points and therefore cannot be estimated. In this case the entries of the corresponding column in J(p) are almost zero. Furthermore, some of the parameters might be linearly dependent, which leads to nearly identical columns in J(p). In both cases the matrix J(p) will be singular or, from a numerical point of view, "nearly" singular. In order to reveal such properties, the linear least squares problem (8) is solved by QR factorization with column pivoting [31]. By a suitable permutation of the columns of the matrix J(p), the diagonal elements of the upper triangular matrix R can be ordered in the form jr 11 j ! jr 22 j ! Á Á Á ! jr qq j ! 0 : As a measure of the term "nearly singular", the sub-condition of parameter p j is defined by Thus, the permutation of matrix columns corresponds to a new ordering of parameters according to increasing sub-condition. The sub-condition indicates whether a parameter can be estimated from the given data or not. Only those parameters can be estimated for which where e is the relative precision of the Jacobian J(p) [32]. Herein, the Jacobian has been computed with a finite difference scheme, resulting in a precision of ffiffiffiffiffiffiffiffiffi mach p , whereby mach = 2 −52 is the machine precision, i.e. the floating-point relative accuracy for double precision. Thus we set e = 10 −7 throughout the numerical computations in this paper.
The iterative scheme (8) can be generalized to the so-called global Gauss-Newton method by introducing damping factors λ k , p kþ1 ¼ p k þ l k Dp k ; 0 < l k 1: The step length λ k is computed successively in each iteration by an adaptive trust-region method [25]. This method for solving a non-linear least squares problem is implemented in the software code NLSCON [25] and part of the software package BioPARKIN [33]. Here, a Matlab-based version of this software package, named POEM 2.0, which is especially adapted to parameter identification in ordinary differential equation models, has been used.
As the data set for parameter identification, we used figures for diagnosed infectious and total population in Cameroon over the period 1994-2010 as published by WHO [1]. The data are listed in Table 1.
Overall, parameter identification is an iterative process because the sub-conditions depend on the current parameter values. At the beginning, all parameter values are fixed to their initial guesses and those parameter values are estimated that have the lowest sub-condition. These parameters are then fixed to their estimated values. Afterwards, parameter estimation is performed with a restricted list of previously non-identifiable parameters. After this phase, the estimation process is repeated with the full set of estimated parameters to check, whether the new values of previously non-identifiable parameters affect values of the identifiable ones. This process is repeated until convergence in the parameter values is achieved.

Results and Discussion
In the following, we show that the model is mathematically well-posed and epidemiologically reasonable [34]. Furthermore, we present the values of model parameters and discuss the results of the sensitivity analysis. In addition, we analyze the effect of increasing access to TB treatment as a result of improved infrastructures and education.

Basic properties of the model
Since model (2) monitors a human population, all its associated parameters and state variables should be non-negative and bounded for all t ! 0. Positivity of the solution. The following result shows that state variables are non-negative and dissipative. with N(t) = S(t)+E(t)+I(t)+J(t)+L(t)+R(t).
The proof of this Lemma follows the ideas from [35,36]. The following steps establish the invariance of the set O r ¼ ðS; E; I; J; L; RÞ 2 R 6 þ ; NðtÞ i.e. solutions remain in O ρ for all t ! 0. This implies that the trajectories of model system (2) are bounded. On one hand, integrating the differential inequality :N L À mN yields NðtÞ Nð0Þe Àmt þ L m ð1 À e Àmt Þ: In particular NðtÞ L m if Nð0Þ L m . On the other hand, if Nð0Þ ! L m , then Λ−μN(0) 0, and _ N ð0Þ L À mNð0Þ 0; i.e. the total population N(t) will decrease until NðtÞ L m : Thus, the simplex O ρ is a compact forward invariant set for model system (2), and for ρ > 0 this set is absorbing. So, we limit our study to this simplex for ρ > 0. The prevalent existence, uniqueness and continuation results hold for model system (2) Basic reproduction number. The global behavior of the TB model crucially depends on the basic reproduction number, i.e., the average number of secondary cases produced by a single infective individual who is introduced into an entirely susceptible population. Obviously, model system (2) has an equilibrium Q 0 = (x 0 , 0) with x 0 = Λ/μ when I = 0. This equilibrium point is the disease-free equilibrium (DFE). On this basis, we are now able to calculate the basic reproduction number R 0 , using the next generation method developed in [37]. For that purpose, let us write system (2) in the form _ x ¼ φðxÞ À nðI; J; LÞx; Vðx; yÞ ¼ nðI; J; LÞ½B 2 he 1 j yi þ B 3 he 5 j yi þ A y: Herein, x = S 2 R + is a state representing the compartment of non transmitting individuals (susceptible), y ¼ ðE; I; J; L; RÞ T 2 R 5 þ is the vector representing the state compartment of different infected individuals, and nðI; J; LÞ ¼ he 1 j yi N þ he 2 j yi; N ¼ x þ y 1 þ y 2 þ y 3 þ y 4 þ y 5 is the force of infection. hÁjÁi is the usual scalar product and A is the constant matrix ; with A 1 , A 2 , A 3 , A 4 and A 5 defined as above in (2). Furthermore, we set e 3 ¼ ð1; 0; 0; 0; 0Þ 2 R 5 ; e 4 ¼ ð0; 0; 0; 0; 1Þ 2 R 5 ; We define the Jacobian matrices at the DFE as F ¼ @F @y ðQ 0 Þ and V ¼ @V @y ðQ 0 Þ: Using the same notations as in [37], the basic reproduction number is given by the spectral radius of FV −1 , For model (14), one has and V ¼ ÀA: Then, according to [38,39], the basic reproduction ratio is given by The following result is established from [37]. Lemma 2: The disease-free equilibrium Q 0 of model system (2) is locally asymptotically stable whenever R 0 < 1, and unstable, if R 0 > 1.
From a biological point of view, Lemma 2 implies that TB can be eliminated from the community (when R 0 1) if the initial size of the population is in the basin of attraction of Q 0 . But if R 0 > 1 the infection will be able to spread through the population. Generally, the larger the value of R 0 , the harder it is to control the epidemic.
Parameter values of the model Fig. 2 shows the column norms kJ(p)(:, j)k 2 of the sensitivity matrix evaluated at the final iterate, i.e. for the final set of parameter values from Table 2. Obviously, the parameter with the largest column norm is p 2 , which means that it has the largest influence on the solution trajectory at the given measurement time-points. In fact, Fig. 2 shows that p 2 is the parameter with lowest subcondition, i.e. it can be best identified from the given measurements. In contrast, parameters σ 1 and p 1 are the two parameters with smallest column norm and largest subcondition. Note, however, that the magnitudes of the sensitivities alone, represented by the column norms in Fig. 2, are only an indicator for identifiability but generally do not provide information on the correct ordering of parameters. Moreover, the sensitivities and subconditions change throughout the modelling process whenever parameter values are changed.
Parameters not contained in the two figures have been estimated by WHO, are well-known for Cameroon [24] or have fixed values according to the National Institute of Statistic of Cameroon (NIS) [40]. In the following, we specify all parameter values.
The natural mortality μ: It is postulated to be equal to the inverse of the life expectancy at birth. We fix μ = 1/53.6 per year according to [40].
The recruitment Λ: According to NIS [40], the average recruitment in the Cameroonian population during the last fifteen year is fixed to Λ = 679685 per year.
TB mortality d 1 , d 2 and d 3 of undiagnosed infectious and lost sight population: Per capita TB-induced mortality rate is 0.193 per year in developed countries, but could be as high as 0.45 per year in some African countries [41]. We fixed the yearly TB-induced mortality rates d 1 = 0.193, d 2 = 0.20 and d 3 = 0.413 TB active cases according to [8,24].
Progression rate parameters p 1 , p 2 and k: For HIV-negative TB people, Bacaer et al. [20,42] estimated that people in a South Africa township have 11% annual risk of developing primary TB disease during five years following the first Mtb infection and a 0.03% annual risk of reactivation after five years. In Cameroon, the estimated average TB prevalence for all forms in HIVpositive is about 431 per 100,000 per year. Starting with this order of magnitude, we estimated that p 1 = 9.50082Á10 −04 , p 2 = 2.38932Á10 −02 and k = 3.31390Á10 −04 per year. Due to the limited access to hospitals, p 2 is expectedly larger than p 1 and k.
Factors σ 1 and σ 2 : Sutherland et al. [43] estimated that a previous Mtb infection reduces the risk of disease after reinfection by 63% for HIV negative males and by 80% for HIV negative females. We estimate that σ 1 = 2.38390Á10 −04 and we use the formula from [20] to set σ 2 = 0.7 Ã (p 1 +p 2 ). Running the simulation with other values for the pre-fractor, e.g. 0.2 or 0.9, leads only to small changes in the simulation result, indicating that σ 1 is indeed a less sensitive parameter that is difficult to estimate from the given data.
Detection rate h: According to WHO data, h 2 [0.5, 0.9] per year. Using POEM, it was estimated h = 0.812279 per year.
Diagnosis rate θ: WHO estimated θ 2 [0.3, 0.6] per year. The results in Fig. 2 show that the model is highly sensitive to θ, which was finally estimated from the data with POEM as θ = 0.497358 per year.
Proportions r 1 and r 2 of successful treatments: Since the chemoprophylaxis is not practiced in Cameroon, we took r 1 = 0 per year and fixed r 2 = 0.758821 per year according to [24].
Rate α at which diagnosed infectious people become lost sight: It has been estimated using POEM as α = 0.215698 per year. Rate δ at which lost sight people return to the hospital: According to the data of TB in Cameroon [24], we fixed δ = 0.309 per year.
Natural recovery rate ρ: In [20], the authors estimated that the natural recovery for HIVnegative TB and HIV-positive TB cases are 0.139 and 0.24 per year, respectively. Herein, we took the average of these values as initial guess for the Gauss-Newton algorithm and estimated ρ = 0.131140 per year.
Recovery rate ω: We fixed ω = 0.5 per year. A reference cannot be given for this value. A decrease (increase) of this value leads to a slight increase (decrease) in the maximum value of population class L but does not qualitatively affect the dynamics of the other population classes because the average population level of L is small.
Relapse rate of recovered individuals: The average relapse rate of recovered individuals is estimated with POEM as γ = 0.0851257 per year.
Numerical values of all parameters are summarized in Table 2. Note that some of the parameters are not specific for Cameroon, but specific for TB or valid for developing countries or sub-Saharan Africa in general. This especially applies to those parameters taken from the studies of Blower et al. [8,41] and Nunn et al. [44]. However, since other parameter values are Cameroon-specific, the model in fact represents a study for Cameroon.

Interpretations in terms of model simulations
In order to illustrate the theoretical results of the foregoing analysis, numerical simulations of model system (2) are carried out using a fourth order Runge-Kutta scheme in the software Matlab, version R2009. The total population of Cameroon in 1994 is given by N = 13,240,337 [40]. The initial values of the other variables were set as in Table 3. We deliberately decided not to fix the population numbers N(t) as equality constraints because the data do not represent true values but rather estimates. Using the data from Table 1, model (2) gives a very good fit to the Cameroonian data for the period 1994-2010 [1], as depicted in Fig. 3. The agreement between data and simulations is not perfect, but there are no evident consistent patterns in the discrepancy. Forward solutions of the deterministic model follow fairly well the observed TB patterns of incidence.
In addition, the obtained dynamics for the other population classes is plausible as well, see    Table 2 and initial values are presented in Table 3. doi:10.1371/journal.pone.0120607.g003 To summarize, with the estimated transmission parameters, the deterministic model appears to capture all the qualitative properties of the observed pattern. Hence, model (2) can be used to gain realistic insight into tuberculosis transmission dynamics at least for a limited period.

Effects of increased access to treatment
Herein, we investigate the impact of the time variation of some specific parameters on the dynamics of model (2). For this purpose, some model parameters are considered as time dependent variables to reflect their possible change within time. However, the variation is assumed to be slow over time.
Effects of increasing the access to TB treatment as a result of improved infrastructures and education are explored by taking into account the following expressions of model parameters, yðtÞ ¼ỹ þ ð1 ÀỹÞt y d þ t ; dðtÞ ¼d þ ð1 ÀdÞt d d þ t ; Herein, t starts at 0, counting the years since the beginning of the strategy, and θ δ , δ δ as well as  Table 2 and initial values are presented in Table 3.
doi:10.1371/journal.pone.0120607.g004 p δ are positive constants to be estimated. The valuesỹ;d;p 1 andp 2 represent the original parameter values from Table 2. The time-dependent functions have been constructed in such a way that they coincide with these values at time t = 0, given that θ δ , δ δ and p δ are larger than zero. Moreover, lim t ! 1 δ(t) = 1 and lim t ! 1 θ(t) = 1, i.e. the rates δ(t) and θ(t) are bounded from above by 1 as necessary. The parameters θ δ and δ δ determine how fast this limit is approached; the smaller their values, the faster the increase in θ(t) and δ(t), respectively. In addition, the functions p 1 (t) and p 2 (t) satisfy p 1 ðtÞ þ p 2 ðtÞ ¼p 1 þp 2 for all t in order to leave the value of parameter σ 2 unchanged. They have been chosen such that lim t!1 p 1 ðtÞ ¼p 1 þp 2 and lim t ! 1 p 2 (t) = 0. Again, the smaller p δ , the faster this limit is approached. Overall, the functions in Equation 18 are the simplest rational functions, i.e. with lowest polynomial order, that satisfy all the conditions discussed above. These functions are assumed to be control functions for the dynamics of TB. Thus, model (2) becomes a non-autonomous controlled system. The goal now is to find two different sets of values for parameters θ δ , δ δ and p δ such that the following two scenarios can be achieved: a) a reduction in the population of undiagnosed infectious (J) and lost sight (L) by 20% until 2025, b) a reduction in the population of undiagnosed infectious (J) and lost sight (L) by 60% until 2025.
To estimate the parameter values, we generated two artificial data sets. Both data sets contain the values of J and L from the previous simulation at the adaptively chosen time-points in the time interval [1994,2035]. However, the values in the time interval [2010,2035] have been reduced by 20% and 60% to obtain the first and second data set, respectively. These data sets were then used in the Gauss-Newton algorithm to obtain the following estimates: a) reduction by 20%: p δ = 8.56043Á10 7 , θ δ = 81.2807, δ δ = 37.2240 b) reduction by 60%: p δ = 85.6660, θ δ = 81.3256, δ δ = 37.1301 The simulation results are presented in Figs. 5 and 6. They show the dynamic of TB inside the population in the presence and absence of continuous effort to diagnose the population. In particular, a relatively small increase in the access to TB treatment could generally result in a decrease in the number of susceptible (S), diagnosed infectious (I), undiagnosed infectious (J), lost sight (L), and recovered (R) individuals, and an increase in the number of latently infected (E) individuals. We also observe that the number of diagnosed infectious people (I) increases at the beginning, but decreases after a few years prior to the beginning of the control strategies. Thus, TB can be reduced within 15 years if some efforts are made to increase the treatment access for rural population, and TB prevention and education for fast and immuno-compromised people.

Conclusion
In this paper, a deterministic model for the transmission dynamics of TB in Cameroon has been presented. The objective is to determine the role of TB diagnosis, treatment, lack of information about the epidemiological status of some people, and the role of traditional medicine and natural recovery on the dynamics of TB. In contrast to other TB models in the literature, the model includes three infective classes emanating from diagnosed infectious, undiagnosed infectious, and lost sight individuals. The undiagnosed and lost sight subclasses are shown to be of particular importance in TB modeling in developing countries like sub-Saharan Africa where public health is under-developed. Model parameters have either been fixed according to  (2) showing the impact of a slow change on parameter valuesθ, δ, p 1 and p 2 with respect to time in order to reduce the TB burden by 20% within 15 years. Solid lines present the model predictions for TB dynamics using parameter values of Table 2 and the dashed lines present the trajectories for parameters θ, δ, p 1 andp 2 set as in Equation (18). Parameter identification with artificial data gavep δ = 8.56043Á10 7 , θ δ = 81.2807, δ δ = 37.2240. All other parameters are defined as in Table 2.   Table 2 and the estimated trajectories (dashed lines) for parameters θ, δ, p 1 and p 2 set as in Equation (18). Parameter identification with artificial data gave p δ = 85.6660, θ δ = 81.3256, δ δ = 37.1301.
doi:10.1371/journal.pone.0120607.g006 data published in literature, or they have been estimated with a Gauss-Newton method using data published by WHO and the National Institute of Statistics of Cameroon (NIS). To the knowledge of the authors, it is the first time that a detailed sensitivity analysis is performed on a TB model in order to infer identifiability of unknown model parameters. In particular, parameters representing the proportion of individuals having access to medical facilities have a large impact on the dynamics of the disease. We showed that a change in these parameters over time can significantly reduce the disease burden in the population within the next 15 year. These parameters can be used to measure the success of educational and diagnosis campaigns that encourage individuals to go for TB screening. In future work, optimal control strategies could be applied to determine the optimal dynamics of these parameters in order to achieve the highest possible reduction of TB in shortest time at low costs. In addition, the model might be extended towards the inclusion of co-infection between TB and HIV, or the resistance to treatment.