A model of auto immune response

Background In this work, we develop a theoretical model of an auto immune response. This is based on modifications of standard second messenger trigger models using both signalling pathways and diffusion and a macro level dynamic systems approximation to the response of a triggering agent such as a virus, bacteria or environmental toxin. Results We show that there, in general, will be self damage effects whenever the triggering agent’s effect on the host can be separated into two distinct classes of cell populations. In each population, the trigger acts differently and this behavior is mediated by the nonlinear interactions between two signalling agents. Conclusion If these interactions satisfy certain critical assumptions this will lead to collateral damage. If the initial triggering agent’s action involves any critical host cell population whose loss can lead to serious host health issues, then there is a much increased probability of host death. Our model also shows that if the nonlinear interaction assumptions are satisfied, there is a reasonable expectation of oscillatory behavior in host health; i.e. periods of remission.


Background
In [1,2] we explore a micro level simulation model of a single host's response to varying levels of West Nile Virus (WNV) infection. In that infection, there is a substantial self damage component and in those papers, we show that this is probably due to the way that the virus infects two cell populations differently. This difference, which involves an larger upregulation of MHC-1 sites on the surface of nondividing infected cells over dividing infected cells, is critical in establishing a self damage or collateral damage response. In [3], we develop a macro level model of the nonlinear interactions between two critical signalling agents that mediate the interaction between these two sets of cell populations. In the case of WNV infection, the two signals are the MHC-1 upregulation level of the cell and the free WNV antigen level. This macro model allowed us to predict a host's health response to varying There are then three nonlinear interaction functions F 1 , F 2 and F 3 because we know C, M and N depend on each other's levels in very complicated ways. Usually, we assume the initial trigger dose V 0 gives rise to some fraction of infected cells and the effect of the trigger will be different in the two cell populations M and N .

Assumption 1
We assume the number of infected cells is p 0 V 0 which is split into p 1 p 0 V 0 in population N and p 2 p 0 V 0 in M, where p 1 + p 2 = 1.
For example, a reasonable choice is p 1 = 0.99 and p 2 = 0.01. Thus, the total amount of trigger that goes into altered cells is p 0 V 0 and the amount of free trigger is therefore (1−p 0 )V 0 . Thus, we could expect C 0 = 0, M 0 = p 2 p 0 V 0 and N 0 = M 0 = p 1 p 0 V 0 . However, we will explicitly assume we are starting from a point of equilibrium prior to the administration of the viral dose V 0 . We could assume there is always some level of collateral damage, C 0 in a host, but we will not do that. We will therefore assume C, M and C have achieved these values C 0 = 0, M 0 = 0 and N 0 = 0 right before the moment of alteration by the trigger. Hence, we don't expect to there to be initial contribution to C (0), M (0) and N (0); i.e. Next, we do a standard tangent plane approximation on the nonlinear dynamics functions F 1 , F 2 and F 3 to derive approximation dynamics. The mathematics behind this approximation come from multivariate calculus and can easily be reviewed if required. We find the approximate dynamics are where we now use a standard subscript scheme to indicate the partials. Now let's add the signals IFN-γ (I), J and K to the mix.

The CDN IJK model
We can think each variable C, M and N as depending on I, J and K . Thus, we have As before assume C, M and C have achieved the same optimal values C 0 = 0, M 0 = 0 and N 0 = 0 prior to the moment of infection with trigger dose V 0 . These correspond to the starting values prior to exposure to the trigger of I 0 , J 0 and K 0 . Initially, we don't expect IFN-γ signals so I 0 = 0. Eventually, we do expect some level of change in J and K due to this initial dose and we will assume this change to be proportional to the level of the dose V 0 applied; that is, we will assume this is a simple scaling factor, i.e. J 0 = q 1 V 0 for some suitable parameter q 1 . Also, once the trigger has been applied dose, we would expect some fraction of it to remain free which will be modeled as K 0 = (1 − p 0 )V 0 . But now, we think of all the initial values as zero; i.e. I 0 = 0, J 0 = 0 and K 0 = 0. We still don't expect to have any contribution to C (0), M (0) and N (0); i.e.
Thus, we have the changes Now we need to estimate i, j and k.

The IJK model
The amount of I, J and K depend on the initial amount of trigger applied when in the equilibrium state; i.e. this is the amount that causes the initial infection. This is V o . We assume the dynamics here are In the model of "The CDN IJK model" section, we assumed C, M and N depended on the perturbations of I, J and K from a zero state. Now, we want to model the I, J and K deviations from a base state I 0 , J 0 and K 0 which is not zero. As previously discussed, we expect K 0 = (1 − p 0 ) V 0 , the initial IFN-γ level I 0 = 0 and the initial upregulation level J 0 = q 1 V 0 . Let the deviations from these equilibrium values be given by i = I − I 0 , J = J − J 0 and k = K − K 0 . The model can then be rewritten as The usual tangent plane approximation on the nonlinear dynamics functions G 1 , G 2 and G 3 then gives the dynamics approximation 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-γ .

Assumption 5
For the signals j and k, we assume Thus, the coefficient matrix above which we call so far looks like Now hold everything constant except j and increase j to j + δj. What happens?
Thus, the coefficient matrix looks like Now hold everything constant except k and increase k to k + δk. What happens?

Assumption 7
We assume an increase in k will not effect IFN-γ levels, i, and hence G V 0 1k = 0. An increase in k must imply G V 0 3k = +. We then assume the signals j and k interact so that G V 0 2k = −.
Thus, the coefficient matrix now looks like

Oscillations in J and K
The eigenvalues of this linearized system are found by setting det(λI − ) = 0. Thus, the coefficient matrix above which we call so far looks like The eigenvalues of the two by two submatrix are the most interesting. We can get complex roots if ⎡ The eigenvalues are then λ 1 = G o 1i and the complex conjugate pair α±β We can solve for j and k to find = e αt ((aV + bW ) cos(βt) We then have Hence, where δ is defined as 1i , β and δ determine a given model. Here, we have J 0 = q 1 V 0 and K 0 = (1 − p 0 )V 0 . Hence, we roughly know at the time of the initial disturbance (infective agent or environmental toxin etc.) Taking a ratio, we find Hence, δ = tan −1 (1−p 0 ) Finally, recall we have α = G V 0 2j and β = G V 0 3j ; thus, the oscillatory solutions for j and k are We do not think the phase shift δ should be a constant; i.e. independent of V 0 . and therefore, we assume that the critical parameters here are proportional to V 0 . Our rough calculation showed us R = V 0 q 2 1 + (1 − p 0 ) 2 , and thus, R should be proportional to V 0 in general. Therefore, we assume for a new parameters r 1 , r 2 , r 3 and r 4 . This leads to our estimate of the dependencies

A health model
Roughly speaking, if the total number of cells is T, the number of healthy cells can be approximated by We know k and so we are looking at deviations from the base values I 0 = 0, J 0 = q 1 V 0 and K 0 = (1 − p 0 )V 0 . It follows we have

k(s)ds
As discussed earlier, we have initially, C 0 = 0, M 0 = p 2 p 0 V 0 and N 0 = p 2 p 0 V 0 . So we have Thus, we have

k(s)ds
Now collect all the terms involving V 0 and set that coefficient to for convenience. Making this replacement, we have This leads to the simplification

k(s)ds
Now we have to compute these integrated transient values. We label them as IT for the transient i integration; JT for the transient J integration; and KT for the transient K integration. We then have Re αt cos(βs − δ) ds

Integration details
The i integration is easy.
The JT is a standard integration by parts.
To evaluate this term, we use integration by parts. We find We can rewrite this is a much better form using our assumptions. First, rewrite as Now we know α, β, δ and R are really dependent of V 0 . For convenience of exposition, let's drop the superscript V 0 in our calculations below Finally, let's define two new parameters, θ 1 and θ 2 as θ 1 = r 1 r 2 2 +r 2 3 and θ 2 = tan −1 r 3 r 2 . Using the above, we can rewrite JT(t) as Using a standard reference triangle for the phase angle θ 2 , we see cos(θ 2 ) = r 2 r 2 2 +r 2 3 and sin(θ 2 ) = r 3 r 2 2 +r 2 3 . We can then rewrite JT(t) again as and using standard trigonometric identities, we then have Next, another standard integration by parts shows We note the same comment on the dependence of R, α, β and δ on V 0 holds still. Now using these values and the terms Q 1 and Q 2 , we we can rewrite KT(t) as follows: Now using the simplifications we obtained for α and β in terms of r 2 and r 3 , we can rewrite this complicated expression as Next, using the phase shift θ 2 , we have This then leads to our final form

Building the health model
Recall the health model is Now plug what we have found for our integrations. We have Then we can rewrite as Now put the e r 2 V 0 t together. We find Let's simplify some more using another phase shift. Define the phase angle θ 3 = tan −1 c j c k ; then, we can rewrite the health like this.
This can be recast as Finally, let s 1 = θ 1 c 2 j + c 2 k . Then, we have the last form of the health estimate: We could also assume the terms ζ V 0 1 and ζ V 0 2 are proportional to V 0 . We would model this by implying ζ V 0 1 = r 5 V 0 and ζ V 0 2 = r 6 V 0 . We then find These parameters depend in complex ways on the initial trigger dose V 0 and it is very difficult to tease out the details.

Collateral damage
We can also work out the functional dependence on collateral damage on initial trigger dose over time. Recall the collateral damage population is given by We can then substitute for IT(t), JT(t) and KT(t) and obtain Now, collect terms as we did in our earlier simplifications. We rewrite as We can also introduce an additional phase shift, phi, as follows. It will be different from the phase shift as here we only use the H 1 partials: . We rewrite as We can then use the the usual cos laws of addition and subtraction of angles to repackage this as 2 and rewrite as Since collateral damage is initially zero, we have as our final form Previously, we used the simplification

Our final collateral damage function is then
Of course, since S 2 , 1 and φ are different from the corresponding values in the health function, it is a bit difficult to compare simulation results, but it is easy to see the qualitative ideas of oscillation in health and collateral. We have done similar experiments in [3] and indeed the graphs we now show were generated using the same MatLab code. The point is that the existence of the oscillation in health, collateral damage and so forth is due to the assumptions we made on the nonlinear interactions between the two populations M and N mediated by the signals J and K . As long as those sorts of interactions are occurring, this kind of interaction behavior is assured; and that is very interesting we feel. We can easily run a quick simulation to see if our predictions of oscillations in the health and collateral damage function are verified. We use the same parameter settings and MatLab code as we used in [3]. The interested reader can look those details up as necessary. We ran the simulation with the chosen parameter values from [3] and plotted both the maximum and minimum collateral values versus the trigger dose in Fig. 1.
Note that there is variation in the collateral damage due to the nonlinear interactions between the J, and K . This is due to the assumptions we have made on the kinds of nonlinear interactions that occur. The remainder of this paper will necessarily discuss why we think these kind of nonlinearities are possible in a variety of auto immune situations.
As noted in [3], it is clear what is happening. The model can be written in terms of decay and push -pull terms as follows: Thus, we have H(t) always decreases unless the pushpull terms counteract that decay. Hence, what is important is the term can oscillate as trigger load increases. To do this, it is important for the two terms in (t) to be out of phase. Hence, roughly speaking cos (r 3 V 0 t − r 4 V 0 − θ 2 + θ 3 ) must be sometimes negative when cos (r 4 V 0 − θ 2 − θ 3 ) is positive. This allows for an increase in health of approximately ξ s 1 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 V 0 t allows growth. The interaction dynamics are determined by Note we can also plot the minimum health obtained over the simulation time for each trigger dose as we did in [3]. We find the percentage minimal values are shown in Fig. 2.
Let's take a moment to reflect on what we are seeing here. If we assume a trigger has an effect on the host's cell population that leads to two distinct cell populations M and N that interact in nonlinear ways following Assumption 1 -Assumption 7 due to the mediating influences of two signals J and K , we inevitably see the host healthy cell population over time have up and down variations due to the size of initial trigger V o . Consider the health plots for four levels of V o shown in Fig. 3.
Note that the minimal values of health do indeed decline for increasing initial trigger dose (although we see up and down variation if we look at all the minimal values versus initial trigger dose as we show in Fig. 1. What we want to focus on here is that the health at many initial trigger doses oscillates over the time of the simulation. For example, the health plot corresponding to V o = 25 has a very low minimum value. These plots are generated by abstractions of health and collateral damage and so it is not clear at all how to relate them to the health of a real host, but it does suggest that the host health can rebound from a low value. Hence, collateral damage can diminish for a given initial trigger which shows a kind of relapse effect.
Also note, the theoretical model we have built so far generates what we call collateral damage and relates it to general health with no discussion of T Cell recognition of infected cells. For us collateral damage is related to IFN-γ signalling which is generated by the lysis of cells in the host. In the West Nile Virus infection studied in [1], the splitting into two separate cell populations due to the WNV antigen causes explicitly changes in the avidity computation of the T Cells that recognize MHC-1 complexes on infected cells. These changes lead to an enhanced probability that T Cells will target healthy cells because self proteins become more visible to the adaptive immune system. In [3], we derive the heath and collateral damage model we also develop in this paper by focusing very explicitly on the splitting phenomena. The possibility of self damage is now redirected much more strongly to the splitting into two cell populations which allows for a much more general treatment of self damage. However, it is now time to think more deeply about what the signals J and K would be in specific cases of auto immune response to a trigger.

General trigger models
We now want to think carefully about the signals J and K and the trigger V . Hence, we consider the transcriptional control of a regulatory molecule which can be co-opted by an external trigger. A good example is the regulatory molecule NFκB whose normal action is changed by the WNV antigen. It always plays a role in immune system response but the virus alters what it does in many ways. We have discussed this sort of modeling in [4] for the purpose of approximating computation in excitable neurons and in that paper, we follow the spirit of the semi-abstract approach in [5]. We are now going to re-task it for our purposes of an auto immune discussion, however, this same way of looking at signalling for the context of altering computation in a cognitive model provides another way of looking at the same idea and we think it is always useful to examine a complicated mechanism from multiple points of view.
Consider a trigger T 0 which activates a cell surface receptor. Inside the cell, there are always protein kinases that can be activated in a variety of ways. Here we denote a protein kinase by the symbol PK. A common mechanism for such an activation is to add to PK another protein subunit U to form the complex PK/U . This chain of events looks like this: T 0 → CSR → PK/U where CSR denotes a cell surface receptor. PK/U then acts to phosphorylate another protein. The cell is filled with large amounts of a transcription factor we will denote by T 1 and an inhibitory protein for T 1 we label as T ∼ 1 . This symbol, T ∼ 1 , denotes the complement or anti version of T 1 . In the cell, T 1 and T ∼ 1 are generally joined together in a complex denoted by T 1 /T ∼ 1 . The addition of T ∼ 1 to T 1 prevents T 1 from being able to access the genome in the nucleus to transcribe its target protein. Using methods similar to those discussed in [4], we find the effects of the trigger T 0 on the change in protein production T 1 , δ T 1 , can be modeled by e >> [PK/U ] e for β >> 1. From this quick analysis, we can clearly see the potentially explosive effect changes in T 0 can have on PK/U . We can use these results for our general immune interaction discussion. We have assumed there are two different cell populations M and N which are effected differently by the trigger. Hence, let's assume there is a pathway involving potentially many steps that leads to an alteration of the fragility of these two cell populations. In a West Nile Virus infection, this fragility is related to the upregulation of MHC-1 complexes but it could be another type of alteration of the overall health of the cell. Let the pathway leading to a change in fragility for M be P M and the one leading to the fragility alteration in N be P N . Then associated such a change in fragility with the alteration of a protein or protein complex we call T M or T N depending on the cell population type. Then we have for parameters μ M , μ N , M and N . Hence, we can estimate a strength level for the trigger effect on each population that leads to increases in fragility.

Computational abstractions
A close look extracellular triggers abstractly helps us understand how to approximate their effects. Let T 0 denote a second messenger trigger which moves though a port P to create a new trigger T 1 some of which binds to B 1 . A schematic of this is shown in Fig. 4. In the figure, r is a number between 0 and 1 which represents the fraction of the trigger T 1 which is free in the cytosol. Hence, 100r% of T 1 is free and 100(1 − r) is bound to B 1 creating a storage complex B 1 /T 1 . For our simple model, we assume rT 1 is transported to the nuclear membrane where some of it binds to the enzyme E 1 . Let s in (0, 1) denote the fraction of rT 1 that binds to E 1 . We illustrate this in Fig. 5.
We denote the complex formed by the binding of E 1 and T 1 by E 1 /T 1 . From Fig. 5, we see that the proportion of T 1 that binds to the genome (DNA) and initiates protein creation P(T 1 ) is thus srT 1 .
The protein created, P(T 1 ), could be many things. Here, let us assume that P(T 1 ) is a protein that binds to the promoter for one of the many proteins that effect MHC-1 creation, peptide binding etc. For our purposes, we will call it Q. Thus, our high level model is sE 1 /rT 1 + DNA → Q. We therefore increase the concentration of MHC-I complexes on the surface of the call, thereby making the cell more visible to the adaptive immune system.
We can model increases in Q as increases in T Cell binding efficiency, where e is a number between 0 and 1. We will not assume all of the sE 1 /rT 1 +DNA to M reaction is completed. It follows that e is similar to a Michaelson -Mentin kinetics constant. Our full schematic is then given in Fig. 6.
We can model the choice process, rT 1 or (1 − r)B 1 /T 1 via a simple sigmoid, f (x) = 0.5(1 + tanh( x−x 0 g )) where the transition rate at x 0 is f (x 0 ) = 1 2g . Hence, the "gain" of the transition can be adjusted by changing the value of g. We assume g is positive. This function can be interpreted as switching from of "low" state 0 to a high state 1 at speed 1 2g . Now the function h = rf provides an output in (r, ∞). If x is larger than the threshold x 0 , h rapidly transitions to a high state r. On the other hand, if x is below threshold, the output remains near the low state 0.
We assume the trigger T 0 does not activate the port P unless its concentrations is past some threshold [ T 0 ] b where [ T 0 ] b denotes the base concentration. Hence, we which has output in [ 0, e). Here, we want to limit how large a change we can achieve in Q. Hence, we assume there is an upper limit which is given by Q = δ Q Q max . Thus, we limit the change in the the expression of Q to some percentage of a baseline value. It follows that h Q (x) is about δ Q if x is sufficiently large and small otherwise. This suggests that x should be [P(T 1 )] and since translation to P(T 1 ) occurs no matter how low [T 1 ] is, we can use a switch point value of x 0 = 0. We conclude Our model of the change in Q expression is therefore We can use these results for our general immune interaction discussion as well. From earlier comments, we

Diffusion trigger models
We will now look at the trigger from a diffusion perspective. This is different as in the earlier models, we focus on how cell populations change in time and derive health and collateral functions that show their dependence on the initial viral dose. We did not consider any spatial relationships between the cellular populations. Now we will do so and the discussion will give us another way to look at the nonlinear interactions between M and N . It is well known that second messenger systems often involve Ca ++ ion movement in and out of the cell. The amount of free Ca ++ ion in the cell is controlled by complicated mechanisms, but some is stored in buffer complexes. The release of calcium ion from these buffers plays a big role in cellular regulatory processes and which protein P(T 1 ) is actually created from a trigger T 0 . The diffusion model is very powerful. Consider some substance u which satisfies where D is called the diffusion constant for this substance. For simplicity this is a one dimensional model where the spatial variable x comes from the line segment [ 0, L]. For example, instead of a three dimensional cell modeled as a sphere, the cell is modeled as a string of finite length. Hence, stuff can only enter the cell from either the right or the left. The condition −D ∂u ∂x | 0,L = J 0,L states that there are conditions on the flux of u through the boundary at x = −0 or x = L that must be satisfied. The term J 0,L can be thought of as an injection of current into the cell. In this model, we think of the substance u as being in equilibrium in the cell and the current injection J 0,L alters that equilibrium and the diffusion model, when solved, tells us what happens to u due to the sudden current injection at the boundary. The critical review on the control of free calcium in cellular processing in [6] notes the concentration of Ca ++ in the cell is controlled by the reversible binding of calcium ion to the buffer complexes, B 1 , B 2 and so forth. There in general are quite a few different buffer complexes which all behave differently. These buffer molecules act as calcium ion sensors that, in a sense, decode the information contained in the calcium ion current injection and then pass on a decision to a target. Many of these targets can be proteins transcribed by accessing the genome. Hence, the P(T 1 ) could be a buffer molecule B j and so the trigger that causes the calcium current injection into the cell could influence the concentration of a buffer B j and therefore influence how it itself is decoded. In this situation, the boundary condition J 0,L plays the role of the entry calcium current. Such a calcium ion input current through the membrane could be due to membrane depolarization causing an influx of calcium ions through the port or via ligand binding to a receptor which in turn indirectly increases free calcium ion in the cytosol. Such mechanisms involve the interplay between the way the endoplasmic reticulum handles release and uptake and also the storage buffers. This boundary current in general therefore determines the u(t, x) solution through the diffusion equation. The exact nature of this solution is determined by the receptor types, buffers and storage sites in the ER. Differences in the period and magnitude of the calcium current u(t, x) resulting from the input J 0,L trigger different second messenger pathways. Hence, there are many possible outcomes due to a given input current J 0,L .
We will now modify discussions in [7,8] and [9] that show us how to model calcium ion movement in the cell to develop a model of trigger movement in and out of the cell populations M and N .
We assume the trigger enters the host and can be sequestered in some form in both M and N cell populations. We also assume the trigger has associated with some sort of diffusion process; letting u(t, x) be the concentration of the trigger in the host at time t and spatial position x, we posit u t = D 0 u xx . Also, note we present our arguments as if the host was one dimensional; i.e., the two cell populations lies along a one dimensional axis measured by the location variable x. This is, of course, very simplistic, but we only want to suggest some functional dependencies here, so it will suffice for our purposes. Let's assume these two populations use or bind the trigger with binding rate constants k + M and k + N and disassociation rate constants k − M and k − N . Let the total number of cells in the host be P. Then the fraction of cells in the M population is M P which we call C M and the fraction in the N population is M P which is denoted by C N . Hence, the number of cells that have not been exposed to the trigger is P − M(t, x) − N (t, x) = F(t, x). Now, let u(t, x) be the concentration of free trigger in the host at (t, x). Some of the trigger has been used to create cells in the populations M and N, but the rest is unused. Hence, if C is a cell which has not been altered by the trigger, we have the reactions where T is the trigger. Of course, the equations above depend on time and space, but we have not written in that dependence to avoid clutter. Note also, in this context, the backward reaction in which trigger is freed from the cellular populations M(t, x) and N (t, x) is typically that of lysis and so the backward rates are part of our nonlinear interaction model. We are just adding low level detail. The corresponding dynamics are and we also know To make the manipulations easier, let B M = 1 − N P and B N = 1 − M P . We can rewrite the equation above as

x)
A similar analysis gives the amount of trigger being added to the N population as Thus, the diffusion dynamics are where D 0 is diffusion coefficient for free trigger. We now assume the spread of cells we collect into the populations M and N satisfy some sort of diffusion law. Certainly, cells are added to these populations as the trigger diffuses throughout the host's body. Hence, this assumption is a good start. We therefore assume the diffusion dynamics for C M and C N are given by Now consider the free trigger plus a correction due to the population fractions C M and C N . We denote this by w(t, x) and note that w = u + C M + C N and Thus, This simplifies to ∂w ∂t Thus, w satisfies where we have not shown the derivation of the boundary terms here as they are less germane to our interests. It seems reasonable to assume that the interaction with the cell populations, determined by k − M and k + M u is fast and reaches equilibrium quickly. Hence, we will assume that Solving, we find We can then rewrite our equations as Plugging in for u, we have Then, another rewrite gives Notice that if we define a new diffusion coefficient, D for the diffusion process that governs w by we obtain

Further approximations
What would it mean if u K M and u K N ? Let's take the first case: u represents the rate that the M cells lose the trigger. This occurs only when the cells are destroyed. These cells are destroyed either because they have exceeded their normal lifespans or because the trigger has made them more fragile. This fragility could mean the cells have caught the attention of the immune system and are being destroyed or the fragility of the cells is such that standard apoptosis strategies are employed to remove the damaged cell. Note losing trigger from the M cells therefore corresponds to increasing trigger concentration. Since k + M is the rate at which M cells are formed, k + M u giving the amount of trigger lost from the formation of the M cells per M cell. We generally assume that in a trigger situation, the trigger is growing inside the M and N cell populations. If the trigger is a virus, replication only occurs after infection and once inside these cells, the virus grows. This additional trigger is then released into the host upon lysis of a M or N cell. It seems reasonable to assume the amount of released trigger is usually quite a bit bigger than the amount of trigger that initiates the formation of these cells. Thus, the inequality k + M u k − M seems reasonable. A similar argument shows that k + N u k − N . Thus, it is reasonable to assume u K m and u K N which leads to the approximations Thus, w has become We can then rewrite C M and C N as We can find the partials with respect to w: We can then redo our calculations for w t .
Letting denote the term 1 + γ M + γ N , then w = u and so we have We also have Thus, Define the new diffusion constantD byD = D 0 +D M γ M +D N γ N . The free trigger dynamics are thus ∂u ∂t =Du xx .

Approximations to trigger modification
Let's examine what might happen if a trigger event T 0 initiated an increase in M. This trigger event initiates a complex pathway culminating in a protein transcription (see Section "General trigger models" for the general trigger discussion). Recall, we let the pathway leading to a change in fragility for M be P M . We associate such a change in fragility with the alteration of a subsidiary signal T M and we derived Let's assume C M is increased to C M + . It is reasonable to assume that both k + M and k − M are independent of the amount of C M that is present. The same comment holds for k + N and k − N . We have B M = 1 − N P stays the same, but Thus, the new value isˆ = − K N . This implieŝ To first order, we know − ξ M ≈ and so D ≈ D N 1 −D D N . We conclude the new diffusion dynamics are on the order of This change in the solution u(t, x) then can then initiate further changes in the distribution of the M and N cells. A similar argument can be used for a change in the N population and it is clear if the signal generates changes in both M and N , to first order we generate an altered diffusion model whose solution gives us clues as to new behavior.
Without boundary conditions, the general solution to a diffusion model with diffusion constant D is given by Hence, our usual trigger solution is 4Dt .

Results and discussion
We have been studying auto immune reactions from a theoretical point of view in this work. We begin by building a model of an auto immune response which is due to nonlinear interactions in two populations of cells M and N which are mediated by the two signals J and K . We do not specify at this time what these two signals are and instead we argue from first principles. There is a third signal also, the IFN-γ level, which is denoted by I. We let the deviations of these signals from base levels be given by i, j and k.
Then the presence of a form of self damage in this model appears to be a consequence of the nonlinear interactions in the i, j and k model: where Note, if the two populations M and N coincide, this model reduces to a two dimensional model, we drop one signal, say k, and we have i k 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 triggering event has an effect on the host that divides into two parts. For a WNV infection, these two cell populations are the dividing and nondividing infected cells, D and N, respectively. But here, we have posited that these two cell populations are given by M and N instead. Hence, any infectious agents or trigger that gives rise to such a split response engenders a similar collateral damage response which we interpret as an autoimmune reaction. We note this give us a significant insight into general autoimmune responses. Note that Fig. 1 shows there is collateral damage that oscillates due to the trigger and we could interpret a downswing in collateral damage as a relapse event.
• We have assumed G V 0 2j = +, G V 0 3j = +, G V 0 2k = −, G V 0 3k = +, which then says the coefficient matrix of the linearized upregulation and free antigen 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. Here, we have posited specific relations that give rise to clearcut oscillations. We have assumed G V 0 2j = G V 0 3k = α V 0 and G V 0 3j = −G V 0 2k = β V 0 which gives rise to the characteristic coefficient matrix The remaining questions are then to try to understand the algebraic sign patterns from a low level analysis of the trigger signal that initiates the potential auto immune reaction. We have provided in this paper, three different ways to look at the trigger response.
• We analyze a general trigger in Section "General trigger models" and show that we can understand how the signal generates protein transcription changes via equations such as which shows how changes in the trigger generate changes in another signal via a cascade of protein transcriptions culminating in a change in fragility for M. • Then, in Section "Computational abstractions", a computational approach using sigmoid activation function h M further found that if we associate the change in fragility in M due to T M to be the alteration of a protein complex P(T M ), then the maximum P(T M ) expression is on the order of

P(T M ) max = h P(T M ) ([ P(T M )] )
where h is a traditional sigmoid response function which switches the protein from a low to a high level.
• We can also model the trigger signal in terms of a diffusion process as we did in Section "Approximations to trigger modification". The trigger effects M and N is complicated ways and in this section, we study what happens if the pathways the trigger operates on generate a change in M itself. We show this in turn alters the diffusion coefficient that controls the trigger dynamics thereby potentially altering all of the trigger pathways.

Conclusion
In conclusion, we have shown that we can build a reasonable model of how a trigger agent such as a virus, a bacteria or an environmental toxin can infect a host's cell and cause an autoimmune reaction. Part of our model is a macro one and we believe it provides clues as to how general auto immune reactions behave. We have posited 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 and collateral damage. Another part of our model is a detailed micro level one which looks at the triggering pathway and examines how the events in that pathway can contribute to the nonlinear interaction assumptions we make. We also discuss in detail how a sudden trigger increase can be modeled in terms of a diffusion based response and how that also can influence the nonlinear interactions we need for an auto immune response. We will finish with a few speculations about how to decide if there is an auto immune response. If one is suspected, we could run the following experiment. Let's assume there are N possible signals ω 1 , . . . , ω N that we think could be important in the analysis of the chosen potential auto immune response. There are many possible cytokines, chemokines and other molecular agents that could be of interest. We then setup a standard M × M well type genomic assay experiment. Each well is designed by growing a three dimensional organoid which will play the role of a host. This is quite possible and the paper [10] provides clear guidelines as to how to do this for mini human brains grown from stem cells. We would first measure cytokine and chemokine etc. levels in all of the wells as well as other expressed proteins for a given trigger level. At this stage, we are trying to find to the populations M and N which handle the trigger level differently as this is part of our assumptions that are needed to cause an auto immune response. If two such populations are discovered, then we need to identify the agents J and K to use in the model. Assume N = 5 for convenience. Then, for each