The flowchart shows how the populations of each of the 6 types of cells effect each other. Because each individual cell has a chance of changing into another type of cell at every unit of time, for the model to be steady the number of cells that stop being type X must also equal the number of cells that become type X.
For example, lets say that there are 5000 Ens cells at the beginning of a unit of time. If during that unit of time 400 cells decide to become Ins (the only type of cell an Ens can become) then 400 cells of some other type have to decide to become Ens for the population of Ens cells to remain at 5000. According to the flowchart those new cells will come from the Ins population

The Deltas

So how does this flowchart become a system of equations? The chart shows us a path of population changes, or the net change in population.
Example: The change in the anticipatory population (As) is equal to the number of new anticipatory cells, minus the number of anticipatory cells that die. For self, the new anticipatory cells is given by Is times the probability of an Is cell becoming anticipatory (or k3), and the number of anticipatory cells that die is given by k6 times the number of anticipatory cells.
In an equation this statement looks like this:
As' = (k3 * Is) - (k6 * As)
Where Is, and As are current populations of initial and anticipatory cells respectively, and As' is the new calculated population of anticipatory cells.

Now we have a system of equations for the net change of each population:

Is'  = (SI * N) + (k4 * Es) - (k5 * Is) - (k3 * Is) - (k2 * Is)
Es'  = (k2 * Is) - (k4 * Es)
Ins' = ((1 - SI) * N) + (k4 * Ens) - (k5 * Ins) - (k2 * Ins)
Ens' = (k2 * Ins) - (k4 * Ens)
As'  = (k3 * Is) - (k6 * As)
Note that N is equal to k1*G in the graph.

We also need an equation to represent the steady size of the entire system. This is accomplished by setting the number of cells entering the system (being created by the stem cells) equal to the number of cells leaving the system (by dying). N is the variable the represents the number of cells produced per unit time.

N  = (k5 * Is) + (k5 * Ins) + (k6 * As)
The N variable is also unknown, but it is different from other cell populations because this value represents the number of cells that are born and that die in the system, for each unit of time, when a steady state is reached.  So we can get a formula for the "change in N per unit time" which is fixed to 0.
N' = (k5 * Is) + (k5 * Ins) + (k6 * As) - N
Thus we have
Is'  = (SI * N) + (k4 * Es) - (k5 * Is) - (k3 * Is) - (k2 * Is)
Es'  = (k2 * Is) - (k4 * Es)
Ins' = ((1 - SI) * N) + (k4 * Ens) - (k5 * Ins) - (k2 * Ins)
Ens' = (k2 * Ins) - (k4 * Ens)
As'  = (k3 * Is) - (k6 * As)
N'   = (k5 * Is) + (k5 * Ins) + (k6 * As) - N
Although we how have 6 equations with 6 variables, we notice that these equations are not linearly independent.  In short,  -Is' =  Es' + Ins' + Ens' + As' + N'.

Another important input of the model is the total system size or Protecton size (T).  This is introduced by the following equation, which effectively scales all the cell groups by the correct amount to have an immune system of size T.

Is + Es + Ins + Ens + As = T
This is now a system of seven equations (six of which are linearly independent) with six unknowns. For simplicy at this point, we will drop the Is' equation, since it can be formed from linear combinations of other equations.

The k values and total system cell population are decided by the user and the model represented by equation, once solved (see "The Steady, Stable State"), can be used to obtain the steady state values of the populations for each cell type. This can be performed quickly by a computer.  The web interface allows for the user to scan ranges of values for any variable, creating a 2D graph. Furthermore it allows the user to scan 2 variables at the same time, to obtain a 3D graph. All the while, raw numbers of the selected values are also returned.

A Steady, Stable State

For the model to have a steady state we would need to have population values for the 6 types of cells that when a system is started at these levels, they will stay at these levels as long as there is no outside interference. The first requirement of a steady, stable state is that the net change of each population equal zero. So our equation for the previous section become...
0 = (k2 * Is) - (k4 * Es)
0 = ((1 - SI) * N) + (k4 * Ens) - (k5 * Ins) - (k2 * Ins)
0 = (k2 * Ins) - (k4 * Ens)
0 = (k3 * Is) - (k6 * As)
0 = (k5 * Is) + (k5 * Ins) + (k6 * As) - N
0 = Is + Es + Ins + Ens + As - T
When the variables are separated from the known values we can form two matrices one known and one unknown.  Then we can use Maple (TM) to solve for the unknown matrix. Like so...
k2*Is -k4*Es 0*Ins 0*Ens 0*As 0*N = 0
0 0 (-k2-k5)*Ins k4*Ens 0 (1-SI)*N 0
0 0 k2*Ins -k4*Ens 0 0 0
k3*Is 0 0 0 -k6*As 0 0
k5*Is 0 k5*Ins 0 k6*As -1*N 0
1*Is 1*Es 1*Ins 1*Ens 1*As 0 T

k2 -k4 0 0 0 0 * Is = 0
0 0 (-k2-k5) k4 0 (1-SI) Es 0
0 0 k2 -k4 0 0 Ins 0
k3 0 0 0 -k6 0 Ens 0
k5 0 k5 0 k6 -1 As 0
1 1 1 1 1 0 N T
When Maple(TM) is given this system of equations and told to solve for Is, Es, Ins, Ens, As & N with T, SI, k2-k6 being an input constant. We get these results:
Is =  -T * SI * k5 * k4 * k6
Ins =  T (-k3 - k5 + SI * k5 + SI * k3)* k6 * k4
Ens =  T * k2 * (-k3 - k5 + SI * k5 + SI * k3) * k6
Es =  -T * k2 * SI * k5 * k6
N T * k6 * k5 * k4 * (-k5 - k3)
As =  -T * K * k5 * k4 * k3
Where D = (-k5*k6*k2 - k6*k4*k3 - k6*k2*k3 - K*k5*k4*k3 + K*k6*k4*k3 + K*k6*k2*k3 - k5*k4*k6).

So the system is solvable as long as D doesn't equal zero and all of the populations are positive, for the given K values. This solution will give us an instant calculation for the steady state, but how do we know that it is a stable steady state? And what's the difference between a steady state and a stable steady state? A stable steady state is one that is resistant to long term change. A pendulum is a stable system because, baring outside interference, the pendulum will always come to the same state, straight down. A coin balanced on it's edge is an unstable state because if an outside force knocks the coin over, the coin will never be able to get up on it's own.

The Eigenvalues

Note: this section requires an understanding of linear algebra.
An eigenvalue is a number, y, that corresponds to a vector, x, and a matrix, A, such that A*x = y*x.
A homogeneous system of equations has non-zero solutions if and only if the coefficient matrix (y*I - A) = 0 is singular; that is, if and only if the determinant is zero.
The system of equations will equal the zero vector if Is = Es = Ens = Ins = As = N = 0, but we are not interested in an empty system. So we have to show that there is a non zero solution.
Because of the nature of this problem (we are applying this matrix of k values an infinite number of times) the eigenvalues cannot be negative, that would result in an alternating pair of values that doesn't converge.
Finding the determinant of this matrix is best done by computer, and even then the solution isn't usable. Maple gives the determinant as a 5th degree polynomial, the simplest coefficient of which is (k6 + k3 + 2*k2 + 1 + 2*k5). Again not solvable by hand, but since the k values are all between 0 and 1 we can run a simulation to see if negative eigenvalues will be produced.
100 randomly picked sets of k values didn't produce a single negative eigenvalue, so we can feel somewhat confident that the model is stable.

The Application

This is how the simulator works.
  1. Get the immune system size (T), number of APCs per peptide, repertoire size (R), fraction of new cells which are anti-Self (SI), and the various rate constants (k2 to k5) from the user.
  2. Start with all cell populations set to 0, except for the stem cells (N).
  3. Use the solutions for the system of questions above, obtained from Maple(TM), to quickly obtain the steady state values for each cell type, given the user's inputs.

  4. i.e. Is = (-T*SI*k5*k4*k6)/D, and N = (T*k6*k5*k4*(-k5-k3))/D, etc.
  5. Display the predicted steady state values to the user.
  6. For variables that are graphed, (1 variable gives a 2D graph, 2 variables results in a 3D graph), when this option is selected, there are a minimum of 100 samples taken per decade. Steps number 2 and 3 above are performed for each new sample and then used to create the graph, which is then returned to the user via the web page.