On a Mathematical Model of a General Autoimmune Disease

Abstract

The proposed paper is devoted to presenting and analyzing a kinetic model describing the development of autoimmune disorders. The proposed model is a nonlinear system of differential equations that considers the biological activity of the interacting populations. The main characteristics of autoimmune diseases are taken into account. Preliminaries to the research area are provided. The modeling problem is discretized and solved approximately. The numerical results illustrate typical outcomes of autoimmune diseases.

1. Introduction

Modeling using mathematical equations is critical in applied and basic science research. Moreover, it often leads to the formulation of new tasks and obtaining new scientific research results. In many areas, using mathematical and computational methods is crucial since their research subjects are often objects making it impossible or very difficult to perform natural experiments. Mathematical modeling of the studied thing allows experimenting with a tangible object to be replaced by purposeful research employing a mathematical model, thus allowing the experiment to be performed on the model using mathematical analysis and computer programs. Examining reality through abstract objects, appearing as models of high accuracy, is the essence of mathematical modeling.

Thus, mathematics, from a means of calculation, becomes a necessary modern research method. At certain stages it can be the only means of revealing the internal properties of the considered objects. Therefore, mathematical modeling applications can be defined as a practical approach to scientific research [1,2].

The current Special Issue can be given as an example of the great interest and activity of scientists in mathematical modeling applications in various scientific fields. It includes, among others, exciting applications in physics, topology, ecology, epidemiology [3], and immunology [4].

In particular, the paper [4] considers a mathematical model describing the development of HIV infection. The central interacting populations are healthy T-cells, infected T-cells, and free HIV viruses. Their dynamics is described by ordinary differential equations (ODE). The authors perform a stability analysis of the model. For the numerical solution, they derive a modified continuous Galerkin–Petrov algorithm. The obtained results are compared with other numerical schemes.

Medicine is among the areas where mathematical modeling has been successfully applied in recent decades. Here, we propose a new model for the development of autoimmune diseases. Our goal is to analyze the role of some critical parameters of the model on the corresponding disorders. The remaining part of our paper is structured in the following sections. Some preliminaries describing the main immunological concepts are given in Section 2. Section 3 reviews some interesting recent applications of mathematical modeling in medicine and other fields. The modeling framework is introduced briefly. The proposed mathematical model is described in Section 4. Section 5 deals with numerical experiments. Section 6 is devoted to conclusions and future scientific plans.

2. Basic Immunological Concepts

Human bodies are complex systems in which their cells and organs interact with foreign microorganisms. These interactions can be mutually beneficial, but they can also cause harm. Examples of valuable functions of various microorganisms, such as digestion, the participation in the production of some vitamins and foods, the detoxification of some harmful chemicals, and the strengthening of some immunological processes can be mentioned [5,6].

On the contrary, various pathogenic microscopic organisms, including viruses, bacteria, parasites, and fungi, can harm the host organisms. Examples of hazardous infectious diseases are Ebola virus hemorrhagic fever, cholera, polio, smallpox, human acquired immunodeficiency syndrome (HIV/AIDS), coronavirus disease (COVID-19), and many others [7].

The defensive mechanisms in human organisms include external barriers, innate and acquired components, and secretions against foreign microorganisms. Collectively, they are called the immune system [8].

One of the basic immunological concepts is that harmful foreign pathogens can directly attack cells, tissues, and organs but also affect the body’s protective functions and mechanisms [5,8]. One of the consequences is that the cause of death may not be an infection but the immune response to it [5,8].

Therefore, an essential property of a correctly functioning immune system must be its ability to discriminate between the cells and molecules of the host organism and foreign pathogens. The immune system should attack only pathogenic substances. Such an immune system is called tolerant to self [9].

The concept of “self-nonself” is related to the clonal selection theory of M. Burnet and D. Talmage, the idea of the immunological homunculus of I. Cohen, and other theoretical concepts in immunology [10].

When an immune system is not tolerant to its substances, it can start working against self components, leading to autoimmune disorders and diseases [11].

The concept of auto-antibodies is related to the idiotype–anti-idiotype network theory of N. Jerne and other theoretical ideas by P. Ehrlich, etc. [10].

The reasons for failures of proper functioning of the protection system and the occurrence of autoimmune disorders can differ: environmental factors, stress, genetics, diets, hygiene conditions (both bad and sound), etc. [12,13,14]. Viruses have also been considered to be among the factors in the occurrence of autoimmune diseases [15,16,17]. Strangely, some observations showed that sometimes viral infections could have a protective role from autoimmunity. To explain this, the so-called “hygiene hypothesis” was invented [9].

Examples of autoimmune diseases, counting more than 80 diseases, are paroxysmal cold hemoglobinuria, lupus, multiple sclerosis, and rheumatoid arthritis [18].

3. Mathematical Modeling in Medicine

Scientists use various research methods in treatment and health care. Clinical studies, laboratory investigations, qualitative health research, and quantitative methods are among them. Quantitative approaches apply statistical analyses, economic experiments, and computational and mathematical modeling methods to perform systematic and rigorous empirical studies [19].

One of the most popular frameworks for the mathematical modeling of diseases is using systems of ordinary differential equations. Some recent models describe the development of autoimmune conditions [20,21,22,23,24].

Another interesting approach is the modeling using fractional differential equations allowing the consideration of memory effects [25,26,27,28,29,30,31]. During the last decade this approach has been successfully used for modeling the development of such diseases as HIV, COVID-19 and toxoplasmosis etc. [28,29,31].

Long-distance travel has significantly increased for professional and tourist purposes in recent years. This change also leads to the possibility of the rapid spread of infections between people whose places of residence are located at great distances. A mathematical model describing such long-distance spread of diseases is proposed in the paper [27]. The concept of superdiffusion is used to model such spreading. The authors present the basic notions and properties of fractional derivatives of noninteger order. Several ways of representing superdiffusion using fractional calculus are reviewed and compared concerning their ability to model the spread of infections over large distances.

The paper [29] is devoted to a new mathematical model describing the development of toxoplasmosis disease in cats and people. The condition is caused by parasites called Toxoplasma gondii. They can infect cats, and humans can also be infected after interactions with cats. The proposed mathematical model is described by a system of fractional differential equations. The authors provide essential information related to calculus using derivatives of fractional order. They describe the mathematical model in detail. The non-integer order paradigm is reviewed and applied to the analysis of the properties of the model. Qualitative analysis of the model is carried out, particularly stability analysis and sensitivity analysis concerning the primary threshold parameter of the model. The qualitative results are supported by the results of simulations performed with the use of a modified Adams–Bashforth–Moulton approach.

Another application of the fractional order calculus can be found in [30], which studies the development of tuberculosis. The disease is caused by Mycobacterium tuberculosis. Detailed information about the disease and its relation to other chronic conditions such as diabetes mellitus, HIV etc. is given. The authors provide a comparison between conventional and non-integer calculus and the possible advantages of each of them. Various areas of successful applications of fractional differentiation and integration are mentioned. The authors describe the basic concepts and properties of fractional derivatives. Further, they introduce the fractional order differential model describing the development of tuberculosis. The stability analysis of the model is performed. The conditions for existence and uniqueness of the model solutions are studied. Numerical simulations are carried out by applying a new predictor–corrector scheme. They confirm the qualitative results.

An application of fractional calculus to modeling the interactions between cancer and the immune system is proposed in the paper [25]. The authors recall the primary factors and populations involved in tumor development and the host organism’s attempts to eradicate the disease. The authors mention the notion of the Caputo–Fabrizio derivative and describe its main properties. The mathematical model describing the interaction is presented and studied numerically applying a modified Adams–Bashforth approach. The role of the main parameters of the model is checked for the occurrence of chaos and oscillations.

The paper [26] is devoted to applying fractional calculus to study the development of autoimmune diseases. The authors extend the model proposed in [32] by applying the Caputo fractal-fractional operator. The problem of existence and uniqueness of the solution is studied. The researchers apply a numerical scheme using Lagrange interpolation to obtain approximate solutions. The role of various parameters of the model is studied numerically. The results of the simulations illustrate typical outcomes of autoimmune diseases.

Our paper uses a different approach: we apply the methodology used in statistical mechanics and generalized kinetic theory for active particles. This approach also considers the functional features of the interacting cells and particles. Additional variables perform the description of these states of activity. Examples of recent results using this methodology can be found in the papers [20,21,22,23,33].

The considered system is divided into several of the most important populations of interacting individuals. Using concepts from the corresponding field of immunology, evolution equations for the concentrations of these populations are described including related gain and loss terms. They are integro-differential equations of the Boltzmann type. Further, qualitative and/or quantitative analysis can be performed.

In Section 4, we introduce our mathematical model describing autoimmune diseases.

4. The Mathematical Model

By our study, we would like to investigate the meaning of the viral particles for developing autoimmune diseases.

Various mechanisms are known to cause autoimmune diseases, which are counted around 100 conditions today. Among these mechanisms is molecular mimicry, which can be caused by multiple viral infections [5,9]. Molecular mimicry can occur when peptides of some agents look like the host’s. This phenomenon may elicit immunological reactions against the foreign pathogens but also against components of the host possessing a similar protein structure or shape. In such a way a cross-reaction of the defense system with host antigens occurs, which may destroy the host’s constituents. The destroyed compartments may release additional isolated antigens that can be further attacked by the immunological cells. In this way we can observe continuous production of antigens that may become targets for the immune cells. Additionally, as a result, the autoimmune response may become solid and autodestructive.

In our mathematical model we consider the competition between four important populations:

(i)

The first population of target cells;

(ii)

The second population of damaged cells;

(iii)

The third population of immune cells;

(iv)

The fourth population of viral particles causing mimicry.

According to the mechanism of molecular mimicry, we suppose that some viruses activate immune responses destroying some healthy (target) cells.

We suppose that the organism can produce healthy cells that become target cells. We use a logistic-like gain term to model the rate of proliferation of the target cells.

Considering the most important features of the interactions between viruses and immune cells, we formulate four equations (one for the dynamics of each population) governing the system’s time evolution under consideration.

As a simplification of the reality, we suppose in our modeling process that the first two subsystems denoted by the labels $i=1$ and $i=2$ (corresponding to healthy (target) and destroyed cells) are unstructured; i.e. homogeneous concerning the activity. The populations denoted by $i=3$ and $i=4$ (corresponding to immunological cells and viruses) can function specifically in a non-homogeneous manner; i.e. they are structured populations. Their activities are represented by the variable $u\in [0;1]$. The meaning of the activities is as follows:

• Immune cells: The activation state refers to their ability to destroy target cells.
• Viruses: The activation state refers to their ability to trigger the production of immunological cells. Increasing values of the activation state denote higher amounts of newly produced immunological cells.

The states of the unstructured populations $i=1$ and $i=2$ are characterized, at time t, by their concentrations

$$n_i\left(t\right),i=1,2.$$

(1)

The states of the structured populations $i=3$ and $i=4$ are described by the probability densities

$$f_i(t,u),i=3,4.$$

(2)

From the values of the probability densities it is possible to determine the concentrations of the structured populations as follows:

$$n_i\left(t\right)={\int }_0^1{f}_i(t,u)du,i=3,4.$$

(3)

Our model is a generalization of previous models [32,34,35] describing the development of autoimmune conditions. These models consider the populations of healthy (target) cells, damaged cells, immunological cells and viral particles.

The model [32] is described using only ordinary differential equations. Therefore, the activation states of the populations are ignored.

The models [32,34] include the state of activity of immune cells. However, these models study the role of different model parameters.

In the present model, the activity of the virus particles is additionally introduced to study more aspects of the role of the pathogenicity of the virus infections for the autoimmune processes.

The modeling of the microscopic interactions between active particles of the various subsystems is performed by taking into account the variation rates of the amounts of active particles with corresponding inlet flux and outlet flux rates caused respectively by proliferative interactions (gain terms, GT) and outlet flux rates caused by destructive interactions (loss terms, LT) included in the general form of kinetic Equation (4).

Thus, we obtain a specific representative of the general form of a kinetic type model as proposed in many foundational papers (e.g., [36]) in the following way:

$$\begin{array}{ccc}\frac{\partial f_i}{\partial t}(t,u)\hfill & =\hfill & AI(t,u)+GT(t,u)-LT(t,u),\hfill \end{array}$$

(4)

where $AI(t,u)$ stands for artificial influx of cells with activity u at time t.

As a result, the following system of ordinary and partial integro-differential equations is obtained:

$$\begin{array}{ccc}\frac{d}{dt}{n}_1\left(t\right)\hfill & =\hfill & S_1\left(t\right)+n_1\left(t\right)\left(p_{11}-\frac{p_{11}}{T_{Max}}{n}_1\left(t\right)\right)\hfill \\ \hfill & \hfill & -d_{11}{n}_1\left(t\right)-d_{13}{n}_1\left(t\right){\displaystyle \underset{0}{\overset{1}{\int }}}u{f}_3(t,u)du,\hfill \end{array}$$

(5)

$$\begin{array}{ccc}\frac{d}{dt}{n}_2\left(t\right)\hfill & =\hfill & d_{13}{n}_1\left(t\right){\displaystyle \underset{0}{\overset{1}{\int }}u{f}_3(t,u)du}-d_{22}{n}_2\left(t\right),\hfill \end{array}$$

(6)

$$\begin{array}{cc}\frac{\partial f_3}{\partial t}(t,u)\hfill & =\hfill \end{array}(1-u)\left(p_{32}{n}_2\left(t\right)+p_{34}{\displaystyle \underset{0}{\overset{1}{\int }}v{f}_4(t,u)dv}\right)-d_{33}{f}_3(t,u),$$

(7)

$$\begin{array}{ccc}\frac{\partial f_4}{\partial t}(t,u)\hfill & =\hfill & p_4{f}_4(t,u)-d_{44}{f}_4(t,u)-d_{43}{n}_3\left(t\right)f_4(t,u)\hfill \end{array}$$

(8)

with non-negative initial conditions

$$n_i\left(0\right)=n_i^{\left(0\right)},i=1,2,f_i(0,u)=f_i^{\left(0\right)}\left(u\right)i=3,4.$$

All parameters of the system (5)–(8) are assumed to be non-negative.

Equation (5) describes the dynamics of the concentration $n_1\left(t\right)$ of the target cells. The meaning of the participating parameters is the following:

• $S_1\left(t\right)$ models the production rate of target cells from sources within the organism;
• $p_{11}$ characterizes the rate of proliferation of the target cells;
• $T_{Max}$ refers to the concentration of the target cells at which their proliferation turns off;
• $d_{11}$ describes the rate of the natural death of the target cells;
• $d_{13}$ describes the damaging rate of the target cells attacked by the immune cells.

Equation (6) models the time evolution of the concentration $n_2\left(t\right)$ of the damaged cells. The parameter $d_{22}$ is the natural death rate of the damaged cells.

Equation (7) describes the evolution of the distribution density $f_3(t,u)$ of the immune cells. The meaning of the participating parameters is the following:

• $p_{32}$ is the rate of production of immune cells due to the presence of self-antigens presented by damaged cells;
• $p_{34}$ is the rate of production of immune cells due to the presence of viruses;
• $d_{33}$ is the rate of the natural death of the immune cells;
• The factor $1-u$ is related to the assumption that the state of activity of the newly produced immune cells is low and they need time for activation.

Equation (8) describes the evolution of the distribution density $f_4(t,u)$ of the viral particles. The meaning of the participating parameters is the following:

• $p_4$ is the rate of production of viruses;
• $d_{43}$ is the rate of destruction of the viruses due to the immune response;
• $d_{44}$ is the rate of natural death of the viruses.

Let us introduce the following notation:

$$X=\{f=(n_1,n_2,f_3,f_4):|n_1|<\infty ,|n_2|<\infty ,f_3\in L_1(0,1)f_4\in L_1(0,1)\},$$

$$X^{+}=\{f=(n_1,n_2,f_3,f_4)\in X:n_1\ge 0,n_2\ge 0,f_3\ge 0,f_4\ge 0,a.e.\}.$$

It can be shown that the following theorem holds.

Theorem 1.

Let $S_1\in C^0([0,\infty );R_{+}^1).$ For each $T>0$, there exists a unique solution

$$f\in C^0([0,T];X)\cap C^1((0,T);X)$$

to the system (5)–(8) with the initial datum $f^{\left(0\right)}=\left(n_1^{\left(0\right)},n_2^{\left(0\right)},f_3^{\left(0\right)},f_4^{\left(0\right)}\right),f^{\left(0\right)}\in X^{+}$. The solution satisfies $f\left(t\right)\in X^{+},\forall t\in [0,T].$

Proof. 

It is easy to check that the right-hand sides of the system (5)–(8) possess the property of Lipschitz continuity in the space denoted by X. This guarantees that the solution exists and is unique locally. Further, by using the method of successive approximations, it can be shown that the solution is non-negative. Appropriate a priori estimations of the solution can be easily found. They are sufficient to have a unique solution on every finite temporary interval.    □

5. Results of Simulations

The numerical algorithm for obtaining the approximate solution to the Cauchy problem corresponding to (5)–(8) includes the discretization of Equations (7) and (8) with respect to the activation states $u\in [0,1]$ by uniform meshes. The included integrals are approximated by using the composite Simpson’s rule.

The obtained system of ordinary differential equations has been solved for various parameter sets by the use of the solver ode15s from the Matlab ODE suite [37].

Due to the complexity of the considered nonlinear systems of ordinary differential equations obtained after the performed discretization, it is not possible to find analytically exact solutions of the models. For this purpose, Matlab’s ode15s solver based on the backward differentiation formulae proposed by Gear [38] was used. The specificity of using such types of solvers is that the solver itself controls the accuracy of each step of the integration. It chooses the integration step so that the user-specified accuracy can be achieved. The values of relative tolerance RelTol and absolute tolerance AbsTol that we used are $RelTol=1.{10}^{-3}$ and $AbsTol=1.{10}^{-4}$. They provide an accuracy of 0.001 at any point in time. Such accuracy guarantees a good applicability and usability of the numerical results presented here.

It should be noted that in the course of the work a comparison was made between the results obtained by the program in the Matlab environment and those obtained by the multistep predictor–corrector method of Adams combined with the one-step Euler and Runge–Kutta methods, guaranteeing an accuracy of up to 0.001.

The assumption for the initial conditions has been that there is a certain amount of target cells, a small amount of damaged cells, and small amounts of immune cells with equally distributed states of activity and viruses with high states of activity:

$$n_1\left(0\right)=100,n_2\left(0\right)=0.5,$$

$$f_3(0,u)=0.1,\forall u\in [0,1],$$

$$f_4(0,1)=0.1,f_4(0,u)=0.0,\forall u\in [0,1).$$

The numerical experiments have been performed with the following values of parameters:

$$S_1\left(t\right)=10,\forall t\ge 0,T_{Max}=100000,$$

$$p_1=0.5,d_{11}=0.2,d_{22}=1.1,$$

$$p_{32}=0.001,p_{34}=0.5,d_{33}=0.1,$$

$$p_4=10.0,d_{43}=4.0,d_{44}=4.0,$$

and various values of the parameter $d_{13}$.

The numerical experiments have been aimed at studying the effects of the damaging ability of the immune cells activated against target self-cells by viral infection.

In Figure 1 and Figure 2 the concentrations of target cells and damaged cells for small values of the parameters $d_{13}=0.011$ and $d_{13}=0.015$ are shown. This situation corresponds to weak immune system self-aggressivity when target cells remain at high levels and autoimmune disorders are not observed or are very mild.

In Figure 3, Figure 4 and Figure 5, respectively, results for higher values of the parameters $d_{13}=0.10$, $d_{13}=0.25$, and $d_{13}=0.30$ are shown. This situation corresponds to strong self-aggressivity of the immune system when damping, periodic or increasing oscillations of the concentrations of target cells and damaged cells are observed. The cases of damping oscillations can be considered autoimmune disorders tending to some equilibrium states representing chronic diseases. Examples of chronic autoimmune diseases are type 1 diabetes and chronic immune thrombocytopenic purpura related to infection with the hepatitis C virus, which has been shown to be able to evade the immune response and to evolve into a chronic state in many infected patients [24].

The cases with periodic and significantly increasing flare-ups can be considered autoimmune diseases characterized by repeated periods of dormancy of the autoaggressive immune cells and their flare-ups triggered by spreading pathogenic proteins, possible failure of regulatory immune cells to control the destructive reaction of effector immune cells against host constituents etc.; an example of repeated flare-ups is MS [24]. The numerous repeated periods when large amounts of target cells are damaged can lead to severe tissue or organ destruction and even fatal results for the patients.

6. Concluding Remarks

This paper presented and studied a mathematical model formulated within the the kinetic theory framework for active particles describing autoimmune diseases. We aimed to present some quantitative and qualitative features of these common diseases. The obtained results from our simulations represent quite characteristic temporary courses of autoimmune diseases that occurred with viruses’ active participation (among them flare-ups, chronic disorders, healthy state etc.).

In our future work we will investigate the role of additional parameters of the model. Further, we will try to develop the modeling approach considering other vital populations’ activity. A significant problem is validating the model by comparing the results with actual clinical data.