Using a Collateral Damage Model to Explain Survival Data in West Nile Virus Infections

Simulation code for a model of the adaptive immune response seen in flavivirus infections is used to explain the immunopathological consequences seen in West Nile Virus virus (WNV) infections. We use a model that specifically handles the differences in how the virus infects resting cells, the G0 state, versus dividing cells, the G1 state, which includes vastly increased MHC-I upregulation for resting cells over dividing cells. The simulation suggests how the infection progresses in a one host model and the results shed insight into the unusual survival curve data obtained for this infection: there is an increase in health even though viral load has increased.


Introduction
The code for the general simulation framework has been discussed in [1] with results in [2], and the reader is referred to those works for the background details.Further, in that work, we outline how a flavivirus such as West Nile Virus infects its host, and so we will only go over some of that information here.
The family Flaviviridae are single-stranded, plus sense RNA viruses transmitted by arthropods such as mosquitoes and ticks.Flaviviruses induce an increase in expression of major histocompatibility complex, MHC-I and II, as well as other immune recognition molecules, including intracellular adhesion molecule-1 (ICAM-1, CD54), vascular adhesion molecule-1 (VCAM-1, CD106) and E-selectin (CD62E) variously and on a wide variety of cells [3].WNV-induced increases in cell surface expression of these molecules therefore result in increased efficiency of recognition and killing of infected cells by WNV-specific cytotoxic T lymphocytes, CTL [4] [5].Increases in MHC-I concentration on infected cells enhance the avidity of interaction of virus-specific CTL with the infected target cell.The increased avidity enables the interaction of infected cells with CTL clones previously below the recognition threshold because of their low affinity for MHC-virus peptide ligand.Furthermore, a polyclonal anti-viral CTL population can include self-reactive low-affinity clones that recognize MHC-I-self peptide configurations as discussed in [5]- [7].The previous simulation results in [2] offered support for this point of view.Due to the increased avidity of their interaction, these low-affinity, self-reactive clones can now lyse both infected and uninfected target cells that express high cell surface MHC-I concentrations [8].
Because the avidity of interaction between T cells and their targets can be increased via MHC associated upregulation, any increases in MHC and/or ICAM-1 expression further enhance CTL recognition of target cells.This in turn increases the chances of low-affinity, self-reactive CTL recognizing these target cells.It is also known that the position of the cell cycle is important in flavivirus infection.Non-dividing cells when infected with WNV increase their MHC-I expression 6-10-fold, while infected cycling cells increase MHC-I expression by only 2-3-fold [9].This results in some 10-fold more lysis of infected target cells than cycling cells by the same CTL [4].Thus, there is enhanced avidity of interaction between CTL and WNV-infected cells in the nondividing state while infected cycling cells are less easily recognized.Further, WNV replicates significantly better in dividing cells than in non-dividing cells.Thus, virus may be eradicated with relatively poor efficiency in a population of infected cycling cells.In vivo, since most cells are not dividing, the decoy hypothesis asserts it is the small population of infected dividing cells that supports virus replication while maintaining a low immunological profile, thus increasing the probability of virus transmission [5].High MHC-I expressing (uninfected) targets are therefore susceptible to lysis by low affinity self-reactive CTL clones.A review summarizing the case for an autoimmune component to the West Nile virus infections is in [10].
In most virus infections, the probability of death of the host varies predictably with viral dose in a dose response curve.Thus, a lower dose results in viral clearance by an efficient immune response that runs its complete course in an infection with no mortality.In contrast, the highest doses produce 100% mortality, which interrupts the immune response before it has managed to clear the virus.In flavivirus infection, however, the mortality associated with intervening doses is unpredictable over several log10 dilutions of virus dose.This results in a ragged dose-response curve depicting population mortality [7] whose data we show in Figure 1.The work in this paper builds on the model of the host-pathogen interactions for flaviviruses detailed in [2] and is used to explain the ragged dose survival curves seen in these infections.

Simple Survival Model Calculations
We can do a simplified calculation to help us see that the presence of MHC-I upregulation should lead to local oscillations in the survival rate.We will do another proof of concept calculation in Section 3 once we have built a model of collateral damage which also suggests that there is a potential for auto immune interaction with the WNV infections.Let H denote the population of uninfected dividing and non-dividing cells and I , the population of infected non-dividing cells (this is the population of interest).The size of these populations is dependent on the viral dose in pfu/cell (plaque-forming units), ρ and the population of non dividing cells that are infected, respectively.Let r denote the overall growth rate of the uninfected dividing and non-dividing cell population H .We assume the rate of misrecognition (i.e., death of an uninfected cell) due to MHC-I upregulation λ is given by βλ I for some fraction β .
For our purposes here, the only parameter of interest is the viral dose ρ .We thus have the system ( ) where ( ) ( ) Let's assume 0 ρ is a point where the survival curve oscillates.If you look at Figure 1, you can see several such candidates.
At such a point, approximate the system by the constant term ( ) 0 A ρ .Near 0 ρ , we expect the dynamics of H to be captured nicely by this approximation.This is a linear system of first order differential equations and the general solution is determined by the eigenvalues of the matrix ( ) 0 A ρ .Our data on survival curves for WNV infections (see Section 4 and Figure 1) tells us that there should be oscillation in H versus ρ plot and hence, the characteristic equation for this system should exhibit complex roots with a negative real part.Then we would see a exponentially damped phase shifted cosine solution which would match the data.The characteristic equation is given by This has roots ( ) ( ) (this is easy to achieve as the viral load increases) and if 0 D < .The discriminant function is therefore a complicated function of the parameters α , λ , β , r and 0 ρ , but it is easy to see that under some combinations of these parameters, D is indeed negative giving rise to damped oscillations near a value of 0 ρ .For example, we can show how for a large enough 0 ρ , these oscilla- tions can arise.From [2], we know 0.01 r ≈ .Let's assume we have high upregulation, 9 λ = .We don't expect the misrecognition rate β to be very large, so let's assume 0.02 β = . If we plot the discriminant function versus the parameter 0 αρ for these values, we can see the results in here and there are many choices of 0 ρ leading to oscillation as we have 0 0.01 ρ α = .Thus 0 ρ can vary enormously in value and give oscillations.We know α is the fraction of viral load that is translated into infection and thus, if .In a typical survival plot, we use the logarithm of 0 ρ , hence we could see oscillations near the value of ( ) . This simple model therefore suggests that we can have survival oscillations as we see in the data at a variety of viral loads 0 ρ .On the other hand, if if we decrease the upregulation λ to a smaller value such as 3 λ = , we see from Figure 3 that the range of 0 αρ giving rise to oscillation is substan- tially diminished.Hence, as the amount of upregulation decreases, there is a lesser chance of oscillation.Finally, we note that is upregulation further decreases to 1 λ = , all chance of oscillation is lost as seen in Figure 4.This calculation with a very simplified model of healthy cell versus viral load dynamics therefore suggests it is possible for there to be combinations of these parameters that lead to oscillations in the survival curve even at high viral dose.This implies that upregulation is the key factor leading to collateral damage and the unusual survival data we see.

Healthy Cell Approximations
When the immune system response causes healthy cells to be lysed, this is called collateral damage and in the case of West Nile virus infections, it has been discussed in [2].We can use the damage model outlined there to shed insight in yet another way.We can rewrite the healthy cell update equation as  where at time t, ( ) the presence of virus that hides in the dividing cells and hence is shrouded from immune system notice.Let's try to make this idea more precise.There needs to be some sort of interaction between the new infections and the collateral damage that leads to these kinds of survival oscillations but it is hard to see it from the dynamics above.Taking a second derivative, we see 1  ( )

y t y t p y t p y t y t y t
In Equation 2, N represents how we discretize our simulation universe of possible cells.We use about 10 8 cells in the simulation (coding details are in [1] and the general derivations in [2]) and use a three dimensional box of size 217 on a side to give us 10 8 possible cells.We divide this universe of cells into 3 N larger cells or coarse cubes for computational ease in our simulations with 512 N = and we let i C denote the collateral damage estimate in course cube i.We assume 7 p cells in a given coarse cube can upregulate due to the infection.In the simulation, the number of new infected cells that are not dividing and are removed in cube i is denoted by 0  n ri , where the subscript r indicates removal.Similarly, the number of new infected cells that are dividing and are removed in cube i is denoted by 1  n ri .In the simulation, we also keep track of the number of upregulation signals sent to the th i coarse cube in the variable i S .As the simulation progresses, the number of removed cells due to upregulation is dependent on how many times a cell receives an upregulation signal.The sum of the removals due to these events in stored in the variable i θ ∑ for a given MHC-I upregulation level θ and coarse cube i.Finally, the term i ε is the amount of random perturbation we use around the base values i θ ∑ .
We let the fraction 0 1 ri ri n n N + and the fraction i S N be denoted by 1 ξ and 2 ξ , respectively.These are probably small but for the moment think of 1 2 ξ ξ as simply another parameter Z .Replace ( ) + Θ and ignore the random perturbation term i ε .Thus, we have ( ) Further, ignore the slight modification due to the intercept term b giving Next, assume this is true for all coarse cubes, giving 14 7 512 .
We think it is reasonable to assume that if the number of healthy cells increases, the amount of collateral damage will increase also.We also believe that the amount of collateral damage will depend on the level of the initial viral dose which we will call I. Hence, let's assume the parameter Z is proportional to interaction between the initial dose and the cumulative amount of healthy cells present from time 0 to the current time t, y , should decrease if the proportion of healthy cells increases; we are assuming qualitatively that a healthier system resists infection.Over the course of an infection cycle of time T (for us 1440 15 minute time units or 15 days), the maximum amount of healthy cell impact is 6 p T as 6 p is the resource allocation limit given by the homeostatic logistics growth law used to control cell population in the simulation [2].This is similar to the idea of a life work: i.e. force times length of time the force is active, T, is a work done point of view.We think of this as the amount of life energy available for use by the infectious agent.Hence, the total amount of life energy possible minus the amount actually used currently is an estimate of what is available for infection.This amount is ( )  − − Θ is going to change at each time step t rather than be constant as we have here and hence, the le- vels of 1 y will move up and down depending on the contribution of ( ) ˆβ α − .Our aim with this graph is to show the positive and negative correction terms move the healthy level 1 y up and down which implies the possibility of the balance of new infection and collateral damage to shift so as to allow the healthy level 1 y to increase.Note, we can't really see this possibility of an increase in 1 y from the original first order dynamics given by 1 y′ .It was important to pass to the second order dynamics 1 y′′ which requires making additional assumptions.The survival data tells us as I Θ increases, sometimes the term ( ) becomes negative, allowing 1 y to increase.The exact dependence of this correction term on Θ and I is, of course, quite complicated and difficult to tease out.We can get a little more insight by modifying our 5 y and 14 y assumptions a bit.As before, we still think it is reasonable to assume that if the number of healthy cells increase, the amount of collateral damage will increase also.However, let's now assume the amount of collateral damage will depend on both the level of the initial viral dose which we will call I and the upregulation level Θ .We now assume the parameter Z is proportional to interaction between the initial dose of virus, I, the upregulation level, Θ to some positive power p and the cumulative amount of healthy cells present from time 0 to the current time t.
Again, modeling the interaction as a product, we have ( ) ( ) Thus, an approximate dynamical model for the healthy cell level 1 y is given by the nonlinear second order model ( ) β α Θ = .This lends further credence to the idea that there is a critical value of I for which we can increase 1 y implying a return to health even though the initial viral dose is large.Also, note it is the presence of the upregulation that plays a major role here.
So we can see that various assumptions about the dependencies of the new infection levels and collateral damage levels on Θ and I lead to predictions of varying 1 y levels that can lead to the survival data we see.Determining these functional dependencies should be amenable to experimental protocols.However, here we will show that by constructing our simulation very carefully, we can see this behavior.Hence, the remainder of our discussion is to buttress our claims that our simulation results provide some validation for the decoy hypothesis.

Survival Models
In a survival experiment, a population of mice of size N is all infected with the same amount of virus, and the number of animals that die is determined.The amount of virus used is measured in plaque forming units (pfu).This determines the concentration of infectious virus by virtue of the number of areas of cell death in a cell monolayer in vitro.The experiment is then repeated, using N mice per group for a range of virus concentrations.If we did this for various pfu levels, we could then graph the number of mice that survive versus the virus concentration.Such a survival experiment is expensive to undertake, requiring many mice and hence, rarely done.With most viruses, a traditional survival experiment in which the virus is titrated in this way, gives a classical dose-response curve which progressively and smoothly decays down to a survival of 0 at high virus concentrations.WNV has a peculiar survival curve which was shown in Figure 1.Our simulations show similar results, although we use an artificial parameter in our simulations for the pfu level and in Figure 6.In this simulation, we modeled host infections in 18 groups using initial viral concentrations ranging from 100 to 3.6 06 e + pfu.The raw results are shown in Table 1.In the table, you will note that there is an increase in all the measures of survival (Normalized Health, Survived, Survived 2 and Survived 3) where each measure attempts to estimate  as it is a sort of clipping of the actual plot.We can also calculate the area of the plot that lies below this threshold and call it a health window and denote it by plot W .We then set another notion of damage to be the normalized value ( ) . If this value is greater than 0.65, we count this plot as a survival and update this measure of survival.
We show normalized health versus the logarithm of the infection level in Figure 6.In Figure 7, we plot all   three different survival measures versus the logarithm of the virus concentration but only for the range of viral dose where there is an increase in survival.We see the characteristic increase in survival with increasing viral load in all three measures as well as in the normalized health plot.We see in both simulation plots, survival sometimes increases with increasing virus concentration, in contrast to the data seen in most other viral infections.

Conclusion
We show the decoy hypothesis is a reasonable way to explain collateral damage and the existing survival data measured in WNV infections.In addition to the simulation results, we have also included arguments based on approximations that also show we can expect some sorts of oscillatory behavior in the level of healthy cells versus viral load.In Section 2, we let H denote the population of uninfected dividing and non-dividing cells and I , the population of infected non-dividing cells and assumed nonlinear interactions between them mediated by the viral load, the variable ρ , and the MHC upregulation level, the variable λ .We showed there was reason to believe this dynamic model allowed for oscillations in the health level, such as we saw in the survival data.Then in Section 3, we use the simulation dynamics from [2] to develop an approximation to collateral damage in the model.This allowed us to develop a second order model for the size of the uninfected cell population, i.e., the healthy cells ( ) y t given by ( ) from which we could also reason that oscillations in the survival rate could be expected.We then presented simulation results that supported these suggested behaviors.A more general model, [11], can also be developed where we assume a large population of cells T consisting of cells which are dividing and are infected, D , cells which are not dividing but are still infected, N , non infected cells, H and non infected cells which will be removed due to auto immune action which we call C , for collateral damage.We then assume all of these variables depend on the three parameters, I , the IFN-γ signal emitted by cells being lysed; the MHC-I upregu- lation factor U , here the variable λ , and the free antigen level A which here is the variable ρ .Our work in this paper does not use the IFN-γ signal, but with the IFN-γ signal a more detailed analysis is possible and we can derive analytical equations predicting health and survival.Numerical results then also support the oscillations in survival we see in the data.The models we present here are based on the micro level one presented in [2] but the new model in [11] essentially uses a macro level approach.These WNV studies then suggest a more general model of autoimmune response which uses both micro and macro level reasoning [12].We show there are self damage effects whenever the triggering agent's effect on the host can be separated into two distinct classes of cell populations such as we see in the WNV infections where the two classes are the infected cells that are dividing and the infected cells that are not.As long the trigger acts differently in each population and this behavior is mediated by nonlinear interactions between two signaling agents which satisfy certain critical assumptions, there is collateral damage.All of these results posit the existence of collateral damage and it is also important to develop a model of how self damage occurs.This has been done in work that uses a new model of the immuno synapse and cytokine/chemokine signalling grammars to understand how a self damage can occur [13].

Figure 2 .
The survival model discriminant function can be negative leading to oscillations: us that the discriminant is negative in the range [ ] 0, 75 or so.Thus, for example,

Figure 2 .Figure 3 .
Figure 2. The survival model discriminant function can be negative leading to oscillations: 0.18 βλ =

Figure 4 .
Figure 4.The survival model discriminant function does not allow oscillations: 0.18 βλ = y t is the uninfected cell population, ( ) 5 y t is the population of newly infected cells and ( ) 14 y t is the population of uninfected cells lyzed by T cells.The parameters 5p and 6 p are the growth factor for uninfected cells and the resource limit for uninfected cells.Equation 1 suggests that 1 y satisfies the differential equation know there is collateral damage that depends on the level of upregulation and the initial dose of virus and that it leads to unusual survival curves.We suspect this is due to

∫
and the initial viral dose I. Interactions are often modeled as products; hence, we assume this is true for all coarse cubes, giving the final model We let β denote the proportionality constant.The number of new infections, 5 Θ .Thus, an approximate dynamical model for the healthy cell level 1 y is given by the nonlinear second order model we subtract from 1 y .This term should be close to a balance, hence, we will solve two problems with positive and negative values of 3.0 5 e − .When − , we call this a negative correction and the graph in Figure5shows 1 y increases past the usual resource allocation limit, 6 p .Of the other hand, if ( ) call this a positive correction and the graph shows the 1 y level decreases below the usual resource allocation limit.Now the term
sign of the correction is determined by the function ( ) ( ) Hence, for this simple extension of our modeling of 5 y and 14 y , we see there is a change in the sign of the reinforcement term around ˆp r I

Figure 7 .
Figure 7.The Percentage of hosts surviving.
This is not necessarily true, but we want to show how assuming the dependence is not linear allows for 1 y behavior that is consistent with the measured survival data.Again, let this new proportionality constant be denoted by α .Combining, we have

Table 1 .
Sample raw simulation results.This is the area under the health curve.Since the plot may terminate early, we also calculate the area , we count this as a survival for this level of virus dose.Since we have said a number of plots per viral dose, this allows us to count how many hosts at this viral dose survived.Survived 2: We can also calculate another measure of survival.If the lowest health value, which is measuring how long the health stays above the mini- mum level over the full time of the simulation.If the normalized smallest health level, H HPT exceeds some threshold, here 0.25, we count this plot as a survival.Survived 3: The third way we can use to estimate survival is to calculate the time interval that the value 2 Y lies below a threshold of P, here 0.9P.We can compute this time and call it plot c t