A singular perturbation approach to epidemics of vector-transmitted diseases.

In vector-borne epidemic models there is often a substantial difference between the vector and host time scales. This makes it possible to use the quasi-steady-state to obtain final size relations.

Introduction

Many diseases are transmitted from human to human indirectly, through a vector. Many vectors are bloodsucking insects that ingest disease-producing microorganisms during blood meals from an infected (human) host, and then inject it into a new host during a subsequent blood meal. The best known vectors are mosquitoes for diseases including malaria, dengue fever, chikungunya, Zika virus, Rift Valley fever, yellow fever, Japanese encephalitis, lymphatic filariasis, and West Nile fever.

Every year there are more than a billion cases of vector-borne diseases and more than a million deaths. Vector-borne diseases account for over 17% of all infectious diseases worldwide. Malaria is the most deadly vector-borne diseases, causing an estimated 627,000 deaths in 2012. The most rapidly growing vector-borne disease is dengue, for which the number of cases has multiplied by 30 in the last 50 years. These diseases are found more commonly in tropical and sub-tropical regions where mosquitoes flourish and in places where access to safe drinking water and sanitation systems is uncertain.

Vector transmitted diseases require models that include both vectors and hosts. For most diseases transmitted by vectors, the vectors are insects, with a much shorter life span than the hosts, who may be humans as for malaria or animals as for West Nile virus.

The compartmental structure of the disease may be different in host and vector species. For many diseases with insects as vectors an infected vector remains infected for life so that the disease may have an or structure in the vectors and an or structure in the hosts. We will describe vector models with structure in the host species and structure in the vector species, but the analysis of other types of vector transmitted diseases is similar.

A basic vector transmission model

We consider a basic model for a vector transmitted disease epidemic, in which we assume throughout that hosts satisfy a simple model and vectors satisfy a simple model, with no recovery from infection. This model is a template for vector disease transmission models for many specific diseases. We are thinking of mosquitoes as vectors, and because a mosquito lifetime is much shorter than that of the human hosts we must include demographics in the vector population.

We assume a constant total population size of hosts (humans), divided into susceptibles, infectives, and removed members. Infected hosts recover at rate γ.

There is a constant birth rate of vectors in unit time and a proportional vector death rate in each class, so that the total vector population size is constant. The vector population is divided into susceptibles and infectives. Infected vectors do not recover from infection.

We assume that the mosquito biting rate is proportional to the size of the host population. Thus an average mosquito makes bites in unit time. The total number of mosquito bites in unit time is and the number of bites received by an average host in unit time is . We assume that is the probability that a bite transmits infection from vector to host and is the probability that a bite transmits infection from host to vector. We define

If we eliminate b from these two equations, we obtain a balance relation

In many models for vector-transmitted diseases, the vector population size is much larger than the host population size but , suggesting that .

A susceptible human receives effective mosquito bites in unit time, of which a fraction is with an infective mosquito. Thus the number of new infective humans in unit time is

A similar argument shows that the number of new mosquito infections is

The model is

Because is the constant , we may replace by and reduce the system (1) to a three-dimensional system,

The basic reproduction number

The basic reproduction number is defined as the number of secondary disease cases caused by introducing a single infective human into a wholly susceptible population of both hosts (humans) and vectors (mosquitoes). For the model (1) this may be calculated directly. There are two stages. First, the infective human infects mosquitoes, at a rate for a time . This produces infected mosquitoes.

The second stage is that these infective mosquitoes infect humans at a rate for a time , producing infected humans per mosquito. The net result of these two stages is infected humans, and this is the basic reproduction number .

We could also calculate the basic reproduction number by using the next generation matrix approach (van den Driessche & Watmough, 2002). This would give the next generation matrix

The basic reproduction number is the positive eigenvalue of this matrix,

In this calculation, the transition from host to vector to host is considered as two generations. In studying vector transmitted diseases it is common to consider this as one generations and use the value that we obtained by our direct approach

This choice is made in (Chowell et al., 2007) and (Kucharski et al., 2016) and is the choice that we make because it conforms to the result obtained directly. However, other references, including (Pinho et al., 2010), use the square root form, and it is important to be aware of which form is being used in any study. The two choices have the same threshold value.

In fact, different expressions are possible for the next generation matrix. This is shown in (Cushing & Diekmann, 2016). Using the next generation matrix approach but considering only host infections as new infections and vector infections as transitions, we would obtain the form (3).

Fast and slow dynamics

In vector-borne disease transmission models in which the vector is an insect, the vector time scale is often much faster than the host time scale. In such models it is possible to consider two separate time scales (Brauer & Kribs, 2016). We can make a quasi-steady-state hypothesis, assuming that the vector population sizes remain almost constant (Segel & Slemrod, 1989). We treat the vector population sizes as constants which depend on the host population sizes and approximate the system by a host model, which has fewer equations than the full system but may be more complicated in form.

A model for dengue fever of the form (2) has been studied in (Rocha, Aguiar, Souza, & Stollenwerk, 2013). In terms of the above model, with time measured in days, this model has parameter values with much smaller than and much smaller than , expressing the fact that the vector dynamics are much faster than the host dynamics, and this suggests a modelwith ε a small positive constant.

Other models describing two time scales may be found in (Souza, 2014; Takeuchi, Ma, & Beretta, 2000). In fact, such a form is valid for all vector disease transmission models in which the vector population time scale is much faster than the host population time scale.

For differential equations or systems of differential equations which depend on a parameter there is a general Theorem to the effect that solutions are continuous functions of the parameter on any finite interval. However, if a derivative is multiplied by a parameter which may be allowed to tend to zero, this is not necessarily true. Such a situation is called a singular perturbation. Many problems in the biological sciences involve actions on very different time scales, and these may lead to a rapid change in some of the variables on a very short initial time interval while other variables act more slowly.

Singular perturbations

Singular perturbation problems arise in models (systems of differential equations) containing a small parameter ε, of the formwith solution . There is a corresponding reduced system obtained by setting ,with solution

Since ε is assumed to be small, the form (4) suggests that the y reaction time is much faster than the z reaction time. Thus y goes to its equilibrium value rapidly, and at its equilibrium . Then we might expect that the reduced problem (5) is a good approximation to the full problem (4) after a short initial time interval near during which y moves to its equilibrium value. In applications one often makes a quasi-steady-state hypothesis, that y remains almost constant, so that . This hypothesis is expressed as ; in singular perturbation language the hypothesis is just that the full problem is approximated by the reduced problem.

Because (5) is a first order differential equation (which requires one initial condition to identify a unique solution) and (4) is a two dimensional system (requiring two initial conditions for a unique solution) we must expect to lose an initial condition in the reduction, and this suggests that the solutions of (5) and (4) (each derived on a different time scale) may not agree close to . Because (the solution to the full problem) and (the solution to the reduced problem) do not match at , we should expect that changes rapidly for t close to 0. It is possible to analyze this by making a change of independent variable to change from the slow time scale to the fast time scale. However, since our interest is mainly in the long term behavior of the system we do not explore this further here.

If the partial derivative we may solve the equation for y as a function of z, . Thus the reduced system (5) is equivalent to the first order initial value problem

Then we have the solution of (5) with the solution of (6) and , and . If it is not possible for the solution of the reduced problem (5) to satisfy the two initial conditions of the full problem (4). The solution of the reduced problem is called the outer solution. In order to use the solution of the reduced problem (5) as an approximation to the solution of the full problem (4), we would need a result to the effect that for each t away from the solution of the reduced problem (5) is the limit as of the solution of the full problem (4).

A change of time scale transforms the system (4) to the system

The second equation of (7) says that the second variable z remains almost constant on a large t-interval corresponding to a small τ-interval. This suggests considering the initial value problem called the boundary layer system, as an approximation valid on a small t-interval called the boundary layer.

The mathematical treatment of singular perturbations began in the 1940's from the perspective of asymptotic expansions. A few years later the qualitative result which justifies the use of the reduced system as an approximation to the full system was obtained independently in the U.S.A. and the Soviet Union (Levinson, 1950; Tihonov, 1948).

Theorem (Levinson-Tihonov): Suppose that

are smooth functions,

the equation can be solved for y as a smooth function of z, ,

the reduced system (5) has a solution on an interval ,

the boundary layer system (8) has an asymptotically stable equilibrium.

Then as for . The convergence of y is non-uniform at

There is an extension of this result to infinite time intervals (Hoppensteadt, 1966).

Theorem (Hoppensteadt): Suppose, in addition to the hypotheses of the Levinson-Tihonov theorem, that the reduced system (5) has a solution which is asymptotically stable and that the boundary layer system (8) has a solution which is asymptotically stable uniformly in . Then the convergence is uniform on closed subsets of .

The essential content of these results is that if ε is sufficiently small the solution of the reduced system is a good approximation to the solution of the singularly perturbed system except very close to . The relation is called the quasi-steady-state. Close to the solution of the boundary layer system (8) describes the behavior of solutions. Thus the analysis of a singular perturbation problem can be decomposed into the analysis of two simpler problems, namely the boundary layer system and the reduced problem. Curiously, in fluid dynamic applications the focus of attention has been on the boundary layer system and matching of the two solutions, whereas in most biological applications the primary interest has been in the long-term behavior, that is, the reduced problem.

The underlying idea in a singular perturbation problem is that there are two different time scales inherent in the problem, and this makes it possible to analyze the problem separately on each time scale. The reduction in dimension because of this separation simplifies the analysis.

A vector-borne epidemic model

We assume that the vector dynamics operate on a much faster time scale than the host dynamics,and we define

We require also

We put the epidemic model (2) into singular perturbation form

The quasi-steady-state, which is the case of the equation for in (9), is given by the equation

This expresses as a function of the other variables,

We substitute this form for into the equations for and to give the reduced equation

This approach gives a way to formulate models for vector-transmitted diseases which consist of a system involving only the host variables. While this system is two-dimensional and the original system is three-dimensional, it does have a more complicated form.

A final size relation

Addition of the first two equations of (9) (or (2)) gives

Thus is a decreasing non-negative function and tends to a limit as . Also, its derivative tends to zero, so that as Integration of this equation with respect to t from 0 to gives

Division of the first equation of (10) by and integration gives

Sincewe obtainand

Combining (12) and (13) we obtain the two-sided estimate final size relation

This final size estimate gives upper and lower bounds for .

Some extensions of the model

We have shown how to view the vector-borne epidemic model (1) as a singular perturbation problem and use this approach to derive a final size relation. This is the simplest possible form for a vector-borne disease model. In fact, many vector-borne epidemic models are of the form form. Such a model is

Since is a constant , we may eliminate the equation from the model to obtain

To calculate the basic reproduction number of the model (15) we proceed in two stages much as we did for the model (9) in the previous section. The net result of these two stages isand this is the basic reproduction number .

As for the model (9) we could also calculate the basic reproduction number by using the next generation matrix approach (van den Driessche & Watmough, 2002). This would give the next generation matrix

The basic reproduction number is the positive eigenvalue of this matrix,

In this calculation, the transition from host to vector to host is considered as two generations. Just as for the model (1) it is common to consider this as one generations and use the value that we obtained by our direct approach

As in the previous section, we assume that the vector dynamics operate on a much faster time scale than the host dynamics,and we define

We require also

We put the model (15) into singular perturbation form

The quasi-steady-state is

We may solve the system (17) for , obtaining

Addition of the first three equations of (16) gives

Thus is a decreasing non-negative function and tends to a limit as . Also, its derivative tends to zero, so that as Integration of this equation with respect to t from 0 to gives

Integration of the first equation of (5) gives

We may now use the estimatesto giveand

Combination of (20) and (21) gives the same final size estimates (14) as those obtained for the model (9). This suggests that there should be a generalization to more complicated vector-borne epidemic models. In fact this has been obtained for an age of infection epidemic model in (Brauer, 2017) but the argument assumes a quasi-steady-state for a more general class of singular perturbation problems.

The models considered here have assumed only vector transmission, but the inclusion of direct transmission, such as for the Zika virus, should be straightforward.

We have considered only epidemic models here, but the methods are applicable to models for endemic diseases as well.