Abstract
In this paper, we construct a linear differential system in both continuous time and discrete time to model HIV transmission on the population level. The main question is the determination of parameters based on the posterior information obtained from statistical analysis of the HIV population. We call these parameters dynamic constants in the sense that these constants determine the behavior of the system in various models. There is a long history of using linear or nonlinear dynamic systems to study the HIV population dynamics or other infectious diseases. Nevertheless, the question of determining the dynamic constants in the system has not received much attention. In this paper, we take some initial steps to bridge such a gap. We study the dynamic constants that appear in the linear differential system model in both continuous and discrete time. Our computations are mostly carried out in Matlab.
1. Introduction
Patients infected with Human Immunodeficiency Virus (HIV) are very likely to develop Acquired Immunodeficiency Disease Syndrome- (AIDS-) related diseases that are usually fatal if not treated with effective antiretroviral therapies. Since the discovery of HIV in 1983, an efficacious vaccine is yet to be developed to fight the deadly virus. Although highly active antiretroviral therapies (HAART) invented in mid-1990s have saved millions of lives and deterred the disease progression of those infected, HIV infection remains a public health threat. Reducing the risk of HIV transmission is of top priority.
One particular challenge in HIV prevention is its long period of latency period. The average time of an HIV infected patient to become symptomatic with AIDS-related diseases can be more than 10 years [1]. In the sexual transmission of HIV, many of the HIV infected patients may not be aware of their HIV infection status, and the virus continues spreading to their HIV negative partners. Therefore, an in-depth understanding of HIV transmission is the key to successful HIV prevention.
HIV dynamics have long been studied in the field of mathematical epidemiology using linear and nonlinear models [2, 3]. The classic model in epidemiology is the SIR model, which considers the dynamics of the susceptible, infected, and recovered populations [4]. This model is not useful for HIV dynamics, as there is no recovered population. An extension of this is the SEIR model, which includes the population of individuals who are exposed but not yet infected. The period between exposure and infectiousness in HIV lasts about two to four weeks [1]. Since a recovered population does not exist, we can consider this period to have a negligible effect on population dynamics.
Hierarchical models are common in HIV modeling due to the high correlation between risky behavior and HIV incidence [5]. In this paper we will incorporate risk indirectly by considering diagnosed and undiagnosed populations. Intuitively, diagnosed individuals would modify their behavior relative to their behavior prior to the diagnosis.
In this paper, we shall form two models: a continuous time linear differential model and a discrete time differential model. These models are the most fundamental among their kinds. The focus will then be given to determination of the parameter estimates, the dynamic constants in these models. As we will show in this paper, the estimates of the dynamic constants depend on the type of model as well as the qualitative properties of the models.
There are two important dynamic constants in our model, namely, the transmission rates for diagnosed HIV population and for undiagnosed HIV population. One important finding in our study is that the transmission rates for the diagnosed and undiagnosed infected populations are comparable. This leads to our conclusion that the transmission rates should be attached to different groups of susceptibles based on their risk level.
2. General Nonlinear Differential Model
One of the frequently used mathematical models for HIV population dynamics can be described as follows. Let be the susceptibles. We divide the HIV positive population into two groups: is the populations that are unaware of the infection; is the populations that are aware of the infection. Let be the mortality rate for the group . Let be the growth rate of the susceptibles. Let be the transmission rate of group and let be the transmission rate of group. Then we have the following nonlinear differential equations: Here counts for those who are infected by group (per unit time), and among them is the proportion of those who are aware of their infection. The constant denotes the rate of the HIV positive population in group who become aware of their infection (per unit time). So there is a flow of from group to once a member from finds out his/her infection through HIV testing.
Many variations of this nonlinear dynamic model have been considered and appeared in the literature to study the HIV population dynamics. For example, in [6], mortality rate of the susceptibles is considered and appears in the differential equation of . In addition, the parameters are allowed to change but are piecewise constant.
In our differential equation model, we have a few constants: , , , , , and . These constants essentially determine the qualitative and quantitative properties of the mathematical model. We shall call these constants the dynamic constants of the model. Notice that some of the constants, like , may have prior estimates, based on the data collected directly from the groups and . Some of the constants, like , will have posterior estimates. The constants , may have prior estimates. Our main focus here is to give posterior estimates of these constants.
We shall remark here that the dynamic constants are model-dependent. This might not be obvious. Even though many of them can be estimated statistically without reference to any models, applying these estimates directly to the model may be problematic, as we shall see in the next section. In this paper, we take some initial steps to estimate the model-based dynamic constants.
3. Linear Differential Model and Preliminary Discussions
We shall now build a simpler linear model. The main assumption is that the susceptible population is a lot larger than and . The change of susceptible population, due to HIV infection, is quite small, comparing with the overall size of susceptible. Therefore, we may ignore the dynamics of susceptible population, by assuming that the susceptible population is a constant. This more or less justifies the use of linear system only involving and .
Let us start with the HIV transmission rate estimates by Pinkerton [7]. The estimates of transmission rates are are estimated in terms of infection transmitted per person per year. Since the overall susceptible population is a lot larger that , we can assume that HIV transmission events are proportional to the size of and . Based on this hypothesis, we may model HIV transmission by linear differential equations: The dynamic constants and remain unchanged. It is also known that [8]. There is no statistics done on . So we can assume as well.
Next, we shall apply the known estimates and study our linear differential model. Notice that remain unknown at this moment. According to [8], is somewhere around . We may tentatively set . Utilizing the estimates of dynamic constants directly from [7, 8], let us consider several cases.
3.1.
We start by assuming that takes the value of the overall portion of those who are aware of their infection. Now we have the following linear equations: We found that the two linear independent solutions have growth rate of However, we know that the growth of is about . Hence our assumption is not valid. Even if we ignore the mortality rate, we have This is still far below the estimated growth rate.
3.2. or
One extreme is that , meaning that the population infected by gets tested and becomes aware of their infection (within the first year). We have Under this assumption will decrease at the rate of , which means that the population will gradually vanish in a few years. This cannot be true.
Another extreme is that , meaning that the population infected by will be initially unaware of their infection (within the first year). We have We have The overall HIV population growth will be less than . This is quite small comparing with the estimate that the growth rate is about .
3.3. , Not Fixed
One might conclude that must be a much smaller number than , what we have initially assumed. We let and be unfixed. In this case, we have We have the matrix We know the growth rates are controlled by the eigenvalues of . In particular, we might assume that with . This will guarantee that the dominant term of the solution will grow at the rate of (per year). Hence we obtain Simplifying it, we have Since , we find that . This suggests that there are between to of getting tested. This percentage seems to be too low comparing with the CDC estimate of about .
We shall remark that our discussion is based on the estimates that and [7]. As we have seen, directly using these estimates as dynamic constants in differential equation modeling will be inadequate to produce the right kind of outcomes and trend. In this paper, we shall discuss posterior estimate of parameters and hope to find some remedy.
4. Posterior Estimate of Parameters
In our earlier discussion, we directly insert the transmission rates from the statistical analysis into the linear differential system. The result is not satisfactory. It is desirable to estimate the transmission rates that will produce the right kind of outcome from the linear differential system model. Let us recall the CDC data from 2007 to 2013 (in thousands) [8].
We first simplify our notation. Let . We rewrite our linear system as where The general solution to this system is Here , are eigenvalues of . They can be both real or complex. There is also a degenerate case that we do not treat here. The behavior of the linear differential system is quite different in these two cases. It is not surprising that we need to use two different methods to estimate the matrix .
4.1. , Real: Simple Curve Fitting
We try a global optimization curve fitting using Matlab. We have Let . Then Notice that the dominant term suggested the overall rate of growth of HIV infected population grows at the rate close to . This seems to be reasonable. But , the rate of flow of population from to , is estimated at . This is completely off the mark. One remedy is that we first estimate the dominant term and then estimate the remainder.
4.2. , Real: Dominant Term Estimate
Suppose that . Then is the dominant term. We shall have Now Using curve to fit this data, we obtain
4.3. Dominant, Real
Now we can assume and use curve fitting to find , , and . We have and , . It follows that We derive that and . These parameters seem to be reasonable. However, . Hence will be a negative number which is not possible.
4.4. Complex with Fixed Real Part
Suppose that and are complex. Then and are conjugate to each other. In particular, the real part of , should be approximately . Write . We should have where is sometimes called a phase constant. A simple curve fitting shows that and and . Hence Let us see what this tells us. We have This roughly says that there are about of that become aware of their infection every year. The annual transmission rate for is . The annual transmission rate for is .
4.5. Complex
We finally use Matlab global optimization to fit the data in the curve We obtain , , Hence we obtain the estimate Now we have It follows that So is neglectable and is about . This again makes the model invalid.
4.6. Discussion
In this section, we choose dynamics constants to fit the temporal data. We have found that these dynamic constants depend on the qualitative properties of the model. Yet, none of the dynamic constants we choose match perfectly with the existing estimates. One reason is that yearly data is not suitable for a continuous time model. Therefore, we shall explore the discrete time model.
5. Discrete Dynamic Model
We may regard () as a discrete time dynamical system. Let us assume that this discrete dynamics is defined by a transition matrix : In principle, based on our earlier discussion, Now we would like to estimate .
5.1. Basic Estimates
The easiest way to find is by considering the following matrix equations: For example, for , we will have Then we find the following estimate of : We can see some consistency among these transition matrices. For example, the -th entry has been around . This translates into This is the rate of transmission for group . It seems to be consistent with the estimate of [7].
5.2. (Arithmetic) Average Estimate of
Now we may average all ’s and obtain Hence Our estimate yields that ; in other words, about of those unaware of their infection will become aware of their infection next year. We also have If the mortality rate is set to be , then we have . Similarly, we have If the mortality rate is set to be , then we have This suggests that the transmission rate of group is twice as large as the transmission rate of group. There may be some truth to it. However, we believe that this estimate is off the mark due to the reason that are correlated with each other. Hence each estimate will be biased. We shall correct this and give a more robust estimate later.
5.3. Least Square Estimate of
Perhaps a good way to estimate is the least square method. We write Applying the least square method, we find that the least square solution to is This estimate seems to be better than the arithmetic average, in the sense that, irregularities will have smaller effect on the least square solution. Because we can reorder ’s and the least square solution does not change, we also avoid the pitfall that and are correlated. We have our posterior estimates: This estimate is similar to the arithmetic average we just computed. The dynamic constant estimates will be very similar. We shall then look for a solution that is more robust. One particular reason that the least square estimate is not satisfactory is that there are additional relations like that least square method does not take into consideration. In other words, , can also be estimated and shall be taken into consideration when we estimate . We shall offer one remedy that avoids this issue.
5.4. A More Robust Estimate
One of the problems with our estimate is that and are correlated to each other. As a remedy, we pick and as far from each other as possible. We observe that
We compute Then We have
5.5. Discussion
This estimate of is quite consistent with the least square estimate. Our estimate seems to suggest that the transmission rates for and may be in the similar range. By [8], assume that . We have Every year about of those unaware of their HIV positiveness become aware of their infection due to testing. About of those infected by group become aware of their infection. This seems to be consistent with some of the observations in [7], with one exception; namely, in our estimates, the transmission rates for and are very close. Figures 1 and 2 show the difference between and .
5.6. Arithmetic Average versus Geometric Average
Now we may state our problem in greater generality. Given a temporal vector , suppose that with transitional matrix . How should one estimate the matrix ?
As we discussed earlier, we can use least squares with the equations The least square estimate of , in some sense, is very similar to the arithmetic mean of the transitional matrix . But what makes better sense is a geometric mean. More precisely, we have to take into consideration that Suppose that is the transitional matrix at time . Then a good estimate of should be the “geometric average” of . For scalars, one can define the geometric average of to be the th root of . But matrix multiplications are not commutative and one cannot define the geometric average of matrices. It remains a challenging problem to define computationally a geometric mean of .
5.7. Roots Estimate
Tentatively, we may define the geometric mean by taking roots. For example, we may now consider Then
We can also considerWe have
Both estimates of are consistent with the least square estimate and the estimates in the previous section. Above all, all our estimates point to the same range of transmission rates for both and .
6. Concluding Remarks
Now we shall compare our dynamic constant estimates in the linear differential model in continuous time and in discrete time.
In the continuous time model, we obtain the transmission rate of about for the group, those who were aware of their HIV infection. Nevertheless, the rate of flow from to due to HIV testing turned out to be too low and often came out to be negative, which cannot be the case. The best results are obtained when we assume the two eigenvalues are complex. In this case Yet seems to be quite low. According to the CDC report [8], is estimated at about .
There are various reasons why our dynamic constants are inconsistent with known estimates. Firstly, the CDC data we use tends to underestimate the and population sizes, particularly in more recent years. The CDC estimates the sizes of the populations infected with HIV by back calculation. The estimates for any given year will increase as new diagnoses are obtained. HIV may go undiagnosed for up to 10 years without causing the death of the patient (Table 1) [5]. Depending on the stage of the disease, the individual will be counted as undiagnosed for a number of years prior to the diagnosis. This causes the estimates of the population sizes to be smaller than the actual size of the population. A new estimate of the HIV prevalence agrees with this conclusion [9]. Although this new estimate is more conservative than the back calculation method, it may still underestimate the population. Both estimates show a downward trend in the data, but this is likely to be erased as more individuals are diagnosed in the later stages of the disease.
Secondly, our computations assume that the susceptible population is much larger than the infected population. However, failing to obtain the right estimates suggests that opposite may be true—the susceptible population could be much smaller. HIV infection is overrepresented in some subpopulations, such as men who have sex with men (MSM) [1]. This subpopulation is only about 5% of the US population, or approximately 15 million individuals. Not all MSM can be considered to have equal risk of contracting the disease, and since over 1 million individuals are currently living with the disease, the susceptible population may be comparable to the size of the infected population. For this reason, any differential system model of the HIV transmission must include the susceptible population.
The discrete dynamics model seems to be most robust against the bias caused by back calculation and the size of the susceptibles. By ignoring the 2013 year’s data, we are able to determine the transmission rates as Also , not too different from the CDC estimate . We see that the dynamic constants in the discrete time model are less affected by the underestimation caused by back calculation. It is also true that the relative size of susceptibles has less effect on the discrete time model than on the continuous time model.
Finally, our estimates of the transmission rates for diagnosed and undiagnosed HIV population are relatively close. This is very different from the previous estimate, where the transmission rate of the undiagnosed population is about 4 times as high as the diagnosed population [7]. This implies that the transmission rates should be attached to the susceptible population. It makes sense to divide the susceptible population into groups depending on the possible transmission rates for them. We shall pursue this in a different paper.
Conflicts of Interest
The authors declare no conflicts of interest regarding the publication of this paper.
Acknowledgments
This research is partially supported by the NIH R01 Grant AI12125903.