Mathematical modeling of trypanosomiasis control strategies in communities where human, cattle and wildlife interact

Abstract

Spillover of trypanosomiasis parasites from wildlife to domestic livestock and humans remains a major challenge world over. With the disease targeted for elimination by 2030, assessing the impact of control strategies in communities where there are human-cattle-wildlife interactions is therefore essential. A compartmental framework incorporating tsetse flies, humans, cattle, wildlife and various disease control strategies is developed and analyzed. The reproduction is derived and its sensitivity to different model parameters is investigated. Meanwhile, the optimal control theory is used to identify a combination of control strategies capable of minimizing the infected human and cattle population over time at minimal costs of implementation. The results indicates that tsetse fly mortality rate is strongly and negatively correlated to the reproduction number. It is also established that tsetse fly feeding rate in strongly and positively correlated to the reproduction number. Simulation results indicates that time dependent control strategies can significantly reduce the infections. Overall, the study shows that screening and treatment of humans may not lead to disease elimination. Combining this strategy with other strategies such as screening and treatment of cattle and vector control strategies will result in maximum reduction of tsetse fly population and disease elimination.

Introduction

Human African trypanosomiasis (HAT), also known as sleeping sickness continues to present a major human and animal public health challenge in developing countries such as Republic of Congo and Uganda (Kasozi et al. 2021). Worldover, it is estimated that approximately 70 million people in sub-Saharan Africa are at risk to HAT infections (Papagni et al. 2023). HAT is targeted for elimination by 2030 (Pepin 2023). The proclamation and desire to eliminate HAT by 2030 was made by the World Health Organization (WHO) in 2014 following a 73% reduction of HAT cases between 2000 and 2012 due to the implementation of sustainable control efforts (Franco et al. 2022).

As we draw towards 2030, the feasibility of eliminating HAT remains unclear, particularly in areas where inhabitants reside in communities close to game reserves (Kasozi et al. 2021). This is mainly because of cross transmission of parasites between livestock and wildlife (Kasozi et al. 2021). Issues of cross-transmission of the disease from animals to humans have prompted WHO to urge modellers to include animal populations in theoretical studies aimed at modelling HAT transmission dynamics (World Health Organization 2013, 2020).

Although there are plenty of mathematical models for HAT in literature (see, for example, Artzrouni and Gouteux (1996); Stone and Chitnis (2015); Hargrove et al. (2012); Medlock et al. (2013); Helikumi et al. (2021); Helikumi et al. (2020, 2019, 2020); Ndondo et al. (2016); Gilbert et al. (2016); Rock et al. (2015b); Rock et al. (2015a)), the inclusion of wildlife in these studies has largely been ignored. As we draw towards 2030, it is imperative to evaluate the feasibility of disease elimination based on current disease intervention strategies, especially in areas where there are human-cattle and wildlife interaction. Motivated by the aforementioned facts, in this study, we use a mathematical model to evaluate the impact of different HAT control strategies on minimizing the disease burden in communities where there are human-cattle and wildlife interactions.

Mathematical models are essential tools that have been used to understand the transmission dynamics of HAT since the last two decades (for example, Artzrouni and Gouteux (1996); Stone and Chitnis (2015); Hargrove et al. (2012); Medlock et al. (2013); Ndondo et al. (2016); Gilbert et al. (2016); Rock et al. (2015b); Rock et al. (2015a)). Aforementioned studies and those cited therein produced many useful results and improved existing knowledge on trypanosomiasis transmission and control. For example., Stone and Chitnis (2015) developed a model to explore the implications of heterogeneous biting exposure of the human population on HAT transmission dynamics. Their findings showed that heterogeneous biting exposure in the human porpulation has a strong influence on disease dynamics. Through a mathematical model, Hargrove et al. (2012) demonstrated that treatment of livestck using trypanocides or insecticides could significantly reduce annual HAT cases. Medlock et al. (2013) used a mathematical model to evaluate the potential impact of Wolbachia-colonized tsetse flies on reducing annual HAT cases.

Despite these efforts, however, several challenges remain in the mathematical modeling of HAT. One limitations of the aforementioned studeis the non-inclusion of wildlife animals in their frameworks. Modeling studies that neglect wildlife population may not clearly account for HAT dynamics in communities where there is human, cattle and wildlife interaction. Implementing disease control strategies in wildlife population is known to be challenging. Current HAT control strategies are mainly a combination of active and passive case-finding, and the treatment of individuals with detected infections (Franco et al. 2022). However, it is imperative to evaluate the feasibility of disease eradication in such communities based on the aforementioned intervention strategies. To evaluate the strength of the intervention strategies to control the disease, we consider time and non-time dependent control efforts. For time dependent controls, optimal control theory is utilized to formulate an optimal control problem in which the intervention efforts are optimized to minimize the number of infected humans and cattle at minimal costs of implementation. In the last two decades, optimal control theory has proved to be an essential tool that can be used to identify the best strategies for the control of infectious diseases.

The outline of the paper is as follows: in the next section, we present the methods and analytical results. In particular, we present the proposed model and perform dynamical analysis and present the proposed optimal control problem. The Methods, results and discussions section is followed by numerical simulations and discussions. The paper rounds up with a discussion of results and future outlook.

Methods, results and discussions

In this section, a mathematical model is presented to study the transmission of trypanosomiasis among humans, cattle, wildlife and tsetse flies, based on Ordinary differential equations (ODE).

Model formulation

We propose a deterministic model that incorporates the interplay of tsetse flies and multiple hosts-humans, cattle and wildlife. Throughout the study, we will use the subscript h, c and w to represent the human, cattle and wildlife populations. Let the variables, $S_{j}(t),$ $E_{j}(t)$, $I_{j}(t)$ and $R_{j}(t),$ $j=c,h,w,$ represent the number of susceptible, exposed, infectious and recovered host at time t,  such that the total population of each host is equivalent to $N_{_j}(t)=S_{j}(t)+E_{j}(t)+I_{j}(t)+R_{j}(t).$ The vector population is divided into three classes such that at time t, there are susceptible $S_{v}(t)$, exposed $E_{v}(t)$, and infected $I_{v}(t)$ vectors. Thus, the total vector population is $N_v(t)=S_{v}(t)+E_{v}(t)+I_{v}(t).$ We assume that susceptible hosts contact the disease upon being bitten by an infectious vector. Similarly, susceptible vectors become infected when they bite an infectious host. Hence, in this model, we consider the following forces of infection $\lambda _c(t)$, $\lambda _w(t),$ $\lambda _h(t)$ and $\lambda _v(t)$ for cattle, wildlife, human and vector populations, respectively,:

$$\begin{aligned} \left. \begin{array}{lll} \lambda _{c}&{}=&{} p_{c} au \frac{I_v}{N_c},\qquad \lambda _{w}=p_{w}au\frac{I_v}{N_w},\qquad \lambda _{h}=p_hab\frac{I_v}{N_h},\\[10pt] \lambda _{v}&{}=&{} p_{h}ac\frac{I_h}{N_h} + p_{c}av\frac{I_h}{N_c} + p_{w}av\frac{I_w}{N_w}. \end{array} \right\} \end{aligned}$$

(1)

In (1), $p_j$, $j=c,h,w$ correspond to the proportion of bloods meals on j host, such that $\sum _{j=c,h,w}p_j=1$, a is the vector biting rate, u is the probability of disease transmission from an infectious vector to a susceptible host (cattle and wildlife) given that a contact between the two occurs and b is the probability of transmission of infection from an infectious vector to a susceptible human given that a contact between the two occurs. Similarly, c and v accounts for the probability that a susceptible vector becomes infected following a blood meal from an infectious human and cattle, respectively.

Upon infection, the susceptible hosts progress to the exposed state where they remain for an average period of $q_j^{-1}$ days, after which they progress to the infectious state (stage 1 of the disease) where they remain for an average period of $\gamma ^{-1}_j$ days. Untreated hosts progress to stage II. Infected animals in stage II are non-infectious (Ndondo et al. 2016). Stage II infected animals recover with temporary immunity, which wanes out after $\kappa _{j}^{-1}$ days.

Similarly, upon infection vectors progress to the exposed state where they incubate the disease. The incubation period lasts $q_v^{-1}$ days and there after they become infectious for their entire life span $\mu _v^{-1}$ days. Natural mortality rate is assumed to be constant in all epidemiological classes for both the hosts and vectors. Mortality in humans, wildlife and cattle is modeled by $b_h$, $b_w$ and $b_c,$ respectively. Furthermore, let $b_hN_h$, $b_wN_w$, $b_cN_c$ and $b_vN_v,$ represent the recruitment rate into the human, wildlife, cattle and vector birth rate, respectively. All new recruits are assumed to be susceptible to infection.

Based on these assumptions, we summarize the transmission dynamics of the disease by a system of ordinary differential equations (ODEs) (2)-(16) while the model transmission diagram is depicted in Fig. 1:

$$\begin{aligned} S'_h(t)= & {} b_hN_h-p_hab\frac{I_v}{N_h}S_h-b_hS_h+\kappa _hR_h,\end{aligned}$$

(2)

$$\begin{aligned} E'_h(t)= & {} p_hab\frac{I_v}{N_h}S_h-(b_h+q_h)E_h,\end{aligned}$$

(3)

$$\begin{aligned} I'_h(t)= & {} q_hE_h-(b_h+\gamma _h)I_h,\end{aligned}$$

(4)

$$\begin{aligned} R'_h(t)= & {} \gamma _hI_h-(b_h+\kappa _h)R_h,\end{aligned}$$

(5)

$$\begin{aligned} S'_c(t)= & {} b_cN_c-p_{c}au\frac{I_v}{N_c}S_c-b_cS_c+\kappa _cR_c,\end{aligned}$$

(6)

$$\begin{aligned} E'_c(t)= & {} p_{c}au\frac{I_v}{N_c}S_c-(b_c+q_c)E_c,\end{aligned}$$

(7)

$$\begin{aligned} I'_c(t)= & {} q_cE_c-(b_c+\gamma _c)I_c,\end{aligned}$$

(8)

$$\begin{aligned} R'_c(t)= & {} \gamma _cI_c-(b_c+\kappa _c)R_c,\end{aligned}$$

(9)

$$\begin{aligned} S'_w(t)= & {} b_wN_w-p_{w}au\frac{I_v}{N_w}S_w-b_wS_w+\kappa _wR_w,\end{aligned}$$

(10)

$$\begin{aligned} E'_w(t)= & {} p_{w}au\frac{I_v}{N_w}S_w-(b_w+q_w)E_w,\end{aligned}$$

(11)

$$\begin{aligned} I'_w(t)= & {} q_wE_w-(b_w+\gamma _w)I_w,\end{aligned}$$

(12)

$$\begin{aligned} R'_w(t)= & {} \gamma _wI_w-(b_w+\kappa _w)R_w,\end{aligned}$$

(13)

$$\begin{aligned} S'_v(t)= & {} b_vN_v-p_hac\frac{I_h}{N_h}S_v-p_{c}av\frac{I_c}{N_c}S_v-p_{w}av\frac{I_w}{N_w}S_v-b_vS_v,\end{aligned}$$

(14)

$$\begin{aligned} E'_v(t)= & {} p_hac\frac{I_h}{N_h}S_v+p_{c}av\frac{I_h}{N_c}S_v+ p_{w}av\frac{I_w}{N_w}S_v-(b_v+q_v)E_v,\end{aligned}$$

(15)

$$\begin{aligned} I'_v(t)= & {} q_vE_v-b_vI_v. \end{aligned}$$

(16)

Fig. 1

Model flow diagram summarizing the dynamics of HAT infection presented in system (2)- (16)

All variables and model parameters of system (2)-(16) are considered to be positive and it can easily be established that the model is epidemiologically and mathematically well-posed in the domain (17):

$$\begin{aligned} \Theta =\bigg \{N_k(t)\le N_k\bigg \},\quad \text {for}, \qquad k = h,c,w,v. \end{aligned}$$

(17)

One can easily verify that this region, $\Theta$ is positively invariant with respect to system (2)-(16). Thus we will study the dynamics of our model in the closed set $\Theta .$ Based on (17), we one can simplify the analysis of model (1) by considering a dimensionless model since the total population is constant. Let;

$$\begin{aligned} s_k=\frac{S_k}{N_k},\ e_k=\frac{E_k}{N_k},\ i_k=\frac{I_k}{N_k},\ r_j=\frac{R_j}{N_j},\ \theta _{vh}=\frac{N_v}{N_h},\ \theta _{vc}=\frac{N_v}{N_c},\ \text {and}\quad \theta _{vw}=\frac{N_v}{N_w}. \end{aligned}$$

(18)

By direct calculation, one can easily verify that the dimensionless form of system (2)-(16) is equivalent to (19):

$$\begin{aligned} \left. \begin{array}{lll} s'_h(t)&{}=&{} b_h-p_hab\theta _{vh}i_vs_h-b_hs_h+\kappa _hr_h,\\ e'_h(t)&{}=&{} p_h\theta _{vh} abi_vs_h-(b_h+q_h)e_h,\\ i'_h(t)&{}=&{}q_he_h-(b_h+\gamma _h)i_h,\\ r'_h(t)&{}=&{}\gamma _hi_h-(b_h+\kappa _h)r_h,\\ s'_c(t)&{}=&{} b_c-p_{c}au\theta _{vc}i_vs_c-b_cs_c+\kappa _cr_c,\\ e'_c(t)&{}=&{} p_{c}au\theta _{vc}i_vs_c-(b_c+q_c)e_c,\\ i'_c(t)&{}=&{}q_ce_c-(b_c+\gamma _c)i_c,\\ r'_c(t)&{}=&{}\gamma _ci_c-(b_c+\kappa _c)r_c,\\ s'_w(t)&{}=&{} b_w-p_{w}au\theta _{vw}i_vs_w-b_ws_w+\kappa _wr_w,\\ e'_w(t)&{}=&{} p_{w}au\theta _{vw}i_vs_w-(b_w+q_w)e_w,\\ i'_w(t)&{}=&{}q_we_w-(b_w+\gamma _w)i_w,\\ r'_w(t)&{}=&{}\gamma _wi_w-(b_w+\kappa _w)r_w,\\ s'_v(t)&{}=&{} b_v-p_haci_hs_v- p_{c}avi_cs_v- p_{w}avi_ws_v-b_vs_v,\\ e'_v(t)&{}=&{} p_haci_hs_v+p_{c}avi_cs_v+p_{w}avi_ws_v-(b_v+q_v)e_v,\\ i'_v(t)&{}=&{}q_ve_v-b_vi_v. \end{array} \right\} \end{aligned}$$

(19)

It can easily establish that model (19) is epidemiologically and mathematically well-posed in the domain:

$$\begin{aligned} \Omega =\left\{ \begin{array}{c} \left( \begin{array}{c} s_j,e_j, i_j,r_j \\ s_v,e_v, i_v \end{array}\right) \in \mathbb {R}^{15}_{+}\left| \begin{array}{llll} 0\le s_j+e_j+ i_j+r_j\le 1,\\ 0\le s_v+e_v+ i_v\le 1. \end{array} \right. \end{array} \right\} , \qquad \text {for}\qquad j = c, h, w. \end{aligned}$$

(20)

Further, direct calculations show that system (19) has a disease-free equilibrium (DFE) given by (21):

$$\begin{aligned} \mathcal {E}^0= & {} \bigg (s_h^0,e^0_h,i^0_h,r^0_h,s_c^0,e^0_c,i^0_c,r^0_c,s_w^0,e^0_w,i^0_w,r^0_w,s_v^0,e^0_v,i^0_v\bigg )\nonumber \\= & {} \bigg (1,0,0,0,1,0,0,0,1,0,0,0,1,0,0\bigg ). \end{aligned}$$

(21)

In Eq. (21), we can observe that the total human, cattle, wild and vector populations are susceptible to the disease.

Disease transmission potential

The basic reproduction number $\mathcal {R}_0$ remains one of the integral threshold quantities in disease modeling. It demonstrates the disease transmission potential. In vector borne disease models, $\mathcal {R}_0$ represents the average number of secondary infections that single host can generate in a totally susceptible population of hosts and vectors. To compute the basic reproduction number associated with model (19) we will use the next generation matrix method presented in Van den Driessche and Watmough (2002). Let, the nonnegative matrix $\mathcal {F}$ represent the generation of new infection terms and the non-singular matrix $\mathcal {V}$ be the remaining transfer terms evaluated at DFE., that is.,:

$$\begin{aligned} \mathcal {F}= & {} \left[ \begin{array}{ccccccccc} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} p_h\theta _{vh} ab\\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0\\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} p_{w}au\theta _{vw}\\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0\\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} p_{w}au\theta _{vw}\\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0\\ 0 &{} p_hac &{} 0 &{} p_{c}av &{} 0 &{} p_{w}av &{} 0 &{} 0\\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \end{array}\right] ,\nonumber \\ \mathcal {V}= & {} \left[ \begin{array}{cccccccc} \tilde{q}_h &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0\\ -q_h &{} \tilde{\gamma }_h &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0\\ 0 &{} 0 &{} \tilde{q}_c &{} 0 &{} 0 &{} 0 &{} 0 &{} 0\\ 0 &{} 0 &{} -q_c &{} \tilde{\gamma }_c &{} 0 &{} 0 &{} 0 &{} 0\\ 0 &{} 0 &{} 0 &{} 0 &{} \tilde{q}_w &{} 0 &{} 0 &{} 0\\ 0 &{} 0 &{} 0 &{} 0 &{} -q_w &{} \tilde{\gamma }_w &{} 0 &{} 0\\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} \tilde{q}_v &{} 0\\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} -q_v &{} b_v \end{array}\right] , \end{aligned}$$

(22)

with;

$$\begin{aligned} \tilde{q}_h= & {} b_h+q_h,\qquad \tilde{\gamma }_h=b_h+\gamma _h,\qquad \tilde{q}_c=b_c+q_c,\qquad \tilde{\gamma }_c=b_c+\gamma _c,\\ \tilde{q}_w= & {} b_w+q_w,\qquad \tilde{\gamma }_w=b_w+\gamma _w,\qquad \tilde{q}_v=b_v+q_v. \end{aligned}$$

It follows from (22) that the next-generation matrix of system (19) is (23):

$$\begin{aligned} \textbf{K}=\mathcal {F}\mathcal {V}^{-1}= \left[ \begin{array}{cccccccc} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} \frac{p_h\theta _{vh} ab}{b_v\tilde{q}_h} &{} \frac{p_h\theta _{vh} ab}{b_v}\\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0\\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} \frac{p_{w}au\theta _{vw}}{b_v\tilde{q}_c} &{} \frac{p_{w}au\theta _{vw}}{b_v}\\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0\\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} p_{w}au\theta _{vw}\\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0\\ \frac{p_hac q_h}{\tilde{q}_h\tilde{\gamma }_h} &{} \frac{p_hac}{\tilde{\gamma }_h} &{} \frac{p_{c}av q_c}{\tilde{q}_c\tilde{\gamma }_c} &{} \frac{p_{c}av q_c}{\tilde{\gamma }_c} &{} \frac{p_{w}av q_w}{\tilde{q}_w\tilde{\gamma }_w} &{} \frac{p_{w}av q_w}{\tilde{\gamma }_w} &{} 0 &{} 0\\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \end{array}\right] . \end{aligned}$$

(23)

To determine the role of different hosts in maintaining transmission, we analyse the host-specific type-reproduction numbers. They are obtained from the spectral radii of the next-generation matrices with or without one or more host types (Roberts and Heesterbeek 2003). With regard to trypanosomiasis, human infection is regarded as accidental, hence we can explore the dynamics of model (22) in the absence and presence of human host. Hence, we considered the follow cases:

$({\textbf {a}})$:

No human host, then the spectral radii of model (19) becomes (24):

$$\begin{aligned} \mathcal {R}_{0cw}=\sqrt{\frac{p^2_ca^2uv\theta _{vc}}{b_v(b_c+\gamma _c)}\frac{q_v}{(b_v+q_v)}\frac{q_c}{(b_c+q_c)}+\frac{p^2_wa^2uv\theta _{vw}}{b_v(b_w+\gamma _w)}\frac{q_v}{(b_v+q_v)}\frac{q_w}{(b_w+q_w)}}. \end{aligned}$$

(24)

From (24) we have (25):

$$\begin{aligned} \mathcal {R}^2_{0cw}= & {} \Bigg (\frac{p^2_ca^2uv\theta _{vc}}{b_v(b_c+\gamma _c)}\frac{q_v}{(b_v+q_v)} \frac{q_c}{(b_c+q_c)}\Bigg ) + \Bigg (\frac{p^2_wa^2uv\theta _{vw}}{b_v(b_w+\gamma _w)}\frac{q_v}{(b_v+q_v)} \frac{q_w}{(b_w+q_w)}\Bigg )\nonumber \\= & {} \mathcal {R}_{0c} + \mathcal {R}_{0cw}, \end{aligned}$$

(25)

where $\mathcal {R}_{0c}$ and $\mathcal {R}_{0w}$ represents the expected number of cattle and wildlife infections, respectively, occurring via implicit stepping stones of infected vectors from a single infected host. In the absence of human host, $p_c+p_w=1.$ Further, $\theta _{vc}$ and $\theta _{vw}$ denotes the ratio of vectors to cattle and wildlife, respectively, $q_c/(b_c+q_c)$ is the probability of an infected cattle surviving the latent stage and becoming infectious for an average period of $1/(b_c+\gamma _c).$ Similarly, $1/(b_w+q_w)$ is the probability of an infected wildlife surviving the exposed stage and becoming infectious for an average period of $1/(b_w+\gamma _w).$ Infected vectors, have a probability $q_v/(b_v+q_v)$ of surviving the latent stage to become infectious for their entire lifespan of $1/b_v$ days.

$(\textbf{b})$:

When humans and other hosts coexist, the spectral radii of model (19) is (26):

$$\begin{aligned} \mathcal {R}_{0}= & {} \sqrt{\frac{p^2_ha^2bc\theta _{vh}}{b_v\tilde{\gamma }_h}\frac{q_v}{\tilde{q}_v}\frac{q_h}{\tilde{q}_h}+\frac{p^2_ca^2uv\theta _{vc}}{b_v\tilde{\gamma }_c}\frac{q_v}{\tilde{q}_v} \frac{q_c}{\tilde{q}_c}+\frac{p^2_wa^2uv\theta _{vw}}{b_v\tilde{\gamma }_w} \frac{q_v}{\tilde{q}_v} \frac{q_w}{\tilde{q}_w}}. \end{aligned}$$

(26)

Eqn. (26) can also be written as (27):

$$\begin{aligned} \mathcal {R}^2_{0cw}= & {} \Bigg (\frac{p^2_ha^2bc\theta _{vh}}{b_v\tilde{\gamma }_h}\frac{q_v}{\tilde{q}_v}\frac{q_h}{\tilde{q}_h}\Bigg )+\Bigg (\frac{p^2_ca^2uv\theta _{vc}}{b_v\tilde{\gamma }_c}\frac{q_v}{\tilde{q}_v}\frac{q_c}{\tilde{q}_c}\Bigg )+\Bigg (\frac{p^2_wa^2uv\theta _{vw}}{b_v\tilde{\gamma }_w}\frac{q_v}{\tilde{q}_v} \frac{q_w}{\tilde{q}_w}\Bigg )\nonumber \\= & {} \mathcal {R}_{0c}+\mathcal {R}_{0cw}, \end{aligned}$$

(27)

where, $\mathcal {R}_{0h}$ is the average number of secondary human infections generated by one infectious tsetse fly introduced in a population wholly of susceptible humans, $\theta _{vh}$ is the ratio of vectors to humans, $q_h/(b_h+q_h)$ is the probability of an infected human to survive the exposed stage and become infectious for an average period of $1/(b_h+\gamma _h).$

Sensitivity analysis

To investigate the relationship between the reproduction number and model parameters, we wil perform sensivity analysis using the partial rank correlation coefficients (PRCC) approach presented in Marino et al. (2008) and the results are in Fig. 2. Baseline values for model parameters utilized in this study were taken from published literature and presented in Table 1. Since model parameters have plausible range of values performing global sensitivity analysis enables us to understand which combinations of parameters maximize disease transmission potential.

Fig. 2

Simulation results on sensitivity analysis of $\mathcal {R}_0$ with respect to key model parameters

Table 1 Model parameters and their biological definitions

The output in Fig. 2 shows that increasing the vector mortality rate strongly reduces the disease transmission potential. One can also note that reducing the average infectious period of the hosts ($\gamma _i$, $j=h,c,w$) has the potential to reduce the disease transmission potential. Hence, it can be concluded that control measures such as vector control, screening and treatment of humans and cattle may strongly reduce disease transmission potential. However, the output also shows that an increase in vector feeding rate, disease transmission probabilities and the ratios of vector to host have significant influence on increasing the disease transmission potential. Hence, minimizing vector-host contact where possible could significantly reduce disease transmission potential.

Figure 3 shows the effects that varying two sample parameters will yield on $\mathcal {R}_0.$ If tsetse fly mortality is sufficiently high, then disease transmission potential reduces significantly. In contrast, if tsetse fly feeding rate is sufficiently high, then disease transmission potential increases.

Fig. 3

Monte Carlo simulations of 1000 sample values for two illustrative parameters (tsetse fly mortality rate and tsetse fly feeding rate) chosen via Latin Hypercube Sampling

Stability analysis of the model steady states

In epidemiological models, it is essential to investigate the global stability of the model steady states, so that one can understand the evolution of the disease over time. As highlighted earlier, in the absence of the disease, model (19) admits a disease-free equilibrium, defined by Eq. (21). In what follows, we now investigate the global stability of the disease-free equilibrium. We now claim the following result.

Theorem 1

If $\mathcal {R}_0<1,$ the disease-free equilibrium point is globally asymptotically stable, otherwise it is unstable.

Proof

Let $\mathcal {Y}(t)=(e_h,i_h,e_c,i_c,e_w,i_w,e_v,i_v)^T$. From (19), by considering infected populations only, that is (36):

$$\begin{aligned} e'_h(t)= & {} p_h\theta _{vh} abi_vs_h-(b_h+q_h)e_h,\end{aligned}$$

(28)

$$\begin{aligned} i'_h(t)= & {} q_he_h-(b_h+\gamma _h)i_h,\end{aligned}$$

(29)

$$\begin{aligned} e'_c(t)= & {} p_{c}au\theta _{vc}i_vs_c-(b_c+q_c)e_c,\end{aligned}$$

(30)

$$\begin{aligned} i'_c(t)= & {} q_ce_c-(b_c+\gamma _c)i_c,\end{aligned}$$

(31)

$$\begin{aligned} e'_w(t)= & {} p_{w}au\theta _{vw}i_vs_w-(b_w+q_w)e_w,\end{aligned}$$

(32)

$$\begin{aligned} i'_w(t)= & {} q_we_w-(b_w+\gamma _w)i_w,\end{aligned}$$

(33)

$$\begin{aligned} r'_w(t)= & {} \gamma _wi_w-(b_w+\kappa _w)r_w,\end{aligned}$$

(34)

$$\begin{aligned} e'_v(t)= & {} p_haci_hs_v+p_{c}avi_cs_v+ p_{w}avi_ws_v-(b_v+q_v)e_v,\end{aligned}$$

(35)

$$\begin{aligned} i'_v(t)= & {} q_ve_v-b_vi_v. \end{aligned}$$

(36)

System (28)-(36), can be written in the compact form (37):

$$\begin{aligned} \mathcal {Y}'(t) \le (\mathcal {F}-\mathcal {V})\mathcal {Y}(t), \end{aligned}$$

(37)

where $\mathcal {F}$ and $\mathcal {V}$ are defined in (22). It follows from the Peron-Frobenius theorem that $\mathcal {V}^{-1}\mathcal {F}$ has a positive left eigenvector $\textbf{v}$ associated with $\mathcal {R}_0,$ that is: $\textbf{v}\mathcal {V}^{-1}\mathcal {F}=\mathcal {R}_{0} {\textbf {v}}.$ Since ${\textbf {v}}V^{-1}$ is a positive vector, we propose the following Lyapunov function $\mathcal {L}= {\textbf {v}}\mathcal {V}^{-1}\mathcal {Y}.$ Taking, the fractional derivative of $\mathcal {L}$ along the solutions of system (19) one gets (38):

$$\begin{aligned} \mathcal {L}'(t)= & {} {\textbf {v}}V^{-1} \mathcal {Y}'(t)\le {\textbf {v}}V^{-1} (\mathcal {F}-\mathcal {V})\mathcal {Y}(t)\nonumber \\= & {} (\mathcal {R}_0-1){\textbf {v}}\mathcal {Y}\le 0,\qquad \qquad \text {if}~~\mathcal {R}_0 \le 1. \end{aligned}$$

(38)

It can easily be verified that the largest invariant subset of $\Gamma$ where $\mathcal {L}'(t)=0$ is the singleton $\{\mathcal {E}^0\}$. Therefore, by LaSalle’s invariance principle (LaSalle 1976), $\mathcal {E}^0$ is globally asymptotically stable in $\Gamma$ when $\mathcal {R}_0 \le 1$. This completes the proof. $\square$

If $\mathcal {R}_0>1,$ then by continuity, $\dot{\mathcal {L}}>0$ in a neighborhood of $\mathcal {E}^0$ in $\mathring{\Gamma }$. Solutions in $\mathring{\Gamma }$ sufficiently close to $\mathcal {E}^0$ move away from the DFE, implying that the DFE is unstable. Using a uniform persistence result from Freedman et al. (1994) and an argument as in the proof of Proposition 3.3 of Li et al. (1999), it can be shown that when $\mathcal {R}_0>1$, the instability of the DFE implies the uniform persistence of the model (19). Uniform persistence implies that there exists an endemic equilibrium point. To study the global stability of the endemic equilibrium point, we shall consider a case $\kappa _j=0$ as in Ndondo et al. (2016). This implies that $r_k^*=\gamma _ki^*_k/b_k.$ Based on this assertion, it implies that one can drop the variables $r_k$ from further consideration. We now turn our attention to investigate the global stability of the endemic equilibrium based on the aforementioned conditions. We begin by claiming the following result.

Theorem 2

If $\mathcal {R}_0>1,$ and $\kappa _k=0$, for $k=h,c,w$, then there exists an endemic equilibrium point which is globally asymptotically stable.

Proof

To investigate global stability of the endemic equilibrium, we consider the following Lyapunov function (39):

$$\begin{aligned} U(t)= & {} \sum _{j= h, c,w,v}\Bigg [\Bigg (s_j-s^*_j-s^*_j\ln \frac{s_j}{s^*_j}\Bigg ) + \Bigg (e_j-e^*_j-e^*_j\ln \frac{e_j}{e^*_j}\Bigg )\nonumber \\{} & {} + \frac{\tilde{q}_j}{q_j}\Bigg (i_j-i^*_j-i^*_j\ln \frac{i_j}{i^*_j}\Bigg )\Bigg ]. \end{aligned}$$

(39)

By taking the fractional derivative of U(t) along the solutions of model (19) leads to (40):

$$\begin{aligned} U'(t)\le \sum _{j= h, c,w,v}\Bigg [\Bigg (1-\frac{s^*_j}{s_j}\Bigg )D^{q}_{t_0} s_{j} +\Bigg (1-\frac{e^*_j}{e_j}\Bigg )D^{q}_{t_0} e_{j} + \frac{\tilde{q}_j}{q_j}\Bigg (1-\frac{i^*_j}{i_j}\Bigg )D^{q}_{t_0} i_j\Bigg ]. \end{aligned}$$

(40)

At endemic equilibrium we have the following identities (41):

$$\begin{aligned} \left. \begin{array}{llll} b_h &{} = &{} (1-p_c-p_{w})ab\theta _{vh}i^*_vs^*_h+b_hs^*_h,\quad \tilde{q}_h e^*_h=(1-p_c-p_{w})\theta _{vh} abi^*_vs^*_h,\\ q_he^*_h &{} = &{} \tilde{\gamma }_hi^*_h,\quad b_c=p_{c}au\theta _{vc}i^*_vs^*_c+b_cs^*_c,\quad \tilde{q}_ce^*_c= p_{c}au\theta _{vc}i^*_vs^*_c,\quad q_ce^*_c=\tilde{\gamma }_ci^*_c,\\ b_w &{} = &{} p_{w}au\theta _{vw}i^*_vs^*_w+b_ws^*_w,\quad \tilde{q}_w e^*_w= p_{w}au\theta _{vw}i^*_vs^*_w,\quad q_we^*_w=\tilde{\gamma }_wi^*_w,\\ b_v &{} = &{} (1-p_c-p_{w})aci^*_hs^{*}_{v} - p_{c}avi^*_cs^*_v + p_{w}avi^*_ws^*_v+b_vs^*_v,\\ \tilde{q}_vE^*_v &{} = &{} (1-p_c-p_{w})aci^*_hs^*_v + p_{c}avi^*_cs^*_v + p_{w}avi^*_ws^*_v,\qquad q_ve^*_v=b_vI^*_v. \end{array}\right\} \end{aligned}$$

(41)

After some algebraic manipulations one gets

$$\begin{aligned} U'(t)\le & {} \Bigg (2-\frac{s_h}{s^*_h}-\frac{s^*_h}{s_h}\Bigg ) + \Bigg (2-\frac{s_c}{s^*_c}-\frac{s^*_c}{s_c}\Bigg ) + \Bigg (2-\frac{s_v}{s^*_v}-\frac{s^*_v}{s_v}\Bigg ) + \Bigg (2-\frac{s_w}{s^*_w}-\frac{s^*_w}{s_w}\Bigg )\nonumber \\{} & {} +(1-p_c-p_w)ab\theta _{vh}i^*_vs^*_h\Bigg (3-\frac{s^*_h}{s_h}-\frac{s_he^*_hi_v}{s^*_he_hi_v^*}+\frac{i_v}{i^*_v} - \frac{i_h}{i^*_h}-\frac{e_hi^*_h}{e^*_hi_h}\Bigg )\nonumber \\{} & {} +(1-p_c-p_w)aci^*_hs^*_v\Bigg (3-\frac{s^*_v}{s_v}-\frac{s_ve^*_vi_h}{s^*_ve_vi_h^*}-\frac{e_vi^*_v}{e^*_vi_v}-\frac{i_v}{i^*_v}+\frac{i_h}{i^*_h}\Bigg )\nonumber \\{} & {} +p_cau\theta _{vc}i^*_vs^*_c\Bigg (3-\frac{s^*_c}{s_c}-\frac{s_ce^*_ci_v}{s^*_ce_ci_v^*}+\frac{i_v}{i^*_v} - \frac{i_c}{i^*_c} - \frac{e_ci^*_c}{e^*_{c}i_{c}}\Bigg )\nonumber \\{} & {} +p_cavi^*_cs^*_v \Bigg (3-\frac{s^*_v}{s_v} - \frac{s_ve^*_vi_c}{s^*_ve_vi_c^*} - \frac{e_vi^*_v}{e^*_vi_v} - \frac{i_v}{i^*_v} + \frac{i_{c}}{i^{*}_{c}}\Bigg )\nonumber \\{} & {} + p_wau\theta _{vw}i^*_vs^*_w\Bigg (3-\frac{s^*_w}{s_w}-\frac{s_we^*_wi_v}{s^*_we_wi_v^*}+\frac{i_v}{i^*_v} - \frac{i_w}{i^*_w} - \frac{e_wi^*_w}{e^*_wi_w}\Bigg )\nonumber \\{} & {} + p_wavi^*_ws^*_{v} \Bigg (3-\frac{s^*_v}{s_v}-\frac{s_ve^*_vi_w}{s^*_ve_vi_w^*}-\frac{e_vi^*_v}{e^*_vi_v} - \frac{i_v}{i^*_v}+\frac{i_w}{i^*_w}\Bigg ). \end{aligned}$$

(42)

From (42) it follows that, if $s_j=s^*_j$, $e_j=e^*_h$, $i_j=i^*_j$ and $r_k=r^*_k,$ for $j=c,h,w,v$ and $k=c,h,w,$ we have $U'(t)=0.$ Since the arithmetic mean is greater or equal to the geometric mean, that is (43):

$$\begin{aligned} \frac{s^{*}_{j}}{s_{j}}+\frac{s_{j}}{s^{*}_{j}}\ge 2\sqrt{\frac{s^{*}_{j}}{s_{j}} \cdot \frac{s_{j}}{s^{*}_{j}}}. \end{aligned}$$

(43)

In addition, we know that $\Phi (x)=1-x+\ln x$, $x>0$. Note that the function $\Phi (x)$ is non positive for $x>0$ and $\Phi (x)=0$ if and only if $x=0$. Based on the properties of $\Phi (x)$ one can easily verify that the terms in the remain brackets are less or equal to zero, therefore one can conclude that $U'(t)\le 0.$ This completes the proof. $\square$

Effects of vector’s host preference on disease dynamics

Prior studies on HAT suggest that vectors prefer animal blood meal and humans are bitten accidentally. This phenomena is known as host preference. it is worth to investigate the effects of host preference on disease dynamics. To accomplish this task we simulated model (19) using parameter values in Table 1 together with the following assumed initial conditions: $s_h=0.9995$, $e_h=0$, $r_h=0$, $i_h=s_h-e_h-r_h$, $s_c=0.9970$, $e_c=0$, $r_c=0$, $i_c=s_c-e_c-r_c$, $s_w=0.98$, $e_w=0$, $r_w=0$, $i_w=s_w-e_w-r_w.$ We simulated solutions for the variables $e_h$ and $i_h$ with the proportion of blood meal taken by tsetse fly on humans $p_h$ fixed to 0.1 while varying $p_c$ and $p_w,$ such that $p_c+p_w=1-p_h.$ The output is shown in Fig. 4. The results indicates that if the proportion of blood feeding on humans and cattle is around 10% and the rest occurs in wildlife, then the number of people infected over time will be higher compared to a scenario when proportion of blood meal taken from cattle is 30%. For $p_c=p_h=0.1$, we found $\mathcal {R}_0=3.0217$ and for $p_c=0.1,$ $p_h=0.3$, $\mathcal {R}_0=2.4612.$ Thus, the availability of wildlife as reservoirs and taking more bites than other hosts have a strong impact on the persistence and extinction of the disease.

Fig. 4

The relationship between host selection and disease dynamics in human population over time, with $p_w = 1 - p_c - p_{h}$

Control of tsetse flies and treatment of humans and cattle

The sensitivity analysis results in Fig. 2 indicates that increasing the tsetse fly mortality rate and exit rates of the hosts from the infectious stage ($\gamma _j$, $j=c,h,w$) may significantly reduce the disease transmission potential. Hence, it is imperative to investigate their influence on disease dynamics. Although the screening and treatment of wildlife is impractical, humans and cattle can be screened and treated. The mortality rate of tsetse flies can also be enhanced insecticide use. To evaluate the effects of these strategies, we redefine that rate at which humans and cattle progress from stage 1 to recovered state. Artzrouni and Gouteux (1996) proposed that the rate of exit of infected hosts to the recovery stage can be modelled as a composite of the intrinsic underlying disease progression (say, $\gamma ^h_{\text {int}})$ and the removal rate by treatment (extrinsic , say $\gamma ^h_{\text {ext}}$), that is (44):

$$\begin{aligned} \gamma _h=\gamma ^h_{\text {int}}+\gamma ^h_{\text {ext}}, \end{aligned}$$

(44)

then, the monthly percent detection is given by (45):

$$\begin{aligned} C_h=100[1-\exp (-30\gamma ^h_{\text {ext}})], \end{aligned}$$

(45)

Consequently, the exit rate from the infected class for the human host is given by (46):

$$\begin{aligned} \gamma _h=\gamma ^h_{\text {int}}-\frac{1}{30}\ln \Bigg (1-\frac{C_h}{100}\Bigg ). \end{aligned}$$

(46)

From (46) one can observe that linear screening and treatment of infected host is not linearly related to $\gamma _h.$ Rock et al. (2015a) opines that modelling recovery rate using this approach allows an extensive analysis of the reproductive number. Similarly, the modified exit rate of cattle is given by (47):

$$\begin{aligned} \gamma _c=\gamma ^c_{\text {int}}-\frac{1}{30}\ln \Bigg (1-\frac{C_c}{100}\Bigg ), \end{aligned}$$

(47)

where $C_c$ is the monthly percent detection of cattle and is given by (48):

$$\begin{aligned} C_c=100[1-\exp (-30\gamma ^c_{\text {ext}})]. \end{aligned}$$

(48)

Furthermore, we model the tsetse fly mortality rate as follows (49):

$$\begin{aligned} b_v=b_{v,\text {int}}+b_{v,\text {ext}}, \end{aligned}$$

(49)

where $b_{v,\text {int}}$ represents natural mortality of tsetse flies and and $b_{v,\text {ext}}$ accounts for additional death rate which may occur as a result of control strategies. Thus, the total mortality rate of vectors due to ‘natural’ and control measures is given by (50):

$$\begin{aligned} b_v=b_{v,\text {int}}-\ln \Bigg (1-\frac{C_v}{100}\Bigg ), \end{aligned}$$

(50)

with $C_v$ representing daily percentage of tsetse killed and is given by (51):

$$\begin{aligned} C_v=100[1-\exp (-30\mu _{v,\text {ext}})]. \end{aligned}$$

(51)

Simulation results in Fig. 5 shows the effects of screening and treatment of humans and cattle on disease dynamics over time, in the absence of insecticide use. The results indicates that increasing monthly detection and treatment of humans may significantly reduce the number of infectious humans and cattle over time and will have a slight impact of reducing the proportion of of infectious wildlife. One can also observe that detection as high as 80% may not sufficiently lead to disease eradication. Thus elimination of hat by focusing on screening and treatment of humans and cattle may not be sufficient for the disease to become extinct. However, in Fig. 6, one can observe that when insecticides are in use, killing 5% of tsetse flies everyday coupled with detection and treatment of humans and cattle at 2.5% success rate each, then the disease may become extinct in the community in a period less than 300 days.

Fig. 5

Effects varying intervention strategies on disease dynamics in human population over time, $C_v$ representing daily percentage of tsetse killed, $C_h$ denotes the monthly percent detection of infected humans and $C_c$ denotes the monthly percent detection of infected cattle

Fig. 6

Effects varying intervention strategies on disease dynamics in human population over time, $C_v$ representing daily percentage of tsetse killed, $C_h$ denotes the monthly percent detection of infected humans and $C_c$ denotes the monthly percent detection of infected cattle

The optimal control problem

Since HAT is prevalent in resource limited settings, it is crucial to desgn control strategies that are cost effective, i.e. that allow minimization of infected population at minimal costs of implementing the efforts. While disease control strategies to mitigate the spread of HAT in humans and domestic animals are available and feasible to implement the same cannot be said with regards to wildlife due to their large population and mobility. Since tsetse flies are more often found in the forests, reducing their population requires aerial spraying. Due to environmental effects associated aerial spraying this approach is no longer recommended. Furthermore, aerial spraying is also associated with higher costs compared to screening and treatment of humans and animals (Okello et al. 2021). Hence in this section we will consider intervention strategies that are applicable to humans and cattle. Thus we seek to develop an optimal control problem and identify an integrated control strategy capable of minimizing the number of infected humans and cattle over time at minimal costs.

We therefore transform model (19) into an optimal control problem by incorporating time dependent screening and treatment of infectious humans and cattle, modeled by $u_h(t)$ and $u_c(t),$ respectively. Mathematically, the objective functional to be minimized is defined as:

$$\begin{aligned} J\Bigl (\gamma _h(t), \gamma _c(t)\Bigr )=\int ^T_0 \Bigg (B_hi_h(t)+ B_c i_c(t)+C_1u_h^2(t)+C_2u_c^2(t)\Bigg )\, dt , \end{aligned}$$

(52)

where the $B_h$ and $B_c$ represent, respectively, the weight constants of infectious human and cattle populations; $C_1$ and $C_2$ are balancing coefficients transforming the integral into cost expended over a finite time period of T time units. Note that $C_1u_h^2(t)$ and $C_2u_c^2(t)$ describe the total costs associated with human and cattle treatment respectively. Note that we have proposed an objective function with quadratic controls since a quadratic structure in the control has mathematical advantages. This implies that the Hamiltonian attains its minimum over the control set at a unique point.

The weights, constant over the prescribed time frame, are a measure of the relative costs of the interventions over a finite time horizon. The optimal control problem hence becomes that we seek optimal functions, $u_h^*(t)$ and $u_c(t)$, such that:

$$\begin{aligned} J\Bigl (u_h^*(t), u_c^*(t)\Bigr )=\min J\Bigl (u_c(t), u_h(t)\Bigr ) \end{aligned}$$

(53)

subject to the state equations in system (16) with initial conditions.

The existence of optimal control follows from standard results in optimal control theory in Lenhart and Workman (2007). The necessary conditions that optimal controls must satisfy are derived using Pontryagin’s Maximum Principle (Pontryagin et al. 1962). Thus, system (19), with $u_h(t)$ and $u_c(t)$, is converted into an equivalent problem, namely the problem of minimizing the Hamiltonian H(t) given by (54):

$$\begin{aligned} H(t) = B_hi_h(t)+ B_c i_c(t)+C_1u_h^2(t)+C_2u_c^2(t)+\sum ^{15}_{i=1}\lambda _i(t)g_i, \end{aligned}$$

(54)

where $g_i$, $i=1,\cdots ,15$ denote the right hand side of system (19), which is the $i^{th}$ state variable equation, $\lambda _i$, $i=1,\cdots ,15,$ are called adjoint variables satisfying the following costate equations, with the terminal condition $\lambda _T=0$, for $i=1,\cdots ,15$. From this Hamiltonian and Pontryagin’s Maximum Principle (Pontryagin et al. 1962), one can demonstrate that given an optimal control pair $(u_h^{*}, \ u_c^{*})$ and solutions $(s_j^*,e^*_j,i^*_j,r_k^*)$, with $j=c,h,w,v$ and $k = c,h,w$ of the corresponding states system (19) there exist adjoint functions $\lambda _i(t)$, satisfying equations (55)-(69):

$$\begin{aligned} \lambda '_1(t)= & {} b_h\lambda _1-p_hab\theta _{vh}i_v(\lambda _2-\lambda _1),,\end{aligned}$$

(55)

$$\begin{aligned} \lambda '_2(t)= & {} b_h \lambda _2-q_h (\lambda _3-\lambda _2),\end{aligned}$$

(56)

$$\begin{aligned} \lambda '_3(t)= & {} -B_1+b_h\lambda _3-(\gamma _h+u_h(t))(\lambda _4-\lambda _3)-p_hacs_v(\lambda _{14}-\lambda _{13}),\end{aligned}$$

(57)

$$\begin{aligned} \lambda '_4(t)= & {} b_h\lambda _4-\kappa _h(\lambda _1-\lambda _4),\end{aligned}$$

(58)

$$\begin{aligned} \lambda '_5(t)= & {} b_c\lambda _5-p_cau\theta _{vc}i_v(\lambda _6-\lambda _5),\end{aligned}$$

(59)

$$\begin{aligned} \lambda '_6(t)= & {} b_c \lambda _6-q_c (\lambda _7-\lambda _6),\end{aligned}$$

(60)

$$\begin{aligned} \lambda '_7(t)= & {} -B_2+b_c\lambda _7-(\gamma _c+u_c(t))(\lambda _8-\lambda _7)-p_{c}avs_v(\lambda _{14}-\lambda _{13}),\end{aligned}$$

(61)

$$\begin{aligned} \lambda '_8(t)= & {} b_c\lambda _8-\kappa _c(\lambda _5-\lambda _8),\end{aligned}$$

(62)

$$\begin{aligned} \lambda '_9(t)= & {} b_w\lambda _9-p_wau\theta _{vw}i_v(\lambda _{10}-\lambda _9),\end{aligned}$$

(63)

$$\begin{aligned} \lambda '_{10}(t)= & {} b_w \lambda _{10}-q_w (\lambda _{11}-\lambda _{10}),\end{aligned}$$

(64)

$$\begin{aligned} \lambda '_{11}(t)= & {} b_w\lambda _{11}-\gamma _w(\lambda _{12}-\lambda _{11})-p_{w}avs_v(\lambda _{14}-\lambda _{13}),\end{aligned}$$

(65)

$$\begin{aligned} \lambda '_{12}(t)= & {} b_w\lambda _{12}-\kappa _w(\lambda _{9}-\lambda _{12}),\end{aligned}$$

(66)

$$\begin{aligned} \lambda '_{13}(t)&= {} b_v \lambda _{13}-p_hac(\lambda _{14}-\lambda _{13})i_{h}-p_{c}av(\lambda _{14}-\lambda _{13})i_{c}\\& \quad- p_wav(\lambda _{14}-\lambda _{13})i_w,\end{aligned}$$

(67)

$$\begin{aligned} \lambda '_{14}(t)= & {} b_v \lambda _{14}-q_v (\lambda _{15}-\lambda _{14}),\end{aligned}$$

(68)

$$\begin{aligned} \lambda '_{15}(t)= & {} b_v \lambda _{15}-p_hab\theta _{vh}(\lambda _2-\lambda _1)s_{h}-p_{c}au\theta _{vc}(\lambda _6-\lambda _5)s_c \\& \quad- p_wau\theta _{vw}(\lambda _{10}-\lambda _9)s_w,\end{aligned}$$

(69)

with transversality conditions $\lambda _j(T)=0$ for $j=1,...,15$. Furthermore, the optimal controls are characterized by the optimality conditions (70):

$$\begin{aligned} u_h (t)=\min \bigg \{1, \max \bigg ( \frac{(\lambda _3-\lambda _4)i_h}{2C_1},0\bigg )\bigg \},\quad \text {and},\quad u_c(t)=\min \bigg \{1,\max \bigg (\frac{(\lambda _7-\lambda _8)i_c}{2C_2},0\bigg )\bigg \}. \end{aligned}$$

(70)

Optimal control results

The results in the preceding section demonstrated that the constant detection and treatment of humans and cattle alone may not be sufficient to eradicate that disease when wildlife are also reservoirs. Hence it is imperative to explore how time dependent screening and treatment of humans and cattle may alter the disease dynamics over time. To investigate the impact of time-dependent screening and detection of humans and cattle, we transformed model (19) into an optimal control problem with an objective function given by Eq. (52). Further, we assume that screening and treatment of humans is given more preference compared to that of cattle , thus we set $B_h>2B_c.$ As a result, the values of $C_1$ and $C_2$ represent the relative costs of their respective controls. We further assume that the screening and treatment of humans incurs higher costs than those for the decontamination so that $C_1 > C_2$. In addition, we have set the upper bounds of our controls to 0.6 since it is highly unlikely for any any control effort to be implemented at 100%.

The simulation results in Fig. 7 show the impact of time-dependent screening and treatment of humans on disease dynamics of over time. Overall, the output shows that time dependent controls have the potential to reduce the disease prevalence to low levels over time. In particular, we can observe the population of infectious humans and cattle will be reduced to values close to zero. In addition, we can also note that this intervention strategy will also reduce exposed and infectious population in wildlife.

Fig. 7

Proportions of exposed and infectious host population over time, with $p_h=0.1,$ $p_c=0.2$, $B_c=1,$ $C_1 = 10^{-3}$ and $C_2 = 10^{-4}$

Figure 8 shows the associated control profiles with the output in Fig. 4, one can observe that all the control profiles starts at their respective maximum and remain there for a long period of time. Precisely, the control profile for control $u_c(t)$ starts at its maximum and stays remains there for the entire time horizon. The control profile for $u_h(t)$ starts at its maximum and remains there for approximately 550 days from the start, then it drops to its minimum. From these results, it can be concluded that to attain the desired results (presented in Fig. 7), all control efforts needs to be maintained at their respective maximum for the entire time horizon. However, the screening and treatment of humans may be reduced after 550 days of implementation.

Fig. 8

Control profiles for the optimal control problem, with $p_h=0.1,$ $p_c=0.2$, $B_c=1,$ $C_1 = 10^{-3}$ and $C_2 = 10^{-4}$

To investigate the impact of the costs on the implementation of control strategies, we also varied the values of the cost parameters ($C_1$ and $C_2$) and the results are in Fig. 9. Here, we set $C_1 = C_2 = 10^{-6}$, implying that the costs are low as compared to those considered earlier (Fig. 8). The results show that when the costs of implementation are low both control efforts can be maintained at their respective maxima for the entire horizon.

Fig. 9

Control profiles for the optimal control problem, with $p_h=0.1,$ $p_c=0.2$, $B_c=1,$ and $C_1 = C_2= 10^{-6}$

Concluding remarks

In this paper, we proposed and analyzed a mathematical model for trypanosomiasis that incorporates multiple hosts, humans, cattle and wildlife. We derived the reproduction number and utilized it to investigate global stability of the model steady states. We performed a sensitivity analysis, to identify model parameters that have positive and negative relationships with the reproduction number. The results showed that increasing the tsetse fly mortality rate strongly reduces the reproduction number. In addition, we also observed that increasing the exit rate of the host from stage 1 of infection can significantly reduce the reproduction number. Results also showed that: tsetse fly feeding rate and the probability of disease transmission from the vector to host and vice-versa are positively correlated with the basic reproduction number. Thus increasing these parameters will increase the reproduction number.

We extended the basic model into an optimal control problem by incorporating time-dependent screening and treatment of humans and cattle, to minimize the proportion of infected humans and cattle over a defined time horizon at minimal costs. The results showed that time dependent control can significantly reduce the disease burden in all the hosts, and to obtain such a favourably outcome, the control needs to be maintained at their maximum intensities for a greater part of the time horizon. Overall, the study showed that time-dependent screening treatment of humans and cattle may significantly reduce trypanosomiasis in communities where the disease is endemic and in proximity to game reserves. However, these control strategies may not lead to disease eradication in areas where there are wildlife reservoirs. Therefore, there may be a need to couple these intervention strategies with others such as use of insecticides. The study presented in this article is not exhaustive, one can extend the work by incorporating host and vector migration as well as temperature variations.