1 Introduction

Historically, mathematical modeling of mosquito-borne disease (VBD) has focused on either within-host pathogen dynamics, or population-level (i.e., between-host) transmission dynamics, with few models linking the two [1, 2, 4, 6, 7, 13, 17, 24, 25, 38]. Only in very recent years has the topic of immuno-epidemiological models begun to grow [1, 12, 20, 22, 23, 26, 31, 36, 37], addressing a critical gap in infectious disease modeling. An underdeveloped topic is the effect of within-human host drug pharmacokinetics (PK) and pharmacodynamics (PD) on population-level transmission dynamics and, in particular, their impact on the spread of drug resistant pathogens [27, 33]. PK is the study of drug dynamics and movement through an organism; PD is the study of the effects of the drug on an organism. Within-human host mathematical models integrated with PK/PD data have been used to evaluate the efficacy of different drugs and dosing protocols [32]. However, the role of PK/PD in population-level disease transmission has not been well studied and seems poorly understood. To date and to the best of our knowledge, only one such model linking PK/PD dynamics to within-human and between-hosts transmission interactions exists [28].

While interest in linking within-human host pathogen dynamics to between human-vector hosts transmission dynamics has increased over the last decade, there are opportunities for exploration of new approaches in creating this link naturally, and as dictated by the natural progression of the disease within a human. Considering more realistic distributions of sojourn times (compared with the usual assumption of exponentially distributed waiting times) is particularly important for these immuno-epidemiological models. For example, how long a human host has been infected is likely to influence their pathogen load and, subsequently, their ability to transmit the pathogen to another host, especially among the immunologically naive individuals in a population and for a disease like malaria. Additionally, the onset of clinical symptoms, if any, is likely to dictate when treatment commences. Likewise, how long an individual has been undergoing treatment for an infection will dictate their current blood-drug concentration and pathogen load, and therefore, their susceptibility to re-infection or superinfection with likely a different pathogen strain, along with their ability to transmit. Extending the ideas developed in [9], incorporating general distributions of sojourn times, we introduce a mosquito-borne disease framework that can be easily extended to incorporate the link between within-human host PK/PD and pathogen load dynamics to population-level disease transmission.

In our previous modeling studies [21, 34] for similar biological questions related to the effect of treatment on malaria dynamics, we used relatively simpler systems of ordinary differential equations (ODEs), simple in the context of this manuscript. Like in most ODE models for infectious diseases, we made the assumption that the disease stages (e.g., latent and infectious stages) were exponentially distributed. Although such an assumption may not be critical when the model is used to address certain biological questions, it is unrealistic and may not be appropriate for use when answering other types of questions. For example, for a PK/PD model, if the timing of treatment after infection and/or if the interaction between drug concentration and pathogen load during treatment are important factors for consideration, which typically are, then more realistic distributions for the waiting times in the different epidemiological stages must be considered in order to provide a more informative intervention strategy. In fact, it was shown in [8] that models which assume exponentially distributed waiting times can produce misleading assessment of control or intervention strategies when compared to those from corresponding models which assume gamma distributed waiting times.

In this paper, we formulate a model which allows arbitrary distributions for several disease stages. We present formulations using multiple approaches that can lead to either integral equations or partial differential equations (PDEs). We illustrate that the models derived from these two approaches are equivalent. We also show that, when these arbitrary distributions are replaced by gamma or exponential distributions, the model equations reduce to ordinary differential equations. One of the important observations is that the reduced ODE system can be very different from that when writing the ODEs without going through the derivation based on arbitrary stage distributions. Non-exponential distributions have been considered in epidemiological models (see, for example, [8,9,10, 15, 16, 18, 19, 30, 38]). However, none of these studies focus on mosquito-borne diseases. For mosquito-borne diseases like malaria in which the use of antimalarial drugs is a critical component of malaria control, particularly for reducing childhood mortality, understanding the role of arbitrary stage distributions on population-level disease dynamics is critical.

This paper is organized as follows. A model based on general survival functions for various disease stages is formulated in Sect. 2. The model consists of ordinary differential and integral equations. The reduction of this system to a system of ODEs is also presented in this section. In Sect. 3, we present an alternative formulation of the model using PDEs. In Appendix A, we show that this PDE model is equivalent to the integral equation model derived in Sect. 2. Analysis of the model is included in Sect. 4, and a discussion of the results included in Sect. 5.

2 Model Formulation with General Disease Stage Distributions

As will be illustrated later in this section, when arbitrarily distributed waiting times are considered for the infected, treated and recovered stages, the reduction of the resulting system of integral equations to ODEs can be very complicated. For demonstration purposes, the model derived in this section includes only drug-sensitive strains of a generic mosquito-borne disease.

2.1 Assumptions on Disease Stage Distributions

We are considering a Susceptible-Infected-Treated-Recovered-Susceptible (SITRS) model for the host population, where \(S=S(t)\) denotes the number of susceptible individuals at time t, who can be infected if bitten by infectious mosquitoes. \(I=I(t)\) denotes the number of infected individuals. This class contains humans who are newly infected or those that have progressed in their infection status but not yet infectious and could be either symptomatic or asymptomatic, as well as those that are infectious. Individuals in this class could die as a result of the disease, or due to natural causes, or may survive to progress to one of two possible classes. They may either be treated with a drug at some time point during their infection, or recover to join the partially immune class. \(T=T(t)\) denotes the number of individuals being treated, and \(R=R(t)\) denotes the number of recovered individuals who may become susceptible again. The definition of these variables is also listed in Table 1.

Table 1 Symbol and definition of state variables used in the model

For the mosquito population, we consider a simple case in which the total vector population size (denoted by M) remains constant, and the model is an SI type with \(S_v\) and \(I_v\) representing the numbers of susceptible and infectious mosquitoes, respectively.

Because our focus in this paper is on the effect of treatment and the interaction between drug concentration and parasite load during treatment, we consider arbitrary distributions for only the transitions associated with the I and T stages. For other processes, including mortality and immunity loss, exponential distribution will be used. Let \(P_X(x)\) denote a general survival function for a probability distribution, which represents the probability that an individuals is still in the class x time units after entering the class. We will use the following definition for the corresponding survival functions:

  • \(P_{IR}(x)\): Probability of an infected human not recovering without treatment x time units after infection

  • \(P_{TR}(x)\): Probability of an infected human not recovering with treatment x time units after receiving treatment.

  • \( P_{IT}(x)\): Probability of an infected individual not being treated x time units after becoming infected.

  • \(P_{RS}(x)\): Probability of a recovered individual not becoming susceptible x time units after recovery.

  • \(P_{D}(x)\): Probability of an individual not dying from natural cause at stage age x.

  • \(P_{ID}(x)\): Probability of a human that is infected not dying from natural or disease related causes x time units after become infected.

Particularly, we assume exponentially distributed deaths and loss of partial immunity so that

$$\begin{aligned} \displaystyle P_{RS}(x)=e^{-w x}, \quad P_{D}(x)=e^{-\mu x}, \quad P_{ID}(x)=e^{-(\mu +\delta ) x}, \end{aligned}$$
(1)

where w is the per capita rate of immunity loss, and \(\mu \) and \(\delta \) are the per capita rates of natural and disease-induced death, respectively. These probabilities define the transition rate from each category of hosts as shown in Fig. 1.

Fig. 1
figure 1

Transfer diagram for the human host

The equation for susceptible humans is given by

$$\begin{aligned} \displaystyle \frac{dS}{dt} = \varLambda -\lambda (t)S-\mu S+w R, \end{aligned}$$
(2)

where \(\lambda (t)\) denotes the rate of infection from mosquitoes to humans given by

$$\begin{aligned} \lambda =\displaystyle \beta r \frac{I_v}{M}= \displaystyle \beta \frac{I_v}{N}. \end{aligned}$$
(3)

The parameter \(\beta \) denotes the transmission rate from infectious mosquitoes to susceptible humans S, \(I_v\) is the number of infected mosquitoes, M is the total number of mosquitoes, and \(r=M/N\) is the ratio of mosquitoes to humans.

For host individuals in other epidemiological classes, we will keep track of their status changes with stage age by using integral equations that allow us to incorporate arbitrary probability distributions for the time spent in each of these epidemiological stages, especially in the infected (I) and treatment (T) stages.

All of the survival functions mentioned above have the following properties:

$$\begin{aligned} P_X(0)=1, \quad \quad P_X'(x) \le 0, \quad \lim _{x\rightarrow \infty } P_X(x)=0, \quad \int _0^{\infty } P_X(x)<\infty . \end{aligned}$$
(4)

2.2 Derivation of Equations with General Distributions

Given initial conditions at \(t=0\), we can keep track of the changes of the number I(t) at time t as follows. For ease of presentation, we use either the notation \(P_X(x)\) for transition probabilities (including death and immunity loss) or the rates \(\mu \), \(\delta \) and w. Assume that an individual was infected at a time u where \(0< u< t\). Then the individual is still in the I class at time t only if she/he has not recovered, nor been treated, nor died, for which the probabilities are each \(P_{IR}(t-u)\), \(P_{IT}(t-u)\), and \(P_{ID}(t-u)\), respectively. Thus, the total number of infected people at time t is given by

$$\begin{aligned} I(t)= \int _0^t \lambda (u) S(u) P_{IR}(t-u) P_{IT}(t-u)P_{ID}(t-u) \, du + I_0(t), \end{aligned}$$
(5)

where

$$\begin{aligned} I_0(t)=I_0 P_{IR}(t) P_{IT}(t)P_{ID}(t) \end{aligned}$$
(6)

denotes the remaining number of people in the I stage who were in I at \(t=0\) with initial number \(I_0\). The initial age distribution is ignored here as \(I_0(t) \rightarrow 0\) as \(t \rightarrow \infty \).

Fig. 2
figure 2

Depiction of transition times from infection to recovery. The top and bottom diagrams are for the cases with and without treatment, respectively

Assume that an individual who was infected at time u received treatment at time \(\tau \), where \(u< \tau < t\) (see the time line illustrated in Fig. 2), then the rate of becoming treated at \(\tau \) is \(-\frac{dP_{IT}}{d\tau }(\tau -u)\), and the probabilities that the individual has not recovered and still alive at time t are \(P_{TR}(t-\tau )\) and \(P_{D}(t-\tau )\), respectively. Then the total number of people in the T class at time t is

$$\begin{aligned} T(t)= & {} \displaystyle \int _0^t \int _0^\tau \lambda (u) S(u) P_{IR}(\tau -u) \Big [\frac{-d P_{IT}}{d \tau }(\tau -u)\Big ] P_{ID} (\tau -u) \nonumber \\&\displaystyle \times P_{TR}(t-\tau ) P_{D}(t-\tau ) du \, d\tau + T_0(t), \end{aligned}$$
(7)

where \( T_0(t)\) denotes the number of people in T class at time t due to the infected people who were present at time \(t=0\). For example,

$$\begin{aligned} T_0(t)= \int _0^t I_0 P_{IR}(u) \Big [-\frac{dP_{IT}}{du}(u)\Big ]P_{ID}(u) P_{TR}(t-u) P_D(t-u) du. \end{aligned}$$

For the recovered class R, there are two inflows, one from the I class and the other from T (see the time line diagram in Fig. 2). For the inflow from I, assume that an individual who was infected at time u recovered at time \(\sigma \), where \(u<\sigma <t\)) before receiving treatment. Then the rate of recovery is \(-\frac{dP_{IR}}{d\sigma }(\sigma -u)\) and the probability the individual is still in R at time t is \(P_{RS}(t-\sigma )\). Thus, the total number in R from I is

$$\begin{aligned} \displaystyle {\int _{0}^{t}}{\int _{0}^{\sigma }}\lambda (u) S(u) \Big [ -\frac{dP_{IR}}{d\sigma }( \sigma -u)\Big ] P_{IT}( \sigma -u) P_{ID}(\sigma -u) P_{RS}(t-\sigma ) P_{D}( t-\sigma ) du \, d\sigma . \end{aligned}$$

For the flow from T, assume that an individual who received treatment at time \(\tau \) recovers at time \(\rho \) where \(\tau< \rho < t\), then the recovery rate is \(-\frac{dP_{TR}}{d\rho }(\rho -\tau )\). The probabilities that the person has not lost immunity (i.e. has not transitioned from R to S) and is still alive at time t are \(P_{RS}(t-\rho )\) and \(P_{D}(t-\rho )\), respectively. Then the number of recovered from T at time t is

$$\begin{aligned}&\displaystyle { \int _{0}^{t}} {\int _{0}^{\rho }}{\int _{0}^{\tau }} \lambda (u) S(u) P_{IR}( \tau -u) \Big [-\frac{dP_{IT}}{d\tau }( \tau -u)\Big ]P_{ID}(\tau -u) \\&\qquad \times \displaystyle \Big [\frac{-dP_{TR}}{d\rho }(\rho -\tau )\Big ] P_{D}(\rho -\tau ) P_{RS}( t-\rho ) P_{D}(t-\rho ) du \, d\tau \, d\rho . \end{aligned}$$

Thus, the total number of people in the R class at time t is

$$\begin{aligned} R(t)= & {} \displaystyle {\int _{0}^{t}}{\int _{0}^{\sigma }}\lambda (u) S(u) \Big [ -\frac{dP_{IR}}{d\sigma }( \sigma -u)\Big ] P_{IT}( \sigma -u) P_{ID}(\sigma -u) P_{RS}(t-\sigma ) \nonumber \\&\times P_{D}( t-\sigma ) du \, d\sigma \nonumber \\&+ \displaystyle { \int _{0}^{t}} {\int _{0}^{\rho }}{\int _{0}^{\tau }} \lambda (u) S(u) P_{IR}( \tau -u) \Big [-\frac{dP_{IT}}{d\tau }( \tau -u)\Big ]P_{ID}(\tau -u) \nonumber \\&\times \displaystyle \Big [\frac{-dP_{TR}}{d\rho }(\rho -\tau )\Big ] P_{D}(\rho -\tau ) P_{RS}( t-\rho ) P_{D}(t-\rho ) du \, d\tau \, d\rho + R_0(t), \end{aligned}$$
(8)

where \( R_0(t)\) denotes the number of people in the R class at time t due to the infected people who were present at time \(t=0\). The detailed expression for \(R_0(t)\) can be written but is much longer than that for \(T_0(t)\) so we ignored it here (as it is not essential for the main results).

If we differentiate the R(t) equation, the term involving \(dP_{RS}(t)/dt\) will provide the flow going from R back to S, which is wR(t) when \(P_{RS}\) is exponential. From the equations in (5)–(8), we obtain the following system for the human population with general distributions:

$$\begin{aligned} \displaystyle \frac{dS}{dt}= & {} \varLambda -\lambda (t)S-\mu S+wR, \quad S(0)=S_0,\nonumber \\ \displaystyle I(t)= & {} \int _0^t \lambda (u) S(u) P_{IR}(t-u) P_{IT}(t-u)P_{ID}(t-u) \, du + I_0(t), \nonumber \\ \displaystyle T(t)= & {} \displaystyle \int _0^t \int _0^\tau \lambda (u) S(u) P_{IR}(\tau -u) \Big [\frac{-d P_{IT}}{d \tau }(\tau -u)\Big ] P_{ID} (\tau -u) \nonumber \\&\times P_{TR}(t-\tau ) P_{D}(t-\tau ) du \, d\tau + T_0(t), \nonumber \\ \displaystyle R(t)= & {} \displaystyle {\int _{0}^{t}}{\int _{0}^{\sigma }}\lambda (u) S(u) \Big [ -\frac{dP_{IR}}{d\sigma }( \sigma -u)\Big ] P_{IT}( \sigma -u) P_{ID}(\sigma -u) P_{RS}(t-\sigma ) \nonumber \\&\times P_{D}( t-\sigma ) du \, d\sigma \nonumber \\&+ \displaystyle { \int _{0}^{t}} {\int _{0}^{\rho }}{\int _{0}^{\tau }} \lambda (u) S(u) P_{IR}( \tau -u) \Big [-\frac{dP_{IT}}{d\tau }( \tau -u)\Big ]P_{ID}(\tau -u) \nonumber \\&\times \displaystyle \Big [\frac{-dP_{TR}}{d\rho }(\rho -\tau )\Big ] P_{D}(\rho -\tau ) P_{RS}( t-\rho ) P_{D}(t-\rho ) du \,d\tau \, d\rho + R_0(t), \end{aligned}$$
(9)

where \(\lambda (t)\) is given in (3) and \(S_0>0\) is the initial number of susceptible humans.

For the mosquito population, under the assumption that \(M=S_v+I_v\) is a constant, we can replace \(S_v\) by \(M-I_v\) and include only the following \(I_v\) equation in the model:

$$\begin{aligned} \displaystyle \frac{dI_v}{dt}= \beta _v \frac{I+\theta T}{N}(M-I_v) -\nu I_v, \quad I_v(0)=I_{v0}, \end{aligned}$$
(10)

where \(\beta _v\) is the constant transmission rate from infected hosts to mosquitoes, \(\theta \in [0,1]\) represents a constant reduction of infectivity by individuals in the T class, and \(\nu \) is the per capita birth and death rate of mosquitoes \(I_{v0}\) is the initial number of infected mosquitoes. We can assume, for simplicity, that \(I_{v0}=0\) (recall that \(I_0>0\)). The full model includes the equations in systems (9) and (10). To incorporate different pathogen loads based on bites from infectious mosquitoes at different time points, or to consider drug concentrations that depend on the time since the onset of an infection or treatment, equation (10) may be replaced by an integral equation for the vector population \(I_v(t)\); this extension is the subject of future work and will be discussed in Sect. 5.

2.3 Reduction of the Integral Equations in (9) to ODEs

In this section, we demonstrate that the integral equations can be reduced to ODEs when the arbitrary distributions are replaced by gamma distributions with the following survival functions:

$$\begin{aligned} \displaystyle P_{IR}(x)= & {} \sum _{j=1}^{J} \frac{(J a_Jx)^{(j-1)}}{(j-1)!} e^{-J a_J x}, \nonumber \\ \displaystyle P_{IT} (x)= & {} \sum _{k=1}^{K} \frac{(K a_K x)^{(k-1)}}{(k-1)!} e^{-K a_K x}, \nonumber \\ \displaystyle P_{TR} (x)= & {} \sum _{\ell =1}^{L} \frac{(L a_L x)^{(\ell -1)}}{(\ell -1)!} e^{-L a_L x}. \end{aligned}$$
(11)

In the survival functions \(P_{IR}, P_{IT}\), and \(P_{TR}\), the integers JKL are the shape parameters and \(1/a_m\) (\(m=J,K,L\)) represents the mean of the corresponding distribution. Other survival functions are exponential as given in (1).

Assume that the survival functions satisfy the following natural conditions:

$$\begin{aligned} \displaystyle P_{X(m)}(0) = 1, \quad \frac{dP_{X(m)}(x)}{dx} \le 0, \quad \lim _{x \rightarrow \infty } P_{X(m)}(x) =0, \end{aligned}$$
(12)

for \(m=J, K, L\) and correspondingly \(X(J) = IR\), \(X(K) = IT\), and \(X(L) = TR\).

The following expression will be used frequently:

$$\begin{aligned} \displaystyle -\frac{d P_{X(m)}(x)}{d x}= & {} \displaystyle m a_m \left( \sum _{i=1}^{m}\frac{(m a_m x)^{(i-1)}e^{-m a_m x}}{(i-1)!} - \sum _{i=1}^{m-1}\frac{(m a_m x)^{(i-1)} e^{-m a_m x}}{(i-1)!} \right) \nonumber \\= & {} \displaystyle \frac{m a_m (m a_m x)^{(m-1)}}{(m-1)!} e^{-m a_m x}. \end{aligned}$$

2.3.1 Derivation of Equations for Infected Sub-classes

When a gamma distribution with mean \(\alpha \) and shape parameter n is considered for a disease stage, a common understanding is that it is equivalent to dividing the stage into n sub-stages with \(\alpha /n\) being the mean residence time within each sub-stage. This does not work in the current model. Because treatment may alter the remaining sojourn in the infected status, and the recovery process is described by a different gamma distribution (\(P_{TR}\)) after receiving treatment, the number of sub-stages for the I class depends not only on \(P_{IR}\) but also on the shape parameter of \(P_{IT}\), as shown below.

Note that

$$\begin{aligned} I(t)= & {} \displaystyle \int _0^t \lambda (u) S(u) P_{IR}(t-u) P_{IT}(t-u)P_{ID}(t-u)du+ {I(0)P_{IR}(t) P_{IT}(t)P_{ID}(t)}, \\= & {} \displaystyle \int _0^t \lambda (u) S(u) \sum _{j=1}^{J} \frac{(J a_J (t-u))^{(j-1)} e^{-Ja_J(t-u)}}{(j-1)!} \\&\times \displaystyle \sum _{k=1}^{K} \frac{(K a_K (t-u))^{(k-1)}e^{-Ka_K(t-u)}}{(k-1)!} e^{- (\mu + \delta ) (t-u)} \, du \\&\displaystyle + {I(0) \sum _{j=1}^{J} \frac{(J a_J t)^{(j-1)} e^{-Ja_J t}}{(j-1)!} \cdot \sum _{k=1}^{K} \frac{(K a_K t)^{(k-1)}e^{-Ka_K t}}{(k-1)!} e^{- (\mu + \delta )t}} \\&\doteq \displaystyle \sum _{j=1}^{J} \sum _{k=1}^{K} I_{j,k}(t), \end{aligned}$$

where

$$\begin{aligned} I_{j,k}(t)= & {} \displaystyle \int _0^t \lambda (u) S(u) \frac{(J a_J (t-u))^{(j-1)} e^{-Ja_J(t-u)}}{(j-1)!} \times \frac{(K a_K (t-u))^{(k-1)}e^{-Ka_K(t-u)}}{(k-1)!} \\&\times \displaystyle e^{- (\mu + \delta ) (t-u)} \, du \\&+ \displaystyle I(0) \frac{(J a_J t)^{(j-1)} e^{-Ja_J t}}{(j-1)!} \\&\cdot \frac{(K a_K t)^{(k-1)}e^{-Ka_K t}}{(k-1)!} e^{- (\mu + \delta )t}, \quad 1\le j\le J, 1\le k \le K. \end{aligned}$$

Because the terms corresponding to \(j=1\) and \(k=1\) in the functions in (11) have much simpler and have very different properties than other terms with larger j or k, we present the derivation of these equations separately. For \(j=1\) and \(k\le 2\), or \(j \le 2\) and \(k=1\), note that

$$\begin{aligned} \displaystyle I_{1,1}= & {} \int _0^t \lambda (u) S(u) e^{-(J a_J+K a_K+\mu +\delta ) (t-u)} \, du+I(0) e^{- ({J a_J+K a_K+} \mu + \delta )t}, \\ \displaystyle I_{1,2}= & {} \int _0^t \lambda (u) S(u) K a_K (t-u) e^{-(J a_J+K a_K+\mu +\delta ) (t-u)} \, du+I(0) Ka_K t e^{- ({J a_J+} Ka_K+\mu + \delta )t}, \nonumber \\ \displaystyle I_{2,1}= & {} \int _0^t \lambda (u) S(u) J a_J (t-u) e^{-(J a_J+K a_K+\mu +\delta ) (t-u)} \, du+I(0) Ja_J t e^{- (Ja_J+{K a_K+}\mu + \delta )t}. \end{aligned}$$

Differentiating these equations we obtain:

$$\begin{aligned} \displaystyle \frac{d I_{1,1}}{dt }= & {} \lambda (t) S(t) - (J a_J + K a_K+ \mu + \delta ) I_{1,1}, \nonumber \\ \displaystyle \frac{dI_{1,2}}{dt}= & {} K a_K I_{1,1} - (J a_J + K a_K+ \mu + \delta ) I_{1,2},\nonumber \\ \displaystyle \frac{dI_{2,1}}{dt}= & {} J a_J I_{1,1} - (J a_J + K a_K+ \mu + \delta ) I_{2,1}. \end{aligned}$$
(13)

For all \(j,k \ne 1\), the equations have similar patterns and are given by

$$\begin{aligned} \frac{dI_{j,k}}{dt}= & {} J a_J I_{j-1,k} + K a_K I_{j,k-1} - (J a_J + K a_K+ \mu + \delta ) I_{j,k}. \end{aligned}$$
(14)

If we adopt the Kronecker delta notation where \(\delta _{ij} = 1\) if \(i=j\) and \(\delta _{ij}=0\) otherwise, then we can write the dI/dt equation succinctly as follows:

$$\begin{aligned} \displaystyle \frac{dI}{dt}= & {} \sum _{j=1}^{J} \sum _{k=1}^{K} \frac{dI_{j,k}}{dt}= \sum _{j=1}^{J} \sum _{k=1}^{K}\bigg [ \delta _{j,1} \delta _{k,1} \lambda (t) S(t) + (1-\delta _{j,1}) J a_J I_{j-1,k} \\&+(1-\delta _{k,1}) K a_K I_{j,k-1} - (J a_J + K a_K+ \mu + \delta ) I_{j,k}\bigg ]. \end{aligned}$$

These sub-classes of the infected stage have different biological interpretations. Figure 3 is a depiction of the transitions between sub-classes \(I_{j,k}\). Particularly, two paths can be observed, one leads to the treated stage and the other leads to recovery. Apparently, among those individuals who were infected at the same time, some may recover faster before receiving treatment while others may receive treatment sooner before recovery.

Fig. 3
figure 3

Depiction of the transitions between classes for the ODE system (18). In this example the values of J, K, and L have been arbitrarily chosen for illustrative purposes. The green dashed arrow indicates transition from an infection state to a treated state while the blue illustrate transitions from an infection state to the recovered state and the red, transition from the treated state to the recovered state. Only individuals in class \(I_{blue(J),k}\), where \(J=6\), can transition from the infected state to the recovered state. Individuals in class \(I_{j,green(K)}\) where, \(K=4\), can transition from the infected state to the treated state meanwhile individuals in class \(T_{j,red(L)}\) (for \(L=3\)) can transition from the treated state to the recovered state. Notice that individuals in class \(I_{blue(6),green(4)}\) have two possible pathways to follow, either directly to the recovered state or to the treated state, depending on the symptomatic prompt requiring treatment (Color figure online)

2.3.2 Derivation of Equations for Treated Sub-classes

For the derivation of equations for sub-classes of T, we need to divide the class into \(J \times L\) sub-classes, where L is the shape parameter of the gamma distribution described by the survival function \(P_{IT}\). Note that

$$\begin{aligned} T(t)= & {} \displaystyle \int _0^t \bigg [\int _0^\tau \lambda (u) S(u) P_{IR}(\tau -u) \Big (\frac{-d P_{IT}}{d \tau }(\tau -u)\Big ) P_{ID} (\tau -u) du \\&+ \displaystyle I(0)P_{IR}(\tau ) \Big (\frac{-d P_{IT}}{d \tau }(\tau )\Big ) P_{ID} (\tau )\bigg ] P_{TR}(t-\tau ) P_{D}(t-\tau ) d\tau \\&\doteq \displaystyle \sum _{j=1}^{J} \sum _{l=1}^{L} T_{j,l}, \end{aligned}$$

where

$$\begin{aligned} T_{j,l}(t)= & {} \displaystyle \int _0^t \bigg [\int _0^{\tau } \lambda (u) S(u) \frac{[J a_J (\tau -u)]^{(j-1)} e^{-Ja_J(\tau -u)}}{(j-1)!} \\&\times \displaystyle \frac{Ka_K[K a_K (\tau -u)]^{(K-1)}e^{-Ka_K(\tau -u)}}{(K-1)!} e^{- (\mu + \delta ) (\tau -u)} du \\&+ \displaystyle I(0) \frac{[J a_J \tau ]^{(j-1)} e^{-(Ja_J+\mu +\delta )\tau }}{(j-1)!} \cdot \frac{Ka_K[K a_K \tau ]^{(K-1)}e^{-Ka_K\tau }}{(K-1)!} \bigg ] \\&\times \displaystyle \frac{[L a_L (t-\tau )]^{(l-1)} e^{-La_L(t-\tau )}}{(l-1)!} e^{- \mu (t-\tau )} d\tau , \quad 1\le j \le J, \ 1\le l \le L. \end{aligned}$$

We first derive equations for the case of \(l=1\). For \((j,l)=(1,1)\),

$$\begin{aligned} \displaystyle T_{1,1}= & {} \int _0^t \bigg [\int _0^\tau \lambda (u) S(u) e^{-(J a_J+K a_K+\mu +\delta ) (\tau -u)} \displaystyle \frac{K a_K [K a_K (\tau -u)]^{(K -1)}}{(K -1)!} du + \displaystyle T^0_{1,1}(\tau ) \bigg ] \\&\displaystyle \times e^{-(La_L+\mu )(t-\tau )} d\tau , \end{aligned}$$

where \(T^{0}_{1,1}(\tau )\) represents the treated individuals from those infected who were presented at \(t=0\) given by

$$\begin{aligned} T^0_{1,1}(\tau )=I(0) e^{-(J a_J+Ka_K+\mu +\delta ) \tau } \frac{Ka_K[K a_K \tau ]^{(K-1)}}{(K-1)!}. \end{aligned}$$

For \((j,l)=(2,1)\),

$$\begin{aligned} \displaystyle T_{2,1}= & {} \int _0^t \bigg [\int _0^\tau \lambda (u) S(u) J a_J(\tau -u) e^{-(J a_J+K a_K+\mu +\delta ) (\tau -u)} \displaystyle \frac{K a_K [K a_K (\tau -u)]^{(K -1)}}{(K -1)!} du \\&+ T^0_{2,1}(\tau ) \bigg ] \displaystyle e^{-(La_L+\mu )(t-\tau )} d\tau , \end{aligned}$$

where

$$\begin{aligned} T^0_{2,1}(\tau )= \displaystyle I(0) J a_J \tau e^{-(J a_J+K a_K+\mu +\delta ) \tau } \frac{Ka_K[K a_K \tau ]^{(K-1)}}{(K-1)!}. \end{aligned}$$

Thus,

$$\begin{aligned} \displaystyle \frac{dT_{1,1}}{dt}= & {} \displaystyle \int _0^t \lambda (u) S(u) e^{-(J a_J +K a_K +\mu + \delta )(t-u)} \frac{K a_K [K a_K (t-u)]^{(K -1)}}{(K -1)!} du +T^0_{1,1}(t) \\&\displaystyle - \int _0^t \left[ \int _0^\tau \lambda (u) S(u) e^{-(J a_J+K a_K +\mu + \delta )(\tau -u)} \frac{K a_K [K a_K (\tau -u)]^{(K -1)}}{(K -1)!} du\right. \\&\quad \left. +T^0_{1,1}(\tau ) \right] \\&\displaystyle \times (L a_L +\mu ) e^{-(L a_L + \mu )(t-\tau )} d \tau , \\= & {} K a_K I_{1,K} - (L a_L + \mu ) T_{1,1}, \\ \displaystyle \frac{dT_{2,1}}{dt}= & {} \displaystyle \int _0^t \lambda (u) S(u) Ja_J(t-u) e^{-(J a_J +K a_K +\mu + \delta )(t-u)} \frac{K a_K [K a_K (t-u)]^{(K -1)}}{(K -1)!} du \\&+T^0_{2,1}(t) \\&\displaystyle - \int _0^t \bigg [\int _0^\tau \lambda (u) S(u) J a_J (\tau -u) e^{-(J a_J+K a_K +\mu + \delta )(\tau -u)} \\&\displaystyle \times \frac{K a_K [K a_K (\tau -u)]^{(K -1)}}{(K -1)!} du+T^0_{2,1}(\tau ) \bigg ] (L a_L +\mu ) e^{-(L a_L + \mu )(t-\tau )} d \tau , \\= & {} K a_K I_{2,K} - (L a_L +\mu ) T_{2,1}. \end{aligned}$$

Following similar patterns, we can derive equations for all \(1 \le j\le J\) and \(l=1\) and obtain:

$$\begin{aligned} \displaystyle \frac{dT_{j,1}}{dt}=K a_K I_{j,K} - (L a_L +\mu ) T_{j,1}, \quad 1\le j \le J. \end{aligned}$$

For \(l \ge 2\), from

$$\begin{aligned} T_{j,l}(t)= & {} \displaystyle \int _0^t \bigg [ \int _0^{\tau } \lambda (u) S(u) \frac{[J a_J (\tau -u)]^{(j-1)} e^{-Ja_J(\tau -u)}}{(j-1)!} \\&\times \displaystyle \frac{Ka_K[K a_K (\tau -u)]^{(K-1)}e^{-Ka_K(\tau -u)}}{(K-1)!} e^{- (\mu + \delta ) (\tau -u)} du+T^0_{j,l}(\tau ) \bigg ] \\&\times \displaystyle \frac{[L a_L (t-\tau )]^{(l-1)} }{(l-1)!} e^{- (L a_L+\mu ) (t-\tau )} d\tau , \quad 1\le j \le J, \ 2 \le l \le L, \end{aligned}$$

we know that, in \(dT_{j,l}(t)/dt\), the term corresponding to the derivative of \(T_{j,l}(t)\) with respect to the upper bound of the integral is zero. Thus,

$$\begin{aligned} \frac{dT_{j,l}}{dt} =\displaystyle L a_L T_{j,l-1} - (L a_L+\mu ) T_{j,l}, \quad 1 \le j \le J, \ 2\le l \le L. \end{aligned}$$

The transitions between sub-classes of \(I_{j,k}\) and \(T_{j,l}\) as well as with treated classes are illustrated in Fig. 3. We observe that some infected people can move more quickly through the infected and treated stages to recover, while others may take much longer time or even recover before being treated. Such detailed division of sub-classes provides a possible approach to investigate the interaction between parasite load and drug concentration during the treatment period. Model applications to such problems will be considered in future studies.

2.3.3 Derivation of the Equations for the Recovered Class

From the derivations for \(dI_{j,k}/dt\) and \(dT_{j,l}/dt\) equations we observe that the extra terms due to initial conditions do not affect the results. We will ignore the terms \( X_0(t)\) (\(X=I,T,R\)) in the following derivation. Recall that

$$\begin{aligned} I_{J,k}(t)= & {} \displaystyle \int _0^t \lambda (u) S(u) \frac{(J a_J (t-u))^{(J-1)} e^{-Ja_J(t-u)}}{(J-1)!} \cdot \frac{(K a_K (t-u))^{(k-1)}e^{-Ka_K(t-u)}}{(k-1)!} \\&\times e^{- (\mu + \delta ) (t-u)} \, du, \end{aligned}$$

and

$$\begin{aligned} \displaystyle -\frac{dP_{IR}(x)}{dx}=Ja_J\frac{(Ja_Jx)^{J-1}}{(J-1)!} e^{-Ja_Jx}. \end{aligned}$$

Thus,

$$\begin{aligned} {\displaystyle \int _{0}^{t}}\lambda (u) S(u) \Big [ -\frac{dP_{IR}}{dt}( t-u) \Big ] P_{IT}(t-u) P_{ID}(t-u) du =Ja_J\sum _{k=1}^K I_{J,k}. \end{aligned}$$
(15)

Note that

$$\begin{aligned} T_{j,l}(t)= & {} \displaystyle \int _0^t \int _0^{\tau } \lambda (u) S(u) \frac{[J a_J (\tau -u)]^{(j-1)} e^{-Ja_J(\tau -u)}}{(j-1)!} \\&\times \displaystyle \frac{Ka_K[K a_K (\tau -u)]^{(K-1)}e^{-Ka_K(\tau -u)}}{(K-1)!} e^{- (\mu + \delta ) (\tau -u)} \\&\times \displaystyle \frac{[L a_L (t-\tau )]^{(l-1)} e^{-La_L(t-\tau )}}{(l-1)!} e^{- \mu (t-\tau )} du d\tau , \quad 1\le j \le J, \ 1\le l \le L, \end{aligned}$$

and that

$$\begin{aligned} \displaystyle -\frac{dP_{IT}(x)}{dx}= & {} Ka_K\frac{(Ka_Kx)^{K-1}}{(K-1)!} e^{-Ka_Kx}, \\ \displaystyle -\frac{dP_{TR}(x)}{dx}= & {} La_L\frac{(La_Lx)^{L-1}}{(L-1)!} e^{-La_Lx}. \end{aligned}$$

Thus,

$$\begin{aligned}&\displaystyle \int _0^t \int _0^\tau \lambda (u) S(u) P_{IR}(\tau -u) \Big [\frac{-d P_{IT}}{d \tau }(\tau -u)\Big ] P_{ID} (\tau -u) \Big [\frac{-dP_{TR}}{dt}(t-\tau )\Big ] \nonumber \\&\qquad \times \displaystyle P_{D}(t-\tau ) du d\tau =La_L\sum _{j=1}^J T_{j,L}. \end{aligned}$$
(16)

From the R equation in (9) and using (15) and (16), we obtain the following differential equation:

$$\begin{aligned} \displaystyle \frac{dR}{dt}= Ja_J\sum _{k=1}^K I_{J,k}+L a_L\sum _{j=1}^J T_{j,L} -(w+\mu ) R. \end{aligned}$$
(17)

The transitions from sub-classes of I and T are shown in Fig. 3. Particularly, we observe that not all sub-classes of I or T can enter the recovered class directly. The reason that the R stage does not need to be divided into sub-classes is because of the assumption that immunity loss follows an exponential distribution, i.e., \(P_{RS}(x)=e^{-wx}\). If the immunity loss is also dependent on time since recovery, then the equations for recovered sub-classes will be much more complicated.

2.3.4 The System of ODEs

Using the ODE equations for sub-classes of I, T, and R derived in the previous sections, we obtained the reduced ODE system of the integral equations in (9) when the stage distributions are either gamma or exponential (which is the special case when the shape parameter of the gamma distribution is 1). The ODE system of the full model reads:

$$\begin{aligned} \displaystyle \frac{dS}{dt}= & {} \varLambda -\lambda (t)S-\mu S+wR, \nonumber \\ \displaystyle \frac{dI_{1,1}}{dt}= & {} \lambda (t) S- (J a_J + K a_K+ \mu + \delta ) I_{1,1}, \nonumber \\ \displaystyle \frac{dI_{1,2}}{dt}= & {} K a_K I_{1,1} - (J a_J + K a_K+ \mu + \delta ) I_{1,2}, \nonumber \\ \displaystyle \frac{dI_{2,1}}{dt}= & {} J a_J I_{1,1} - (J a_J + K a_K+ \mu + \delta ) I_{2,1}, \nonumber \\ \displaystyle \frac{dI_{j,k}}{dt}= & {} J a_J I_{j-1,k} + K a_K I_{j,k-1} - (J a_J + K a_K+ \mu + \delta ) I_{j,k},\nonumber \\&\qquad \displaystyle 2\le j \le J, \ 2\le k \le K, \nonumber \\ \displaystyle \frac{dT_{j,1}}{dt}= & {} K a_K I_{j,K} - (L a_L +\mu ) T_{j,1}, \quad 1\le j \le J, \nonumber \\ \displaystyle \frac{dT_{j,l}}{dt}= & {} \displaystyle L a_L T_{j,l-1} - (L a_L+\mu ) T_{j,l}, \quad 1 \le j \le J, \ 2\le l \le L, \nonumber \\ \displaystyle \frac{dR}{dt}= & {} Ja_J\sum _{k=1}^K I_{J,k}+L a_L\sum _{j=1}^J T_{j,L} -(w+\mu ) R, \nonumber \\ \displaystyle \frac{dI_v}{dt}= & {} \beta _v \frac{I+\theta T}{N}(M-I_v) -\nu I_v, \end{aligned}$$
(18)

with the appropriate initial conditions given as \(S(0)=S^0\), \(I_{j,k} (0)= (I_{j,k})^0\) for \(1 \le j \le J \) and \(1 \le k \le K\), \(T_{j,l} (0)= (T_{j,l})^0\) for \(1 \le j \le J\) and \(1 \le l \le L\), \(R(0)=R^0=0\) and \(I_v(0)=(I_v)^0\), with \(S^0 + \sum _{j ,k} (I_{j,k})^0 + \sum _{j, l} (T_{j,l})^0 = 1\) and \(0 \le (I_v)^0 \le 1\).

Remark

We observe from the system (18) and more clearly from the transition diagram in Fig. 3 some unusual characteristics of this system. Particularly, it is not quite straightforward to capture the intermediary transition patterns between sub-classes of \(I_{j,k}\), \(T_{j,l}\), and R, since the transition rate from \(I_{j,k}\) to \(I_{j+1,k}\) is different from the transition rate from \(I_{j,k}\) to \(I_{j,k+1}\) and from \(I_{j,k}\) to \(I_{j,l}\). However, individuals that can transition to R are only those that are in the last infected column compartments (\(I_{J,k}, k = 1,2, \cdots ,K \)) or last treated row compartments (\(T_{j,L}, j = 1,2, \cdots ,J \)). In most ODE models using gamma distribution for a disease stage, e.g., the I stage, the number n of sub-stages \(I_i\) with \(I=\sum _{i=1}^n I_i\) is the same as the shape parameter in the gamma distribution. However, in our model, the number of sub-stages is \(J\times K\) where J and K are the shape parameters for the gamma distributions for I and T stages, respectively. The more detailed separation of sub-stages is beneficial for the need of considering interaction of drug concentration and parasite load during treatment to study the effect of treatment. Results of these applications will be published elsewhere.

3 An Equivalent Model Formulation Using PDEs

To simplify notation, we assume in this section that the survival function for mortality is a negative exponential with \(\mu \) and \(\delta \) being the natural and disease-related death rates, respectively, just as was done in reducing model (9) to ODEs.

Let x(ta), y(ta) and z(ta) denote the age densities of infected, treated, and recovered individuals, respectively. Then

$$\begin{aligned} \displaystyle I(t)=\int _{0}^{\infty } x(t,a)da, \quad \displaystyle T(t)=\int _{0}^{\infty } y(t,a)da, \quad \displaystyle R(t)=\int _{0}^{\infty } z(t,a)da. \end{aligned}$$

Let \(x_0(a)=x(0,a), y_0(a)=y(0,a), z_0(a)=z(0,a)\) denote the age densities at \(t=0\), and assume that \(y_0(a)=z_0(a)=0\). Let \(B_I(t)\) denote the number of new infections in humans at time t, i.e.,

$$\begin{aligned} B_I(t)=\lambda (t)S(t). \end{aligned}$$
(19)

Then the equation for x(ta) is given by

$$\begin{aligned} \displaystyle x(t,a)=\left\{ \begin{array}{ll} B_I(t-a) P_{IR}(a)P_{IT}(a) e^{-(\mu +\delta )a}; &{} 0 \le a<t, \\ x_0(a-t)\frac{P_{IR}(a)P_{IT}(a)}{P_{IR}(a-t)P_{IT}(a-t)}e^{-(\mu +\delta )t}; &{} 0<t \le a, \\ x_0(a); &{} t=0\le a, \end{array} \right. \end{aligned}$$
(20)

where \(x_0(a)\) is the density distribution of infected at time \(t=0\). The fraction \(P_{j}(a)/P_{j}(a-t)\) (\(a>t\), \(j=IR, IT\)) represents the conditional probability that the person is still in the respective stage at stage age a given that the person was in the stage at stage age \(a-t\). From (20) we know that for \(a<t\),

$$\begin{aligned} \displaystyle \frac{\partial x}{\partial t}= B'_I(t-a)P_{IR}(a)P_{IT}(a)e^{-(\mu +\delta )a} \end{aligned}$$

and

$$\begin{aligned} \displaystyle \frac{\partial }{\partial a}\Big (\frac{x(t,a)}{P_{IR}(a)P_{IT}(a)} \Big ) =-\frac{\frac{\partial }{\partial t}x(t,a)}{P_{IR}(a)P_{IT}(a)}- (\mu +\delta )\frac{x(t,a)}{P_{IR}(a)P_{IT}(a)}. \end{aligned}$$

It follows that

$$\begin{aligned} \displaystyle \frac{\partial }{\partial t}x(t,a) +P_{IR}(a)P_{IT}(a)\displaystyle \frac{\partial }{\partial a}\Big (\frac{x(t,a)}{P_{IR}(a)P_{IT}(a)} \Big ) =-(\mu +\delta ) x(t,a). \end{aligned}$$
(21)

For \(a \ge t>0\),

$$\begin{aligned} \displaystyle \frac{\partial }{\partial t}x(t,a)= & {} \Bigg [-\frac{x'_0(a-t)}{P_{IR}(a-t)P_{IT}(a-t)} +x_0(a-t) \frac{\Big (P_{IR}(a-t)P_{IT}(a-t)\Big )'}{P^2_{IR}(a-t)P^2_{IT}(a-t)}\bigg ] \\&\times P_{IR}(a)P_{IT}(a)e^{-(\mu +\delta )t} -(\mu +\delta ) x(t,a) \end{aligned}$$

and

$$\begin{aligned} P_{IR}(a)P_{IT}(a) \frac{\partial }{\partial a}\Big (\frac{x(t,a)}{P_{IR}(a)P_{IT}(a)} \Big ) =-\frac{\partial }{\partial t}x(t,a)-(\mu +\delta ) x(t,a), \end{aligned}$$

which is identical to (21). Therefore, the PDE for x(ta) is given by (21).

Note that the instantaneous per capita rate of leaving the I stage at stage age a due to treatment is \(P'_{IT}(a)/P_{IT}(a).\) Let \(B_T(t)\) denote the flux from I to T at time t, then

$$\begin{aligned} \displaystyle B_T(t)= - \int _0^{\infty } x(t,a) \frac{P'_{IT}(a)}{P_{IT}(a)}da, \end{aligned}$$
(22)

and

$$\begin{aligned} \displaystyle y(t,a)=\left\{ \begin{array}{ll} B_T(t-a) P_{TR}(a) e^{-\mu a}; &{} 0 \le a <t, \\ 0; &{} 0\le t\le a. \end{array} \right. \end{aligned}$$
(23)

Similarly to the derivation of the PDE for x(ta) we have

$$\begin{aligned} \displaystyle \frac{\partial }{\partial t}y(t,a) +P_{TR}(a)\displaystyle \frac{\partial }{\partial a}\Big (\frac{y(t,a)}{P_{TR}(a)} \Big ) =-\mu y(t,a). \end{aligned}$$
(24)

Let \(B_{IR}(t)\) and \(B_{TR}(t)\) denote the flux from I to R and from T to R, respectively, at time t, then

$$\begin{aligned} \displaystyle B_{IR}(t)=- \int _0^{\infty } x(t,a) \frac{P'_{IR}(a)}{P_{IR}(a)}da, \quad \displaystyle B_{TR}(t)= -\int _0^{\infty } y(t,a) \frac{P'_{TR}(a)}{P_{TR}(a)}da. \end{aligned}$$
(25)

Thus,

$$\begin{aligned} \displaystyle z(t,a)=\left\{ \begin{array}{ll} \Big [B_{IR}(t-a) +B_{TR}(t-a)\Big ] P_{RS}(a) e^{-\mu a}; &{} 0 \le a <t, \\ 0; &{} 0 \le t \le a. \end{array} \right. \end{aligned}$$
(26)

and the PDE for z(ta) is

$$\begin{aligned} \displaystyle \frac{\partial }{\partial t}z(t,a) +P_{RS}(a)\displaystyle \frac{\partial }{\partial a}\Big (\frac{z(t,a)}{P_{RS}(a)} \Big ) =-\mu z(t,a). \end{aligned}$$
(27)

Therefore, we obtain the following PDE equations for the I, T and R classes:

$$\begin{aligned}&\displaystyle \frac{\partial }{\partial t}x(t,a) +P_{IR}(a)P_{IT}(a)\displaystyle \frac{\partial }{\partial a}\Big (\frac{x(t,a)}{P_{IR}(a)P_{IT}(a)} \Big ) =-(\mu +\delta ) x(t,a), \nonumber \\&\displaystyle \frac{\partial }{\partial t}y(t,a) +P_{TR}(a)\displaystyle \frac{\partial }{\partial a}\Big (\frac{y(t,a)}{P_{TR}(a)} \Big ) =-\mu y(t,a), \nonumber \\&\displaystyle \frac{\partial }{\partial t}z(t,a) +P_{RS}(a)\displaystyle \frac{\partial }{\partial a}\Big (\frac{z(t,a)}{P_{RS}(a)} \Big ) =-\mu z(t,a), \end{aligned}$$
(28)

with initial conditions

$$\begin{aligned} x(0,a)=x_0(a), \ \ y(0,a)=z(0,a)=0, \end{aligned}$$
(29)

and boundary conditions

$$\begin{aligned} x(t,0)=B_I(t), \ \ y(t,0)=B_T(t), \ \ z(t,0)=B_{IR}(t)+B_{TR}(t), \end{aligned}$$
(30)

where (see (19), (22) and (25))

$$\begin{aligned} B_I(t)= & {} \lambda (t)S(t), \nonumber \\ \displaystyle B_T(t)= & {} - \int _0^{\infty } x(t,a) \frac{P'_{IT}(a)}{P_{IT}(a)}da, \nonumber \\ \displaystyle B_{IR}(t)= & {} - \int _0^{\infty } x(t,a) \frac{P'_{IR}(a)}{P_{IR}(a)}da, \nonumber \\ \displaystyle B_{TR}(t)= & {} -\int _0^{\infty } y(t,a) \frac{P'_{TR}(a)}{P_{TR}(a)}da. \end{aligned}$$
(31)

Most of the symbols used in the PDE model (28)–(31) are listed in Table 2.

Table 2 List of symbols used in the PDE systems (28) and (32)

3.1 An Alternative Notation for the PDE Formulation

Let

$$\begin{aligned} \displaystyle {\tilde{x}}(t,a)=\frac{x(t,a)}{P_{IR}(a) P_{IT}(a)}, \quad \displaystyle {\tilde{y}}(t,a)=\frac{y(t,a)}{P_{TR}(a)}, \quad \displaystyle {\tilde{z}}(t,a)=\frac{z(t,a)}{P_{RS}(a)} \end{aligned}$$

and \(t>a\). Then

$$\begin{aligned} {\tilde{x}}(t,a)= & {} B_I(t-a) e^{-(\mu +\delta )a}, \\ {\tilde{y}}(t,a)= & {} B_T(t-a) e^{-\mu a},\\ {\tilde{z}}(t,a)= & {} \Big [B_{IR}(t-a)+B_{TR}(t-a)\Big ] e^{-\mu a}, \end{aligned}$$

and the equations in (28) can be written as

$$\begin{aligned}&\displaystyle \bigg (\frac{\partial }{\partial t}+\frac{\partial }{\partial a}\bigg ) {\tilde{x}}(t,a) =-(\mu +\delta ) {\tilde{x}}(t,a), \nonumber \\&\displaystyle \bigg (\frac{\partial }{\partial t}+\frac{\partial }{\partial a}\bigg ) {\tilde{y}}(t,a) =-\mu {\tilde{y}}(t,a), \nonumber \\&\displaystyle \bigg (\frac{\partial }{\partial t}+\frac{\partial }{\partial a}\bigg ) {\tilde{z}}(t,a) =-\mu {\tilde{z}}(t,a), \end{aligned}$$
(32)

with initial conditions

$$\begin{aligned} {\tilde{x}}(0,a)=\displaystyle \frac{x_0(a)}{P_{IR}(a) P_{IT}(a)}, \ \ {\tilde{y}}(0,a)={\tilde{z}}(0,a)=0, \end{aligned}$$

and boundary conditions

$$\begin{aligned} {\tilde{x}}(t,0)=B_I(t), \ \ {\tilde{y}}(t,0)=B_T(t), \ \ {\tilde{z}}(t,0)=B_{IR}(t)+B_{TR}(t). \end{aligned}$$

To demonstrate that the PDE system is equivalent to the system of integral equations (9), we provide in 1 a derivation of the integral equations from the PDE setting.

4 Model Analysis

The model analysis presented in this section uses the system including equations in (9) and (10), which consists of differential and integral equations. We will focus on the case when \( P_{ID}, P_{D}\) and \(P_{RS}\) are given in (1). That is, only \(P_{IR}(\tau )\), \(P_{TR}(\tau )\) and \(P_{IT}(\tau )\) are arbitrary. To simplify the analysis, we consider the case of no disease-induced host mortality (i.e., \(\delta =0\)). In this case, by differentiating the IT and R equations we can show that \(dN/dt=\varLambda -\mu N\) (a proof of this can be found in Appendix B). Thus, \(N(t) \rightarrow \varLambda /\mu \) as \(t \rightarrow \infty \). We assume that the total population has already reached the equilibrium and replace the birth rate \(\varLambda \) by the constant \(\mu N\). We include a numerical comparison of these results for the case with disease-induced death, \(\delta \ne 0\) in Sect. 4.4.

4.1 The Basic and Control Reproduction Numbers

For ease of presentation, we introduce the following notation:

$$\begin{aligned} \displaystyle {\mathcal {T}}_{IR}= & {} \int _0^{\infty } \big [-\frac{d}{du}P_{IR}(u) \big ]P_{IT}(u) e^{-\mu u}du,\nonumber \\ \displaystyle {\mathcal {T}}_{IT}= & {} \int _0^{\infty } P_{IR}(u) \big [-\frac{d}{du}P_{IT}(u) \big ]e^{-\mu u}du,\nonumber \\ \displaystyle {\mathcal {T}}_{TR}= & {} \int _0^{\infty } \big [-\frac{d}{du}P_{TR}(u) \big ]e^{-\mu u}du, \end{aligned}$$
(33)

and

$$\begin{aligned} \displaystyle {\mathcal {D}}_I= & {} \int _0^{\infty }P_{IR}(u)P_{IT}(u) P_{ID}(u)du = \int _0^{\infty }P_{IR}(u)P_{IT}(u) e^{-\mu u}du,\nonumber \\ \displaystyle {\mathcal {D}}_{T}= & {} \int _0^{\infty }P_{TR}(u) P_D(u)du = \int _0^{\infty }P_{TR}(u) e^{-\mu u}du,\nonumber \\ \displaystyle {\mathcal {D}}_{R}= & {} \int _0^{\infty }P_{RS}(u)P_{D}(u) du = \int _0^{\infty }e^{-w u}e^{-\mu u} du =\frac{1}{w+\mu }. \end{aligned}$$
(34)

These quantities have clear biological interpretations: \({\mathcal {T}}_{IR}\) is the treatment- and death-adjusted probability of recovery (i.e., probability of recovery without being treated or dead), \({\mathcal {T}}_{TR}\) is the death-adjusted probability of recovery after being treated, \({\mathcal {T}}_{IT}\) is the recovery- and death-adjusted probability of being treated, \({\mathcal {D}}_X\) represents the mean residence time in stage X with \(X=I, T, R\).

Let \({\mathcal {R}}_0\) and \({\mathcal {R}}_c\) denote the basic and control reproduction numbers, respectively. Then

$$\begin{aligned} {\mathcal {R}}_0= \Big (\beta _v {\mathcal {D}}_I\Big )\Big ( \frac{\beta r }{\nu }\Big ), \quad {\mathcal {R}}_c= \Big [\beta _v \big ({\mathcal {D}}_I+\theta {\mathcal {T}}_{IT} {\mathcal {D}}_T\big )\Big ]\Big (\frac{\beta r }{\nu }\Big ), \end{aligned}$$
(35)

In the expression of \({\mathcal {R}}_0\), the first factor (\(\beta _v {\mathcal {D}}_I\)) describes the average number of mosquitoes infected by one infected host who was not treated during the whole infectious period (\({\mathcal {D}}_I\)), and the second factor (\(\beta r /\nu \)) describes the average number of hosts infected by one infected mosquito during its entire period of infection (\(1/\nu \)). Similarly, in the expression of \({\mathcal {R}}_c\), the term \(\beta _v \theta {\mathcal {T}}_{IT} {\mathcal {D}}_T\) gives the average number of mosquitoes infected by an infected individual who received treatment (with probability \({\mathcal {T}}_{IT}\)) during the time before recovery (\({\mathcal {D}}_T\)).

4.2 Equilibria

Consider the order of variables \((S,I,T,R, I_v)\). The system always has the disease-free equilibrium

$$\begin{aligned} U_0=(N,0,0,0,0). \end{aligned}$$

Let \(U^*=(S^*,I^*,T^*,R^*, I^*_v)\) denote a positive equilibrium (i.e., \(I^*>0\)). Note that \(U^*\) can be obtained by solving the following equations:

$$\begin{aligned}&\displaystyle \mu N-\lambda ^* S^*-\mu S^*+w R^*=0, \nonumber \\&\displaystyle I^*=\lambda ^* S^* {\mathcal {D}}_I, \nonumber \\&\displaystyle T^*=\lambda ^* S^* {\mathcal {T}}_{IT} {\mathcal {D}}_T,\nonumber \\&\displaystyle R^*=\lambda ^* S^* \Big ({\mathcal {T}}_{IR}+ {\mathcal {T}}_{IT}{\mathcal {T}}_{TR} \Big ) {\mathcal {D}}_R, \end{aligned}$$
(36)

and

$$\begin{aligned} \displaystyle \beta _v \frac{I^*+\theta T^*}{N}(M-I^*_v) -\nu I^*_v=0, \end{aligned}$$
(37)

where \(\lambda ^*=\beta r I^*_v/M\). Note that

$$\begin{aligned} \displaystyle T^*=\frac{{\mathcal {T}}_{IT} {\mathcal {D}}_T}{{\mathcal {D}}_I} I^*, \quad R^*=\frac{({\mathcal {T}}_{IR}+ {\mathcal {T}}_{IT}{\mathcal {T}}_{TR}) {\mathcal {D}}_R}{{\mathcal {D}}_I} I^*. \end{aligned}$$
(38)

From (37) we have:

$$\begin{aligned} \displaystyle \frac{I^*_v}{M}=\displaystyle \frac{\beta _v C I^*/N}{\nu +\beta _v C I^*/N}, \end{aligned}$$

where \(C=1+\theta {\mathcal {T}}_{IT} {\mathcal {D}}_T/{\mathcal {D}}_I\). Using the second equation in (36) and noting that \(I^* \ne 0\) we have:

$$\begin{aligned} \displaystyle \frac{\beta \beta _v r ({\mathcal {D}}_I+ \theta {\mathcal {T}}_{IT}{\mathcal {D}}_T) S^*/N}{\nu +\beta _v C I^*/N}=1, \end{aligned}$$

from which we obtain

$$\begin{aligned} \displaystyle \frac{S^*}{N}=\frac{\nu +\beta _v C I^*/N}{\beta \beta _vr ({\mathcal {D}}_I+ \sigma {\mathcal {T}}_{IT}{\mathcal {D}}_T) }. \end{aligned}$$
(39)

Using the first equation in (36) and (39), and noticing that \(\lambda ^*S^*=I^*/{\mathcal {D}}_I\), we have

$$\begin{aligned} \displaystyle \frac{I^*}{N}=\displaystyle \frac{\mu \beta r{\mathcal {D}}_I}{\mu + \beta r\big [1-w({\mathcal {T}}_{IR}+ {\mathcal {T}}_{IT}{\mathcal {T}}_{TR}) {\mathcal {D}}_R \big ]}\Big (1-\frac{1}{{\mathcal {R}}_c}\Big ). \end{aligned}$$
(40)

Substituting (40) into (39) we also get

$$\begin{aligned} \displaystyle \frac{S^*}{N}= 1- \frac{\beta r \big [1-w({\mathcal {T}}_{IR}+ {\mathcal {T}}_{IT}{\mathcal {T}}_{TR}) {\mathcal {D}}_R \big ]}{\mu +\beta r \big [1-w({\mathcal {T}}_{IR}+ {\mathcal {T}}_{IT}{\mathcal {T}}_{TR}) {\mathcal {D}}_R \big ]} \bigg (1-\frac{1}{{\mathcal {R}}_c}\bigg ). \end{aligned}$$

In the case of no immunity loss (i.e., \(w=0\)), it is clear that \(I^*>0\) if and only if \({\mathcal {R}}_c>1\). When \(w>0\), notice that \(w{\mathcal {D}}_R=w/(w+\mu )\le 1\) and that

$$\begin{aligned} {\mathcal {T}}_{IR}+ {\mathcal {T}}_{IT}{\mathcal {T}}_{TR}\le & {} {\mathcal {T}}_{IR}+ {\mathcal {T}}_{IT} \\= & {} \displaystyle \int _0^{\infty } \big [-\frac{d}{du}P_{IR}(u) \big ]P_{IT}(u)e^{-(\mu +\delta )u}du \\&+ \int _0^{\infty } P_{IR}(u) \big [-\frac{d}{du}P_{IT}(u) \big ]e^{-(\mu +\delta ) u}du \\\le & {} -\displaystyle \int _0^{\infty } \bigg (\frac{d}{du}P_{IR}(u) P_{IT}(u) +P_{IR}(u) \frac{d}{du}P_{IT}(u)\bigg ) du\\= & {} - \displaystyle \int _0^{\infty } \frac{d}{du}\bigg (P_{IR}(u) P_{IT}(u)\bigg )du =1. \end{aligned}$$

Thus, \(1-w({\mathcal {T}}_{IR}+ {\mathcal {T}}_{IT}{\mathcal {T}}_{TR}) {\mathcal {D}}_R\ge 0. \) It follows that \(I^*>0\) when \({\mathcal {R}}_c>1\), in which case, \(S^*, T^*\), \(R^*\) and \(I_v^*\) are also positive. We can also show that solutions of (36) satisfy \(S^*+I^*+T^*+R^*=N\) (see Appendix C). It follows that \(0<S^*, I^*, T^*, R^* <N\). We have proved that, \(U^*\) exists and is unique if and only if \({\mathcal {R}}_c>1\).

4.3 Stability

For ease of presentation, introduce the following notation:

$$\begin{aligned} a_1(\tau )= & {} \displaystyle P_{IR}(\tau ) P_{IT}(\tau ) P_{ID}(\tau ), \nonumber \\ a_2(\tau )= & {} \displaystyle \int _0^{\tau } P_{IR}(\tau -s)\big [-\frac{d}{d{\tau }} P_{IT}(\tau -s)\big ] P_{ID}(\tau -s) P_{TR}(s)P_{D}(s)ds. \end{aligned}$$
(41)

Then,

$$\begin{aligned} {\mathcal {R}}_c=\frac{\beta r \beta _v }{\nu } \int _0^{\infty }\big [a_1(u)+\theta a_2(u)\big ]du. \end{aligned}$$

Consider the I and T equations in system (9). When t is sufficiently large, the terms \({\tilde{I}}_0(t)\) and \({\tilde{T}}_0(t)\) are close to zero and can be ignored when considering solution behaviors as t tends to infinity. Then, we can rewrite the I equation in system (9) as

$$\begin{aligned} I(t)=\int _0^t \lambda (u)S(u) a_1(t-u) du. \end{aligned}$$
(42)

By changing the order of integration we can rewrite the T equation as

$$\begin{aligned} T(t)=\int _0^t \lambda (u)S(u)a_2(t-u) du \end{aligned}$$
(43)

Because \(I_v(t)\ge 0\) for \(t\ge 0\) and all parameters are non-negative, from (10) we have

$$\begin{aligned} \displaystyle \frac{dI_v}{dt} \le \beta _v \frac{I+\theta T}{N}M -\nu I_v =\beta _v r (I+\theta T)-\nu I_v. \end{aligned}$$

Assume that \(I_v(0)=0\). Then

$$\begin{aligned} I_v(t)\displaystyle\le & {} e^{-\nu t} \int _0^t e^{\nu s} \beta _v r \big [I(s)+\theta T(s)\big ]ds \\ \displaystyle\le & {} \beta _v r e^{-\nu t} \int _0^t e^{\nu s} \int _0^s \lambda (u)S(u) \big [a_1(s-u)+\theta a_2(s-u)\big ] du ds. \end{aligned}$$

Let \(A(\tau )=a_1(\tau )+\theta a_2(\tau )\). Multiplying both sides of the the above inequality by \(\beta r/M\) and noticing that \(S\le N\), we have

$$\begin{aligned} \displaystyle \lambda (t) \le \beta \beta _v r e^{-\nu t}\int _0^t e^{\nu s} \int _0^s \lambda (u)A(s-u)du ds. \end{aligned}$$

Integrating by parts:

$$\begin{aligned} \lambda (t) \le \displaystyle \frac{\beta _v \beta r}{\nu }\int _0^t \lambda (u) A(t-u) du. \end{aligned}$$
(44)

Let \(t_n\) be a sequence with the property that

$$\begin{aligned} t_{n+1}-t_n \rightarrow \infty , \quad \lambda (t_n) \rightarrow \lambda ^{\infty }, \quad \hbox {as} \ \ n\rightarrow \infty . \end{aligned}$$

Then,

$$\begin{aligned} \lambda (t_{n+1})\le \frac{\beta _v \beta r}{\nu } \int _0^{t_n} \lambda (u)A(t_{n+1}-u)du + \frac{\beta _v \beta r}{\nu } \int _{t_n}^{t_{n+1}} \lambda (u)A(t_{n+1}-u)du. \end{aligned}$$

Note that \(\lambda (u)=\beta r I_v(u)/M \le \beta r\). Then

$$\begin{aligned} \int _0^{t_n} \lambda (u)A(t_{n+1}-u)du \le \beta r \int _{t_n}^{t_{n+1}}A(\tau )d\tau \rightarrow 0 \quad \hbox {as} \ n \rightarrow \infty , \end{aligned}$$

and

$$\begin{aligned} \int _{t_n}^{t_{n+1}} \lambda (u)A(t_{n+1}-u)du \le \sup _{t\ge t_n}\lambda (t) \int _{0}^{\infty }A(\tau )d\tau \rightarrow \lambda ^{\infty }\int _{0}^{\infty }A(\tau )d\tau \quad \hbox {as} \ n \rightarrow \infty . \end{aligned}$$

It follows that

$$\begin{aligned} \lambda ^{\infty } \le \displaystyle \frac{ \beta _v\beta r}{\nu } \lambda ^{\infty } \int _{0}^{\infty }A(\tau )d\tau =\lambda ^{\infty } {\mathcal {R}}_c. \end{aligned}$$

Therefore, if \({\mathcal {R}}_c<1\) then \(\lambda ^{\infty }=0\), i.e.,

$$\begin{aligned} \lim _{t\rightarrow \infty } \lambda (t)=0, \end{aligned}$$
(45)

and equivalently, \(\lim _{t\rightarrow \infty } I_v(t)=0.\) From this and the fact that \(I_v(t) \ge 0\), we know that \(I_v'(t) \rightarrow 0\) as \(t \rightarrow 0\). From the \(I_v\) equation (10) we have

$$\begin{aligned} \lim _{t \rightarrow \infty }\beta _v \frac{I(t)+\theta T(t)}{N}\big (M-I_v(t) \big )-\nu I_v(t) =0, \end{aligned}$$

from which we obtain

$$\begin{aligned} \lim _{t \rightarrow \infty }\big [ I(t)+\theta T(t)\big ] =0. \end{aligned}$$

It follows that

$$\begin{aligned} \quad \lim _{t\rightarrow \infty } I(t) =\lim _{t\rightarrow \infty } T(t) =\lim _{t\rightarrow \infty } R(t)=0, \quad \lim _{t\rightarrow \infty } S(t) = \frac{\varLambda }{\mu }. \end{aligned}$$

This shows that the disease-free equilibrium is globally asymptotically stable if \({\mathcal {R}}_c<1\). Using a similar approach as in [8] we can show that \(U^*\) is locally asymptotically stable if \({\mathcal {R}}_c>1\).

4.4 Numerical Comparison with Nonzero Disease-Induced Death

To simplify our model analysis in the previous subsections, we assumed a zero disease-induced death, \(\delta =0\). In this section we present numerical comparisons between the control reproduction number, \(\mathcal {R}_c\) calculated in Eq. 35 to the value calculated for the ODE system given by Eq. 18 using the next generation method. The parameter values, shown in Table 3, used in the numerical comparison are for malaria in a high transmission region. The shape parameters J, K, and L are chosen to match the illustration of stages in Fig. 3. To model a particular situation these shape parameters would be determined by a fitting a gamma distribution to data. In Table 3, we see that the malaria induced death rate for adults is small compared to the natural death rate, \(\delta \approx \mu /40\). Figure 4 illustrates that the difference between the control reproduction numbers as a function of \(\beta \beta _v\) is very small. Note that if \(\beta \beta _v = 1\) every mosquito bite would transmit the malaria parasite from an infected human to a mosquito. In high transmission regions with endemic malaria, \(\beta \beta _v < 0.2\).

Table 3 Parameter values for numerical comparison: malaria high transmission regions, adult host
Fig. 4
figure 4

Difference between \(\mathcal {R}_c\) from PDE Analysis with \(\delta =0 \) in Eq. 35 and \(\mathcal {R}_c\) from next generation matrix from ODE system Eq. 18 with \(\delta \ne 0\), using parameters in Table 3

In Fig. 5 we use the host integral equation system described by Eq. 9 and mosquito differential equation in Eq. 10 to plot the infection status of the hosts and vectors versus time for \(\mathcal {R}_c>1\) and \(\mathcal {R}_c<1\). The numerical algorithm used to solve Eqs. 9, 10 and to produce Fig. 5 is described in [14]. The parameter values chosen are from Table 3, including disease-induced death. The values for all the parameters, including \(\beta _V\) and \(\beta \) at endemic equilibrium (EE) are for a high transmission region with endemic malaria. The values for \(\beta _V\) and \(\beta \) at disease-free equilibrium (DFE) are artificial. We see that as our analysis shows for \(\mathcal {R}_c>1\), shown in Fig. 5a, the system reaches an endemic equilibrium and for \(\mathcal {R}_c<1\) shown in Fig. 5b, the system reaches the disease-free equilibrium after 200 days. In Fig. 5 both the EE and DFE situations include disease-induced death, the equilibrium values will change whether or not disease-induced death is included. In both situations, \(\delta = 0\) or \(\delta \ne 0\), for \(\mathcal {R}_c>1\) the host and vector populations move to an endemic equilibrium and for \(\mathcal {R}_c<1\) the host and vector populations move to the disease-free equilibrium. We do not include the figures for \(\delta =0\), as they are not markedly different than Fig. 5.

Fig. 5
figure 5

Infection status of hosts and vectors versus time in days for the integral and differential equation system in Eqs. 9 and 10 showing the approach to equilibrium. a \(\mathcal {R}_c>1\) with endemic equilibrium values from Table 3. b \(\mathcal {R}_c<1\) with DFE values from Table 3

5 Discussion

The main contribution of this modeling study is the incorporation of general distributions of disease stages into models for host-vector diseases such as malaria to investigate effects of drug treatment. Because of the realistic description of the stage duration from time since infection and from treatment, such models are capable of describing the dynamic changes in drug concentration and parasite load during treatment, which is impossible if stage distributions are assumed to be exponential as commonly done in most ODE models for mosquito-borne diseases.

We present two approaches for the derivation of the model with general distributions: one uses integral equations and the other uses partial differential equations. We showed that the models derived from these two approaches lead to the same system. Depending on goals of model analysis and available mathematical tools, one form might be easier to use than the other. In this paper, the model analysis is based on the formulation of integral equations. This formulation is also used to obtain the corresponding system of ODEs when the distributions are gamma or exponential.

Another important finding in this study is that it highlights the shortcoming of writing models of ODE equations (assuming exponentially distributed stage durations) without using the model with arbitrary sojourn distributions (see Sect. 2.3). For example, when a gamma distribution is used for a disease stage, the resulting ODE equations usually includes n sub-stages, where n is the shape parameter in the gamma distribution. However, this method does not work for the model studied in this paper. As illustrated in Fig. 3, the I class in the integral equation (9) actually needs to be divided into \(J \times K\) sub-classes (denoted by \(I_{j,k}\), \(1 \le j \le J, \ 1\le k \le K\)), where J and K are the shape parameters of the gamma distributions for the I and T classes, the two classes representing the two pathways from the infected human class. This, in fact, highlights the heterogeneous nature of the pathways infected humans take and how a pathogen manifests within each individual, even when the humans are infected at the same time. This is even more true and important for a disease like malaria where because of the role of acquired (adaptive) immunity (which is positively correlated to how often an individual may have been exposed to infectious mosquito bites), the parasite manifestation within each individual will differ. Some individuals may be asymptomatic with a lower severity of the diseases and likely to recover without treatment, while others may exhibit symptoms early and may seek treatment early. The progression then of a malaria-infected individual to either the recovered class without treatment or to the treated class, based on the severity of symptoms will differ among individuals. Some infected individuals may recover faster before seeking treatment, if they do seek treatment (likely as a result of having a better developed adaptive immune response), while others may receive treatment sooner before recovery due to early symptoms or the severity of the symptoms and disease. The latter is likely to occur among the immunologically naive populations. Thus, the path to recovery without treatment for an individual will require a minimum of J steps traversed horizontally on I as in Fig. 3, or \(K+L\) traversed vertically on I and T or \(K+J\) traversed vertically and then horizontally on I. Our modeling approach sets the stage for us to model these different cases as the waiting times for an infected human before recovery or before treatment are no longer assumed to be exponentially distributed.

Since the goal of this manuscript was to demonstrate the modeling approach, the model studied includes only drug-sensitive strains. One of the serious public health concerns is the development of drug-resistant strains due to inappropriate treatment schedules. We have considered ODE models with both sensitive- and resistant strains in previous studies (see [21, 34]). We did not include resistance strains in the current model because, as seen in Sect. 2.3, the derivation for the reduction of equations from the integral equations to ODEs is already very complicated. We will extend the general model in this paper to include a drug-resistant strain in future studies. Furthermore, to extend the general model in this paper to capture the interaction between vectors and humans for a mosquito-borne disease, and address the feedback between within-human host PK/PD and between-hosts transmission, we must replace the ODE vector equation \(dI_v(t)/dt\) with an integral equation for \(I_v(t)\) of the following form:

$$\begin{aligned} I_v(t)&= I_v(0)e^{-\nu t} + \int _0^t \int _0^\eta \beta _v(\eta -u) (M - I_v(\eta )) \\&\qquad \frac{\lambda (u)S(u) P_{IR}(\eta -u) P_{IT}(\eta -u) P_{ID}(\eta -u)}{N(u)}e^{-\nu (t-\eta )} du d\eta \\&\quad + \int _0^t \int _0^\tau \int _0^\eta \beta _v(\tau -u) \theta (\tau -\eta )(M- I_v(\tau )) \\&\quad \times \frac{\lambda (u)S(u) P_{IR}(\eta -u) \left[ -\frac{dP_{IT}(\eta -u)}{d\eta } \right] P_{ID}(\eta -u)P_{TR}(\tau -\eta ) P_{TD}(\tau -\eta )}{N(u)} \\&\quad \times e^{-\nu (t-\eta )} du d\eta d\tau , \end{aligned}$$

where \(\eta \) denotes the time at which a susceptible mosquito takes a blood meal from an infected human, and u the time the human was infected by an infectious vector. This formulation incorporates the transmission potential as a function that depends on how long the infected human has been in stage I or T at the time of transmission. This can be linked to the pathogen load and drug concentration in the infected human. For example, a newly infected individual will have a \(\beta _v\) value that is close to zero, but increases with time, with a sharp rise, as the time the transmissible forms of the pathogens to appear approaches (thus capturing the incubation period of the infection). We remark that in the above integral equation for \(I_v(t)\), the infection of a mosquito by a human host in the treatment class T depends on the transmission rate \(\beta _v (\tau -u)\) (which depends on the age-since-infection \((\tau -u)\) and the reduction factor \(\theta (\tau -\eta )\), which depends on the age-since-treatment \((\tau -\eta )\)). Because the age-since-infection can be associated to parasite load whereas the age-since-treatment is associated with drug concentration, the product \(\beta _v (\tau -u) \theta (\tau -\eta )\) allows for the possibility to describe and model the interaction between parasite load and drug concentration.

The differential equation that corresponds to the above integral equation is

$$\begin{aligned} \frac{dI_v(t)}{dt}&= \int _0^t \beta _v(t-u) (M - I_v(t)) \frac{\lambda (t)S(t) P_{IR}(t-u) P_{IT}(t-u) P_{ID}(t-u) }{N(t)} du\\&\quad + \int _0^t \int _0^\eta \beta _v(t-u) \theta (t-\eta )(M - I_v(t)) \\&\quad \times \frac{\lambda (u)S(u) P_{IR}(\eta -u) \left[ -\frac{dP_{IT}(\eta -u)}{d\eta } \right] P_{ID}(\eta -u)P_{TR}(t-\eta ) P_{TD}(t-\eta )}{N(u)} du d\eta \\&\quad - \nu I_v(t), \end{aligned}$$

Notice that when we assume constant transmission rate \(\beta _v\), and constant reduction parameter \(\theta \in [0,1]\) , we can show using the forms of I(t) and T(t) defined in (9), that the above differential equation reverts to the form

$$\begin{aligned} \displaystyle \frac{dI_v(t)}{dt}= \beta _v \frac{I(t)+\theta T(t)}{N}(M-I_v(t)) -\nu I_v(t), \end{aligned}$$
(46)

as was obtained in (10).