A threshold delay model of HIV infection of newborn infants through breastfeeding.

The breast milk of HIV infected women contains HIV virus particles, therefore children can become infected through breastfeeding. We develop a mathematical epidemiological model of HIV infection in infants, infected children and infected women that represents infection of an infant/child as a series of exposures, by incorporating within-host virus dynamics in the individuals exposed to the virus through breastfeeding. We show that repeated exposures of the infant/child via breastfeeding can cause bi-stability dynamics and, subsequently, infection persistence even when the epidemiological basic reproduction number R 0 is less than unity. This feature of the model, due to a backward bifurcation, gives new insight into the control mechanisms of HIV disease through breastfeeding.

Introduction

In recognition of the benefits of breastfeeding on the ultimate health of an infant/child, the World Health Organization (WHO) recommends early initiation of breastfeeding, exclusive breastfeeding during the first 6 months of life, and continued breastfeeding until 24 months of age (see Chasela et al. (2010); WHO (2017a)). However, women infected with HIV can transmit the HIV virus to an infant through breastfeeding Kourtis, Lee, Abrams, Jamieson, and Bulterys (2006); Van de Perre et al. (1991). It has been shown that this risk is greatly reduced when infected breastfeeding women adhere to antiretroviral therapy regimens Mofenson (2010). However, when HIV status is unknown, or when antiretroviral therapy doses are missed, the benefits of breastfeeding in infant/child health will not be optimal.

The effects of individual level exposure to HIV through breastfeeding on its spread on the population level are not well understood. The scenario is also more complex when antivirals are introduced into the system, as these affect the population level transmission of HIV, and also the dose size of HIV virus that infants are subjected to. Mathematical models that integrate the in-host and epidemiological levels of HIV pathogenesis and persistence can provide useful tools needed to study the effects of HIV exposure to infants in a population. Our goal is to determine the effects that individual-level parameters have on the population level persistence of the disease, and vice-versa. Here, we develop an immuno-epidemiological model that integrates the in-host pathogenesis of HIV infection in infants that are exposed to the virus through breastfeeding, and embed this model into a population level model of HIV persistence in a sub-population of breastfeeding women. This represents a first step towards the development of general immuno-epidemiological model frameworks of HIV exposure through breastfeeding and HIV persistence at the population level.

There is strong evidence that the longer the duration of breastfeeding, the greater the risk of HIV transmission, i.e. the risk is cumulative (see Leroy et al. (1998); Miotti et al. (1999); Breastfeeding and HIV International Transmission Study Group et al. (2004)) with the cumulative risk of transmission through breastfeeding being 10–14% for infants aged 4–6 weeks and 22% for infants aged 18–24 months (Leroy et al. (2002); Nduati et al. (2000); The Petra study team (2002)). Several research teams have successfully designed experiments which confirm this hypothesis (see, for example, Hessell et al. (2009); Kersh et al. (2009); Van Rompay et al. (2005)). Despite great success of mathematical models having been used to investigate the impact of a single-exposure on infection diseases (see Anderson and May (1992); Brauer (1990); Cooke and van den Driessche (1996)) and references therein), the body of literature on the consequences of multiple-exposures for disease dynamics is rather scant, with only an abstract representation given in Qesmi, Heffernan, and Wu (2015). In this paper, we present an immuno-epidemiological model with threshold delay that captures the dynamics of HIV transmission to infants through repeated low-dose exposures through breastfeeding. The model presented here can be viewed as an extension of that given in Qesmi et al. (2015) with particular focus on an identified route of HIV transmission from mothers to infants.

The modelling study Qesmi et al. (2015) depicts an infection process as a series of exposures to a pathogen, whereby, upon the first encounter with the pathogen, an individual is considered to be exposed but not yet infected. The pathogen load in-host is affected by this first exposure to the pathogen and any subsequent additional exposures until a viral load threshold is reached, upon which the individual is considered to be fully infected and infectious. Susceptible-Exposed-Infectious (SEI) framework was used to capture infection dynamics at the population level. This approach gave rise to a threshold-delay system that was subsequently reduced to a state-dependent delay system. The study showed the possibility of a backward bifurcation, and therefore, of a regime of bi-stability when the basic reproductive number, , is less than 1. The study also established the existence of a threshold value for , such that if , the disease could be eradicated from a population of hosts.

Here, we have adapted this multiple-exposure approach to model HIV dynamics in a population of newborn infants and HIV-positive females of child-bearing age. Inherently, such a population will have different population- and in-host-level dynamics than the more general one considered in Qesmi et al. (2015). To wit, infants upon entering the fully-infected stage lack the capacity to infect. Only a fraction of infected children will eventually become infectious adult females (subject to both gender, as well as to the increased probability of not surviving to reach maturity.) Thus it is necessary to distinguish between the two types of infected individuals. On the other hand, infection via breastfeeding is not a sole or the main route of HIV infection in a population. Since this model focuses on the target population of breast-fed infants and infected adult females of child-bearing age and that are breastfeeding (henceforth, named infectious adult females or infectious women), we assume that viral growth and decay are functions of the total population size of infectious women , thus forgoing HIV dynamics in the entire population of healthy and infected children, women and men. This means that we assume that the prevalence of HIV over the entire population is reflected in the population of , and we assume that this prevalence is relatively stable i.e., we are modelling a large population that is stable and sustains a large population of infectious adult females. We note that the case considered in Qesmi et al. (2015), with the viral growth rate being a constant, is a particular case of the more general setting used here. As a consequence, the resulting system is more complex and general.

Analysis of the new model also shows the possibility of the coexistence of multiple stable equilibria when the epidemiological basic reproduction number is less than unity but greater than some positive constant . Here, we are able to derive an expression for the constant in terms of the model parameters. Finally, we implement very intensive numerical simulations to examine the impact of in-host parameters on the occurrence of bi-stability, and the domains of attraction of co-existing stable equilibria. Our findings show that the in-host dynamics parameters are key determinants whether bi-stability is present. This has important implications for the disease control strategies (such as antiretroviral therapy programs for lactating women) that need to be placed.

The article is organized as follows. Section 2 presents the model, that originally takes the form of an age-structured system involving the age-since exposure variable. This system is then reduced to a system of differential equations with threshold delay. Using the Global Implicit Function Theorem, in Section 3, we transform the system with threshold delay to a state-dependent delay system of differential equations. Properties of the resulting system, such as positivity and existence of steady states, are established in Section 4. Then, in Section 5, using a Lyapunov-Razumikhin function we give a sufficient condition for global asymptotic stability of the trivial steady state. In Section 6, we investigate system's dynamics using bifurcation analysis and show that the system can go through a transcritical and a backward bifurcation. We follow up by the sensitivity analysis with respect to parameters and initial data in the region of bi-stability of solutions, and we end with a discussion of the model results.

The model

In-host model

HIV attacks CD4 T-cells by integrating the viral DNA into the cell DNA. Subsequently virus particles bud from the infected T-cells and infect new cells Zhang et al. (2015). As a result of the infection, infected cells have a shorter life span than their non-infected counterparts, and their population starts to deplete (McCune (2001); Perelson and Nelson (1999)). The course of HIV infection includes an acute phase in the beginning of infection (characterized by an increase to a peak viral load followed by a decrease to a low level of pathogen), and a latent phase, when the T-cell count and viral load reach a quasi-steady state at a low pathogen load, which can last for many years. The final stage of HIV infection, called AIDS, begins when the CD4 T-cell count in the blood declines below a certain threshold (200 cells per ), and the immune system becomes compromised, leaving the individual susceptible to opportunistic infections and diseases.

For the sake of modelling viral dynamics within an infant/child, we consider an in-host process at the start of infection, at the beginning of the acute stage.

In the early stage of the breastfeeding, the virus population dynamics is simply governed by the growth rate, denoted by θ. However, If the mother is infected, breastfeeding can cause transfer of virus particles through lactation. This was corroborated by an experiment Van Rompay et al. (2005) where repeated low-dose exposures of infant macaques led to an increase of the virus load until an infection threshold was reached. Along with the HIV-infected cells, neutralizing HIV-specific antibodies are passed to the infant in the breast milk Overbaugh (2014). This means that infants have a clearance mechanism for virus that is budded from cells. Therefore, the growth rate should be written aswhere is the intrinsic rate of the virus growth and h is the clearance rate of the virus due to neutralizing antibodies transferred through the breast milk of infected females.

Each of the above intrinsic growth rate and clearance rate can depend on the level of infectious woman , which reflects on the level of infection in the entire populations. In particular, the larger the , the greater the burden on the health resources, and thus the smaller resources available to control the virus growth and increase the clearance rate. Therefore, the growth rate is an increasing function of . While the passive antibodies conferred to infants through breastfeeding contribute to protection by reducing virus levels in breast milk and thus, diminishing its infection force (see John-Stewart et al. (2004); Shapiro et al. (2007)), it has been shown that these antibodies are not completely effective in blocking HIV transmission (see Overbaugh (2014); WHO (2007)). Also, neutralizing antibodies do not contribute to the inhibition of virus production from infected cells in the infant's body (see Klimpel (1996)). We thus assume that the intrinsic growth rate dominates over the killing rate of the virus

During the early stage of the breastfeeding, the viral load of the child is changed due to the internal viral dynamics and, at the same time, due to the repeated exposure to the infection. Hence, it is important to track the viral load with respect to the time a since the first exposure of an uninfected infant. This leads to the following structured virus dynamics model:where is the intrinsic growth rate of the virus in the infant and expresses the additive viral load in the infant as a function of infectious adults females, , due to repeated exposure through breastfeeding. We assume that the function F is positive and non-decreasing. This, again, reflects that fact that, when the number of infectious females is low (high), the health resources available are high (low), and the added amount of virus through breastfeeding would be low (high). Finally, we assume that there is a maximum amount of virus that can be added through each exposure to the virus through breastfeeding. A Holling type II functional response (see Holling (1965)) is a suitable choice as the formulation for the additive viral load due to multiple exposures

Here, b is the number of effective contacts with infectious adult females, is the viral load introduced at each exposure, and k is an adjustable parameter which controls how fast the saturation occurs.

Similar to Qesmi et al. (2015) we assume that once a certain viral load threshold, A, is achieved, the exposed infant will become infected. Therefore, Eq. (1) governs the viral load when only. Finally, since A is the viral load needed for infection, we assume that the viral load from a single exposure is less than A, i.e. .

In the following, we will state some properties of the viral load function which will be used in the next sections.where is the history function of the infectious population defined for by

There exists a positive function , such that

for all ;

is equivalent to ;

is bounded;

for all , the following is true:

The proof of Lemma 1 is similar to the one given in Qesmi et al. (2015), and hence is omitted here.

Between-host model

The general population dynamics, captured by Figure H.1, is modelled as follows. First, the population under consideration is comprised of susceptible infants that can be exposed to infection through breastfeeding (S), exposed infants (E), infected children/infants () and infectious adult females ().

We note that this joint subpopulation is a part of a general population that resides at an endemic equilibrium. Contribution to the general population level infection spread by the individuals infected through breastfeeding is deemed negligible.

For simplicity, we assume that all infants enter into the model through the susceptible infant class (S) with rate π (individual year−1). We assume that the population is sufficiently large such that the birth rate can be approximated as a constant. Individuals leave this compartment with rate d (year−1) subject to death and termination of breastfeeding. While in the compartment, infants are exposed to the infection with rate year−1. Although an infant is likely to be breastfed by their mother exclusively, we model the probability of infection as being directly correlated to the overall infection prevalence in the population, and thus, the number of infected women who breastfeed reflects the overall level of infection in the population. This assumption allows us to integrate the availability of health resources into the model. Since a small (large) size of reduces (increases) the probability a mother is on effective drug therapy, and thus reduces (increases) the probability that an infant will be exposed to the HIV virus through breastfeeding from its mother. It is assumed that infants born to infected women are exposed to the virus with rate (year−1) at birth, but that infants exposed to the virus through breastfeeding are exposed with rate (year−1). Thus, where N constant is the total population of healthy and infected women.

Once an infant has been exposed to the virus, they move to the exposed infants compartment (E). The length of sojourn in this compartment is subject to an average in-host dynamics of viral growth and a death rate, , a positive, continuous, non-decreasing function of age-since-exposure. We assume that the duration of stay in the exposed compartment is relatively short (see Van Rompay et al. (2005)). The dynamics of this compartment is traced in terms of time, t and age since the first exposure, a.

When the internal cumulative viral load exceeds the threshold value A, or equivalently, when the age since exposure passes the period , an exposed infant becomes infected and moves to the compartment. We track individuals in this compartment in terms of time, t and age since reaching viral threshold, a. Observe, that since the sojourn in the exposed compartment is relatively short (several weeks) the latter can also be used to track the population dynamics in this compartment by age. Individuals leave it either due to dying or due to reaching maturity, with rates and , respectively. Both are subject to the age of infection and are positive, continuous non-decreasing functions of age.

A proportion p of infected children in class proceed to the mature infectious female class . However, individuals can be drafted into this compartment of infectious female breastfeeding mothers in two ways: infectious females that were infected with HIV by routes different than breastfeeding and infectious females that acquired the infection through breastfeeding, matured and upon maturation proceeded to become breastfeeding mothers. The former term is given by reflects the fact that the incoming rate of individuals into the breastfeeding compartment ( is indicative of the overall prevalence in the population. We track the dynamics of the compartment in terms of time, t and age since initial entrance to the compartment, a. Individuals leave the compartment due to death or termination of breastfeeding with rates α (year−1) and , respectively.

The model takes the following form:subject to the following boundary conditions

Reduction to a state-dependent delay system

First, observe that using the method of characteristics (see Webb (1985)), the solution of Eq. (3b) with boundary condition given by Eq. (4a), is given, for and , by

Let be the maximal delay, i.e. . Then, following a sequential set of procedures similar to Qesmi et al. (2015), for , the threshold delay model (3) reduces tosubject to the following boundary conditions

Set μ (year−1) given by to be an average death rate of an infected child and γ given by - an average rate of leaving the compartment due to the maturation. We reduce system (6) to the following form:for . For the details see Appendix A.

Note that using the process we applied to solve Eq. (3b), we can reduce system (6) to a nested threshold-dependent delay system. While the resulting model is more accurate than its non-nested analog, the two share many dynamical properties (attributable to the threshold delay) with the latter being more analytically and numerically tractable.

Next, we will show that the threshold-delay system (8) can be reduced to a system of differential equations with a state-dependent delay. Since the exposure period is relatively short, we can assume that the average virus growth rate is relatively constant, i.e.,where r is a positive constant. Thus, Eq. (2) can be read as

Let C denote the space of continuous functions defined on . A function is said to be decreasing if, for all , (for all ) implies that .

There exists a unique continuously differentiable function such that, for , . Furthermore, function σ is decreasing and

The proof of Proposition 1 can be found in Appendix B. Thus, for , the threshold delay system (8) is equivalent to the following state-dependent delay differential equation system:where , a history function of .

Set , and . Then, for , system (9) can be written as a state-dependent delay differential systemwhere , a history function of . System (10) writeswhere is defined by for , and function is given, for , the space of continuous functions on , by

Since function σ is continuously differentiable on C, g is also continuously differentiable on . Therefore, existence and uniqueness in forward time of a solution of system (10) with initial data in (space of continuously differentiable functions on ) follow from Walther (2003) and Hartung, Krisztin, Walther, and Wu (2006).

Properties of the model and existence of steady states

Here we focus on basic properties of solutions of system (10), such as positivity and boundedness. The expression for the disease-free equilibrium is stated and conditions for the existence of the endemic equilibrium are derived, as well as, the threshold condition for the eradication of the disease in a population.

We consider solutions of system (10) with initial data in . Therefore unique solutions of system (10) exist for . The following result is straightforward and can be obtained, for instance, using a method similar to the one presented in Adimy, Crauste, Hbid, and Qesmi (2010).

Solutions of system (10) with non-negative initial data are non-negative and bounded for .

System (10) has a disease-free equilibrium (DFE) given by. Next, we consider the existence of endemic steady states of system (10). Let and be positive real numbers given byandsuch that is the global maximum of . Note that, here, is the basic reproductive ratio (or basic reproduction number) of system (10).holds. Then, system (10) has:

Assume that, for all z positive,

at least one endemic equilibrium if ;

no endemic equilibrium if .

If inequality (14) does not hold, then system (10) has:

at least one endemic equilibrium if ;

at least two endemic equilibria if ;

no endemic equilibrium if .

Moreover, every endemic equilibrium, satisfies

The proof of Proposition 3 is omitted since it is somewhat similar to the one given in Qesmi et al. (2015).

In the case when and δ is sufficiently close to 0, the system has a unique endemic equilibrium. Also, in the case of being a constant function and , the system has either no or exactly two equilibria.

Asymptotic stability of the trivial steady state

In the following theorem, we give a sufficient condition for global asymptotic stability of the trivial steady state of system (10), provided some conditions are satisfied. Let ρ be the positive real value given by such that holds, then the disease-free equilibrium of system (10) is globally asymptotically stable with respect to solutions with non-negative initial data.

If

The proof of the theorem can be found in Appendix C.

Next, we consider the local asymptotic stability of the disease-free equilibrium of system (10), .

Let be the map defined by

The derivative of J is given by

The linearization of system (10) at an equilibrium is given by

Substitute the Anstaz , where to obtain

Set and divide by , to obtain and where λ is a root of the characteristic equation

From the analysis of Eq. (18) we obtain the following result (more details can be found in Appendix D.)

If , then the disease-free equilibrium (DFE) of system (10) is unstable. If , then is locally asymptotically stable.

Bifurcation analysis

Transcritical bifurcation

This section is devoted to the analysis of local asymptotic stability of an endemic equilibrium of system (10) when condition (14) holds. To show existence of a transcritical bifurcation, we investigate the sign of real parts of roots of Eq. (18).

Throughout this section, we assume that function σ is given by , where ν is a positive constant and is a positive, decreasing, bounded and differentiable function. Here, is a space of positive continuous functions. Then, system (10) becomes

We assume that condition (14) holds. Define

Note that h is an increasing continuous function on and thus, is invertible. Therefore, the condition is equivalent to

Intuitive explanation of the above equation is that a long period of exposure renders an endemic equilibrium impossible. Therefore, existence of a positive steady state depends on parameter ν. On the other hand, according to Proposition 3, satisfies

Thus, using the Implicit Function Theorem, is a decreasing continuously differentiable function of ν. Furthermore, using Eq. (20) and Eq. (21), we obtain

Define for ,

Let denote the equilibrium for which is the first solution of equation such that . Treating ν as a control parameter, we set out to prove the existence of the transcritical bifurcation of system (19), as ν moves through (respectively, moving through 1).

From Eq. (18), the characteristic equation associated with iswhere

The next result states the existence of a transcritical bifurcation of the positive steady state at .

If condition (14) holds at , one of the positive steady states of system (19) undergoes a transcritical bifurcation. That is for , ν close to , the positive steady state is locally asymptotically stable and the disease-free steady state is unstable, and for the DFE is locally asymptotically stable.

A proof of the above theorem can be found in Appendix E.

Backward bifurcation

If condition (14) does not hold and if holds, then according to Proposition 3, there exist at least two endemic equilibria, and . Therefore, a backward bifurcation may occur for values of close to 1.

Denote the two equilibria by and , such that and are the first two solutions of equation given by Eq. (23) with and , i.e.respectively. We state the following theorem.

If condition (14) does not hold at , then system (19) undergoes a backward bifurcation. That is, for , ν close to , there is a positive steady state that is locally asymptotically stable and the disease-free equilibrium is unstable; and for , ν close to , both the disease-free equilibrium and one of the endemic equilibria are locally asymptotically stable, whereas a second endemic equilibrium exists and is unstable.

A proof of the above theorem can be found in Appendix F.

Note that, similar to the model in , in the case of single exposure, i.e, if , system (10) is a system of constant-delay differential equations, which does not go through a backward bifurcation, since in such a case condition (14) holds true for all viable parameter values.

The remark can be generalized as follows.

There exists a threshold value of b, , above which the backward bifurcation appears, and a threshold value of b, , below which the backward bifurcation disappears.

A proof of the above proposition can be found in Appendix G.

Numerical simulations and sensitivity analysis

The system

In this section we provide illustrations for Theorems 3 and 4. Observe that it is challenging to explicitly solve for σ and therefore, it is difficult to locate the value of . However, since the condition of in Theorems 3 and 4 is equivalent to , we will use π to demonstrate bifurcation dynamics.

In order to generate time trajectories of solutions we consider system (8). By substituting and we obtain

Differentiating Eq. (26) with respect to t and solving for yields

Set

Then

We replace with its definition and therefore obtain a system in S, , , τ and with a state-dependent delay. A number of available numerical solutions software can be applied to generate the simulation results presented here. We have used the MATLAB (2015) function, ddesd.

To demonstrate a setting where both bifurcations are possible, we consider a linear function approximation for θ:

Parameter value calculations

We commence from the discussion of susceptible, exposed and infected infants mortality rates. Observe that the equations capturing dynamics of the infected infants and infectious adults females are given by

In the age-structured model for infected infants and infected females, the dynamics are described as followssubject to the following boundary conditions

Here, in system (32) is the death rate, a non-decreasing positive function of infection age such that .

Observe that in age-structured formulation of compartment, using the characteristics method, the number of individuals at time t with the age of infection a is given by:

Therefore the probability of surviving until the age of T (after becoming fully infected) is given by . On the other hand, from system (31) it follows that the probability of surviving to maturity age is given by . Therefore, γ in system (31) is given by

We estimate μ using the data from Congo. The disease is endemic there, and, according to WHO data, the prevalence of HIV in adults aged 15 to 49 was at roughly 3% from 2005 to 2016. Thus it is accurate to say that the country is witnessing the endemic equilibrium (see WHO (2017b)). We obtain the infected infants death rate of year−1 and transfer rate out of infected infants compartment year−1.

The leaving rate of a healthy susceptible infant, d, is a sum of a death rate and growing out of breastfeeding (presumed to be 2 years) rate. The death rate for healthy infants was taken from a country whose HIV-prevalence rate was low (Tunisia, in 2016, WHO (2017b)) by considering the death rates in populations of infants less than one years old and 1–4 years. Therefore, we obtain year −1.

Similarly, we assume that women leave the infected compartment either at the end of the child-bearing process (50 years of age) or due to the death (we have considered weighted average of death rates for females aged 20–49 years in Congo, 0.004784 year−1). Thus, . We assume that half of the infants are female, this . The incoming rate, π is taken in a range where backward or transcritical bifurcation take place and infection rate was adjusted so that the necessary condition for endemic equilibrium viability, , holds. Parameter values are summarized in Table 1.

Table 1

Parameter values used to simulate the dynamics of system (8).

Parameter	Description	Units	Value
d	leaving rate of susceptible infants	year−1	0.512
μ	death rate of infected infants	year−1	0.0099
γ	transition rate from infected infant to infected adult	year−1	0.995
p	fraction of infected infants that grow up to become infected adult females		0.5
α	exit rate of infected females	year−1	0.0381
β	transmission rate	(individualyear)^−1	7.4e−3
δ	death rate of exposed infants	year−1	1e-2
r	average virus growth rate	year−1	0.1, 1
m	adjustable constant controlling the intrinsic virus growth rate	year−1 individual−1	1e-11, 1e-5
c	viral load introduced at each exposure	virions	20
b	number of effective contacts with infectious adult females	year−1	16
k	adjustable saturation constant		1e-6, 1
A	viral load needed for infection	virions	1e3, 1e5

Numerical investigation of Eq. (26) yields that the switch between backward bifurcation and transcritical bifurcation happens due to the variations of the in-host parameters. In particular, if the virus growth rate constant r is sufficiently small, the exposure delay causes appearance of a regime of bi-stability for a range of values of π, such that when , solutions can get attracted to either an endemic equilibrium or the disease-free equilibrium. We also observed that the smaller is r, the larger is the interval where bi-stability takes place. The results are shown in Figure H.2, plotting vs π.

Next we investigate the dependence of components of an attracting endemic equilibrium on parameters using sensitivity analysis. It is also used to evaluate the basin of bistability (in the case of a backward bifurcation) for different values of π in the region of bistability. The method is adapted from Marino (2008) and is based on the method outlined in Marino, Hogue, Ray, and Kirschner (2008), which works as follows. First, we identify inputs: parameters and their ranges, and generate a large number of unique parameter sets (bins) using Latin Hypercube Sampling (LHS) Mckay, Beckman, and Conover (2000). Then the magnitude and strength of the relationship between the output metric and each of the input parameters is assessed using Partial Rank Correlation Coefficients (PRCC). The coefficient assumes values in the interval and the magnitude of the coefficient equal to 0.5 and higher signifies a strong correlation between an input parameter and the output variable.

We start by considering the sensitivity of the system's attractor to the perturbations of parameters. We selected the parameter values corresponding to the attractor landscape where a backward bifurcation takes place, as π increases from zero. The initial data was taken to be a function mapping to , set to DFE for and to a perturbed stable endemic equilibrium in the region of bi-stability at time . Each parameter was sampled within 80–120% interval of its base value assuming the uniform distribution. Altogether 8000 sets of parameters were generated. For each set of parameters, a condition for the existence of an endemic equilibrium was evaluated and the sets with no viable endemic equilibrium were discarded (4202 sets were accepted).

The results are captured on Figure H.3. We observe that while the virus dynamic parameters determine the attractor landscape, they do not influence the location of the equilibrium itself. The most important parameters for all four components of the system are the infection rate β, the proportion of infected infants that will proceed to be infectious adults p, and the average lifetime of an infectious female . S does not have a dependence on any other parameters. However, the remaining three components have a pronounced relationship with both the recruitment and death rate of susceptible infants. In short, an increase of the average life span of a susceptible infant causes an increase in the equilibrium values of both infected compartments. Not surprisingly, the average death rate of an infected infant, γ, has a correlation with the number of infected infants at an endemic equilibrium.

Next, we consider sensitivity of the system to the initial data in the region of bi-stability with three different values of π sampled and the rest of parameters fixed. The base for the initial data was a perturbed stable endemic equilibrium in the region of bi-stability, at time and DFE for . The components of the state at time , , were sampled on the following intervals:using LHS with 3000 bins. The results are captured in Figure H.4. We observe that in the region of bi-stability, provided that there is a non-zero number of infected individuals (infants or adults) at time , the initial value of the delay determines whether the solution will converge to a DFE or an endemic equilibrium. If it is significantly large, the solution will converge to the DFE.

Next we investigate the role of the initial size of the infected populations on the attractor of the system. From Figures H.5 and H.6 it follows (for the interval that we have considered) that the initial value of τ is the key determinant of the shape of the basin of attraction. However, for a certain range of , the total number of infected individuals in the system has a significant correlation with the attractor. We observe (Figure H.5) that joined initial population of infected adults and children determine the basin of attraction, although they fail to do so separately. For example, it is evident that there are values of which cause a trajectory to converge to an endemic equilibrium and (depending on the remaining components of the initial data) can also cause the trajectory to converge to the disease-free equilibrium. We also observe that as the incoming rate of susceptible infants, π, increases, the basin of attraction of the DFE narrows, i.e., as π increases, the influence of total number of infected people on the division of basins of attraction becomes less pronounced. Finally, we observe (Figure H.7) that the size of does not affect the basin of attraction.

Discussion

Breastfeeding is an important route of mother-to-child HIV transmission. Immuno-epidemiological models, that embed individual levels of pathogenesis into population level models of disease, provide the tools necessary to study the effects of individual-level effects of infectious diseases on the population, and vice versa.

In this paper we proposed an immuno-epidemiological model with threshold delay to investigate the effects of the duration of breastfeeding on the disease dynamics in a population of infants, infectious women, and infected children. This model represents an important first step towards the development of a larger immuno-epidemiological framework that also includes the disease dynamics of HIV in the entire population (including of healthy and infected men, women, children, and infants), not only given the fact that it is the first of its kind, but that the model analysis shows a clear connection between population level results and individual-level parameters. We now summarize the results and provide a discussion on the connection between the two scales.

The basic reproduction number, , as given by Eq. (12), is a key parameter in the model. The non-delayed analog model has a globally asymptotically stable DFE when the associated , and it is unstable . However, an may not be sufficient to eradicate the disease. When condition (14) holds achieving an will eradicate the disease (the DFE is locally asymptotically stable), but when condition (14) does not hold and the system undergoes a backward bifurcation. The existence of the backward bifurcation depends on individual-level parameters. In Proposition 4 we showed that the existence of the backward bifurcation depends on the effective number of contacts b. Therefore, the analysis of the proposed model indicates that while the population level processes determine the level of the disease prevalence in the population being studied, an in-host process parameter determines the possibility of bi-stability. Given the connection between the population level outcomes and the in-host level parameters found here, a direct conclusion is to allocate resources towards the determination of accurate approximations of these parameters, and a thorough knowledge of the in-host process and response functions.

There are several limitations to the model presented here. Important extensions include: a density dependent birth rate so that the effects of mother-to-child breastfeeding can be studied in smaller populations; an in-host model of HIV viral load in infected women, to allow for different infectious doses over the entire course of HIV infection; a stochastic in-host model of HIV pathogenesis in infants to study the probability of viral clearance after exposure; the explicit consideration of health resources available to the population (i.e., therapeutics, community programs), and adherence regimens to any drug therapies being studied; a comprehensive model of HIV transmission in an entire population, including healthy and infected men, women, children and infants.

Threshold models for infectious disease transmission dynamics have been formulated and analyzed before, and some of the earlier relevant developments can be found in the monograph Walthman (1974). A case study in Smith (1993) suggested that in many instances, as is in this paper, structured population models can be reduced to threshold delay equations through integration along characteristics. In Barbarossa, Hadeler, and Kuttler (2014), it was shown that a class of state-dependent delay equations can be derived from some fundamental biological principles governing the birth, death, and maturation processes in which the length of the juvenile phase depends on the total adult population size. The breastfeeding infection model we studied describes a similar epidemiological process, where the length of time, during which infants remain exposed but not yet infected, depends on the total infected female population.

It remains to be seen whether the novel class of neutral equations similar to that derived by Barbarossa et al. Barbarossa et al. (2014) can also occur in the epidemiological setting. While our present study focuses on steady state dynamics and backward bifurcations, an interesting issue to be addressed in a future study is whether nonlinear oscillations can be observed when the assumption of the birth rate being a constant is relaxed and replaced by a more nonlinear feedback depending on the total female population. This will require the applications of some type of Hopf bifurcation theory or fixed-point theory. As the pioneering work Arino, Hadeler, and Hbid (1998a, b) shows, this will be a highly non-trivial task since the state-dependent duration from initial exposure to final stage of being infection of infants due to breastfeeding adds additional multiplicity of zero eigenvalues for the linearized systems.