AN AGE-STRUCTURED MODEL OF HIV INFECTION THAT ALLOWS FOR VARIATIONS IN THE PRODUCTION RATE OF VIRAL PARTICLES AND THE DEATH RATE OF PRODUCTIVELY INFECTED CELLS

Mathematical models of HIV-1 infection can help interpret drug treatment experiments and improve our understanding of the interplay between HIV-1 and the immune system. We develop and analyze an agestructured model of HIV-1 infection that allows for variations in the death rate of productively infected T cells and the production rate of viral particles as a function of the length of time a T cell has been infected. We show that this model is a generalization of the standard differential equation and of delay models previously used to describe HIV-1 infection, and provides a means for exploring fundamental issues of viral production and death. We show that the model has uninfected and infected steady states, linked by a transcritical bifurcation. We perform a local stability analysis of the nontrivial equilibrium solution and provide a general stability condition for models with age structure. We then use numerical methods to study solutions of our model focusing on the analysis of primary HIV infection. We show that the time to reach peak viral levels in the blood depends not only on initial conditions but also on the way in which viral production ramps up. If viral production ramps up slowly, we find that the time to peak viral load is delayed compared to results obtained using the standard (constant viral production) model of HIV infection. We find that data on viral load changing over time is insufficient to identify the functions specifying the dependence of the viral production rate or infected cell death rate on infected cell age. These functions must be determined through new quantitative experiments. 2000 Mathematics Subject Classification. 92B05, 35L99, 45D05.

1. Introduction.We develop an age-structured model for HIV-1 infection dynamics that tracks the length of time a T cell has been infected.Age-structured models have been widely used to study the epidemiology of infectious diseases, such as HIV [1], hepatitis C [2] and tuberculosis [3,4].These works explore issues of global stability analyses [1,2] and numerical analyses [4] with the focus on the analysis.The model by Castillo-Chavez and Feng did explore the host dynamics of TB [3].However, age-structured models have not been widely considered for within-host dynamics for HIV; the only case to date dealing with the context of treatment strategies is [5].We pursue a different application, were we focus on the analysis of the age-structured equations via the Jacobian matrix as well as focus on the application of these models to the study of HIV patient data.The work by Zheng et al. [4] utilizes the Jacobian matrix to calculate the characteristic equation and we follow a similar line in our analysis.Age-structure allows more realistic representations of the biology of HIV-1 infection.In particular, it allows us to account for the fact that the production of new virus particles (virions) by an infected cell does not occur at a constant rate, but rather ramps up as viral proteins and unspliced viral RNA accumulate within the cytoplasm of an infected cell.Our model also allows the rate of death of an infected cell to vary according to the time a cell has been infected.Special cases of this model are the standard ordinary differential equation model of HIV-infection [6,7,8], delay models that account for the fact that production of new HIV virions from an infected cell cannot occur instantaneously after infection [9,10,11,12,13], and models that assume cell death due to HIV infection may be delayed [14].The age-structured model developed here has now been used to examine optimal viral production rates [15] and optimal viral fitness [16] under certain simplifying assumptions.
2. Age-structured model of HIV-1 infection.Our model considers a population of uninfected target cells, T (t), infected T cells structured by the age, a, of their infection, T ⋆ (t, a), and virus, V (t).The equations defining the model are dT dt = s − dT (t) − kV (t)T (t), In this model target cells, T , are assumed to be produced at a constant rate, s, and die at a fixed rate, d, per cell.The infection of target cells is assumed to occur via the law of mass-action as a second-order kinetic process between uninfected cells, T , and virus particles, V , with an interaction-infection rate constant, k.The virion production rate, P (a), and the death rate, δ(a), of infected cells of age a, T ⋆ (a, t), are assumed to be functions of the age of cellular infection, a. Virions, V , produced by infected cells are assumed to be cleared at a fixed rate per virion, c.We also assume da dt = 1; that is, that the time unit for age of infection is the same as that for clock time.
Since this model contains a first-order hyperbolic equation, we need to introduce the appropriate boundary and initial conditions.First, infected T cells of age zero are created by infection; that is, Because infected cells of any age class have units of 1/age or equivalently 1/time, both the left and right-hand sides of (2) have the same units, 1/time.
To study primary infection we impose the initial conditions T (0) = T 0 , T ⋆ (a, 0) = 0 and V (0) = V 0 , where T 0 is the level of target cells prior to infection and V 0 is the initial virus concentration.If we wish to study the effects of drug or monoclonal antibody therapy on patients with established infections, we assume the patients are initially at steady state and choose T (0) = T ss , T ⋆ (a, 0) = f (a), and V (0) = V ss , where T ss and V ss are the steady state levels of target cells and virus, respectively, and f (a) is the steady state distribution of infected T cells of different ages.Under these conditions, time t no longer represents the time since the infection of the host, but instead represents time since the administration of therapy.
With the above boundary and initial conditions, we note that by standard methods, it is possible to prove existence and uniqueness of the solutions for (1) (see [17,18]) and to show that the solutions all remain bounded and non-negative for t > 0.

2.1.
Varying viral production, P (a).The functional form of the viral production kernel, P (a), is unknown and remains to be determined experimentally [19].We consider two possible kernels that capture features of the biology.Both have a maximum production rate, P max , because cellular resources will ultimately limit how rapidly virions can be produced, but they differ in how they approach the maximum.First, we study a delayed exponential function where β controls how rapidly the saturation level, P max , is reached.We have also included a term a 1 to represent a delay in viral production; that is, it takes time a 1 after initial infection for the first viral particles to be produced.This kernel can mimic either a very rapid increase to maximal production or a slow increase to maximal production depending on the value of β (Fig. 1a).
Second, we consider a Hill type function, which also allows for saturation of viral production, where K a is the half-saturation level and n is a constant called the Hill coefficient.This function allows for quick growth to a maximal level, depending upon the value of K a , but can also approximate the delayed effect seen in ( 3) for large values of n (see Fig. 1b).Hence, both functions can account for the fact that immediately after infection a number of biological processes need to occur before the first virus particle is produced.The Hill function does this without an explicit delay, although one can be added (Fig. 1c).
To contrast the current approach with previous ones, note that in the earliest models of within-host HIV dynamics [8,20] the viral production rate was chosen as a constant; that is, P (a) = constant.Later, delay models were introduced that chose P (a) = 0 for a ≤ a 1 and P (a) = constant for a > a 1 , where a 1 was either fixed or given by a probability distribution [11,12,13].Thus, the current model is a generalization of both the fixed production rate and delay models.Recent work [15,16] examines the possibility that the viral production rate may evolve to maximize viral production.Also, it is possible that δ is a function of the virion production rate P (a) or of total virion production to date.These ideas are considered elsewhere [15,16].

2.2.
Varying the death rate, δ(a).Another feature of our model is that the death rate, δ(a), of infected cells can vary with the age of an infected cell.The correct form for this age-dependent death rate distribution is unknown but experiments that examine changes in nef and rev gene expression [21,22] allow conjecture about its behavior.Once a target cell becomes infected, it takes some time before viral epitopes begin to appear on the cell surface, leaving initial cell death to be dominated by background, or natural death.Further, for a strong CTL (cytotoxic T cell) response to occur, large numbers of epitopes must be expressed on the cell surface [23].Thus, an infected cell is expected to become susceptible to CTL-mediated killing only at some as yet unknown age of infection.Hence, for δ(a) we will consider an increasing function, where the rate of cell death increases with the age of infection (Fig. 2).
In addition, we assume a minimal infection time is needed before epitopes are expressed and the cell becomes susceptible to cell-mediated killing or for enough viral products to be made that the cell may die from the infection itself; that is, from viral cytopathic effects.One possible choice is where δ 0 + δ m is the maximal death rate, γ controls the time to saturation and a 2 is the delay between infection and the onset of cell-mediated killing.The term δ 0 represents a background death rate.given by ( 5) for γ = 0.1, 0.5, and 1.0.The background death rate δ 0 = 0.5 day −1 .The additional death rate due to both cytotoxic T cells and viral cytopathic effects, reaches a maximum of δ m = 1 day −1 , at a rate determined by the parameter.
Changes in γ can mimic killing of infected cells by CD8 + T cells because a higher γ could relate to a more rapid expression of viral peptides on the cell surface and hence create a better recognition by CD8 + T cells.The model can also account for a weak CD8 + T cell response by allowing γ to remain very low, so that most of the early killing would be due to viral cytopathic effects and background death.However, as the model stands, we cannot isolate the killing effects by CD8 + T cells or other methods of infected cell loss, and so we only consider a combination of all effects.
2.3.Burst size.The total number of viral particles produced over the lifespan of an infected cell is called the burst size, N .For cells with age-dependent viral production and death rates we have where σ(a) = e − a 0 δ(s)ds (7) is the probability an infected cell survives to age a.
3. Model Analysis.To determine the steady states of the age-structured model we set the time derivatives in (1) equal to zero and solve with It is easily seen that one steady state is the trivial or noninfected steady state There is also a non-trivial or infected steady state, which we will denote with overbars.Solving the second equation in (8) with initial condition (9), we get Equations ( 8) also imply Substituting ( 11) and ( 12) into ( 13) and rearranging yields and substituting ( 14) into (12), we obtain with Lemma 1.The infected steady state exists if and only if or for the case when P (a) is given by (3) and δ(a or written another way, highlighting the infectivity rate k, if and only if If condition (ii) is violated, then the only steady state is the noninfected steady state.
Proof of (i) At the infected steady state The existence of a biologically realistic infected state, requires V to be greater then zero.Hence 16) is positive.
Proof of (ii) Calculating N = ∞ 0 P (a)σ(a)da for P (a) given by ( 3), and δ = δ(a) = constant, gives which when substituted into ( 17) yields (18).In addition, When the inequality ( 18) is violated, V < 0, since Hence from the proof of (i) we see this to be violated and the only non-negative steady state is the noninfected steady state.Equation ( 24) can easily be re-written to get (19).QED The condition for the existence of an infected state state can also be seen by direct examination of the model equations.In the standard ordinary differential equation model of HIV-1 infection presented in [6,7,24] a transcritical bifurcation occurs at c = πkT 0 /δ, where π is the constant virion production rate and T 0 = s/d.The viral production rate can also be written as π = N δ, where N is the burst size.Thus, if viral clearance is faster then viral production; that is, c > πkT 0 /δ, or N kT 0 , then the noninfected steady state is stable [7].Alternatively, if c < N kT 0 , then the infected steady state is stable [7].
In our model with age-structure, we find the same bifurcation.At steady state, viral clearance equals viral production, or from (15) which is equivalent to c = N k T .If we use (3) for P (a), we get If we use T = T 0 = s/d, so that the infected and uninfected steady states merge, and rearrange (26) we get which is the bifurcation point corresponding to (19) (see Fig. 3).
The bifurcation point at N = dc sk is equivalent to the classical formulation involving R 0 , the basic reproductive number.Nowak and May [6] show that for the standard model without age-structure, R 0 = kπs δdc .With π = N δ, R 0 = kN s dc ; thus R 0 > 1, the criterion for a stable infected steady state (17), is equivalent to N > dc sk .
If we consider the Hill function ( 4) for the production kernel, it is possible to determine a bifurcation analytically for the case n = 1 and δ(a) = δ.In this case we can integrate (6) to get where s ds.From ( 17) the bifurcation occurs at c = N ks d .If we set K a = 1 and δ = 0.4, day −1 for example, we find the bifurcation occurs at c = 1.45P max ks d (see Fig. 3).For other values for the Hill coefficient, n, we would have to resort to numerical methods to determine the bifurcation point.Note from Fig. 3 that we find almost identical results using the two different production kernels, (3) and (4).

Stability results.
The following Lemma shows how the stability of a steady state of a model with age-structure can be determined directly from the application of the Jacobian matrix about the steady state [4].Jacobian matrices are used to determine the characteristic equations for finite dimensional systems, but we will show that they can be used for our age-structured model after we apply the method of characteristics.This process will require a Laplace transform of each age-dependent term.
We present this generalization of the standard Jacobian matrix method for determining stability in full detail here.This derivation using the Laplace transform is essentially the same as presented in earlier examinations of age-structured models by Pruess [26], Thieme and Castillo-Chavez [1], Webb [18], and Hethcote [27].However, we believe that expressing the criteria in matrix terms even for this simple model is instructive and worthwhile, since the framework will extend directly to stability analysis of more complex age-structured models., a max = 10 days, and T (0) = 10/µl and V (0) = 0.02 /ml, as in [25].Other parameters come from Stafford et al. [24].With these parameters, the bifurcation value for c is 11.27 day −1 (see ( 26)) (a) and c = 12.54 day −1 (b).For c = 10 day −1 , the infected state is stable and V ss = 3569/ml (a) and 7177 (b).The bottom picture shows the transcritical bifurcation diagram for the case when (3) is used for the production kernal.The bifurcation diagram is similar for the case where (4) is used.
Proof of Lemma 2: Consider the model (1) with the defined boundary (2) and initial conditions.First, note that the general solution for T ⋆ can be immediately determined using characteristics as where T ⋆ 0 (a − t) = T ⋆ (a, 0), is the initial value and B(t − a) is an unknown function that is determined by the boundary condition T ⋆ (0, t) = kV (t)T (t).From ( 29) with a = 0, we find B(t) = kV (t)T (t).Hence the trajectories of T ⋆ are controlled by the value of δ and the concentration of V since V is given in the boundary condition.If V = 0 (that is, a noninfected steady state), then the only valid solution of ( 29) is T ⋆ ≡ 0. Only for V > 0 can we find a nonzero solution of T ⋆ and hence an infected steady state.Substituting ( 29) into (1), we obtain where da and since we assume P (a) ≤ P max < ∞, Equation (30) has two differential equations and one algebraic equation because the dynamics of T ⋆ are controlled by the boundary condition.Once the stability of ( 30) is found, it then gives the stability of (1).
We can write (30) as where K 1 (a) = P (a)σ(a) and is the convolution, in a, of the functions K 1 and B. The solutions of s − dT ss − kV ss T ss = 0, give us the points about which to linearize.Linearizing about any given fixed point Taking the Laplace transform of (35) we get, denoting transforms with hats, where λ is the Laplace variable.We begin the solution of these equations by rewriting the double integral.Using Fubini's theorem we switch the order of integration and set α = t − a to obtain Since B = 0 for α < 0 (that is, t < a) the lower bound of integration of the interior integral can be set to zero without loss of generality and we get The resulting interior integral ∞ 0 B(α)e −λα dα is the Laplace transform of B, i.e., B(λ).Substituting this result back in gives (39) Equations (36) now form a simple linear system.We can solve the equations using Cramer's rule to find for some functions f T , f V and f B , and using the Jacobian matrix for (32) with a Laplace transform for the production kernel, The absence of the λ in the (3, 3) element of the matrix is a consequence of the equation for B(t) being algebraic (not a differential equation) [27].On taking the inverse Laplace transform, the growth rates λ are found at the poles of the solutions (40).After checking that there are no common factors in f T , f V or f B , these poles are at the zeroes of det(A), that is, the solutions of We now determine the characteristic equation for our model using the Jacobian matrix and show that it is the same as (42).Evaluating the Jacobian matrix at the noninfected steady state [(T ss , V ss , B ss ) = (s/d, 0, 0)] gives det with the characteristic equation Note that for V ss = 0 and T ss = s/d, (44) is the same as (42).The roots of this equation are λ 1 = −d, and solutions of where if we use and set σ(a) = e −δa we find the roots to be solutions of Rearranging, we get where we immediately see that one condition for stability, provided by the Routh-Hurwitz condition, is that the constant term of the polynomial must be positive; that is, cδβ + cδ 2 − kT ss P max (β + δ(1 − e βa1 )) > 0, which is satisfied by Lemma 1.
Since the coefficient of λ 2 is also positive, the final condition for stability is We can rewrite (49) to get Notice the two terms with the negatives in front.The first with the exponential is always positive as 1−e βa1 ≤ 0 and the second, cδ(β+δ), is included in the product of the first positive term and hence will cancel out.Hence, the noninfected steady state is locally stable, provided (24) is satisfied (that is, if c − P max (β+δ−δe βa 1 )ks δ(β+δ)d > 0), and therefore, the infected steady state does not exist.
To study the infected steady state, when it exists, we linearize (30) Solving gives If we now consider the case where P (a) is defined by (3), δ(a) = δ and a 1 = 0 we have where The conditions for stability are obtained by employing the Routh-Hurwitz conditions for a fourth order polynomial, as a 1 > 0, a 4 > 0, B 1 ≡ a 1 a 2 − a 3 > 0, and By inspection we see that a 1 > 0 for all parameter values.Lemma 1 guarantees that a 4 > 0, provided the infected steady state exists.Using Maple, one can multiply out B 1 and see that all the negative terms cancel out and hence B 1 > 0 for all parameters.The final condition for stability is C 1 > 0 and again using Maple one can find that all the negative terms can be matched with a positive term that is bigger in magnitude, given Lemma 1.Hence, the infected steady state is stable provided the condition in Lemma 1 is satisfied; that is, or in our case with a 1 = 0, Numerical results.The continuous age-structured model, (1), was solved numerically by using MatLab 6.0 (Mathworks, Natick MA, USA).We first converted the partial differential equation in (1) to a series of coupled ordinary differential equations by discretizing the age classes of T ⋆ into an array of equally sized age classes between 0 and a max .Here a max does not necessarily correspond to the maximum age of infection of a cell, but instead corresponds to the age at which the production rate has closely approached the asymptotic value P max and, therefore, is essentially constant with age.Thus, the age class a max includes all cells whose age of infection was greater than or equal to a max and was chosen such that the viral production rate in the last class satisfies P (a max ) = (1 − 10 −6 )P max .The flux of individuals between age classes was calculated using a fourth-order finite difference method [28], except near the boundary for the penultimate age class a max −1 .At this upper age boundary we calculated the flux from age class a max −1 into a max using a first order upwind finite difference method.Using this method near the upper age boundary damped any instabilities introduced by defining the T ⋆ density in age class a max as ∞ amax T ⋆ (a, t)da.We verified that the relative error for this approach was less than 10 −5 for special cases where the equilibrium densities are known analytically.The set of coupled ordinary differential equations and the corresponding equations for T (t) and V (t) were numerically solved using a variable order Adams-Bashforth-Moulton method with a variable time step to maintain the temporal error below 10 −6 .4.1.Primary infection.As a test case to study the effects of age-structure we examine the kinetics of primary HIV infection.After an individual is infected with HIV-1, the viral load measured in plasma typically rises to a peak within a few weeks after infection and then declines, reaching a quasi-steady state or set-point value.The time to the peak, the amplitude of the peak, and the set-point viral load are all important characteristics of primary infection that any model needs to match.Data on the kinetics of primary infection for ten patients are given in Stafford et al. [24], as well as fits and parameter estimates based on the standard ordinary differential equation model of HIV infection.This corresponds to the non-age-structured version of model (1).We reasoned that if the virion production rate varied with infected cell age then the time to reach the peak viral load and the subsequent dynamics might be affected.Thus, we compared predictions of our model to those in Stafford et al. [24].The initial conditions we used were T (0) = T 0 , T ⋆ (0, a) = 0 and V (0) = V 0 , where T 0 and V 0 were chosen to be the same as in [24].
In Stafford et al. [24], the rate of viral production, denoted π in that paper, was held constant and independent of infected cell age.In their model the total number of virions produced per infected cell during its lifetime, the burst size N , is given by π/δ, where 1/δ is the average infected cell lifetime.To ensure that the effects we see with an age-structured model are not due to choosing a different burst size, we held N constant in all of our comparisons between models.In Fig. 4 we compare the results in Stafford et al. to simulations with the age-structured model for two representative patients.For the age-structured model we choose P (a) given by (3), with δ(a) = δ, and a 1 = 0, so that N = Pmaxβ δ(δ+β) .Notice, in Fig. 4, that if viral production ramps up slowly to P max then the time to peak viral load is delayed compared to the standard model with constant production.For both patients examined, the delay is about ten days when the characteristic time for the viral production to increase is one day; that is, β = 1 day −1 .However, when β is increased to 10 day −1 so that viral production ramp up rapidly (see Fig. 1a) and P max is decreased so as to keep N constant, a closer fit to the dynamics seen in Stafford et al. [24] is obtained.
If we use P (a) given by (4), we obtain similar results to those for P (a) given by (3).However, in this case the closest fit to the solutions in Stafford et al. is obtained when we increase P max and decrease n (Fig. 5).In this limit the production rate is approximately constant (Fig. 1). Figure 4. Numerical solutions of (1) using the parameters from Stafford et al. [24], with P (a) given by ( 3), with a 1 = 0 and P max and β constrained so that N = Pmaxβ δ(δ+β) had the value estimated in Stafford [24].Top panel (Patient 1): The value of s = 0.13 cells/(µl-day) was chosen from the steady state conditions, given that d = 0.013 day −1 , k = .46x10−6 ml/(virus-day), δ(a) = δ = 0.4 day −1 , a max = 10 days and T (0) = 10/µl, and V (0) = 10 −6 /ml.The solid line corresponds to the case studied in Stafford et al. with constant virion production, π = N δ.The other two lines show results for the age dependent model.Notice, when β = 1/day and hence ramping up to full viral production is slow, the time to reach the peak viral load is delayed.However, for β = 10/day the solution approximates the constant production case and as β → ∞ the solutions coexist (not shown).Bottom panel (Patient 8): We used s = 0.085 cells/(µl-day) to get the steady state conditions given by d = 0.0085 day −1 , k = .66x10−6 ml/(virus-day), δ(a) = δ = 0.17 day −1 , a max = 10 days and T (0) = 10 µl and V (0) = 10 −6 /ml.
Numerical simulations of the model with a varying infected cell death rate, where δ(a) is given by ( 5), also show that the time to the viral peak and the overall level of virus at times before the steady state is established are sensitive to the choice of δ(a) (Fig. 6).
If we allow the production rate to vary with age of infection we find that different production schedules and different burst sizes can give similar viral dynamics (Fig. 7), suggesting that it may be difficult or impossible to deduce P (a) from data  4), K a = 1 day and the other Hill parameters listed in the graph.The chosen values for P max and n were constrained to maintain consistency with the estimated value of N in Stafford et al. [24].The solid line is the result given in [24].The remaining parameters are given in Fig. 4(top).3), with a 1 = 0 and δ(a) given by ( 5), with a 2 = 0. We compared the results of allowing for a constant δ, δ(a) = 0.4 (solid line) to those with a death rate that varied with age of cellular infection, with δ 0 = 0.05 day −1 and δ m = 0.35 day −1 .We examined cases where there was a slow growth to maximal killing (γ = 0.2 day −1 ) and a faster growth to maximal killing (γ = 1 day −1 ).The remaining parameters were P max = 1019 virions/cell-day, β = 10 day −1 , s = 0.13 cells/day, d = 0.013 day −1 , k = .46x10−6 ml/(virion-day), a max = 15 days, T (0) = 10/µl, and V (0) = 0.02/ml.on virus concentration changes with time.However, when we consider these production rate changes coupled with a varying death rate, the virus dynamics curves that were nearly identical become less similar; for example the time to achieve the viral peak changes (Fig. 7).

4.2.
Sensitivity of model to initial values.One of the difficulties of estimating parameters that characterize primary infection, with this model as well as with non-age-structured models [29,25,24], is the model sensitivity to initial data.For later stages of infection this sensitivity is less severe.There is no way to measure the amount of virus that initially infects a patient, and in general the exact time of infection is unknown.Also, the initial density of target cells is generally unknown, especially if one considers target cells to be the activated fraction of CD4 + T cells, as is done in the Stafford et al. model [24].Hence, we studied changes in both T 0 and V 0 to determine how they affect the time to peak viral load (Fig. 8).
We found that for low initial values of target cells (that is, T 0 = 10 cells/µl the value used in Stafford et al. [24]), changes in the initial viral load led to a greater variation in the time to the viral peak than if we chose a higher value of T (0); say, T 0 = 100 cells/µl (Fig. 9).
More importantly, we found that it is difficult to distinguish between the effects of changes in initial values and the effects of changes in P (a) (Fig. 10).This sensitivity to changes in initial values reduces the significance of parameter fitting to determine P (a).It also implies that one can not use the time to peak viral load to distinguish between different functional forms of P (a).Instead, direct experimental measurements of P (a) are needed.

5.
Conclusions.We developed a model of HIV-1 infection that accounts for variations in the death rate of productively infected T cells and viral production that are due to the age of the cellular infection.Our model considers time, t ≥ 0 to be the time since a person has become infected with HIV-1 and a time or age of infection, a, which keeps track of the time between a T cell becoming infected and its ultimate death.The age of cellular infection plays a key role in determining the rate of viral particle production per productively infected T cell and how long the productively infected T cell lives.This is the first account of the behavior of HIV models with viral production and death rates that vary as functions of age of infection.
We provided a detailed local stability analysis of the fixed points of this system and found as in the non-age structured models that there were two steady states, an uninfected one and an infected steady state, with the infected steady state arising by a transcritical bifurcation.
We considered two separate functions for the production of viral particles as function of the age of cellular infection.In our model the average burst size, N , is given by ∞ 0 P (a)σ(a)da.Not surprisingly, our model was able to simulate dramatically different scenarios for viral production while providing similar overall total viral loads.However, we did find in certain scenarios that variations in the burst size and production schedule (see Fig. 7) did not affect the timing or concentration of viral peaks, while in other scenarios it did (see Fig. 4).
When we examined the effects of a varying death rate, δ(a), we found that predicted viral loads, and in particular the time to viral peak and the amplitude of the damped oscillations that occur before reaching steady state, were sensitive to the choice of δ(a).However, since other effects such as the initial viral load or the initial number of target cells also affect the time to peak, it would be difficult to identify the age-dependent death rate from an examination of viral load data alone.While we currently do not have any experimental justification for the chosen functions, P (a) and δ(a), our work lays the foundation for future studies of HIV pathogenesis that we hope will be, in part, directed to the determination of these functions.In other publications [15,16], we have asked whether natural selection in the host can determine P (a) and have posed the question of finding the P (a) that maximized the burst size under conditions where the death rate of infected cells, δ(a), may be influenced by the viral production rate P (a).The model presented here allows consideration both of primary infection and of the establishment of chronic infection.An important application will be the consideration of the effects of drug therapy on chronically infected patients.Much work has already been done in this area using either ordinary differential equation or delay differential equation models [8,11,30,31,32].Whether the use of age-structured models will change the interpretations of drug-perturbation experiments remains to be determined.

1 Figure 1 .
Figure 1.Virion production kernel, P (a), showing the ability to represent rapid production or delayed production depending on the chosen parameters.Here P max = 850 virions/day.(a): P (a) given by (3) with a 1 = 0, except where noted.(b): P (a) given by (4) with K a = 1.(c): The two functions mimic one another for certain parameters.

1 Figure 2 .
Figure 2. Death rate of productively infected T cells with δ(a)given by (5) for γ = 0.1, 0.5, and 1.0.The background death rate δ 0 = 0.5 day −1 .The additional death rate due to both cytotoxic T cells and viral cytopathic effects, reaches a maximum of δ m = 1 day −1 , at a rate determined by the parameter.

Figure 5 .
Figure 5. Numerical solutions of (1) using parameters characteristic of patient 1 in Stafford et al. with P (a) given by (4), K a = 1 day and the other Hill parameters listed in the graph.The chosen values for P max and n were constrained to maintain consistency with the estimated value of N in Stafford et al.[24].The solid line is the result given in[24].The remaining parameters are given in Fig.4(top).

Figure 7 .
Figure 7. Numerical solutions of (1) with P (a) given by (3), with a 1 = 0 days, δ(a) = constant = .40,and N = 3133 when β = 0.8 and P max = 1880, or N = 2176 when β = 5 and P max = 940 (top figure) or δ(a) given by (5) with δ 0 = 0.05 day −1 , δ m = .35day −1 , γ = 0.5, and a 2 = 0.The top figure shows that similar timing and concentrations at viral peaks can be obtained with greatly different maximum production rates, burst sizes, and rates of gearing up production γ.The bottom figure shows that with the timing of viral production as in the top panel but an age varying death rate that the viral concentration profiles no longer are similar.All other parameters are the same as in Fig.6.