Abstract
The Marchuk model of infectious diseases is considered. Distributed control to make convergence to stationary point faster is proposed. Medically, this means that treatment time can be essentially reduced. Decreasing the concentration of antigen, this control facilitates the patient’s condition and gives a certain new idea of treating the disease. Our approach involves the analysis of integro-differential equations. The idea of reducing the system of integro-differential equations to a system of ordinary differential equations is used. The final results are given in the form of simple inequalities on the parameters. The results of numerical calculations of simulation models and data comparison in the case of using distributive control and in its absence are given.
1. Introduction
Distributed feedback control system is a challenging problem, but only a few papers were devoted to it (see, for example, [1] and literature therein). A noise in the feedback delay control is one of the main reasons for developing mathematical models with distributed inputs. The second one is resulted from the fact that the dynamics of the processes rather depends on the average value of the process than on the value at a corresponding moment.
We propose a method reducing the stability analysis of the system with memory to the analysis of a corresponding finite-dimensional system. This idea was proposed, as well as we know, in [2] and realized, in finite spectrum assignment (see, for example, [3–6]). Our realization of this idea is different. We propose a simple method that allows to reduce the analysis of systems with memory to one of a higher order but finite-dimensional systems. Questions of stability and estimates of solutions can be considered on this base.
The model of infectious diseases, constructed by Marchuk in his well-known book [7], reflects the most significant patterns of the immune system acting during these diseases. This model was studied in many works; note, for example, the recent papers [8–12] and the bibliography therein. The adding control was proposed, for example, in [11–14]. In [10, 15], the basic mathematical model that takes into account the concentrated control of the immune response was proposed. It can be noted that the use of information about behavior of a disease and the immune system for a long time (defined by distributed control, for example, in the form of an integral term) looks very natural in choosing strategy of a possible treatment. Optimal control in the basic model of the infectious diseases was considered in the work [15], where the control function characterizing the realization of an immunotherapy which consists in administration of ready immunoglobulins or donor antibodies is proposed. In [16], the model of influence of an immunotherapy on dynamics of an immune response which represents a generalization of basic model was considered. On the basis of the proposed model, the problem of determination of coefficients on the basis of laboratory data was considered and also the management was proposed in [13, 14, 17]. Such task is called control in uncertain conditions [11]. The control algorithm in the case of uncertain conditions was proposed in the work ([15], see pages 71–73).
The Lyapunov exponents characterize the rate at which solutions approach each other with increasing time. In our case, one of these solutions is a certain stationary state of the system. The stationary state determines the conditions of recovery of the body after exposure to infection. The speed of approaching this state has a very important role. It allows us to estimate the duration of treatment. In some cases, this can determine the choice of treatment strategy. It is enough to imagine a situation when the patients’ body is weakened by some chronic diseases and is not able to endure too long treatment.
Our approach to treatment modeling can be described as follows. We introduce a distributed feedback control in the form of an integral term. This management is based on long-term monitoring of the patient’s condition and comparing the patients’ condition with a certain stationary state over a long period of time. Thus, a closed controlled system, the purpose of which is to improve the patient’s condition, is created. Intuitively, such a treatment strategy is understandable. Further, the scheme of reducing the integro-differential system to a system of ordinary differential equations is proposed. This ordinary differential system is studied. The technique based on properties of the Cauchy matrices is used to estimate possible deviations of solutions in the case of various modifications of this model, such as nonlinear models or delay models, from the solution of this system ordinary differential equation. Note the papers [1, 18–21] where the distributed control was also considered for various problems.
In this paper, we consider the Marchuk model of infectious diseases [7]:where is the antigen concentration rate, is the plasma cell concentration rate, is the antibody concentration rate, and is the relative features of the body. Denote as corresponding coefficients described in [7]. takes into account the destruction of the normal functioning immune system, . Denote as the plasma rate concentration of the healthy body and the antibody concentration that we wish to achieve after the treatment.
Let us discuss every equation in the model (1) in more detail. The first equation presents the block of the virus dynamics. It describes the change in the antigen concentration rate and includes the amount of antigen in the blood. Antigen concentration decreases as a result of interaction with antibodies. The immune process is characterized by the number of antibodies, whose concentration is described by the equation . The value of decreases as a result of the interaction and mutual destruction with antigens. The amount of the antibody cells also decreases as a result of natural destruction. However, the plasma restores antibodies and therefore the plasma state plays an important role in the immune process. Thus, the change in the plasma cell concentration rate is included in several equations of the system. Taking into account the healthy body level of plasma cells and their natural aging, the term is included in the second equation of system (1). It is assumed that the plasma is restored as a result of the interaction of the antigen and the antibody cells. The second and third equations present the block of the humoral immune response dynamics.
Concerning the last equation of system (1), we can note the following. The value of m increases with the antigen’s concentration rate . The maximum value of m is 1 in the case of organ damage and is equal to 0 for a fully healthy organ. The coefficient describes the rate of degeneration of the target organ.
The main goal of this paper is to demonstrate new ideas in the use of distributed control in the model of infectious diseases [7]. Our goal is to make convergence to a stationary point faster. This allows us to decrease the duration of disease’s treatment. We add a distributed control in the equation describing antibody concentration rate. Numerical simulations demonstrate that this approach could open some new possibilities for a treatment.
2. Modeling Distributed Control
Let us consider the following system:
Comparing this system with the model of Marchuk (1), we see the control added in the third equation describing the antibody concentration rate. Note that this control is a reasonable one from the medical point of view [10, 15].
We choose the control in the following form:
It can be noted that influence of a corresponding average value instead of at the point t looks reasonable since the control is rather dependent on an “average” value of the difference on a corresponding time-interval than on this difference at the moment t only. The integral term in (3) increases the influence of the previous moments which are closer to the current moment t.
Following [15], we can pass to the dimensionless case, substituting in (2). In the dimensionless case, system (2) is of the following form:
Consider the system
Denoting in (5)we obtain
Let us find possible stationary points of system (7). Their coordinates satisfy the following algebraic system:
Only in two cases, the first equation is satisfied: or . Let us start with the first case. If , then from the second equation and from the fourth one. We have the system for f and :
From the last equation, we have , and substituting this into the first one, we obtain , and then . Thus, starting with , we can come only to the stationary point . In the second case of , we consider only the case of , since in the case of , we come to the same stationary point. From the fourth equation , we obtain that . This means that there is no possibility to completely disinfect the affected organ. That is why we study below the behavior of solutions only in the neighborhood of the noted above stationary point .
Remark 1. It was obtained in [15] on the basis of the laboratory data that .
It was noted above that . Linearizing systems (7) and (4) in the neighborhood of the stationary point, we obtain the corresponding linear systems:where , respectively. Denote the matrix of the coefficients of system (10):Note that a corresponding linear system for the Marchuk model of infectious diseases without control (see system (1)) can be written in the form:
3. About Comparison of the Lyapunov Exponents
Our approach is based on the following auxiliary assertion.
Lemma 1. The solution-vector of system (4) and 4 first components of the solution-vector of system (5) satisfying the condition coincide.
Proof of Lemma 1. Let us consider the last equation of system (5). Using the representation of solution of the first-order scalar equation, we can write . Note that system (5) is considered with the condition . Substituting now this representation of into the third equation of system (5), we obtain the third equation of system (4). Thus in systems (4) and (5), the first, second, and fourth equations coincide and the third equation of (4) is equivalent to the last equation of system (5) with the initial condition .
This completes the proof of Lemma 1.
Remark 2. It is clear that Lemma 1 can be used also for integro-differential system (11) and system of ordinary differential equations (10). Every solution-vector of integro-differential system (11) coincides with the first 4 components of the solution-vector.
Take of system (10) with the initial condition . Thus, if we take a part of the space of solutions of ordinary differential system (10), satisfying the initial condition and delete then the fifth component of solution-vectors, we come to the 4-space of solutions of integro-differential equation (11). Now, it is clear that the exponential stability of the 5-dimensional ordinary differential system (10) implies the exponential stability of 4-dimensional integro-differential system (11). The roots of the characteristic equation for ordinary differential system (10) can be used in the representation of vector-solution in the space of solutions of integro-differential system (11). Negativity of real parts of the roots of the characteristic equation of system (10) with constant coefficients allows to make the conclusion about the exponential stability of the integro-differential system (11). The root with maximal real part of the characteristic equation for (11) cannot be greater than the one of the characteristic equations for (10).
The characteristic polynomial of system (10) of ordinary differential equations,has 5 roots .
Denote as the roots of the characteristic polynomial of systems (13), and .
Theorem 1. If , and , then integro-differential system (11) is exponentially stable, and if in addition, the inequality is fulfilled then .
The proof is based on the following assertion.
Lemma 2. For the roots of the characteristic polynomials of systems (11) and (13), the following facts are true:and if , then
Proof of Lemma 2. The characteristic polynomial of systems (13) isand the characteristic polynomial of system (10) is (14).
It is clear from the definition of the coefficients (6) thatThe inequality in the conditions of Theorem 1 implies that and consequently . It is clear thatIn the case of complex roots and , we have if . In the case of real roots and , we have . Thus we obtain the exponential stability of system (10). According to Remark 2 mentioned above, system (11) is exponentially stable.
Comparing the values of the roots of the characteristic polynomials (17) and (14) we have to verify only the inequalityin (16). Substituting the values of and to (20), we have to obtain the inequalities:for the case of real roots , andfor the case of the complex roots .
The inequalities (21) and (22) are fulfilled if and .
Substituting the values of in Remark 1, we see that . If the inequalities are fulfilled, then (16) is true.
Note the assertion following from Theorem 1.
Corollary 1. Let the coefficients be defined by the formulas in Remark 1, and then the Lyapunov exponents of system (11) are less than the ones of system (13).
4. Numerical Simulations and Comments
The mathematical model can be represented in the following form:where is the right-hand side of (7). Values of the parameters are defined in Remark 1, and the initial conditions are . We use the second-order Runge–Kutta’s scheme with a constant step h.
In Figures 1–4, the solution of model with the natural flow of data without the control of disease is presented by curves of red color and disease in the case of considered distributed control by curves of green color.
Figure 1 demonstrates the dynamics in antigen concentration during the course of the disease. The insert detailing the process in the first two days was performed on a different scale and demonstrates the fact that the management transfers the disease from the acute form to the subclinical one (the antigen concentration only decreases after injection). Figure 2 demonstrates the dynamics in plasma cell concentration during the disease process. It can be seen from the figure that control leads to a faster increase in the concentration of plasma cells, which in this case ensures a transition to the subclinical form of the disease. In addition, it is necessary to note a fourfold increase in the maximum concentration of plasma cells in the case of control, compared with the option without control. Figure 3 demonstrates the dynamics in antibody concentration during the disease process. The graph shows that the concentration of antibodies in the solution with control practically does not change, because in this case, they are replaced by donor antibodies, which is what the control actually consists of. The dynamics in the proportion of target organ cells destroyed by antigen during the disease process is presented in Figure 4. The values for the variant with control are given with an increase of times. Thus, control allows to reduce the maximum proportion of affected cells of the target organ by more than times. The dynamics of the control (concentration of donor antibodies) during the disease process is presented in Figure 5.
Remark 3. The integral term can accumulate small mistakes made in the process of numerical integration. This explains difficulties in the use of numerical methods for solving integro-differential equations. The idea to reduce a system of integro-differential equations to a corresponding system of ordinary differential ones and the use of the well-developed technique for their solution could be one of the possible ways to develop numerical methods for integro-differential equations.
5. System with Delay in Equation of Plasma Concentration
It was explained in [7] why the delay can appear in the second equation of system (1) in the model of infectious diseases:
We have to define what should be set instead of and when . We consider all delay systems with the initial function for . Let us consider the following system with the control in antibody concentration:where .
We can pass to dimensionless case
Denoting , according to (6), we obtain
For stability studies, only behavior of solutions for sufficiently large t is considered. Below, we assume that . We can write the second equation in the following form:
We can write the expression under the integral in the following form:
Linearizing system (27) in the neighborhood of the stationary pointwe obtain
Consider the systemwhere is a -matrix and is n-vector. The general solution can be represented in the following form:where -matrix is called the Cauchy matrix. Its j-th column () for every fixed s as a function of t is a solution of the corresponding homogeneous system:satisfying the initial conditions , where
Construction of the Cauchy matrix of system with ordinary differential equations can be found, for example, in [22].
Let us denote .
Theorem 2. If andwhere is the Cauchy matrix of system (10), then system (31) is exponentially stable.
Proof of Theorem 2. Consider the nonhomogeneous system corresponding to homogeneous system (31):where are measurable essentially bounded functions on the semiaxis (i.e., they are from the space ).
We can write system (37) in the following form:whereand A is defined by (12). It is clear that the use of the standard formula of solutions’ representation leads to the following system:It follows from Theorem 1 that the Cauchy matrix satisfies the exponential estimate, i.e., there exist such positive γ and N, then all elements of the Cauchy matrix satisfy the following estimate:Now it is clear that the vector-function is bounded for .
Define the operator ( is the space of 5-dimensional vector-functions with continuous components) by the following formula:System (40) can be rewritten in the formwhereThe inequality (36) implies that the norm of the operator is less than one, then there exists the operator and it is bounded. Thus, the solution-vector is bounded for , for every bounded on the semiaxis right-hand side . Now the exponential estimate of the solution of the homogeneous system (31) follows the Bohl–Perron theorem [23, 24]. We have proven the exponential stability of system (31).
Remark 4. It is clear that system (31) is exponentially stable if for sufficiently small delay .
Remark 5. Let us compare the difference of solution of system (10) and solution of system (31) with the same initial conditions . It is clear thatDenote . satisfies the following equation:The solution can be represented in the following form:It follows from the exponential estimates of and that tends to zero exponentially, when .
In order to estimate , we can use the results of the paper [25], where the elements of the Cauchy matrix were constructed.
6. Construction of the Cauchy Matrix and Stability of Model with Delay
In order to use Theorem 2, we have to obtain the estimates of . In this section, we present these estimates.
Let us denote ,
Lemma 3 (see [25]). If then the second column of the Cauchy matrix of (10) is as follows:Let us denote
Lemma 4 (see [25]). f then the second column of the Cauchy matrix of (10) is as follows:Let us denote
Lemma 5 (see [25]). If then the second column of the Cauchy matrix of (10) is as follows:The following assertions are results of substitutions of presented in Lemmas 3–5 into inequality (36) in Theorem 2.
Theorem 3. If , , andthen system (31) is exponentially stable.
Theorem 4. If , , andthen system (31) is exponentially stable.
Theorem 5. If , , andthen system (31) is exponentially stable.
Data Availability
The data used to support the findings of this study are available within the article.
Disclosure
This paper is part of the first author’s Ph.D. thesis which is being carried out in the Department of Mathematics at Ariel University.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this article.