Using Mathematical Modelling For Autoimmunity Biology Essay

Published: Last Edited:

This essay has been submitted by a student. This is not an example of the work written by our professional essay writers.

Mathematics is used a lot in the world that average folk are blissfully unaware of. Everyone goes through the education system learning mathematics in the early years as counting from 1 through to 100, and then later on, in their teens, as trigonometry and random equations. From the education we gain as a youngster, it's hard to imagine mathematics to be useful in many areas of employment. This is untrue. Mathematics is everywhere, working behind the scenes. One of these unidentified areas is Medicine.

Mathematic is used extensively to help gain advances in medical research; from the Pharmacology of new drug concoctions to the understanding of the mechanics of complex diseases. Mathematicians are working with medical researchers every day to understand how different triggers affect the human body, creating some form of mathematical equation or system to show this, and then to create a mathematical system to predict how adverse effects can be counteracted.

Immunology is a complex area of Medicine and Biology. Rather than the molecular and cellular study, it deals with complex non-linear biological systems. Although there are an increasing number of attempts to use and develop computational software, these often don't consider the nature of the interactions between the various cellular and molecular components nor do they usually give much biological insight into how the system works. Mathematical modelling is used to predict tumour growth and cancer spread, where as another branch of mathematics, statistics are used to interpret data collected from clinical trials.

As the title suggests, this paper will be concentrating on the mathematical modelling used to aid research into Autoimmune Diseases. Firstly, there will be a fairly brief overview of what autoimmunity is, followed by the basic modelling. The modelling will start with the basic linear models, where a brief insight will be given to their properties. Towards the end of the paper, non-linear models will be presented with a much deeper insight to the workings of the model.

Biological Background

So, what is an Autoimmune Disease? An Autoimmune disease is when the body's immune system becomes overactive or 'confused' and starts to produce a response against the person's healthy cells. In essence, attacking and damaging them. There are two different responses the immune system can take; it can produce antigens to cause direct damage to a single organ or tissue causing a 'localised' autoimmune disease, or produces a response against multiple tissues and organs resulting in a 'systemic' autoimmune disease.

Basic Biological Terms

Here is the definition of the basic biological terms used throughout this paper.

Target Cell - the body's healthy cells

Damaged Cell - a cell that has been damaged, i.e. by a virus or the body's immune response.

Antigen - a molecule which when released into the body activates the immune response, producing antibodies to attack what it recognises as a potential harmful invader.

Antibody - used by the immune system to recognise and destroy potential invaders.

Key Effector Cell - A cell that performs the body's reaction to a stimulus

APC's - The cell that process the antigens to allow T cells to recognise them

T Cells - A set of cells belonging to a group of immune cells (white blood cells) called lymphocytes

"in vivo" - where experimentation has been done in live isolated cells as opposed to a whole organism.

Vicious Cycle of Autoimmunity

Before continuing, the cycle of autoimmunity needs to be explained, in brief.


(5)Damaged cell


Target cell

Key effector cell




Lymph Vessel




Fig 1. A Cycle of Autoimmunity

Dynamical Properties of autoimmune disease models: Tolerance, flare-up, dormancy

(Iwami et al., 2007)





When autoimmune disease is initiated by some event, an existing cell becomes a key effector cell (1). Then, this key effector cell attacks and damages a healthy cell (2). Here, the protein of the damaged cell (antigen) (3) is captured by an APC (4) and the protein is shown as a self-antigen at the lymph vessel (5). Then the immune cells which are specific to the protein are induced (6), and these specific immune cells again attack and damage target cells, (7) starting the cycle again (8).


Autoimmune diseases generally don't just have singular symptoms, or necessarily a specific test for diagnosis, therefore it can be hard for medical practioners to provide a diagnosis. Also many symptoms can overlap from one disease to another, as the body reacts in similar ways; it can be hard for them to establish which disease the patient is suffering from.

There can be many potential triggers to activate the disease activity, and also to determine the severity. Many common life components; like the weather, stress etc.; can cause a completely non-linear reaction to a sufferers immune response.

Therefore, because of the many varying symptoms and triggers, there ranges a large number of different drug and physical therapies. Through medical literature there exists plenty of papers on the subject of autoimmune disease and it's complications, yet there is limited theoretical work in mathematical literature. However, over the years, mathematical models have been proposed in relation to T cell vaccination (Iwami et al., 2008), the behaviour of the immune network and the activation of immune cell activation, for example. For the latter, there is a paper (Wodarz and Jansen, 2003) investigating the ratio of cross-presentation to direct-presentation of antigen presenting cells assuming their number is variable. However these models are fairly specific to certain autoimmune trends.

Current understanding of Autoimmunity is that tissue injury causes T cell reactions and the activation of some immune cells. Consequently, this leads us to assume that patient's autoimmune symptoms are based on the population size of healthy (target) cells. Hence, the smaller the population size, the more severe the symptoms. In the following model we will assume that the number of antigen presenting cells (henceforth known as APC's) is constant. More specifically, the Dendritic Cells (DC's), a type of APC, are known to carry out almost all antigen presentations to T cells and the number of these are said to be constant in vivo.

Throughout the remainder of this paper, we will investigate a simple model for general autoimmune disease.

Basic model

In the following model we will be showing the effects of;

The immune system cells, C , which naturally die at the rate γ

The healthy (target) cells, T , which naturally die at the rate μ and are

produced at a rate of λ

The damaged cells, D , which naturally die at the rate α ;

where α ˃ μ

The immune system cells, C, damage the target cells, T, at a rate proportional to their population, βTC. Here, β represents the efficacy of the process.

From the above background, a basic dynamical model is obtained, which is a system of non-linear differential equations. In order to simplify the model we assume that damaged cells already exist and disregard the dynamics of the key effector cells and APC's; these could have been caused by an earlier viral or bacterial infection, or tissue injury. We combine the dynamics of immune cells and target cells and obtain a model as follows:

Target cell growth - force of damage by immune cells

Force of damage by immune cells - natural rate of death of damaged cells

Personal immune response - natural rate of death of immune cells

The immune system is a lot more complicated than our simple mathematical model above based on the target cell growth function and the personal immune response function , but this will be adequate for the moment. includes the rates that immune cells find and succeed in attacking target cells. The effects of target cell growth and the personal immune response will be investigated.

Target Cell Growth Function

The above mathematical model is based on a common model used to understand HIV infection. Reasonable functions, to represent target cell growth in humans have been investigated. (Martin A. Nowark and Robert M. May, 2000; Martin A. Nowark et al., 1996; Alan S. Perelson and Patrick W. Nelson, 1998)

The first is a simple equation which is just , the rate at which the new target cells are produced, minus (rate of death of target cells x population size of target cells) which was investigated by Martin A. Nowak et al. (2000; 1996)

The second equation is slightly more complex, (Alan S. Perelson and Patrick W. Nelson, 1998), with the addition of an extra logistical term to take into account natural target cell proliferation. Here, represents the maximum proliferation rate and L; the target cell population density at which proliferation shuts off.

It is noted that the above equation is density dependent. However, an alternative equation has been investigated by Liancheng Wang and Michael Y. Li (2006), in order to include density dependence, and is as follows:

However, for mathematical simplicity and without changing the qualitative behaviour of the system, we will assume that the first equation of is the functional form including density dependence.

Personal Immune Response Function

If the personal immune response function is defined as , the relationship between the immune cell inducement and the symptoms of autoimmune disease can be investigated. Immune response functions can vary from person to person or they may depend on a patient's condition or the kind of immune cells. Immune cell inducement is considered to be a reasonable function shown as:

Fig 2. Personal Immune Response

Dynamical Properties of autoimmune disease models: Tolerance, flare-up, dormancy,

Iwami et al. (2007)

The above diagram shows APCs do not produce immune cells if only a few antigens exist, but then immune cells are gradually induced when relatively many antigens exist.

Henceforth, the paper will investigate the behaviour of two estimations of the personal immune response functions;


And Non-Linear:


Linear Personal Response Function

In this section, the linear personal response function and it's representation of autoimmune disease symptoms will be investigated, i.e. , where represents the average magnitude of activation of Immune Response by APCs per damaged cell.

Here, the number of immune cells induced by APC's is proportional to the number of damaged cells, so D can be considered as the proliferation rate of immune cells by APCs.

Linear Target Cell Growth Function

To start, we select the simplest equation for the population dynamics of target cells; . This gives the following simple model of autoimmune disease dynamics:


The above model is the same as a basic HIV model studied by several researchers (Nowak and May, 2000; Leenheer and Smith, 2003; and others), but only the implications for autoimmune disease will be explored in this paper. Numerical solutions were found for this system of differential equations with parameters (Iwami et. Al, 2007)

Now, if is small then the damaged cells, D, vanish and the immune cells, C, are not activated hence the target cells, T, do not deplete. However, if is large then the number of damaged cells, D, increases activating the immune cells, C, and depleting the number of target cells, T. Hence, the larger is, the more likely a patient is to develop an autoimmune disease.

This system has two equilibria:

The tolerance equilibrium, , where ,

The chronic infection equilibrium,

Next, the basic reproductive number is defined as, which shows how many newly damaged cells produced from one damaged cell in an individual who does not have an autoimmune disease currently. It has been shown that is globally asymptotically stable, if , by De Leenheer and Smith (2003). If , under certain conditions is globally asymptotically stable.

The argument of continues here. The chronic infection equilibrium coordinates, from direct calculation, are given by:

Taking the derivative of and , we get;

As increases there is a transferal from tolerance to autoimmune disease activity from the result of . From acute symptoms the disease transfers to a chronic phase. Here, increasing deteriorates the symptoms shown by the results of , and . Namely, increasing activates the immune cells, and reduces the number of target cells. Hence, , which represents the average magnitude of activation of the immune response, affects the severity of the autoimmune disease attack and it's symptoms. It is a person's responsiveness to a trigger or treatment.

Density-Dependant Target Cell growth Function

In the following section the density-dependant function is used to obtain a similar model, which is also a HIV model that has been previously researched (Wang and Li, 2006; Nevo et al., 2003; and others), once more only the autoimmune disease implications are considered.

Again, (Iwami et al., 2007), numerical solutions of this system where found with parameters . This leads the model to appear as;

Iwami et al (2007) demonstrated the following values of ;

As , target cells reproduce themselves. The different values of determine the state of the disease;

gives a tolerance of the immune response resulting in no autoimmune disease,

gives a slow progression of the disease resulting in mild symptoms for the patient: The target cells decrease gradually, however, in the chronic phase, there is still a relatively high level of target cells.

gives repeated flare-ups of an autoimmune disease.

gives a more severe reaction than (ii). There is a quick progression of the disease resulting in the rapid decrease of target cells, and subsequently, a low level of target cells in the chronic phase.

These results show that the values of in (i) and (ii) are similar to the previous results with the linear target cell growth function ; hence, the higher the value of the more severe the disease symptoms for the patient. The result for (iii) behaves interestingly. It shows a periodic pattern for the disease symptoms, which relate to repeated flare-ups of the autoimmune diseases.

The reason for this flare-up pattern in (iii) can be explained as the following; if is relatively small and the number of target cells is small, then the number of target cells increases because of the nature of (i.e. the logistical term: ) . Then there is an associated increase in damaged cells and immune cells. The immune cells attack target cells can cause a decrease in the number of target cells and the cycle repeats.

In (iv), target cells cannot be increased any more due to the logistic nature of , if is large. Therefore, the effect of density-dependant target cell growth dramatically changes the symptoms of autoimmune disease.

Mathematically, this system is more interesting than the previous one, as a positive equilibrium point could be unstable. This system has been investigated by De Leenheer and Smith (2003) and also by Iwami et al. (2007). It has two equilibria:

The tolerance equilibrium, ,

where ,

The chronic infection equilibrium,

The basic reproductive number is defined as

De Leenheer and Smith (2003) show that if , then is globally asymptotically stable and if then is only globally asymptotically stable under certain conditions but can also be unstable under certain conditions.

Non-Linear Personal Response

In the subsequent section, the non-linear personal immune response function, , and it's relationship with the symptoms of autoimmune disease is to be investigated. In ecological population dynamics, this function is called the "functional response", or more accurately a Holling type III or sigmoid functional response.

Here, the parameters of the Functional Response are defined;

= maximum proliferation rate of immune cells caused by APCs

= number of damaged cells at which the proliferation of immune cells is half of the maximum

Thus, the function can be regarded as the proliferation rate of immune cells by APCs.

This non-linear response function intensely changes the structure of the equation systems. Namely, both the tolerance and chronic infection equilibriums can be simultaneously stable under certain conditions. Specifically, is always stable.

Linear Target Cell Growth Function

As with the linear personal response function, we start with the linear target cell growth function . Using this, the model changes to the following;

Now, with (i.e. the target cell's natural rate of death is higher than the immune cell's natural rate of death), the implications for autoimmune disease suggested by the model in this section are considered. Numerical solutions are investigated with parameters: , (Iwami et. al, 2007)

Substituting the parameter values into the template of the model we get;

Then using the initial conditions:

gives the following results;

gives a representation of dormancy of an autoimmune disease. In the beginning, the target cells (T) remain at a fairly steady large level, while the immune cells (C) remain steady at a low level. At some point, without some evolutionary event (i.e. without changing parameters), the number of target cells suddenly falls to a fairly small level while the amount immune cells increases. This result suggests that for an autoimmune disease to become apparent, evolutionary events are not essential.

gives a representation of tolerance within the immune system, though this is different to the previous results as the parameters stay the same, only changing the initial values. As long as the initial value of damaged cells is small, this system will always show a tolerance.

Here is given a stability analysis for this system. This model has three equilibria:

, where

Hereafter, a detailed explanation of the stability of these equilibria is given. To find the stability of these equilibria, the eigen-values of a Jacobian matrix for each of these equilibria coordinates are investigated.

The formula for the Jacobian matrix for this differential equation system is;

Hence, the Jacobian matrix of the first equilibrium point is as follows;

Consequently, the eigen-values of this matrix are and . And so, is always stable.

Moving onto the Jacobian matrix for , the coordinates of which are found to be;

The characteristic equation of this matrix is found as;

In order to make this equation simpler, the coefficients for the indeterminate of the polynomial, s, are denoted by,

All eigen-values have negative real parts if;

is clearing true, so the final two conditions need to be proven.

By using the coordinates for , is the same as


Then if and only if

which implies that,

which is interpreted as for . However, if it exists, is always unstable.

Finally, expanding the works of Iwami et al. (2007) the last condition is investigated:

So, if , then and is always stable when it exists.

Henceforth is the argument for . Using the coordinates, we get:

Therefore, this is similar to the previous systems with regards to the shift of the severity of symptoms, where m is large. On the other hand, the fact this system is bi-stable gives the possibility for a shift to disease tolerance with no matter how large is. If the number of damaged cells is small, this system will always show tolerance.

Density-Dependant Target Cell Growth Function

In this subsection, the second representation of g from Section 3 is chosen, i.e.

changing the basic model to the form:

From the paper by Iwami et al. (2007), numerical solutions were found with the parameters Hence, the model becomes;

To carry out some simulations of this model, different sets of initial values are used, for different values of:

Mathworks© Matlab is used in order to gain detailed results.

[Matlab Codec in Appendix A]


Time = t

Target Cells , T ,

Damaged Cells , D ,

Immune Cells , C ,











































Target Cells, T,

Damaged Cells, D,

Immune Cells, C,











































Target Cells, T,

Damaged Cells, D,

Immune Cells, C,











































Target Cells, T,

Damaged Cells, D,

Immune Cells, C,









































Using the above tables, and diagrams, it is easy to see the patterns caused by the variants in initial values and in addition the variation of . A.(i) shows a dormancy phase of the disease followed by repeated 'flare-ups'. The level of immune cells starts of at a consistently low rate during the dormancy phase, but as this ends, suddenly there is a significant drop in the steady high rate of target cells, with a then large increase of the number of immune cells. However, this pattern of autoimmunity does not then shift the symptoms to chronic phase, but instead, shows the repeated flare-ups.

A.(ii) shows a tolerance of autoimmune disease. It can be explained as if only a few antigens exist the APC's do not induce the activation of immune cells, and no autoimmune disease develops.

B.(i) shows a dormancy phase to begin but, where the number of immune cells, with extreme slowness, increases. When this finishes, the target cells suddenly decrease and the immune cells suddenly increase. As the number of antigens increases, and the target cells decrease, the APC's induce the activation of the immune cells, causing the shift from tolerance to disease activity, where this model shows the symptoms settling into a chronic state. B.(ii) again shows a consistent tolerance of autoimmune disease.

The above figures show two things;

is connected to the mechanism of repeated flare-ups of autoimmune disease.

is connected to the mechanism of dormancy and attack of autoimmune disease.

Now is given a stability analysis for this model. From reviewing the model, it is apparent that a boundary equilibrium exists, i.e.

Next, the eigen-values of the Jacobian matrix for are investigated. Hence the Jacobian matrix is as follows;

The eigen-values of the above matrix are . Inserting the equation for , and simplifying, . Hence, all eigen-values are negative and so is always stable.

The existence of interior equilibria still needs to be shown. Using the equations for gives;

Now substituting the equation for D into C' gives;

Furthermore, using the equation for , a representation for C is obtained, as follows:

Substituting this equation for C back into the previous equation, gives;

The positive roots of this equation are the T element.

The following graphs from "Dynamical properties of autoimmune disease models, by Iwami et al (2007)" give a visual representation of

and , and


Fig. 4

Dynamical properties of autoimmune disease models

The roots of are shown in fig 4.(ii) by the intersection of the line and the curve . Here, it is shown that two interior equilbiria exist denoted as and on the diagram.

The larger root (as shown in the diagram) corresponds to the equation:


Consequently, an interior equilibria, , exists near to .

Next, allowing to be the maximum value that satisfies the derivative , there exists two equilibria where .

Using the parameters given at the start of section 5.2, the boundary equilibrium is always stable for all four cases. Hence, there is a tolerance to autoimmune disease if only a few antigens exist. However, in the cases for A. , the interior equilibrium and are both unstable. Alternatively, is stable for the cases of B.. The stability of represents whether the autoimmune disease shows repeated flare-ups, or chronic symptoms. From observing the patterns of A.(i) and B.(i), the significant difference the parameter makes is apparent. Repeated flare-ups and chronic symptoms are reliant on the density-dependant growth of target cells. The effect of the density-dependant growth is strong when is relatively small resulting in the repeated flare-ups of disease symptoms, which remains unchanged6 when the parameter is changed within a proper biological range. Otherwise when is relatively large, the effect of density-dependant growth is weak, which therefore results in the autoimmune disease presenting chronic symptoms.


Throughout this paper, there have been multiple differentiations of a basic model for autoimmune disease. Starting off with the linear interpretations and ending with non-linear functions.

The personal immune response function has been shown to encourage tolerance or dormancy. For the linear function , the size of the variable determines whether there is tolerance of autoimmune disease, or whether it will be active. There is always tolerance if is small, otherwise Autoimmune disease becomes active for large. Hence, with a linear personal immune response, only a few antigens cause APC's to present and bring about the activation of immune cells.

The non-linear function determines dormancy and tolerance of auto-immune disease. It is dependent on the initial population of antigens (damaged cells, D). If the level of damaged cells is small in the beginning, then there is a tolerance of autoimmune disease, in dependent of the size of . On the other hand, if the initial condition of damaged cells is not small, then it induces a dormancy of autoimmune disease, which eventually results in autoimmune disease activity. Hence, with a non-linear personal immune response, APC's do not present for few antigens.

It turned out that the personal immune response function is closely related to autoimmune disease. Every patient is different, and has a different immune response, and hence displays different symptoms. Therefore, in order to gain an accurate model for a particular patient for therapy of autoimmune disease, that person's immune response would need to be investigated.

The target cell growth function has been shown to cause repeated flare-ups of autoimmune disease. especially shows repeated flare-ups of the disease. If , then presents no repeated flare-ups, just chronic symptoms.

If the level of target cells is small, the multiplication of the logistic term of , [], can increase their number. Additionally, an increase of damaged and immune cells can be related to an increase in target cells. This then has the effect of the immune cells 'attacking' the target cells, which decreases their number. Hence, the pattern of repeated flare-ups.

Using the function , keeping small then the target cell growth function can increase the number of target cells and hence repeated flare-ups are observed. However if is large, is unable to increase the population of target cells, and no flare-ups of symptoms occur.

With the function , changes to have no effect on whether flare-ups present or not. However, if the parameter is set to a small value, the density-dependant has a strong effect, and hence flare-ups of symptoms occur.


The differences in the discussion suggest that it is the target cell growth function that determine if flare-ups occur. According to Iwami et al. (2007), it is thought that "the difference of target cell growth functions corresponds to the differences between internal organs which initiate autoimmune disease.". So, similar to the personal immune response function, each individual internal organ's target cell growth needs to be investigated in order to gain an accurate model for therapy of autoimmune disease.

The Personal Immune Response, is just that: personal. One specific function will not represent the population. However, empirical studies help to provide typical responses in order to give an idea of an immune response. Data is collected from a multitude of candidates by medical professionals and passed to mathematicians. The mathematicians then interpret it and create viable models in order to represent the phenomena of Autoimmune disease. (Jo Stevenson, 2010)

These models can also be used to help understand the effect of drug therapy used for autoimmune disease patients. For example, immunosuppressant's. These drugs, as suggested in their name, work in order to reduce the immune response from a patient. In the model investigated by this paper, these drugs effect the personal response functions, , by reducing the parameters .

As already suggested, for tolerance with the immune response , needs to be a very small value. Hence, unless the drugs can succeed in this, there would be little effect from them. However, for , drug therapy may be more effective regardless of the size of .


Appendix A


function dy = nlddtcgfA(t,y)

dy = zeros(3,1); % a column vector

dy(1) = 0.1-0.1*y(1)+3*y(1)*(1-y(1)/100)-0.5*y(1)*y(3);

dy(2) = 0.5*y(1)*y(3)-1.1*y(2);

dy(3) = ((10*y(2)^2)/((60^2)+y(2)^2))-0.1*y(3);



function dy = nlddtcgfB(t,y)

dy = zeros(3,1); % a column vector

dy(1) = 5-0.1*y(1)+3*y(1)*(1-y(1)/100)-0.5*y(1)*y(3);

dy(2) = 0.5*y(1)*y(3)-1.1*y(2);

dy(3) = ((10*y(2)^2)/((60^2)+y(2)^2))-0.1*y(3);


Matlab Codec > Command window


>> [T,Y]=ode45(@nlddtcgfA,[0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95],[100 4 0])

>>[T,Y]=ode45(@nlddtcgfA,[0 100],[100 4 0]);



>> [T,Y]=ode45(@nlddtcgfA,[0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95],[100 3 0])

>>[T,Y]=ode45(@nlddtcgfA,[0 100],[100 3 0]);



>> [T,Y]=ode45(@nlddtcgfB,[0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95],[100 4 0])

>>[T,Y]=ode45(@nlddtcgfB,[0 100],[100 4 0]);



>> [T,Y]=ode45(@nlddtcgfB,[0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95],[100 3 0])

>>[T,Y]=ode45(@nlddtcgfB,[0 100],[100 3 0]);


Appendix B