The dynamics of a simple, risk-structured HIV model

1.   Introduction

Prevention can help curb the HIV epidemic. Many studies ^[1,2,3,4,5] have demonstrated the efficacy of pre-exposure prophylaxis (PrEP) in reducing HIV acquisition in at-risk populations. New biomedical tools, including long-acting injections ^[6,7,8], slow-release implants ^[9,10], and multi-dose HIV vaccines ^[11,12], are being tested and may provide additional preventive measures.

In the future, scientists will need to design strategies to account for both the variety of preventive measures and the exact distribution of HIV risk as determined by age, gender, mixing preference, sexual behavior, and other individual characteristics. The goal of these new strategies will be to both increase the proportion of individuals practicing prevention and to decrease HIV transmission using newly approved tools. Mathematical modeling will likely become the tool of choice for evaluating these strategies due to its ability to investigate multiple intervention scenarios with differing coverage and targeting criteria.

Despite an increase in the use of individual-based network models in the last decade, deterministic compartmental models built around systems of ordinary differential equations (ODEs) still dominate HIV modeling. Many investigators find deterministic systems easier to both implement and calibrate. In an extensive survey of models studying the effects of HIV/AIDS intervention ^[13], 86.9% of the models were deterministic. In a recent analysis of antiretroviral therapies ^[14], 9 of 12 models were deterministic compartmental models, while 3 were individual-based microsimulation models.

Compartmental models for isolated populations deal with HIV risk in a variety of ways. Many models simply ignore individual differences and assume that risk is spread equally across the entire population.

Another group of models ^{[15,16,17,18]} divide the population into several (usually 2 to 6) distinct risk groups and assign specific behavioral attributes (number of partners, sexual frequency, condom use) to each group. These models are typically built around either transmission or contact matrices ^[19] or contact fraction or mixing matrices ^[20] with many elements. Because of this complexity, these models are often hard to analyze.

More rarely ^[15,21,22], risk is treated as a continuous variable, leading to integrodifferential equations with bivariate functions (kernels) rather than matrices. In general, continuous-risk models are more complicated than models with a finite number of risk classes, and most analyses of continuous-risk models rely on numerical methods. For separable kernels, however, analyses may simplify, as stressed by ^[23] and ^[24].

There is also a growing literature on heterogeneity and risk in metapopulation models. Metapopulation models are compartmental models with a large number of distinct populations. Quite commonly, each population has its own homogeneous dynamics, and the populations are coupled together linearly, due to migration or movement. Many published models ^{[25,26,27,28]} contain susceptible-infectious-susceptible (SIS) dynamics, and these papers often focus on the spread of flu-like diseases between cities due to air travel. Other authors ^[29,30], however, have used metapopulation models to analyze the effects of migrant workers on HIV transmission.

In this paper, we consider a model in which the internal transmission dynamics are not homogeneous. Rather, we focus on a model for a single, structured population with an arbitrary number of discrete risk classes and direct, nonlinear transmission between all of the risk classes. To use the terminology of ^[31], we focus on a cross-coupled model rather than on a mobility metapopulation model. Interestingly, methods developed for continuous-risk models can sometimes simplify the analysis of models with discrete risk classes. In this study, we borrow methods from the study of integral and integrodifferential equations to analyze our model, and we investigate how the number of risk classes affects the model's dynamics. We consider distinct scenarios in which added groups have lower risk due to increased awareness or in which average risk across classes remains constant.

In the next section, we describe our model and spell out our simplifying assumptions. Then, in section 3, we find the equilibria of our model. In section 4, we outline methods for determining the stability of these equilibria. We complement these general analyses with two examples. In section 5, we add low-risk groups, while in section 6, we add groups of intermediate risk. Finally, in section 7, we summarize our work, discuss its limitations, and suggest future avenues for research.

2.   A structured model

Let $r$ and $s$ be continuous and dimensionless risk variables on a scale from zero to one. Zero is the lowest risk. One is the highest risk.

Suppose that there are $n$ risk classes, each of width $\Delta r = 1/n$ or $\Delta s = 1/n$, and that $i$ and $j$ are indices for these risk classes. Let $S _ i (t)$ be the number of susceptible individuals in risk class $i$ at time $t$, $I _ i (t)$ be the number of sexually-active infectious individuals in risk class $i$ at time $t$, and $A _ i (t)$ be the number of infected individuals in risk class $i$ at time $t$ who have developed AIDS and are no longer sexually active. We now consider the model

where $i = 1, ..., n$.

System (2.1) is a risk-structured susceptible-infectious-removed (SIR) model with demography. We assume that each risk class has its own (constant) birth rate, $\alpha _ i > 0$, and that the natural per-capita death rate ($\mu$), AIDS-induced per-capita death rate ($\delta$), and per-capita AIDS development rate ($\gamma$) are all constant and independent of risk class. The $\beta _ {ij}$ are the transmission coefficients between individuals of class $j$ and class $i$. That is, $\beta _ {ij}$ is the per-capita rate of infection of class-$i$ susceptibles per class-$j$ infective. Model (2.1) should, arguably, contain standard incidence, but we have assumed mass-action incidence to greatly simplify subsequent calculations.

To ease comparison with an integrodifferential model, we now shift from numbers to densities. In particular, let

be the densities of susceptible, infectious, and inactive individuals in class $i$ at time $t$. (We have, in other words, divided $S _ i$, $I _ i$, and $A _ i$ by the width $\Delta r = 1/n$.) Since equations (2.1a) and (2.1b) do not depend on equation (2.1c), we now consider the reduced system

$i = 1, ..., n$, where

As the number of risk classes goes to infinity, we expect system (2.3) to approach the integrodifferential system

for suitably defined birth-rate density $\nu (r)$ and transmission kernel $\beta (r, \, s) \,$.

Students of integral equations ^[32,33] often focus on separable kernels. These are kernels that can be written as a finite sum of products of pairs of univariate functions, with

for some $n$. In epidemiology ^[23,24], the term separable is most commonly used for a single product,

in the transmission kernel.

For models with discrete risk classes, the term separable is again used ^[34,24], when the elements of the transmission matrix can be written as the simple product

The transmission matrix $[ \beta _ {ij} ]$ is now just the dot product of a column vector and a row vector. This matrix is of rank one, and it has a one-dimensional range.

In our analysis of system (2.3), we will focus on symmetric and separable transmission matrices with nonnegative elements

In section 5, we will let

so that

In section 6, we will let

for $n \ge 2$ so that

We thus assume well-mixed populations, with encounters proportional to densities, but with transmission coefficients that vary multiplicatively, as products of simple functions of susceptible and of infective risk. Risk levels, in turn, differ due to differences in behavior or prevention (e.g., vaginal versus anal sex or safe versus unprotected sex). Assumption (2.9) greatly simplifies our analyses, but still allows us to capture a broad range of behavioral interactions, ranging from very low- to very high-risk encounters.

3.   Equilibria

We begin our analysis of system (2.3) by finding the system's equilibria. To do so, we set the time derivatives of system (2.3) equal to zero and look for roots of the system

for $i = 1, ..., n$.

For separable transmission matrices, the above system has two equilibria. There is, first and foremost, a disease-free equilibrium at

for $i = 1, ..., n$.

To find our second, endemic equilibrium, we begin by noting that equation (3.1a) implies that

If

as in equation (2.9), then

where the sum

still needs to be determined. We call $a$ the infectivity potential.

To find the densities of infectious individuals at equilibrium, we sum equations (3.1a) and (3.1b),

and obtain

In light of equation (3.5), it now follows that

Inserting equation (3.9) into the right-hand side of equation (3.6) gives us

This equation can be solved, either analytically or numerically, for the infectivity potential $a$.

The product $a \, b _ i$ is the force of infection, at equilibrium, for class $i$. The infectivity potential $a$, in turn, acts like the amplitude for the endemic equilibrium. Indeed, when $a = 0$, the endemic equilibrium with coordinates (3.5) and (3.9) reduces to disease-free equilibrium (3.2). Equation (3.10), in turn, reduces to the equation

We thus expect the endemic equilibrium to pass through the disease-free equilibrium in a transcritical bifurcation as the basic reproduction number, $R _ 0$, passes through one. For $R _ 0 = 1$, $a = 0$, while for $R _ 0 > 1$, $a > 0$, and we have a unique endemic equilibrium.

4.   Stability

We have seen that we have a disease-free equilibrium with coordinates

and an endemic equilibrium with coordinates

for infectivity potentials $a > 0$ satisfying equation (3.10). As usual, here and throughout, $i = 1, ..., n$.

Let us now analyze the stability of these equilibria by introducing small perturbations, $\xi _ i (t)$ and $\eta _ i (t)$, about the coordinates of the equilibria,

into system (2.3). After simplifying, using system (3.1), and linearizing, we find that

For separable kernel (2.9), system (4.4) simplifies to

where the infectivity potential $a$ is given by equations (3.6) and (3.10).

We now look for solutions of the form

After substituting these solutions into system (4.5) and canceling exponentials, we get the linear algebraic system

Let us now define

Our eigenvalue equations now simplify to

Adding these two equations, we also note that

We must infer the $2n$ (allowing for multiplicity) eigenvalues $\lambda$ for each equilibrium from these equations.

4.1.   Disease-free equilibrium

For disease-free equilibrium (4.1), with $a \, = \, 0$, eigenvalue equations (4.9a) and (4.9b) reduce to

We now have several possibilities for our eigenvalues and eigenvectors.

If $b _ i$ is zero, the above equations reduce to

and we get two eigenvalues, $\lambda = - \, \mu$ and $\lambda = - \, (\gamma + \mu)$. The first eigenvalue has an eigenvector with $u _ i$ nonzero and all other components zero. The second eigenvalue has an eigenvector with $v _ i$ nonzero and all other components zero. These eigenvalues recur for each zero $b _ i$.

More commonly, $b _ i$ is positive. In that case, if the $v _ i$ are all zero, the right-hand sides of equations (4.11a) and (4.11b) still vanish, as does the left-hand side of equation (4.11b). This leaves us with

and with the eigenvalue $\lambda = - \, \mu$. The corresponding eigenvector has $u _ i$ nonzero and all other components zero. This eigenvalue recurs for each positive (or zero) $b _ i$. It thus occurs $n$ times.

For positive $b _ i$, we also have eigenvalues $\lambda = - \, (\gamma + \mu)$. These eigenvalues are consistent with system (4.11) if all of the $u _ i$ are zero and if, in addition,

so that $w$ is zero. In general, there are $n-1$ vectors orthogonal to a non-zero vector with $n$ components. Thus, if any of our $b _ i$ are positive, $\lambda = - \, (\gamma + \mu)$ occurs as an eigenvalue $n-1$ times.

For (at least one) positive $b _ i$, we also have one more eigenvalue, corresponding to $w \neq 0$. Now,

Inserting the above expression for $v _ i$ into our equation (4.8) for $w$, we get

It follows that our last eigenvalue is

This last eigenvalue is positive if

in which case the disease-free equilibrium is unstable. If, instead, $R _ 0 < 1$, the disease-free equilibrium is asymptotically stable. The above formula for $R _ 0$, the basic reproduction number is consistent with our earlier definition of $R _ 0$ in equation (3.11).

4.2.   Endemic equilibrium

For endemic equilibrium (4.2), with $a > 0$, eigenvalue equations (4.9a) and (4.9b) now take the form

As with our disease-free equilibrium, we now have several possibilities for our eigenvalues and eigenvectors.

If $b _ i$ is zero, the above equations again reduce to

and we again get the eigenvalues $\lambda = - \, \mu$ and $\lambda = - \, (\gamma + \mu)$. These eigenvalues recur for each zero $b _ i$.

For positive $b _ i$, we again have eigenvalues $\lambda = - \, (\gamma + \mu)$. These eigenvalues are consistent with system (4.19) if all of the $u _ i$ are zero and if $w = 0$. Since there are $n-1$ vectors orthogonal to a nonzero, $n$-component vector, $\lambda = - \, (\gamma + \mu)$ now occurs as an eigenvalue $n-1$ times.

Finally, for some $b _ i$ positive and $w \neq 0$, we have eigenvectors with components

and, in light of equation (4.10),

Inserting the above expression for $v _ i$ into equation (4.8) for $w$, we find that our remaining eigenvalues satisfy the characteristic equation

If all of our $b _ i$ are positive, we get $n+1$ eigenvalues from this equation.

The left-hand side of equation (4.23) is a nasty expression with many singularities in $\lambda$. Rather than hurting us, however, these singularities instead help us locate eigenvalues of characteristic equation (4.23). We will explore and illustrate this phenomenon more fully in the next section, in the context of an example.

5.   Adding low-risk strata

Let us now consider a simple example in which the birth-rate densities are all constant, $\nu _ i = \nu$, and the transmission coefficients take the form

with

for $i = 1, ..., n$. All of these $b _ i$ are positive.

For $n = 1$, there is only one risk class, and $b _ 1 = c$. As we increase $n$ (see Table 1), we, in effect, add low-risk classes. The mean of the $b _ i$,

decreases towards $c/2$ as $n$ goes to infinity.

5.1.   Basic reproduction number

The basic reproduction number, equation (4.18), is now

and the right-most fraction decreases from $1$, at $n = 1$, to $1/3$, as $n$ approaches infinity. For $R _ 0 > 1$, the disease-free equilibrium is unstable and we have a unique endemic equilibrium.

5.2.   Endemic equilibrium

Endemic equilibrium (4.2) now has coordinates

for infectivity potentials $a$ satisfying

In principle, we could specify the number of classes $n$ and the parameters $\nu$, $\gamma$, $\mu$, and $c$ for any particular problem, use equation (5.6) to solve for $a$, and use equation (5.5) to find the coordinates of the endemic equilibrium. Since, however, $a$ and $c$ often occur as a product in equations (5.5) and (5.6), we instead take a simpler and more efficient approach.

Let

We now treat $\tau > 0$ as a parameter and use equations (5.6) and (5.7) to write $c$ and $a$ as the parametric equations

where

These equations may be used to plot the infectivity potential $a$ as a function of $c$, as in Figure 1.

As an alternative, $c$ and $a$ may each be plotted as a functions of the parameter $\tau$, as in Figure 2. This means that you can pick a value of $\tau > 0$ and simply read off the corresponding values of $c$ and $a$.

For a given $\tau$, the coordinates of the endemic equilibrium in equation (5.5) are now just

We may thus plot the equilibrium densities of susceptibles and infectives, by class, for different $\tau$ values, as in Figure 3. Increasing $\tau$ decreases the equilibrium densities of susceptibles and increases the equilibrium densities of infectives.

5.3.   Stability of the endemic equilibrium

We now wish to determine the stability of the endemic equilibrium. This equilibrium has $2n$ eigenvalues and we know, from section 4.2, that $n-1$ of the eigenvalues have the value $\lambda = - \, (\gamma + \mu)$. Using $\nu _ i = \nu$ and equations (4.23) and (5.2), the remaining $n+1$ eigenvalues satisfy the characteristic equation

where the $\omega _ j$ are defined in equation (5.9) and

Dividing both sides of characteristic equation (5.11) by the first fraction on its left-hand side and substituting $c$ from parametric equations (5.8), we see that the remaining $n+1$ eigenvalues satisfy

If all of the eigenvalues have negative real part, the endemic equilibrium is asymptotically stable.

The right-hand side of equation (5.13) does not depend on $\lambda$ and is a positive constant. The function on the left-hand side approaches zero through positive values as $\lambda$ goes to plus infinity, and it approaches zero through negative values as $\lambda$ approaches minus infinity. (Feel free to look ahead to the next subsection and to Figures 4–7 for graphical illustrations.) The left-hand function has an obvious zero at $\lambda = - \, \mu$, to the right of all singularities, and it has singularities (vertical asymptotes) at $\lambda = - \, (\gamma + \mu)$ and at $\lambda = - \, \rho _ i$\, , for $i = 1, ..., n$. That is, it has $n+1$ singularities. (For ease of exposition, we now assume that all of these singularities are distinct.) These $n+1$ singularities help us locate the eigenvalues of equation (5.13).

As we increase $\lambda$, from left to right, between neighboring singularities that are both to the right of the singularity at $\lambda = - \, (\gamma + \mu)$, the function on the left-hand side of equation (5.13) increases (monotonically) from minus infinity to plus infinity. If we instead increase $\lambda$, from left to right, between neighboring singularities that are both to the left of the singularity at $\lambda = - \, (\gamma + \mu)$, the function on the left-hand side of equation (5.13) decreases (monotonically) from plus infinity to minus infinity. If we look for eigenvalues by plotting each side of equation (5.13) as a function of lambda, checking for intersections, these intervals each give us one real eigenvalue.

The situation is more complicated as we move between neighboring singularities, one of which is $\lambda = - \, (\gamma + \mu)$. We have either one or two such intervals depending on whether $\lambda = - \, (\gamma + \mu)$ is an endpoint or an interior singularity.

On one interval, the function on the left-hand side of equation (5.13) increases to plus infinity at both ends of the interval. We call this interval the parabolic interval.

On the second interval, if it occurs, our left-hand function either increases from minus infinity to plus infinity or decreases from plus infinity to minus infinity, depending on whether the second interval is to the right or to the left of the parabolic interval. This second interval thus gives us another real eigenvalue. We thus have $n-1$ real roots between the $n$ singularities $\lambda = - \, \rho _ i$, not counting real roots from the parabolic interval. For two intervals with $\lambda = - \, (\gamma + \mu)$ as an endpoint, which interval has one real eigenvalue and which is the parabolic interval depends on the exact location of the singularities.

Let us now return to the parabolic interval. If the function on the left-hand side of equation (5.13) dips low enough to intersect the right-hand-side constant, we have two real eigenvalues in the parabolic interval. Otherwise, we have two complex conjugate eigenvalues or, possibly, two real eigenvalues that live elsewhere.

It is tempting to think that the real parts of the aforementioned complex eigenvalues always lie within the parabolic interval, but this is incorrect. We can, however, show that the complex eigenvalues always have negative real part.

To do so, we first rewrite characteristic equation (5.13) in polynomial form,

and note that the leading coefficient, the coefficient in front of $\lambda ^ { n+1 } \,$, is

The next coefficient, the one in front of $\lambda ^ n$, is, in turn,

It is a well-known result, from the theory of equations ^[35], that the negative of the ratio of these two coefficients is just the sum of the roots of the polynomial. Thus,

Since we know that $n-1$ of these roots are real and nested between the $n$ singularities $\lambda = - \, \rho _ i$,

it quickly follows that

Thus, if these last two eigenvalues are a complex-conjugate pair, they must have negative real part with

Alternatively, consider the possibility that the last two eigenvalues, $\lambda _ n$ and $\lambda _ { n+1 }$, are real and lie to the right of all of the singularities. The rightmost singularity is, depending on our parameters, either $\lambda = - \, \rho _ 1$ or $\lambda = - \, (\gamma + \mu)$. In either case, the current assumptions imply that $\lambda _ n > - \, \rho _ 1$ and $\lambda _ { n+1 } > - \, \rho _ 1$. In light of inequality (5.19), it now follows that the last two eigenvalues must lie in the region bounded by the inequalities

This is a triangular region that lies entirely within the interior of the third quadrant of the $(\lambda _ n, \, \lambda _ { n+1 })$ plane. It follows that both $\lambda _ n$ and $\lambda _ { n+1 }$ must be negative.

In summary, our endemic equilibrium has $2n$ eigenvalues. At most, two of these eigenvalues are complex. All of the eigenvalues are negative or have negative real part. Our endemic equilibrium is thus asymptotically stable.

5.4.   Graphical illustrations

To illustrate some possibilities and to show how the loci of the eigenvalues of the endemic equilibrium can change, we start with $\mu = 0.01$ per year, $\gamma = 0.03$ per year, and $n = 4$. In Figure 4, for $\tau = 0.3$ per year, we plot both sides of characteristic equation (5.13) as a function of $\lambda$. We see three real eigenvalues in the intervals bounded by the four singularities at $\lambda = - \, \rho _ i$, $i = 1, ..., 4$, and two real eigenvalues in the parabolic interval bounded by the singularities at $\lambda = - \, \rho _ 1$ and $\lambda = - \, (\mu + \gamma) \, = - \, 0.04$ (per year). In addition, we have $n-1 = 3$ eigenvalues (not shown) at $\lambda = - \, (\mu + \gamma)$. Since all of our eigenvalues are negative real numbers, our endemic equilibrium is asymptotically stable.

As we decrease $\tau$ to $0.2$ per year (see Figure 5), we lose the two real eigenvalues in the parabolic interval. A standard numerical eigenvalue program ^[36] applied to system (4.7) reveals complex eigenvalues at $\lambda \approx - \, 0.0532 \pm 0.0179 \, i$. In this example, the real part of the complex eigenvalues lies within the parabolic interval. Our endemic equilibrium is again asymptotically stable.

In Figure 6, we kept all other parameters as in Figures 4 and 5, but set $\tau = 0.08$ per year. As we decreased $\tau$, $n$ of our $n+1$ singularities shifted to the right, and our parabolic interval, now $(- \, 0.05, \, - \, 0.04)$, shifted to the left relative to the other intervals. Our complex eigenvalues, however, lagged behind. Our numerical eigenvalue program indicates that the complex eigenvalues are now at $\lambda \approx - \, 0.0344 \pm 0.0305 \, i$.

After we decrease $\tau$ even further, to $\tau = 0.02$ per year (see Figure 7), the parabolic interval, $(- \, 0.04, \, - \, 0.03)$, lies to the left, but the complex eigenvalues, $\lambda \approx - \, 0.0133 \pm 0.021 \, i$, lie to the right. Even so, inequalities (5.20) guarantees that the real parts of the two complex eigenvalues lie to the left of $\lambda = - \, { \rho _ 1 } / 2 = - \, 0.0075$.

Finally, we may plot solutions of system (2.3), for birth-rate densities $\nu _ i = \nu$ and transmission coefficients (5.1) and (5.2), by computing these solutions numerically. Figure 8 shows, as an example, the densities of susceptible and infectious individuals for system (2.3) for $\mu = 0.01$ per year, $\gamma = 0.03$ per year, $\tau = 0.2$ per year, and $n = 100$ risk classes. The solutions were computed using a standard, variable-step, Runge-Kutta algorithm ^[37]. The solutions do eventually approach the asymptotically stable endemic equilibrium. At the same time, the eigenvalues for the endemic equilibrium may not tell the whole story, since our solutions quickly overshoot the endemic equilibrium before slowly returning to the equilibrium. We will return to this topic in the discussion.

6.   Filling intermediate strata

Let us now briefly consider a second example. We again assume that the birth-rate densities are constant, $\nu _ i = \nu$, but we now let the transmission coefficients take the form

with

for $n \ge 2$ and $i = 1, ..., n$. The first of these $b _ i$ is always zero. All of the other $b _ i$ are positive.

For $n = 2$, there are two risk classes, with $b _ 1 = 0$ and $b _ 2 = c$. As we increase $n$ (see Table 2), we add intermediate risk classes. The mean of the $b _ i$,

remains constant while the variance of the $b _ i$,

decreases from $3 c ^ 2 / 4$, at $n = 2$, to ${ c ^ 2 } / 4$, as $n$ approaches infinity.

6.1.   Basic reproduction number

If we again use equation (4.18), we see that the basic reproduction number is now

In this instance, the right-most fraction decreases from $1/2$, at $n = 2$, to $1/3$, as $n$ approaches infinity. This reduction occurs despite the fact that the mean of the $b _ i$ remains constant.

6.2.   Endemic equilibrium

Endemic equilibrium (4.2) now has coordinates

for infectivity potentials $a$ satisfying

We now proceed much as we did in section 5. We let

treat $\tau > 0$ as a parameter, and use equations (6.7) and (6.8) to write $c$ and $a$ as the parametric equations

with

For a given $\tau$, the coordinates of the endemic equilibrium are now just

6.3.   Stability of the endemic equilibrium

We now wish to determine the stability of the endemic equilibrium for $\tau > 0$. Since $b _ 1 = 0$, we know, from section 4.2, that $\lambda = - \, \mu$ is an eigenvalue. Since the remaining $b _ i$ are positive, we also know that $n-1$ of the eigenvalues have the value $\lambda = - \, (\gamma + \mu)$.

For the remaining $n$ eigenvalues, we proceed much as we did in section 5. Using $\nu _ i = \nu$ and equations (4.23), (6.2), and (6.9a), we quickly determine that the remaining $n$ eigenvalues satisfy the characteristic equation

with $\omega$ and $\omega _ j$ as in (6.10) and

Because $b _ 1 = 0$, the first term of each series drops out. So, counting really starts with $j = 2$.

The function on the left-hand side of the above equation has singularities (vertical asymptotes) at $\lambda = - \, (\gamma + \mu)$ and at $\lambda = - \, \rho _ i$ for $i = 2, ..., n$. These $n$ vertical asymptotes constrain and help determine the location of the $n$ eigenvalues of characteristic equation (6.11), much as in our previous example. All of the eigenvalues are negative or have negative real part, and the endemic equilibrium is asymptotically stable.

7.   Discussion

In this paper, we used a structured infectious-disease model with demography, mass-action incidence, and an arbitrary number of risk classes to investigate the effects of risk on the spread of HIV. Using a simple, separable transmission matrix, we obtained equations for the disease-free equilibrium, the basic reproduction number, the endemic equilibrium, and the eigenvalues that determine the stability of the endemic equilibrium.

Separable transmission matrices can capture a broad range of behavioral interactions. For real populations, the actual distribution of interactions will also depend on the initial distributions of risk amongst susceptible and infective individuals. (Interactions may also, of course, depend on other details, such as the exact mixing pattern ^[38] and sexual network ^[39], that we did not consider.) Scientists ^[40] have recently made tremendous progress in classifying the risk-structure of real populations. Our model thus holds the promise that it may help us assess the effects of risk for real populations in an analytically tractable way.

We complemented our general analysis with two specific examples. In our first example, we stratified risk so that the mean of the separable components of the transmission coefficients decreased as we added more risk classes. We found a simple, closed-form expression for the basic reproduction number, and this number decreased as we added more risk classes. Thus, initiating HIV-prevention strategies that add low-risk categories lowers average risk and reduces the basic reproduction number, as expected.

We also determined the nontrivial endemic equilibrium and characterized the spectrum of eigenvalues at this equilibrium. For $n$ risk classes, the endemic equilibrium had $2n$ eigenvalues. No more than two of the eigenvalues are complex and all eigenvalues are negative or have negative real part. Thus, when the endemic equilibrium exists, it is asymptotically stable.

Despite the asymptotic stability of our endemic equilibrium, the number of infectives can greatly overshoot endemic levels, as seen in Figure 8. This is not surprising since many mass-action disease models have density thresholds, above which infectives increase and below which infectives decrease ^[41,42,43]. In other words, our model burns through many susceptibles before equilibrating.

It is also important to remember, moreover, that the short-term behavior of a perturbation from an equilibrium may differ from its long-term behavior. Transient amplification can, for example, occur near asymptotically stable equilibria in reactive systems governed by nonnormal matrices ^[44,45,46]. In this paper, we focused solely on the stability of our endemic equilibria. Preliminary analyses (not shown) suggest that our endemic equilibria may be reactive, despite their asymptotic stability. We hope to look at the reactivity of our endemic equilibria in future work.

In our second example, we briefly looked at a model in which the mean of the separable components of the transmission coefficients remained constant but the variance decreased as we added intermediate risk classes. Here, the basic reproduction number again decreased as we added more risk classes.

This second example highlights the value of getting the risk structure right, since different risk structures with comparable mean properties can lead to different predictions. The seemingly unimportant decision of the number of risk strata can influence the dynamic behavior of our system, even if every effort has been made to keep the overall risk identical. Choosing the wrong number of strata can even lead to qualitatively different behavior if the basic reproductive number falls below one. There has recently been great progress in analyzing data-driven structured-population models in other fields ^[47], and we hope to apply some of these new methods to our own model. In a future study, we plan to look at real risk data to explore how feasible the above scenarios are for realistic HIV epidemics.

In formulating our models, we used mass-action rather than standard incidence and a separable, well-mixed transmission matrix rather than assortative mixing. We made these assumptions to increase analytic tractability. We hope to loosen these assumptions in future work. We note, however, that our use of mass-action incidence may, in fact, make our model useful for the study of other diseases, such as coronavirus (COVID-19).

We also introduced but did not analyze integrodifferential system (2.5). Our analysis of discrete system (2.3), in the limit as the number of risk classes goes to infinity, suggests that the characteristic equation at the endemic equilibrium of system (2.5) will have a continuum of singularities and that it would be hard to understand system (2.5) without first understanding discrete system (2.3). At the same time, many investigators find it easier to estimate parameters and to assess goodness of fit for models with continuous traits ^[47]. We hope to analyze integrodifferential system (2.5) more carefully in the future.

Acknowledgments

We thank Shane D. Wilson for useful discussions that helped motivate this study. We also wish to thank the anonymous reviewers for their hard work and helpful comments.

Conflict of interest

All authors declare that there are no conflicts of interest regarding the publication of this paper.