The weakest link: uncertainty and sensitivity analysis of extinction probability estimates for tsetse (Glossina spp) populations

Background A relatively simple life history allows us to derive an expression for the extinction probability of populations of tsetse, vectors of African sleeping sickness. We present the uncertainty and sensitivity analysis of extinction probability for tsetse population, to offer key insights into parameters in the control/eradication of tsetse populations. Methods We represent tsetse population growth as a branching process, and derive closed form estimates of population extinction from that model. Statistical and mathematical techniques are used to analyse the uncertainties in estimating extinction probability, and the sensitivity of the extinction probability to changes in input parameters representing the natural life history and vital dynamics of tsetse populations. Results For fixed values of input parameters, the sensitivity of extinction probability depends on the baseline parameter values. For example, extinction probability is more sensitive to the probability that a female is inseminated by a fertile male when daily pupal mortality is low, whereas the extinction probability is more sensitive to daily mortality rate for adult females when daily pupal mortality, and extinction probabilities, are high. Global uncertainty and sensitivity analysis showed that daily mortality rate for adult females has the highest impact on the extinction probability. Conclusions The strong correlation between extinction probability and daily female adult mortality gives a strong argument that control techniques to increase daily female adult mortality may be the single most effective means of ensuring eradication of tsetse population. Author summary Tsetse flies (Glossina spp) are vectors of Trypanosomiasis, a deadly disease commonly called sleeping sickness in humans and nagana in livestock. The relatively simple life history of tsetse enabled us to model its population growth as a stochastic branching process. We derived a closed-form expression for the probability that a population of tsetse goes extinct, as a function of death, birth, development and insemination rates in female tsetse. We analyzed the sensitivity of the extinction probability to the different input parameters, in a bid to identify parameters with the highest impact on extinction probability. This information can, potentially, inform policy direction for tsetse control/elimination. In all the scenarios we considered, the daily mortality rate for adult females has the greatest impact on the magnitude of extinction probability. Our findings suggest that the mortality rate in the adult females is the weakest link in tsetse life history, and this fact should be exploited in achieving tsetse population control, or even elimination.


Introduction
Tsetse flies (Glossina spp) are biting flies of both public health and economic 16 importance in 36 Sub-Saharan Africa countries. They feed exclusively on the blood of 17 vertebrates -game animals and livestock, and also humans, and provide the link that 18 drives the transmission of African trypanosomiasis, a tropical disease called Sleeping 19 Sickness in humans and nagana in livestock. According to a WHO 2018 factsheet for 20 Human Sleeping Sickness, the disease still occurs in about 36 countries in sub-Saharan 21 Africa, mostly among poor farmers living in rural areas. Due to sustained control 22 efforts, the number of cases of the disease has reduced. In 2015 there were about 2804 23 cases recorded: 97% of these were chronic infections with Trypanosoma brucei 24 gambiense [1]. To sustain the reduction in cases, it is important to continue to improve 25 understanding of the tsetse fly vector, in a bid to develop more effective control 26 techniques: with improved cost effectiveness.

27
A recent study [2] employed the theory of branching processes to derive an expession 28 for the extinction probability for closed populations of tsetse. This equation involves 29 numerous parameters representing death, development and fertility rates during the 30 fly's lifecycle. These results allow us to determine, by sensitivity analysis, the relative 31 importance of changing, through control techniques, the various parameters. Sensitivity 32 analysis is often used to investigate the robustness of model output to parameter 33 values [3][4][5]. In this context, it is important to identify the parameters that have the 34 greatest influence on extinction probabilities of tsetse, since this information will provide 35 insight to the eradication of tsetse, and inform policy on the direction of control efforts. 36 Here we adopt the model developed by Kajunguri et al [2] and Hargrove [6] for the 37 reproductive performance of female tsetse flies inseminated by a fertile male. We then 38 use a framework, developed by Harris [7], to derive a fixed point equation for the 39 extinction probability for a tsetse population. This approach allows us to obtain the 40 same expression for extinction probability as [2], but it is derived with fewer steps and 41 with less mathematical complexity. We compute local sensitivity indices of extinction 42 probability with respect to all input parameters by allowing the value of daily mortality 43 rate of female pupae (χ) to vary between 0.001 per-day and 0.025 per-day. Due to 44 nonlinearities and interdependencies between input parameters, local sensitivity may be 45 highly dependent on the baseline values of the parameters [8]. uncertainty and the sensitivity of the extinction probability to all input parameters.

56
In the next section, we present the branching process model developed in [2] and [6] 57 and present an approach based on a method used in [7] to derive a fixed point equation 58 for the extinction probability for a tsetse population. In section 3, we present the local 59 sensitivity analysis for the extinction probability at two fixed baseline values of the 60 input parameters and the mathematical derivation of the sensitivity indices of 61 extinction probability w.r.t all input parameters. Section 4 presents global uncertainty 62 and sensitivity analysis using LHS/PRCC methods. The results are discussed in detail 63 in section 5.

65
The aim of this section is to develop a stochastic model for tsetse population growth in 66 the form of a branching process and to use the model to obtain a fixed point equation 67 for extinction probability for tsetse populations ( [2], [6], [14], [15]). We develop the 68 branching process focusing only on female tsetse flies [6]. We follow a framework 69 developed in [6], assuming a female tsetse fly is fertilized with probability and survives 70 to deposit her first larva with probability λ ν+τ :ν is days to first ovulation, τ is the 71 inter-larval period, and λ is adult female daily survival probability. She produces a 72 female pupa with probability β, and the pupa survives to adulthood with probability φ g 73 (where g is the pupal duration and φ is the daily survival probability of the pupa). The 74 mother dies before the next pregnancy, having produced a single surviving daughter, 75 with probability (1 − λ τ ). The probability that an adult female tsetse dies after 76 producing a single surviving daughter after surviving one pregnancy is: Equation (1) can be generalized by induction to obtain the probability that a female 78 tsetse produces k surviving female offspring after surviving n pregnancies. Thus where n k = n! (n−k)!k! is the binomial coefficient.

81
surviving female offspring in her lifetime, respectively. Suppose also that p 0 + p 1 < 1, to 82 avoid the trivial case where a tsetse fly only produces 0 or 1 female offspring.

83
Summing equation (2) over all n, gives p k , the probability that a female produces k 84 surviving female offspring in its lifetime.
Equation (3) was used in [2] to obtain the mean and variance of the population size, 86 extinction probability and time to extinction of populations of tsetse. Proofs of 87 equations (1) and (2) are provided in [2] (Supplementary Information).

88
It can be shown easily that p 0 , p 1 , p 2 , ... follow a geometric series, such that Following a framework developed by Harris [7], the generating function g(θ) of p k , is 93 a fractional linear function given by; Extinction probability 95 The extinction probability for tsetse population is the non-negative fixed point of where βφ g λ τ = 0. In practice, 0 < β < 1, 0 < λ < 1 and 0 < φ < 1. This implies that In this section, we perform local sensitivity analysis, otherwise known as elasticity 106 analysis, on the extinction probability for tsetse populations. Given that the extinction 107 probability θ, depends differentiably on each input parameter, the normalized forward 108 sensitivity (elasticity) index of θ w.r.t all input parameters is: where ρ i is the set of all input parameters of the extinction probability. This method 110 has been used extensively in the literature to determine the sensitivity of the 111 reproduction number R o of epidemiological models to model parameters [4,5,16]. When 112 the initial population consists of N female tsetse, the extinction probability is θ N . The 113 sensitivity indices of θ N w.r.t all input parameters is; Notice that, when there are N female flies in the initial population, the sensitivity 115 indices of θ N w.r.t all input parameters is the sensitivity indices of θ multiplied by N . 116 The larger the size of the initial population, the more sensitive extinction probability is 117 to input parameters.

118
Writing equation (6) in terms of daily mortality rate for adult females (ψ), and the 119 daily mortality rate for female pupae (χ), yields: (9) Table 1 shows the derivations of the sensitivity indices of extinction probability with 121 respect to all seven input parameters. These expressions were derived from equations 122 (7) and (9) with a simple code in MAPLE 17 environment.
123 Table 1. Expressions for the sensitivity indices of extinction probability with respect to all seven input parameters.

Parameters
The sensitivity of extinction probability (θ) to input parameters

Results
124 Table 2 shows the sensitivity indices of extinction probability w.r.t all input parameters 125 for different values of extinction probabilities. For instance, the sensitivity indices of θ 126 w.r.t to (probability female is inseminated by a fertile male) decreases by > 60% when 127 θ (extinction probability) approaches 1, implying that, at θ = 0.419, a 10% decrease in 128 will yield a 22% increase in θ, whereas, at θ = 0.96, a 10% decrease in will only yield 129 an 8.7% increase in θ. Table 2. List and description of parameters affecting extinction probabilities for tsetse populations, and the sensitivity indices for these parameters, at two different values of extinction probability. Here we investigate the changes that occur in the sensitivity indices of extinction 133 probability with respect to six input parameters by allowing χ to vary between 0.1% to 134 2.5%. A simple script was written in MAPLE 17 environment to calculate the local 135 sensitivity indices of θ w.r.t to the six remaining input parameters for different values of 136 χ. Figure 1(A and B) show changes in the sensitivity indices of θ w.r.t to each 137 parameter as the daily mortality rate for female pupae (χ) varies from 0.1% to 2.5%,

139
As χ increases from 0.001 to 0.0065, the sensitivity index of θ w.r.t reduces below 140 the sensitivity index of θ w.r.t ψ. At that point extinction probability becomes more 141 sensitive to ψ than . When χ increases further to 0.013, the sensitivity of extinction 142 probability to drops further below the sensitivity of extinction probability to τ (Fig 1 143  (A and B)).

144
Local sensitivity analysis may not be robust enough to capture the actual influence 145 of all input parameter values on the extinction probability since there are 146 interdependencies between input parameters. We, therefore, proceed to carry out global 147 uncertainty and sensitivity analysis of the extinction probability for tsetse population. 148

Global uncertainty and sensitivity analysis of θ 149
The exact values of the input parameters are not known in field suituatins, where many 150 of these parameters depend on temperature and other climatic factors. It is therefore 151 important to quantify the uncertainty involved in estimating the extinction probability. 152 To quantify the uncertainty involved in estimating the extinction probability (θ) and to 153 establish the most important input parameters, we use LHS and PRCC methods for the 154 global uncertainty and sensitivity analysis of the extinction probability. The method 155 follows the approach of Samsuzzoha et al [4].

7/12
Uncertainty analysis 157 We aim to analyse the uncertainty involved in quantifying extinction probability (θ) 158 based on the uncertainties associated with the input parameters. Accordingly, in order 159 to investigate the sensitivity to this uncertainty we sample values from prior  Table 3, where β, N and U denote beta, normal and uniform distributions, respectively. 164 Table 3. List of parameters and their prior probability distributions .
Parameters Prior probability distribution Using LHS, we obtain the uncertainty output for all the input parameters and also 165 for the extinction probability. LHS is used to sample from the stratified probability 166 distribution functions for different parameters. Using 1000 intervals of equal 167 probabilities. Figure 2 shows the uncertainty output for all the input parameters and 168 the shape of their probability distribution together with their summary statistics. The 169 uncertainty output for extinction probability (θ) shows that it is beta distributed with 170 mean = 0.415 and standard deviation = 0.386.

172
To identify key input parameters, we carry out a sensitivity analysis by calculating the 173 PRCC between each input parameter and the extinction probability. The parameter 174 with the highest PRCC has the largest influence on the magnitude of the extinction 175 probability. Figure 3 shows the PRCC outputs for all input parameters, where the 176 probability ( ) that a female fly is inseminated by a fertile male is essentially equal to 1. 177 In the field, males manage to find and mate with females, even when population levels 178 are quite low [19]. For most tsetse populations, the probability of insemination is thus 179 close to 1. Accordingly, we allow to vary between 0.999-1. In figure 5, the prior 180 probability distributions are kept the same save for which is sampled between 0.885 181 and 1 [10,11,20].

182
LHS is used to sample from the prior probability distributions, where is sampled 183 from a uniform distribution U (0.999, 1). Figure 3 shows that daily mortality rate for 184 adult females (ψ) has a strong correlation with the extinction probability with PRCC 185 score 0.91, followed by daily mortality rate for female pupae (χ) and inter-larval period 186 (days) (τ ) having PRCC scores 0.47 and 0.058, respectively.

187
The female tsetse fly generally mates only once in her lifetime, storing the sperm in 188 spermathecae and using small amounts to fertilize her eggs one at a time [21,22]. When 189 sterile males are introduced into tsetse population, the probability ( ) that a female is 190 inseminated by a fertile male falls below unity, by an amount that depends on the ratio 191 of sterile to fertile males in the population.

192
The Sterile Insect Technique (SIT) has been used in attempts to control tsetse flies 193 populations [23,24] and was used to eradicate a small population of G. austeni on 194 Unguja Island, Zanzibar, Tanzania [25]. The probability that a female is inseminated by 195 a sterile male is 1 − . We now allow baseline values of to vary over a wide range, in 196 order to assess the sensitivity of extinction probability to changes in at varying The simple life history of the tsetse fly enabled us to model its population dynamics as 204 a stochastic branching process. We derived an expression for the extinction probability 205 for tsetse populations and performed local and global sensitivity analyses, as well as 206 global uncertainty analysis, on the extinction probability. We calculated all results for 207 two fixed baseline values for χ, corresponding to values that resulted in low or high 208 extinction probabilities. We obtained the sensitivity indices of the extinction probability 209 to seven input parameters. When the extinction probability (θ) is fixed at either low or 210 high levels (0.419 or 0.960) θ is more sensitive to changes in daily adult mortality (ψ)  These techniques can be pooled into two fundamental control philosophies -those which 218 aim, primarily, to increase mortality rates in adult flies and those, like SIT, which aim 219 to reduce tsetse birth rates [24]. Sensitivity analysis will indicate which parameter out 220 of the two has the highest impact on the extinction probability.

221
From Table 2, observe that the sensitivity indices of θ to the input parameters 222 depends on the value of the extinction probability. We allowed the daily mortality rate 223 for pupae (χ) to vary from 0.001 to 0.025. The lower and upper bound values result in 224 low and high extinction probabilities, respectively. We then calculated the sensitivity 225 indices of θ w.r.t. the remaining six parameters. Figure 1(A and B) shows that the 226 sensitivity of θ to each of the input parameters changes as extinction probability 227 increases with increasing values of χ. Observe that for χ ≥ 0.018, the sensitivity indices 228 of all the six parameters converged to zero. This is expected since the set baseline 229 parameters values for all input parameter will correspond to extinction probability . We defined prior probability density functions for the seven input 235 parameters and we sampled from intervals of equal probability using LHS. The PRCC 236 score of all input parameters was obtained for three sets of the probability distribution 237 function, fixed for six parameters and varied only for . In all cases, ψ has the strongest 238 impact on the extinction probability. The PRCC score for increases as we allow for 239 more variability in its prior probability distribution.

240
The SIT is an effective technique used to suppress or eradicate populations of tsetse, 241 but its major drawback is the large number of sterile flies that have to be produced and 242 introduced into the wild [23]. Our results confirm this; the higher the number of sterile 243 males introduced into the wild, the higher the impact of on the extinction probability 244 (Figs 3-6).

246
In all scenarios considered, control techniques which can achieve high mortality rates for 247 adult female flies have the strongest impact on extinction probability. Control 248 techniques such as SIT, which can reduce reproductive rates, without increasing 249 mortality, can also have a strong impact on extinction probability. This happens only 250 when the number of sterile males, introduced into the population, massively outnumber 251 wild males, such that the probability is low that a virgin female will mate with a fertile 252 male.

253
A limitation of our work is the assumption that the tsetse flies experience fixed 254 environmental conditions throughout their life history. This assumption is not true in 255 the wild, where tsetse experience daily and seasonal changes in various climatic effects. 256 We will address this problem in future work. 257 Fig 1. Variation in the sensitivity of extinction probability θ to six input parameters (β, , ν, g, ψ, τ ) as a function of the values of the background daily rate (χ) of pupal female mortality. (A) The sensitivity indices of extinction probability to six input parameters with signs. (B) The sensitivity indices of extinction probability to six input parameters, in absolute value. The arrow through B indicates the point where θ becomes more sensitive to ψ than . Fig 2. The uncertainty output for all input parameters, together with uncertainty output of the extinction probability, obtained from Latin hypercube sampling using a sample size of 1000 for the seven input parameters. Each parameter appears at the top of the corresponding sub-plot.