In this work, we develop a theoretical model that explains the survival data in West Nile Virus infection.
Results
We build a model based on three cell populations in an infected host; the collateral damage cells, the infected dividing cell, and the infected non-dividing cells. T cell-mediated lysis of each of these populations is dependent on the level of MHC-1 upregulation, which is different in the two infected cell populations, interferon-gamma and free virus levels.
Conclusions
The model allows us to plot a measure of host health versus time for a range of initial viral doses and from that infer the dependence of minimal health versus viral dose. This inferred functional relationship between the minimal host health and viral dose is very similar to the data that has been collected for WNV survival curves under experimental conditions.
Viruses in the family Flaviviridae are single-stranded, plus-sense RNA viruses, principally transmitted by mosquitoes and ticks. They are either viscerotropic, causing diseases such as dengue and yellow fever, or neurotropic, causing central nervous system disease, like West Nile virus encephalitis, Japanese encephalitis or Saint Louis encephalitis, all of which may be fatal. These viruses are found worldwide and the diseases they cause pose a significant public health burden. A careful discussion of the many facets of flavivirus infections can be found in [1], and we therefore confine our introductory remarks to the features salient to this paper.
Virus-infected cells process virus proteins into peptide fragments in the proteasome and these are bound to class I major histocompatibility complex molecules (MHC-I) in the endoplasmic reticulum. With the transport of MHC molecules to the cell surface, viral peptides displayed in the context of MHC on the infected cell can be recognized by virus-specific cytotoxic T lymphocytes (CTL). Such CTL kill infected cells by a variety of lytic effector molecules prior to the release of mature virus progeny from infected cells, thereby progressively reducing virus numbers and ultimately eradicating virus from surviving hosts (reviewed in [2]).
Paradoxically, flaviviruses such as West Nile (WNV) and others, directly induce increased expression of MHC-I, as well as MHC-II, and several adhesion molecules involved in immune recognition by CTL. This increased expression results in a marked increase in the efficiency of recognition and killing of infected cells by WNV-specific CTL ([3, 4]), because although the affinity of individual T cell receptor (TcR)-MHC-virus peptide interaction is unchanged, the multiple intermolecular interactions increases the avidity of interaction of virus-specific CTL with the infected cell. This increased avidity also enables the functional interaction with MHC hi, infected cells by CTL clones of low MHC-virus peptide affinity, i.e., clones previously below the recognition threshold. Some of these low-affinity CTL clones are likely to be self-reactive [4] or even able to recognize MHC without peptide specificity [5]. Thus, the increased avidity brought about by high MHC expression enables low-affinity, self-reactive clones, not normally involved in anti-viral immune responses, to lyse both infected and uninfected target cells [6]. Adhesion molecules such as ICAM-1, as non-specific accessory molecules upregulated by WNV infection, also increase the avidity of CTL-target cell interactions, further lowering the affinity threshold for T cell recognition and target cell lysis [7]. In addition, interferon- γ (IFN- γ), released by CTL on recognition of their cognate ligand, strongly increases MHC and ICAM-1 expression on neighboring target cells, further contributing to the progressive increase in avidity of interaction between CTL and target cells [8]. In this context, the stage of the cell cycle in which a cell is infected is also important; cells infected by WNV in G0 (resting) increase MHC-I expression by 6-10-fold, compared to a 2-3-fold increase observed in cells infected during the cell cycle (G1, S, G2+M) [3]. WNV-infected G0 cells are approximately 10-fold more susceptible to CTL lysis than infected cycling cells exposed to the same CTL [3]. Thus, WNV-infected cycling cells are less easily recognized, while the avidity of interaction between CTL and infected G0 cells is significantly enhanced by the higher levels of MHC-I and ICAM-1. Notably, WNV replicates significantly better in the poorly recognized cycling cells than in G0 cells. In vivo, most cells are in G0, presumably presenting an easy CTL target once infected, but a small population of productively infected cycling cells maintaining a low immunological profile, could substantially increase the probability of virus transmission to the next host [4]. As indicated above, the release of IFN- γ upon target cell recognition by CTL would also increase MHC and ICAM-1 expression on uninfected cells in the vicinity of virus- infected cells, making high MHC-I-expressing (uninfected) cells a potential target for lysis by low affinity virus-specific and/or self-reactive (i.e., cross-reactive) CTL clones. The collateral destruction of uninfected cells by low-affinity clones of this kind would cause substantial additional damage to the brain, with corresponding increases in morbidity and mortality [9]. We have developed a model of the collateral damage caused by a West Nile Virus infection, which is supported by simulation results [10]. A discussion of the underlying simulation code can be found in [11]. We have also used these simulation results to explain the unusual ragged survival data seen in West Nile Virus infections [12]. The previous work, in focusing on very low level details, based on first principles, could be regarded as a micro model. Here we develop a macro model, based on a much higher level approach and ideas. Each model has its pros and cons but we believe each brings complementary insights into a very complicated problem. Our focus here is therefore on developing a macro level theoretical model of how the survival of a host depends on the level of initial viral dose. We believe this approach provides an abstract focus which can also be directed towards more general models of immunopathology and we will briefly mention those connections at the end. These more general models of autoimmune response are the focus of additional work we are doing.
In a survival experiment, a population of animals such as mice are infected with a common amount of virus and the number of animals that die are counted. This experiment is repeated for a range of virus. With most viruses, the survival experiment gives a classical dose-response curve which progressively and smoothly decays down to a survival of 0 animals at high virus. However, WNV has a peculiar survival curve as shown in [12] but we can easily simulate what such a curve looks like as we have done in Fig. 1. The purpose of this paper is to find a theoretical model that explains this data. The derivations here use standard ideas from advanced calculus and differential equations. The necessary background can be found in any nonlinear modeling text such as [13]. We assume we have a large population of cells T which consists of cells which are dividing and are infected, D, cells which are not dividing but are 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 assume all of these variables depend on the three parameters, I, the IFN- γ signal; the MHCI upregulation factor U and the free virus level A. The approach we take here is a common one in such a nonlinear interaction environment. An example of how this is used to develop a model of diabetes detection can be found in [13].
There are then three nonlinear interaction functions F1, F2 and F3 because we know C, D and N depend on each other’s levels in very complicated ways. Usually, we assume the initial dose S0 gives rise to some fraction of infected cells which will differ in both the dividing and nondividing cells. In our previous WNV simulations of collateral damage [10], we assume p0=1.0e−3 of the initial dose gives rise to infected cells. This gives the number of infected cells to be p0S0 which is split into p1p0S0 infected non dividing cells - i.e. in N, and p2p0S0 dividing infected cells, i.e. in D, where p1+p2=1. Typically, we let p1=0.99 and p2=0.01 but other choices could be made. Thus, the total amount of virus that goes into infected cells is p0S0 and the amount of free virus is therefore (1−p0) S0. Thus, we could expect C0=0, D0=p2p0S0 and N0=D0=p1p0S0. However, we will explicitly assume we are starting from a point of equilibrium prior to the administration of the viral dose S0. We could assume there is always some level of collateral damage, C0 in a host, but we will not do that. We will therefore assume C, D and N have achieved these values C0=0, D0=0 and N0=0 right before the moment of infection. Hence, we don’t expect there to be initial contribution to C′(0), D′(0) and N′(0); i.e. F1(C0,D0,N0)=0, F2(C0,D0,N0)=0 and F3(C0,D0,N0)=0. We are interested in the deviation of C, D and N from their optimal values C0, D0 and N0, so let c=C−C0, d=D−D0 and n=N−N0. We can then write C=C0+c, D=D0+d and N=N0+n The model can then be rewritten as
Next, we do a standard tangent plane approximation on the nonlinear dynamics functions F1, F2 and F3 to derive approximation dynamics.
Linearization details
To develop a useful approximation to the CDN dynamics, we use an appropriate tangent plane approximation on the nonlinear dynamics functions F1, F2 and F3. The mathematics behind this approximation come from multivariate calculus and can easily be reviewed if required. The standard tangent plane expansions is as follows.
It seems reasonable to assume that since we are so close to ordinary operating conditions, the errors \(E_{1}^{F}\), \(E_{2}^{F}\) and \(E_{3}^{F}\) will be negligible. Also, to save space, we let ()o denote that we are evaluating the partial derivatives at the point (C0,D0,D0). Thus our model approximation is
As before assume C, D and C have achieved the same optimal values C0=0, D0=0 and N0=0 prior to the moment of infection with virus dose S0. These correspond to the starting values of prior to infection for base values I0, U0 and A0. Initially, we don’t expect IFN- γ signals so I0=0. Eventually, we do expect some level of upregulation due to this initial dose and from experimental data, we expect the upregulation to be proportional to the level of the dose S0; we will assume this is a simple scaling factor, i.e. U0=q1S0 for some suitable parameter q1. Also, once the viral dose is administered, we would expect some fraction of it to remain as free virus which as discussed in [10] is modeled as A0=(1−p0)S0. We will deal with these initial values after infection later in the Section “The IFN- γ, upregulation and free virus model”. But now, we think of all the initial values as zero; i.e. I0=0, U0=0 and A0=0. We still don’t expect to have any contribution to C′(0), D′(0) and N′(0); i.e. H1(I0,U0,A0)=0, H2(I0,U0,A0)=0 and H3(I0,U0,A0)=0. We are interested in the deviation of C, D and N from their optimal values C0, D0 and N0 due to the changes i, u and a from the base IFN- γ, upregulation and virus values. So let c=C−C0, d=D−D0 and n=N−N0. We can then write C=C0+c, D=D0+d and N=N0+n The model can then be rewritten as
Next, we again perform a tangent plane approximation on the nonlinear dynamics functions H1, H2 and H3 just as we did in Section “Linearization details” for the F1, F2 and F3 we use for the CDN model.
Linearization details
Once again, it seems reasonable to assume that since we are so close to ordinary operating conditions, the tangent plane errors are negligible. Letting ()o denote that we are evaluating the partial derivatives at the point (I0,U0,A0), we find
where we now use a standard subscript scheme to indicate the partials.
The algebraic signs of the linearization matrix
If we hold everything constant except i which we increase to i+δi, what happens? Increasing the IFN- γ level should increase collateral damage. Hence, \(H_{1\boldsymbol {i}}^{o} = +\). What about the number of dividing cells that are infected? We do not think this should have an effect; hence, \(H_{2\boldsymbol {i}}^{o} = 0\) too. A similar argument suggests \(H_{3\boldsymbol {i}}^{o} = 0\) as well. Thus, the coefficient matrix above which we call Ψ so far looks like
Now hold everything constant except the upregulation u and increase u to u+δu. What happens? Increasing the upregulation level should increase collateral damage. Hence, \(H_{1\boldsymbol {u}}^{o} = +\). Now infected cells are lysed once their IFN- γ upregulation hits a certain level and the increase in upregulation that makes the infected cells more visible certainly effects this as it will take fewer additional IFN- γ signals to trigger lysis. In general, there is a lag of say \(T_{L}^{D}\) and \(T_{L}^{N}\) time steps before this happens. With the increase δu, we expect \(T_{L}^{D}\) and \(T_{L}^{N}\) to both go down. Hence D and N decrease and we have \(H_{2\boldsymbol {u}}^{o} = -\) and \(H_{3\boldsymbol {u}}^{o} = -\) for each type of infected cell. Thus, the coefficient matrix looks like
Next hold everything constant except the virus level a and increase a to a+δa. What happens? Let’s think of the virus increase δa as giving rise to an increase in the amount of virus stored inside a dividing cell or a non dividing cell. Now if the amount of virus in the cell goes up, that means when the cell is lysed, there is more virus available to infect cells which means more cells will be infected in later times. An increase in virus means an increase in collateral damage in general, so \(H_{1\boldsymbol {a}}^{o} = +\). At a given time then, A is the virus level. We can write A=AF+AD+AN where AF is the free virus, AD is the virus inside the dividing cells and AN is the virus inside the nondividing cells. So if AD goes up, we expect the amount of A and AF to remain constant. Hence, if AD goes up, AN goes down.
The increase in AD would mean more virus is released on lysis and hence an increase in D. However, that increase in AD is a concomitant decrease in N. So we should have \(H_{2\boldsymbol {a}}^{o} = +\) and \(H_{3\boldsymbol {a}}^{o} = -\). Thus, the coefficient matrix now looks like
The amount of interferon level I, the upregulation level U and the free virus level A depend on the initial amount of virus applied when in the equilibrium state; i.e. this is the amount that causes the initial infection. This is So. We assume the dynamics here are
In the model of “The CDN IFN- γ, upregulation and virus model” section, we assumed C, D and N depended on the perturbations of I, U and A from a zero state. Now, we want to model the I, U and A deviations from a base state I0, U0 and A0 which is not zero. As previously discussed, we expect A0=(1−p0) S0, the initial IFN- γ level I0=0 and the initial upregulation level U0=q1S0. Let the deviations from these equilibrium values be given by i=I−I0, u=U−U0 and a=A−A0. The model can then be rewritten as
The approximation to the i, u and a model is handled in a way that is very similar to the previous two expansions which were explained in some detail in “Linearization details” section.
Linearization details
The linearization again involves a tangent plane approximation to nonlinear dynamics functions, which are now the functions G1, G2 and G3. It seems reasonable to assume that since we are so close to ordinary operating conditions, the tangent plane errors as usual will be negligible. Thus, using \(()^{\boldsymbol {S_{0}}}\phantom {\dot {i}\!}\) to denote the partials evaluated at the base point the model approximation, we have
The analysis of the signs of these partials is next. This is similar to what we did for the previous model. If we hold everything constant except i which we increase to i+δi, what happens? Increasing the IFN- γ level should increase IFN- γ. Hence, \(G_{1\boldsymbol {i}}^{\boldsymbol {S_{0}}} = +\). What about the upregulation level? We do not think this should have an effect; hence, \(G_{2\boldsymbol {i}}^{\boldsymbol {S_{0}}} = 0\) too. A similar argument suggests \(G_{3\boldsymbol {i}}^{\boldsymbol {S_{0}}} = 0\) as well. Thus, the coefficient matrix above which we call Φ so far looks like
Now hold everything constant except the upregulation u and increase u to u+δu. What happens? Increasing the upregulation level should increase the IFN- γ level. Hence, \(G_{1\boldsymbol {u}}^{\boldsymbol {S_{0}}} = +\). If upregulation increases, the u′ goes up, hence \(G_{2\boldsymbol {u}}^{\boldsymbol {S_{0}}} = +\). What about the level of virus? If upregulation increases, the increase in upregulation makes infected cells more visible as it will take fewer additional IFN- γ signals to trigger lysis. In general, there is a lag of say \(T_{L}^{D}\) and \(T_{L}^{N}\) time steps before this happens. With the increase δu, we expect \(T_{L}^{D}\) and \(T_{L}^{N}\) to both go down. Hence, free virus should increase; i.e. \(G_{3\boldsymbol {u}}^{\boldsymbol {S_{0}}} = +\). Thus, the coefficient matrix looks like
Now hold everything constant except the virus level a and increase a to a+δa. What happens? Let’s think of the virus increase δa as giving rise to an increase in the amount of virus stored inside a dividing cell or a non dividing cell. Now if the amount of virus in the cell goes up, that means when the cell is lysed, there is more virus available to infect cells which means more cells will be infected in later times. An increase in virus should not effect IFN- γ levels, so \(G_{1\boldsymbol {a}}^{\boldsymbol {S_{0}}} = 0\). An increase in virus must imply \(G_{3\boldsymbol {a}}^{\boldsymbol {S_{0}}} = +\) Now reason as we did earlier: at a given time then, A is the free virus level. We can write A=AF+AD+AN where AF is the free virus, AD is the virus inside the dividing cells and AN is the virus inside the nondividing cells. So if AF goes up and the total amount of virus at this time is constant, we expect AD and AN to go down. If the amount of virus inside the cells decreases, we expect the upregulation to decrease. So \(G_{2\boldsymbol {a}}^{\boldsymbol {S_{0}}} = -\). Thus, the coefficient matrix now looks like
Or letting \(G_{2\boldsymbol {u}}^{\boldsymbol {S_{0}}} = a\), \(G_{3\boldsymbol {u}}^{\boldsymbol {S_{0}}} = b\), \(G_{2\boldsymbol {a}}^{\boldsymbol {S_{0}}} = -c\), \(G_{3\boldsymbol {a}}^{\boldsymbol {S_{0}}} = d\), we have the coefficient matrix now looks like
$$\begin{array}{@{}rcl@{}} \boldsymbol{\Phi} &=&\left[\begin{array}{lll} + & + & 0 \\ 0 & a & -c\\ 0 & b & d \end{array}\right] \end{array} $$
We now use standard results from the theory of linear differential equation systems. This analysis relies on the eigenvalues of our model systems. The eigenvalues of this linearized system are found by solving the det(λI−Φ)=0. Thus, for the coefficient matrix Φ, we have
or \(G_{2\boldsymbol {u}}^{\boldsymbol {S_{0}}} = G_{3\boldsymbol {a}}^{\boldsymbol {S_{0}}}\) and \(G_{3\boldsymbol {u}}^{\boldsymbol {S_{0}}} = -G_{2\boldsymbol {a}}^{\boldsymbol {S_{0}}}\). Since an eigenvalue equation of the form
$$\begin{array}{@{}rcl@{}} \det \left[\begin{array}{ll} \lambda -a & c \\ -b & \lambda - d \end{array}\right] = \lambda^{2} - (a + d) \lambda + ad + bc \end{array} $$
we have complex roots if (a−d)2−4bc<0. In particular, if a=d=α and b=β and c=β, then we have −4β2<0 and the roots are complex. We can have complex roots for other choices of a, b, c and D which implies conditions on the partials which is something we can explore, but we will not do that here. For our purposes, we will examine closely what happens when \(a = G_{2\boldsymbol {u}}^{\boldsymbol {S_{0}}} = d = G_{3\boldsymbol {a}}^{\boldsymbol {S_{0}}}\) and \(b = G_{3\boldsymbol {u}}^{\boldsymbol {S_{0}}} = -c = -G_{2\boldsymbol {a}}^{\boldsymbol {S_{0}}}\) The general solution to the system with the complex eigenvalues we have been discussing is this. The eigenvalues are \(\lambda _{1} = G_{1\boldsymbol {i}}^{o}\) and the complex conjugate pair α±βj where \(\alpha = G_{2\boldsymbol {u}}^{\boldsymbol {S_{0}}} = G_{3\boldsymbol {a}}^{\boldsymbol {S_{0}}}\) and \(\beta = G_{3\boldsymbol {u}}^{\boldsymbol {S_{0}}} = -G_{2\boldsymbol {a}}^{\boldsymbol {S_{0}}}\). The eigenvectors here are simple
Finally, recall we have \(\alpha = G_{2\boldsymbol {u}}^{\boldsymbol {S_{0}}}\) and \(\beta = G_{3\boldsymbol {u}}^{\boldsymbol {S_{0}}}\); thus, the oscillatory solutions for upregulation and free virus are
It seems unreasonable to us that the phase shift δ should be a constant; i.e. independent of S0. After all, the reasoning above is approximate and we should not think of this as actually fixed. So we will assume that all critical parameters here are proportional to S0. Our rough calculation showed us \(R = \boldsymbol {S_{0}} \sqrt {q_{1}^{2} + (1-p_{0})^{2}}\), so it seems reasonable that R is proportional to S0 in general. Therefore, we now assume
for a new parameters r1, r2, r3 and r4. This leads to our estimate of the dependence of the upregulation and free virus on the initial dose S0: we have
Now we have to compute these integrated transient values. We label them as IT for the transient i integration; UT for the transient u integration; and AT for the transient a integration. We then have
$$\begin{array}{@{}rcl@{}} \boldsymbol{IT}(t) &=& \int_{0}^{t} \boldsymbol{i}(s) ds = \int_{0}^{t} \: A e^{-G_{1\boldsymbol{i}}^{o} s} \: ds\\ \boldsymbol{UT}(t) &=& \int_{0}^{t} \boldsymbol{u}(s) ds = \int_{0}^{t} \: R e^{\alpha t} \cos(\beta s - \delta) \: ds\\ \boldsymbol{AT}(t) &=& \int_{0}^{t} \boldsymbol{a}(s) ds = - \int_{0}^{t} \: R e^{-\alpha t} \sin(\beta s - \delta) \end{array} $$
To evaluate this term, we use integration by parts. Now we have to compute the integrated transient values required to find the health estimate. We have labeled the integrations as UT for the transient u integration; and AT for the transient a integration. We then have
$$\begin{array}{@{}rcl@{}} \boldsymbol{UT}(t) &=& \int_{0}^{t} \boldsymbol{u}(s) ds = \int_{0}^{t} \: R e^{\alpha t} \cos(\beta s - \delta) \: ds\\ \boldsymbol{AT}(t) &=& \int_{0}^{t} \boldsymbol{a}(s) ds = - \int_{0}^{t} \: R e^{-\alpha t} \sin(\beta s - \delta) \end{array} $$
Integration details
First, let’s calculate UT. To evaluate this term, we use integration by parts. We find
Finally, let’s define two new parameters, θ1 and θ2 as \(\theta _{1} = \frac {r_{1}}{\sqrt {r_{2}^{2}+r_{3}^{2}}}\) and \(\theta _{2} \: = \: \tan ^{-1}\left (\frac {r_{3}}{r_{2}} \right)\). Using the above, we can rewrite UT(t) as
Using a standard reference triangle for the phase angle θ2, we see \(\cos (\theta _{2}) = \frac {r_{2}}{\sqrt {r_{2}^{2}+r_{3}^{2}}}\) and \(\sin (\theta _{2}) = \frac {r_{3}}{\sqrt {r_{2}^{2}+r_{3}^{2}}}\). We can then rewrite UT(t) again as
Another standard integration by parts similar to what was done in “Integration details” allows us to find AT. We note the same comment on the dependence of R, α, β and δ on S0 holds still. Now using these values and the terms Q1 and Q2, we we can rewrite AT(t) as follows:
Let \(c_{u} = H_{1\boldsymbol {u}}^{o} + H_{2\boldsymbol {u}}^{o} + H_{3\boldsymbol {u}}^{o}\) and \(c_{a} = H_{1\boldsymbol {a}}^{o} + H_{2\boldsymbol {a}}^{o} + H_{3\boldsymbol {a}}^{o}\). Then we have
Let’s simplify some more using another phase shift. Define the phase angle \(\theta _{3} = \tan ^{-1}\left (\frac {c_{u}}{c_{a}} \right)\); then, we can rewrite the health like this.
Next, we can combine the ratio \(H_{1\boldsymbol {i}}^{o} \frac {A}{G_{1\boldsymbol {i}}^{\boldsymbol {S_{0}}}}\) into the new parameter \(\zeta _{1}^{\boldsymbol {S_{0}}}\) and rewrite \(G_{1\boldsymbol {i}}^{\boldsymbol {S_{0}}}\) as \(\zeta _{2}^{\boldsymbol {S_{0}}}\) to give
We could also assume the terms \(\zeta _{1}^{\boldsymbol {S_{0}}}\) and \(\zeta _{2}^{\boldsymbol {S_{0}}}\) are proportional to S0. We would model this by implying \(\zeta _{1}^{\boldsymbol {S_{0}}} = r_{5} \boldsymbol {S_{0}}\) and \(\zeta _{2}^{\boldsymbol {S_{0}}} = r_{6} \boldsymbol {S_{0}}\). We then find
These parameters depend in complex ways on the initial virus dose S0 and it is very difficult to tease out the details. The only data we have to help fit this model is the survival data, so the next step of seeing how this model of health gives rise to the experimental survival data requires additional analysis.
We are going to graph Eq. 4 for various values of the parameters. First, we set the parameters as showin in Fig. 2.
Fig. 2
Set the parameters
Then, we define some auxiliary functions in Fig. 3.
Fig. 3
Define needed functions
Then, we generate a plot of health versus time for a range of choices of S0, here represented by S0 as shown in Fig. 4. For each plot, we find the minimum value and store it in the variable Min. We can then plot Min versus initial viral dose S0 and see if the plot looks like a traditional WNV survival curve. We then place all of this code into a MatLab function and use it in the usual way to generate the plots. For these choice of parameters, we can plot the minimal health curve which is shown in Fig. 5. Another choice of parameters, leads to an even better minimum health curve which captures the essence of the WNV survival plot. The parameter choices for this run are shown in Fig. 6. For these parameters values, all of the health plots versus time can be seen in Fig. 7 and show much oscillation. The corresponding minimal health curve is then shown in Fig. 8. We can make this look more like the real survival data by scaling the minimal health values and plotting them as a percentage. When we do this, we find a plot that is a bit easier to compare to the real data. It is not perfect, of course, but it is very interesting that we are capturing the essential quality of the data using our theoretical model. The percentage minimal values are shown in Fig. 9. From these experiments, it is clear what is happening. The model
can oscillate as viral load increases. To do this, it is important for the two terms in Δ(t) to be out of phase. Hence, roughly speaking cos(r3S0t−r4S0−θ2+θ3) must be sometimes negative when cos(r4S0−θ2−θ3) is positive. This allows for an increase in health of approximately ξs1 where ξ is the difference between the two terms. This is possible when the two cos arguments are out of phase by about π radians. Note it is also important that the exponential term \(e^{r_{2} \boldsymbol {S_{0}} t}\phantom {\dot {i}\!}\) allows growth. The interaction dynamics are determined by
but then the real part of the eigenvalues would be negative and we would have to model the exponential term as \(e^{-r_{2} \boldsymbol {S_{0}} t}\phantom {\dot {i}\!}\). The induced oscillations would then be damped and the potential for the WNV survival curve would vanish.
Hence, although we can use some of parameter choices from our health simulations, it is difficult to compare completely as S2, Λ1 and ϕ are different. Our final collateral damage function is then
We can easily run a quick simulation to see if our prediction that the collateral damage will have oscillations is true or not. We use the function Collateral() for this which is listed in Fig. 10. The parameter values we use are as similar as possible to the ones we used in generated the survival curves, although the values of s2, Λ1 and ϕ are necessarily new choices. We ran the simulation with these parameter values and then plotted both the maximum and minimum collateral values versus the viral dose in Fig. 11. Note that there is variation in the collateral damage due to the nonlinear interactions between the upregulation, u, and the free virus, a.
The presence of a form of self damage in our WNV infection model therefore appears to be a consequence of the nonlinear interactions in the i, u and a model:
and the chance of oscillation between the cellular population groups is lost. Hence, we can note some consequences and predictions due to our model.
The crucial assumption here is that the viral infections effect on the host divides into two parts. For a WNV infection, these two cell populations are the dividing and nondividing infected cells, D and N, respectively. We can envision other infectious agents or triggers that give rise to such a split response which then in principle could engender a similar collateral damage response which we interpret as an autoimmune reaction. So there is hope that this approach could perhaps give us insight into more general autoimmune responses. Note that Fig. 11 shows there is collateral damage that oscillates due to the infectious agent which here is WNV. It is clear that other triggering events, another virus or bacteria or even an environment toxin, could give rise to this behavior as well.
Specific to the WNV model, we assume that \(G_{2\boldsymbol {u}}^{\boldsymbol {S_{0}}} = +\), \(G_{3\boldsymbol {u}}^{\boldsymbol {S_{0}}} = +\), \(G_{2\boldsymbol {a}}^{\boldsymbol {S_{0}}} = -\), \(G_{3\boldsymbol {a}}^{\boldsymbol {S_{0}}} = +\), which then says the coefficient matrix of the linearized upregulation and free virus model has the form
This algebraic sign pattern itself can give rise to complex eigenvalues for the linearized nonlinear interaction model and we have not explored this more general problem. We have noted in our discussion in Section “Results and discussions” that if we did not have \(G_{2\boldsymbol {u}}^{\boldsymbol {S_{0}}} = +\), we could still have oscillatory behavior but it would be damped and therefore it would not explain the data we see in the survival experiments. Here, we have posited specific relations that give rise to clearcut oscillations. We have assumed \(G_{2\boldsymbol {u}}^{\boldsymbol {S_{0}}} = G_{3\boldsymbol {a}}^{\boldsymbol {S_{0}}}\) and \(G_{3\boldsymbol {u}}^{\boldsymbol {S_{0}}} = -G_{2\boldsymbol {a}}^{\boldsymbol {S_{0}}}\) which gives rise to the characteristic coefficient matrix
The assumptions above then give rise to a general model of how the minimal health changes with varying initial virus dose, which appears to allow us to approximate the data we can measure.
The exact replication of the data found in the biological situation is unlikely to occur. Indeed, a standard dose response curve of survival does not usually repeat exactly from experiment to experiment, except in the form of it, even when using genetically identical animals. The ragged dose response of mortality seen in WNV infection is similarly subject to biological variability. Clearly, with no virus there will be 100% survival, and with a large amount of virus there will be 100% death, as expected in a standard dose response curve. In between these two doses, however, the response to infection is subject to probability, which affects the outcome (survival or death). Thus, if the experiment were undertaken several times it would show the ragged form on each occasion, but not exactly the same percentage survival at each dose used. This implies that small biological differences at the starting point of infection, albeit in genetically identical mice, may subtend a large range of endpoints. It is of interest to note that bypassing the early initiation of the adaptive immune response, by inoculating virus intracranially, that the standard graded dose response occurs with WNV. This is because the replication of the virus overtakes the animal before an effective immune response can be generated, emphasising the role of the immune system in generating this ragged survival curve [12].
We have shown that we can build a reasonable model of how WNV infects a host’s cell in such a way that damage to the host can decrease, even though the inoculating viral dose increases.
Our model is a macro one and we believe it provides insight as to how we can model more general auto-immune reactions. We propose that for an infectious agent or trigger to cause oscillations in health it is required that the trigger causes alterations in two distinct cell populations. Then, if the nonlinear interactions between these two populations satisfies the conditions for damped oscillatory response we have mentioned here, we should see oscillations in the host health. We consider this work essentially a theoretical model and as we explore what we can do to extend the results to more general auto-immune settings, we hope that we can generate greater mechanistic insight into auto-immune disease.
Publication of this article was funded by National Health and Medical Research Council Project Grants, Numbers, 512413, 1030897, to NJCK and Australian Research Council Discovery Project Number DP0666152 to AK and NJCK.
Availability of data and materials
Not applicable.
Authors’ contributions
All authors were equally responsible for the ideas and development of the work with JP responsible for the mathematical and computational details. All authors have read and approved the manuscript.
Competing interests
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver(http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
Department of Biological Sciences, Department of Mathematical Sciences, Clemson University, Martin Hall O-304, BOX 340975, Clemson, SC, 340975, USA
(2)
Discipline of Child and Adolescent Health and the Marie Bashir Institute for Emerging Infectious Diseases and Biosecurity, University of Sydney, Sydney, NSW, Australia
(3)
Sydney Children’s Hospital Network, Westmead, Sydney, NSW, Australia
(4)
Discipline of Pathology, University of Sydney, Sydney, NSW, Australia
Diamond MS. (ed.): West Nile Encephalitis Virus Infection: Viral Pathogenesis and The Host Immune Response. New York: Springer; 2009.View ArticleGoogle Scholar
Abbas AK, Lichtman AH, Pillai S. Cellular and Molecular Immunology. New York: Saunders Elsevier; 2010.Google Scholar
Douglas DW, Kesson AM, King NJC. CTL recognition of west Nile virus-infected fibroblasts is cell cycle dependent and is associated with virus-induced increases in class I MHC antigen expression. Immunology. 1994; 82:561–70.PubMedPubMed CentralGoogle Scholar
Kesson AM, Cheng Y, King NJC. Regulation of immune recognition molecules by flavivirus, West Nile. Viral Immunol. 2002; 15:273–83.View ArticlePubMedGoogle Scholar
Müllbacher A, Hill AB, Blanden RV, Cowden WB, King NJC, Hla RT. Alloreactive cytotoxic T cells recognize MHC class I antigen without peptide specificity. J Immunol. 1991; 147:1765–72.PubMedGoogle Scholar
King NJC, Müllbacher A, Tian L, Rodger JC, Lidbury B, Hla RT. West Nile virus infection induces susceptibility of in vitro outgrown murine blastocysts to specific lysis by paternally directed allo-immune and virus-immune cytotoxic T cells. J Reprod Immunol. 1993; 23:131–44.View ArticlePubMedGoogle Scholar
Kuhlman P, Moy VT, Lollo BA, Brian AA. The accessory function of murine intercellular adhesion molecule-1 in T lymphocyte activation. Contributions of adhesion and co-activation. J Immunol. 1991; 146:1773–82.PubMedGoogle Scholar
Shen J, T-To SS, Schrieber L, King NJC. Early E-selectin and VCAM-1 and ICAM-1 and and late major histocompatibility complex antigen induction on human endothelial cells by flavivirus and comodulation of adhesion molecule expression by immune cytokines. J Virol. 1997; 71:9323–32.PubMedPubMed CentralGoogle Scholar
Getts DR, Matsumoto I, Muller M, Getts MT, Radford J, Shrestha B, Campbell IL, King NJC. Role of IFN- γ in an experimental murine model of West Nile virus-induced seizures. J Neurochem. 2007; 103:1019–30.View ArticlePubMedGoogle Scholar
Peterson J, Kesson AM, King NJC. A Simulation For Flavivirus Infection Decoy Responses. Adv Microbiol. 2015; 5(2):123–42. doi:10.4236/aim.2015.52013.
