A Mosquito-Borne Disease Model with Non-exponentially Distributed Infection and Treatment Stages

Most epidemiological models for mosquito-borne disease assume exponentially distributed residence times in disease stages in order to simplify the model formulation and analysis. However, these models do not allow for accurate description of the interaction between drug concentration and pathogen load within hosts. To improve this, we formulate a model by considering arbitrarily distributed sojourn for various disease stages. The model formulation is presented using two proven equivalent approaches: integral equations and partial differential equations. Analysis of the model includes the existence of equilibrium solutions and stability, which are shown to be dependent on whether the control reproduction number Rc\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {R}}_c$$\end{document} is less or greater than 1. It is also shown that, when the general distributions are replaced by gamma distributions, the system of integral equations can be reduced to a system of ODEs, which has some non-trivial characteristics which are only captured by non-exponential distributions for disease stages.

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.

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 mosquitoborne disease.

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.
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 I R (x): Probability of an infected human not recovering without treatment x time units after infection  Particularly, we assume exponentially distributed deaths and loss of partial immunity so that where w is the per capita rate of immunity loss, and μ and δ 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.
The equation for susceptible humans is given by where λ(t) denotes the rate of infection from mosquitoes to humans given by The parameter β 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:

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 μ, δ 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 I R (t − u), P I T (t − u), and P I D (t − u), respectively. Thus, the total number of infected people at time t is given by where I 0 (t) = I 0 P I R (t)P I T (t)P I D (t) (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) → 0 as t → ∞. Assume that an individual who was infected at time u received treatment at time τ , where u < τ < t (see the time line illustrated in Fig. 2), then the rate of becoming treated at τ is − d P I T dτ (τ − u), and the probabilities that the individual has not recovered and still alive at time t are P T R (t − τ ) and P D (t − τ ), respectively. Then the total number of people in the T class at time t is 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, 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 σ , where u < σ < t) before receiving treatment. Then the rate of recovery is − d P I R dσ (σ − u) and the probability the individual is still in R at time t is P RS (t − σ ). Thus, the total number in R from I is For the flow from T , assume that an individual who received treatment at time τ recovers at time ρ where τ < ρ < t, then the recovery rate is − d P T R dρ (ρ − τ ). 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 − ρ) and P D (t − ρ), respectively. Then the number of recovered from T at time Thus, the total number of people in the R class at time t is 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 d P RS (t)/dt will provide the flow going from R back to S, which is w R(t) when P RS is exponential. From the equations in (5)-(8), we obtain the following system for the human population with general distributions: where λ(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: where β v is the constant transmission rate from infected hosts to mosquitoes, θ ∈ [0, 1] represents a constant reduction of infectivity by individuals in the T class, and ν 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.

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: In the survival functions P I R , P I T , and P T R , the integers J , K , L 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: for m = J , K , L and correspondingly X (J ) = I R, X (K ) = I T , and X (L) = T R.
The following expression will be used frequently:

Derivation of Equations for Infected Sub-classes
When a gamma distribution with mean α 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 α/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 T R ) after receiving treatment, the number of sub-stages for the I class depends not only on P I R but also on the shape parameter of P I T , as shown below. Note that where 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 ≤ 2, or j ≤ 2 and k = 1, note that Differentiating these equations we obtain: For all j, k = 1, the equations have similar patterns and are given by If we adopt the Kronecker delta notation where δ i j = 1 if i = j and δ i j = 0 otherwise, then we can write the d I /dt equation succinctly as follows: 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.

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 × L sub-classes, where L is the shape parameter of the gamma distribution described by the survival function P I T . Note that We first derive equations for the case of l = 1. For ( j, l) = (1, 1), where T 0 1,1 (τ ) represents the treated individuals from those infected who were presented at t = 0 given by For ( j, l) = (2, 1), Thus, Following similar patterns, we can derive equations for all 1 ≤ j ≤ J and l = 1 and obtain: For l ≥ 2, from 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, 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.

Derivation of the Equations for the Recovered Class
From the derivations for d I 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 Note that and that From the R equation in (9) and using (15) and (16), we obtain the following differential equation: 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.

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: with the appropriate initial conditions given as S(0) = S 0 , I j,k (0) = (I j,k ) 0 for 1 ≤ j ≤ J and 1 ≤ k ≤ K , T j,l (0) = (T j,l ) 0 for 1 ≤ j ≤ J and 1 ≤ l ≤ L, R(0) = R 0 = 0 and I v (0) = (I v ) 0 , with S 0 + j,k (I j,k ) 0 + j,l (T j,l ) 0 = 1 and 0 ≤ (I v ) 0 ≤ 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, · · · , K ) or last treated row compartments (T j,L , j = 1, 2, · · · , 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 = n i=1 I i is the same as the shape parameter in the gamma distribution. However, in our model, the number of sub-stages is J × 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.

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 μ and δ being the natural and disease-related death rates, respectively, just as was done in reducing model (9) to ODEs. Let x(t, a), y(t, a) and z(t, a) denote the age densities of infected, treated, and recovered individuals, respectively. Then Then the equation for x(t, a) is given by 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 = I R, I T ) 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,

It follows that ∂ ∂t x(t, a) + P I R (a)P I T (a) ∂ ∂a x(t, a) P I R (a)P I T (a)
= −(μ + δ)x(t, a).
For a ≥ t > 0, and

P I R (a)P I T (a)
∂ ∂a which is identical to (21). Therefore, the PDE for x(t, a) 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 I T (a)/P I T (a). Let B T (t) denote the flux from I to T at time t, then and Similarly to the derivation of the PDE for x(t, a) we have Let B I R (t) and B T R (t) denote the flux from I to R and from T to R, respectively, at time t, then Thus, and the PDE for z(t, a) is Therefore, we obtain the following PDE equations for the I , T and R classes:   (28) and (32) Symbol Description and boundary conditions where (see (19), (22) and (25) Most of the symbols used in the PDE model (28)-(31) are listed in Table 2.

An Alternative Notation for the PDE Formulation
Letx

(t, a) = x(t, a) P I R (a)P I T (a) ,ỹ(t, a) = y(t, a) P T R (a) ,z(t, a) = z(t, a) P RS (a)
and t > a. Thenx  (28) can be written as with initial conditions ,ỹ(0, a) =z(0, a) = 0, and boundary conditions 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.

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 I D , P D and P RS are given in (1). That is, only P I R (τ ), P T R (τ ) and P I T (τ ) are arbitrary. To simplify the analysis, we consider the case of no disease-induced host mortality (i.e., δ = 0). In this case, by differentiating the I , T and R equations we can show that d N /dt = Λ − μN (a proof of this can be found in Appendix B). Thus, N (t) → Λ/μ as t → ∞. We assume that the total population has already reached the equilibrium and replace the birth rate Λ by the constant μN . We include a numerical comparison of these results for the case with disease-induced death, δ = 0 in Sect. 4.4.

The Basic and Control Reproduction Numbers
For ease of presentation, we introduce the following notation: These quantities have clear biological interpretations: T I R is the treatment-and death-adjusted probability of recovery (i.e., probability of recovery without being treated or dead), T T R is the death-adjusted probability of recovery after being treated, T I T is the recovery-and deathadjusted probability of being treated, D X represents the mean residence time in stage X with X = I , T , R. Let R 0 and R c denote the basic and control reproduction numbers, respectively. Then In the expression of R 0 , the first factor (β v D I ) describes the average number of mosquitoes infected by one infected host who was not treated during the whole infectious period (D I ), and the second factor (βr /ν) describes the average number of hosts infected by one infected mosquito during its entire period of infection (1/ν). Similarly, in the expression of R c , the term β v θ T I T D T gives the average number of mosquitoes infected by an infected individual who received treatment (with probability T I T ) during the time before recovery (D T ).

Equilibria
Consider the order of variables (S, I , T , R, I v ). The system always has the disease-free equilibrium 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: μN − λ * S * − μS * + w R * = 0, and where λ * = βr I * v /M. Note that From (37) we have: where C = 1 + θ T I T D T /D I . Using the second equation in (36) and noting that I * = 0 we have: from which we obtain Using the first equation in (36) and (39), and noticing that λ * S * = I * /D I , we have Substituting (40) into (39) we also get In the case of no immunity loss (i.e., w = 0), it is clear that I * > 0 if and only if R c > 1.
When w > 0, notice that wD R = w/(w + μ) ≤ 1 and that Thus, 1 − w(T I R + T I T T T R )D R ≥ 0. It follows that I * > 0 when 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 R c > 1.

Stability
For ease of presentation, introduce the following notation: Then, Consider the I and T equations in system (9). When t is sufficiently large, the termsĨ 0 (t) andT 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 By changing the order of integration we can rewrite the T equation as Because I v (t) ≥ 0 for t ≥ 0 and all parameters are non-negative, from (10) we have Assume that I v (0) = 0. Then Let A(τ ) = a 1 (τ ) + θa 2 (τ ). Multiplying both sides of the the above inequality by βr /M and noticing that S ≤ N , we have Integrating by parts: Let t n be a sequence with the property that Then, Note that λ(u) = βr I v (u)/M ≤ βr . Then It follows that Therefore, if R c < 1 then λ ∞ = 0, i.e., and equivalently, lim t→∞ I v (t) = 0. From this and the fact that I v (t) ≥ 0, we know that I v (t) → 0 as t → 0. From the I v equation (10) we have Fig. 4 Difference between R c from PDE Analysis with δ = 0 in Eq. 35 and R c from next generation matrix from ODE system Eq. 18 with δ = 0, using parameters in Table 3 from which we obtain It follows that This shows that the disease-free equilibrium is globally asymptotically stable if R c < 1.
Using a similar approach as in [8] we can show that U * is locally asymptotically stable if R c > 1.

Numerical Comparison with Nonzero Disease-Induced Death
To simplify our model analysis in the previous subsections, we assumed a zero diseaseinduced death, δ = 0. In this section we present numerical comparisons between the control reproduction number, 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, δ ≈ μ/40. Figure 4 illustrates that the difference between the control reproduction numbers as a function of ββ v is very small. Note that if ββ v = 1 every mosquito bite would transmit the malaria parasite from an infected human to a mosquito. In high transmission regions with endemic malaria, ββ v < 0.2.
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 R c > 1 and 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 β V and β at endemic equilibrium (EE) are for a high transmission region with endemic malaria. The values for β V and β at disease-free equilibrium (DFE) are artificial. We see that as our analysis shows for R c > 1, shown in Fig. 5a, the system reaches an endemic equilibrium and for 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, δ = 0 or δ = 0, for R c > 1 the host and vector populations move to an endemic equilibrium and for R c < 1  Table 3 the host and vector populations move to the disease-free equilibrium. We do not include the figures for δ = 0, as they are not markedly different than Fig. 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 substages, 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 × K sub-classes (denoted by I j,k , 1 ≤ j ≤ J , 1 ≤ k ≤ 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 d I v (t)/dt with an integral equation for I v (t) of the following form: where η 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 β 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 β v (τ − u) (which depends on the age-since-infection (τ − u) and the reduction factor θ(τ −η), which depends on the age-since-treatment (τ −η)). Because the age-since-infection can be associated to parasite load whereas the age-since-treatment is associated with drug concentration, the product β v (τ − u)θ (τ − η) 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 Notice that when we assume constant transmission rate β v , and constant reduction parameter θ ∈ [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 as was obtained in (10).
Note also from (31) that Thus, which is identical to the I equation in (9) except the term for initially presented infected people. The term involving x 0 (a − t) above includes a more detailed description as it takes into consideration of the age density at time t = 0 and the time already spent in the infected class at time t > 0 (the denominator), whereas the term I 0 (t) in (6) ignores the initial age distribution by assuming a constant number of infected I 0 . This does not affect the analysis for equilibrium and stability (because the term tend to zero as t → ∞) but it makes the derivation of reduction to ODEs much more tractable analytically. From the solution expressions for y(t, a) given in (23) and using the expression for B T given above, we have Changing variables by letting

Appendix B: Proof of dN(t) dt = 3 − N for General Distributions
This proof is for the system (9) with general distributions for P I R , P I T , and P T R in the case of δ = 0. Other distributions are exponential as given in (1). We ignore the X 0 (t) terms (X = I , T , R) as they follow the same argument as for the main part of the corresponding variables.
Differentiating the I equation in (9) Similarly, from the T and R equations in (9), we obtain From the definition of T X and D X in (33) and (34) It follows that From the first equation in (36) and using (47)