A model of auto immune response

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 levels of initial virus dose. Hence, we could begin to understand the oscillations in collateral damage and host health that lead to the survival data we see in WNV infections. The more general musings of [3] will now be extended to the setting of auto immune interactions in general. The derivations here use standard ideas from advanced calculus and differential equations.

We assume we have a large population of cells T which consists of cells which are infected or altered in two distinct ways by a trigger V. based on signals I, J and K. These two distinct populations of cells will be labeled M and N. There are also 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 will be using the same approach to studying nonlinear interactions that was used in [3].

There are then three nonlinear interaction functions F ₁,F ₂ and F ₃ 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 ₀ 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.

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.

where we now use a standard subscript scheme to indicate the partials. 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 other two cell populations, M and N?

Now hold everything constant except j and increase j to j+δ j. What happens?

Next hold everything constant except the signal level k and increase k to k+δ k. What happens?

Now we need to estimate i,j and k.

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- γ.

Now hold everything constant except j and increase j to j+δ j. What happens?

Now hold everything constant except k and increase k to k+δ k. What happens?

These parameters depend in complex ways on the initial trigger dose V ₀ and it is very difficult to tease out the details.

Of course, since S ₂, Λ ¹ 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.

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 {V_{0}} t}\phantom {\dot {i}\!}\). The induced oscillations would then be damped it would not be possible to see oscillatory grown in the collateral damage function.

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.

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 N F κ 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 ₀ 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 to form the complex . This chain of events looks like this: where CSR denotes a cell surface receptor. then acts to phosphorylate another protein. The cell is filled with large amounts of a transcription factor we will denote by T ₁ and an inhibitory protein for T ₁ we label as \(T_{1}^{\sim }\). This symbol, \(T_{1}^{\sim }\), denotes the complement or anti version of T ₁. In the cell, T ₁ and \(T_{1}^{\sim }\) are generally joined together in a complex denoted by \(T_{1}/T_{1}^{\sim }\). The addition of \(T_{1}^{\sim }\) to T ₁ prevents T ₁ 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 ₀ on the change in protein production \(T_{1}, \delta _{T_{1}}\phantom {\dot {i}\!}\), can be modeled by

for β>>1. From this quick analysis, we can clearly see the potentially explosive effect changes in T ₀ can have on .

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 ₀ denote a second messenger trigger which moves though a port P to create a new trigger T ₁ some of which binds to B ₁. 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 ₁ which is free in the cytosol. Hence, 100r% of T ₁ is free and 100(1−r) is bound to B ₁ creating a storage complex B ₁/T ₁. For our simple model, we assume r T ₁ is transported to the nuclear membrane where some of it binds to the enzyme E ₁. Let s in (0,1) denote the fraction of r T ₁ that binds to E ₁. We illustrate this in Fig. 5.

We denote the complex formed by the binding of E ₁ and T ₁ by E ₁/T ₁. From Fig. 5, we see that the proportion of T ₁ that binds to the genome (DNA) and initiates protein creation P(T ₁) is thus s r T ₁.

The protein created, P(T ₁), could be many things. Here, let us assume that P(T ₁) 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 \(\boldsymbol {\mathcal {Q}}\). Thus, our high level model is \(sE_{1}/rT_{1} \: + \: \text {DNA} \rightarrow \boldsymbol {\mathcal {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 \(\boldsymbol {\mathcal {Q}}\) as increases in T Cell binding efficiency, where e is a number between 0 and 1. We will not assume all of the s E ₁/r T ₁ + 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, r T ₁ or (1−r)B ₁/T ₁ via a simple sigmoid, \(f(x) = 0.5 (1 + \tanh (\frac {x-x_{0}}{g}))\) where the transition rate at x ₀ is \(f^{\prime }(x_{0}) = \frac {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 \(\frac {1}{2g}\). Now the function h=r f provides an output in (r,∞). If x is larger than the threshold x ₀, 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 ₀ does not activate the port P unless its concentrations is past some threshold [T ₀] _b where [T ₀] _b denotes the base concentration. Hence, we can model the port activity by \(h_{p}([T_{0}]) = \frac {r}{2} \left ({\vphantom {\frac {[T_{0}]-[T_{0}]_{b}}{g_{p}}}1 + \tanh }\right.\left.\left (\frac {[T_{0}]-[T_{0}]_{b}}{g_{p}}\right)\right)\) where the two shaping parameters g _p (transition rate) and [T ₀] _b (threshold) must be chosen. We can thus model the schematic of Fig. 4 as h _p ([T ₀]) [T ₁] _n where [T ₁] _n is the nominal concentration of the induced trigger T ₁. In a similar way, we let \(h_{e}(x) = \frac {s}{2}(1 + \tanh (\frac {x-x_{0}}{g_{e}}))\) Thus, for x = h _p ([T ₀]) [T ₁] _n , we have h _e is a switch from 0 to s. Note that 0 ≤ x ≤ r[T ₁] _n and so if h _p ([T ₀]) [T ₁] _n is close to r[T ₁] _n ,h _e is approximately s. Further, if h _p ([T ₀]) [T ₁] _n is small, we will have h _e is close to 0. This suggests a threshold value for h _e of \(\frac {r[T_{1}]_{n}}{2}\). We conclude

$${} {{\begin{aligned} h_{e} \left(h_{p}([T_{0}]) \: [T_{1}]_{n}\right) = \frac{s}{2} \left(~ 1 + \tanh\left(\frac{h_{p}([T_{0}])[T_{1}]_{n} -\frac{r[T_{1}]_{n}}{2}}{g_{e}}\right) \right) \end{aligned}}} $$

which lies in [0,s). This is the amount of activated T ₁ which reaches the genome to create the target protein P(T ₁). It follows then that [P(T ₁)]=h _e (h _p ([T ₀])[T ₁] _n )[T ₁] _n . The protein is created with efficiency e and so we model the conversion of [P(T ₁)] into a change in \(\boldsymbol {\mathcal {Q}}\) as follows. Let

$$\begin{array}{@{}rcl@{}} h_{\boldsymbol{\mathcal{Q}}}(x) &=& \frac{e}{2} \left(~1 + \tanh\left(\frac{x-x_{0}}{ g_{\boldsymbol{\mathcal{Q}}}}\right) ~ \right) \end{array} $$

which has output in [0,e). Here, we want to limit how large a change we can achieve in \(\boldsymbol {\mathcal {Q}}\). Hence, we assume there is an upper limit which is given by \(\Delta \: \boldsymbol {\mathcal {Q}} \: = \: \delta _{\boldsymbol {\mathcal {Q}}} \: \boldsymbol {\mathcal {Q}}^{max}\). Thus, we limit the change in the the expression of \(\boldsymbol {\mathcal {Q}}\) to some percentage of a baseline value. It follows that \(h_{\boldsymbol {\mathcal {Q}}}(x)\) is about \(\delta _{\boldsymbol {\mathcal {Q}}}\) if x is sufficiently large and small otherwise. This suggests that x should be [P(T ₁)] and since translation to P(T ₁) occurs no matter how low [T ₁] is, we can use a switch point value of x ₀=0. We conclude

$$\begin{array}{@{}rcl@{}}{} h_{\boldsymbol{\mathcal{Q}}}(\left[P(T_{1})\right]) = \frac{e}{2} \: \delta_{\boldsymbol{\mathcal{Q}}} \: \boldsymbol{\mathcal{Q}}^{max} \left(~1 + \tanh\left(\frac{\left[P(T_{1})\right]}{ g_{\boldsymbol{\mathcal{Q}}}}\right)\right) \end{array} $$

(3)

Our model of the change in \(\boldsymbol {\mathcal {Q}}\) expression is therefore \(\Delta \: \boldsymbol {\mathcal {Q}}^{max} \: = \: h_{\boldsymbol {\mathcal {Q}}}(\left [P(T_{1})\right ])\).

We can use these results for our general immune interaction discussion as well. From earlier comments, we know that associated with a change in fragility in the two cell populations M and N is the alteration of a protein or protein complex called P(T _M ) or P(T _N ) depending on the cell population type. Then we have

$$\begin{array}{@{}rcl@{}} h_{\boldsymbol{M}}(\left[P(T_{M})\right]) &=& \frac{e_{P(T_{M})}}{2} \: \delta_{P(T_{M})} \: \boldsymbol{P(T_{M})}^{max} \\&& \left(1 + \tanh\left(\frac{\left[P(T_{M})\right]}{g_{\boldsymbol{P(T_{M})}}}\right)\right) \end{array} $$

Our model of the change in maximum P(T _M ) expression is therefore \(\Delta \: \boldsymbol {P(T_{M})}^{max} \: = \: h_{\boldsymbol {P(T_{M})}}(\left [P(T_{M})\right ])\). We can do this for the N cells as well and obtain

$$\begin{array}{@{}rcl@{}} h_{\boldsymbol{N}}(\left[P(T_{N})\right]) &=& \frac{e_{P(T_{N})}}{2} \: \delta_{P(T_{N})} \: \boldsymbol{P(T_{N})}^{max}\\&& \left(1 + \tanh\left(\frac{\left[P(T_{N})\right]}{g_{\boldsymbol{P(T_{N})}}}\right) \right) \end{array} $$

Our model of the change in maximum P(T _N ) expression is therefore \(\Delta \: \boldsymbol {P(T_{N})}^{max} \: = \: h_{\boldsymbol {P(T_{N})}}([P(T_{N})])\).

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 \frac {\partial u }{ \partial x } \mid _{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 C a ⁺⁺ in the cell is controlled by the reversible binding of calcium ion to the buffer complexes, B ₁,B ₂ 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 ₁) 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 ₀ 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 \(\frac {\boldsymbol {M}}{P}\) which we call C _M and the fraction in the N population is \(\frac {\boldsymbol {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).

Notice that if we define a new diffusion coefficient, for the diffusion process that governs w by

$$\begin{array}{@{}rcl@{}} \mathcal{D} &=& D_{0} + (D_{M}- D_{0}) \frac{\partial C_{M} }{ \partial w } ++ (D_{N}- D_{0}) \frac{\partial C_{N} }{ \partial w } \end{array} $$

(7)

we obtain

from which can conclude we probably have oscillations if the algebraic sign of θ _3j is opposite to θ _2k and if the algebraic signs of θ _2j and θ _3k match. We also want θ _2j positive so we get undamped oscillations. The more exact equality conditions are probably not actually needed although the analysis was easier when we made them.

How do we use this information? Once we identify signals j and k useful for the dynamical model of interest, we need to experimentally estimate \(\frac {\delta \boldsymbol {j}}{ \boldsymbol {j}}\), \(\frac {\delta \boldsymbol {j} }{\boldsymbol {k}}, \frac {\delta \boldsymbol {k}}{\boldsymbol {j}}\) and \(\frac {\delta \boldsymbol {k}}{\boldsymbol {k}}\). This will give us estimates for the algebraic signs. We believe there will be an autoimmune interaction is the sign patterns we have discussed hold.

We consider this work essentially a theoretical model and we hope that we can generate some insights into the many troubling auto immune problems we face.