Distributed Control and the Lyapunov Characteristic Exponents in the Model of Infectious Diseases

e 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-dierential equations. e idea of reducing the system of integro-dierential equations to a system of ordinary dierential equations is used. e nal results are given in the form of simple inequalities on the parameters. e results of numerical calculations of simulation models and data comparison in the case of using distributive control and in its absence are given.


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.e 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 nite-dimensional system. is idea was proposed, as well as we know, in [2] and realized, in nite spectrum assignment (see, for example, [3][4][5][6]).Our realization of this idea is di erent.We propose a simple method that allows to reduce the analysis of systems with memory to one of a higher order but nite-dimensional systems.Questions of stability and estimates of solutions can be considered on this base.e model of infectious diseases, constructed by Marchuk in his well-known book [7], re ects the most signi cant patterns of the immune system acting during these diseases.is model was studied in many works; note, for example, the recent papers [8][9][10][11][12] and the bibliography therein.e adding control was proposed, for example, in [11][12][13][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 (de ned 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 in uence 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 coe cients 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].e control algorithm in the case of uncertain conditions was proposed in the work ( [15], see pages 71-73).
e 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.e stationary state determines the conditions of recovery of the body after exposure to infection.e 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. is management is based on longterm monitoring of the patient's condition and comparing the patients' condition with a certain stationary state over a long period of time.us, 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 integrodifferential system to a system of ordinary differential equations is proposed.is ordinary differential system is studied.e 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][19][20][21] where the distributed control was also considered for various problems.
In this paper, we consider the Marchuk model of infectious diseases [7]: where V(t) is the antigen concentration rate, C(t) is the plasma cell concentration rate, F(t) is the antibody concentration rate, and m(t) is the relative features of the body.Denote α, β, c, ρ, σ, η, μ c , μ f , and μ m as corresponding coefficients described in [7].ζ(m) takes into account the destruction of the normal functioning immune system, ζ(0) � 1. Denote C * as the plasma rate concentration of the healthy body and F * the antibody concentration that we wish to achieve after the treatment.Let us discuss every equation in the model (1) in more detail.e first equation dV/dt � βV(t) − cF(t)V(t) 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.e immune process is characterized by the number of antibodies, whose concentration is described by the equation dF/dt � ρC(t) − η cF(t)V(t) − μ f F(t).e value of F(t) decreases as a result of the interaction and mutual destruction with antigens.e 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.us, 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 μ c (C(t) − C * ) 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.e second and third equations present the block of the humoral immune response dynamics.
Concerning the last equation dm/dt � σV(t) − μ m m(t) of system (1), we can note the following.e value of m increases with the antigen's concentration rate V(t).e maximum value of m is 1 in the case of 100% organ damage and is equal to 0 for a fully healthy organ.e coefficient μ m describes the rate of degeneration of the target organ.
e 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.
is 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.

Modeling Distributed Control
Let us consider the following system: Comparing this system with the model of Marchuk (1), we see the control u(t) 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 u(t) in the following form: It can be noted that influence of a corresponding average value instead of F(t) − F * at the point t looks reasonable since the control u(t) is rather dependent on an "average" value of the difference F(s) − F * on a corresponding timeinterval than on this difference at the moment t only.e 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 (2).In the dimensionless case, system (2) is of the following form: Denoting in ( 5) Let us find possible stationary points of system (7).eir coordinates v, s, f, m, and  u satisfy the following algebraic system: Only in two cases, the first equation is satisfied: v � 0 or f � α 1 /α 2 .Let us start with the first case.If v � 0, then s � 1 from the second equation and m � 0 from the fourth one.We have the system for f and  u: From the last equation, we have  u � (f − 1)/k, and substituting this into the first one, we obtain f � 1, and then  u � 0. us, starting with v � 0, we can come only to the In the second case of f � α 1 /α 2 , we consider only the case of v ≠ 0, since in the case of v � 0, we come to the same stationary point.From the fourth equation α 6 v � α 7 m, we obtain that m ≠ 0. is means that there is no possibility to completely disinfect the affected organ.at is why we study below the behavior of solutions only in the neighborhood of the noted above stationary Remark 1.It was obtained in [15] on the basis of the laboratory data that α 1 � 0.25;α 2 � 8.5000332;α 3 � 1792175675; α 4 � 1.95992344 • 10 − 7 ;α 5 � 0.5;α 6 � 10;α 7 � 0.4;α 8 � 1.7 • 10 − 3 .It was noted above that ζ(0) � 1. Linearizing systems ( 7) and ( 4) in the neighborhood of the stationary point, we obtain the corresponding linear systems:

Complexity
where 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:

About Comparison of the Lyapunov Exponents
Our approach is based on the following auxiliary assertion.
Proof of Lemma 1.Let us consider the last equation of system (5).Using the representation of solution of the firstorder scalar equation, we can write Note that system ( 5) is considered with the condition  u(0) � 0. Substituting now this representation of  u(t) into the third equation of system (5), we obtain the third equation of system (4).us 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  u(0) � 0. is 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 (x 1 (t), x 2 (t), x 3 (t), x 4 (t)) of integro-differential system (11) coincides with the first 4 components of the solution-vector.
Take (x 1 (t), x 2 (t), x 3 (t), x 4 (t), x 5 (t)) of system (10) with the initial condition x 5 (0) � 0. us, if we take a part of the space of solutions of ordinary differential system (10), satisfying the initial condition x 5 (0) � 0 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).
e 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).e root with maximal real part of the characteristic equation for (11) cannot be greater than the one of the characteristic equations for (10).
e characteristic polynomial of system (10) of ordinary differential equations, , and λ * 5 .Denote λ i , i � 1, 4 as the roots of the characteristic polynomial of systems (13), and (11) is exponentially stable, and if in addition, the inequality k e 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: Proof of Lemma 2. e characteristic polynomial of systems ( 13) is and the characteristic polynomial of system ( 10) is (14).
Comparing the values of the roots of the characteristic polynomials ( 17) and ( 14) we have to verify only the inequality in (16).Substituting the values of λ * 5 and λ 4 to (20), we have to obtain the inequalities: for the case of real roots λ * 4 and λ * 5 , and for the case of the complex roots λ * 4 and λ * 5 .e inequalities ( 21) and ( 22) are fulfilled if k > α 4 and b > 0.
Substituting the values of α i in Remark 1, we see that λ 4 > max 1≤i≤3 λ i .If the inequalities k > α 4 and b > 0 are fulfilled, then ( 16) is true.
Note the assertion following from eorem 1.

□
Corollary 1.Let the coefficients α i be defined by the formulas in Remark 1, and k > 1.95992344 • 10 − 7 , b > 0, then the Lyapunov exponents of system (11) are less than the ones of system (13).

Numerical Simulations and Comments
e mathematical model can be represented in the following form: where Φ(x(t)) is the right-hand side of (7).Values of the parameters are defined in Remark 1, k � 4 and b � 1 and the initial conditions are v(0) � 10 − 6 , f(0) � 1, s(0) � 1,m(0) � 0 and u(0) � 0. 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.
e 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.e 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.e dynamics in the proportion of target organ cells destroyed by antigen during the disease process is presented in Figure 4. e values for the variant with control are given with an increase of 10 4 times.us, control allows to reduce the maximum proportion of affected cells of the target organ by more than 2.5 × 10 4 times.e dynamics of the control (concentration of donor antibodies) during the disease process is presented in Figure 5.

Remark 3.
e integral term can accumulate small mistakes made in the process of numerical integration.is explains difficulties in the use of numerical methods for solving integrodifferential equations.e idea to reduce a system of integrodifferential 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.

System with Delay in Equation of Plasma Concentration
It was explained in [7] why the delay τ(t) 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 V(t − τ(t)) and F(t − τ(t)) when t − τ(t) < 0. We consider all delay systems with the initial function V(ζ) � 0 for ζ < 0. Let us consider the following system with the control in antibody concentration: 6 Complexity where Denoting α i (i � 1, . . ., 8), according to (6), we obtain For stability studies, only behavior of solutions for sufficiently large t is considered.Below, we assume that t − τ(t) ≥ 0. We can write the second equation in the following form: We can write the expression under the integral in the following form: (29) Linearizing system (27) in the neighborhood of the stationary point

Complexity
Consider the system where P(t) is a (n × n)-matrix and G(t) is n-vector.e general solution X(t) � col x 1 (t), . . ., x n (t)   can be represented in the following form: where n × n-matrix C(t, s) � C ij (t, s)   n i,j�1 is called the Cauchy matrix.Its j-th column (j � 1, . . ., n) for every fixed s as a function of t is a solution of the corresponding homogeneous system: satisfying the initial conditions x i (s) � δ ij , where Construction of the Cauchy matrix of system with ordinary differential equations can be found, for example, in [22].
Proof of eorem 2. Consider the nonhomogeneous system corresponding to homogeneous system (31): where g 1 (t), . . ., g 5 (t) are measurable essentially bounded functions on the semiaxis (i.e., they are from the space L ∞ ).We can write system (37) in the following form: where and 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 eorem 1 that the Cauchy matrix C(t, s) satisfies the exponential estimate, i.e., there exist such positive c and N, then all elements C ij (t, s) of the Cauchy matrix C(t, s) satisfy the following estimate: Now it is clear that the vector-function Define the operator K : C 5 ⟶ C 5 (C 5 is the space of 5dimensional vector-functions x : [0, ∞) with continuous components) by the following formula: System (40) can be rewritten in the form where e inequality (36) implies that the norm of the operator K : C 5 ⟶ C 5 is less than one, then there exists the operator (I − K) − 1 : C 5 ⟶ C 5 and it is bounded.us, the solutionvector x(t) is bounded for t ∈ [0, ∞), for every bounded on the semiaxis right-hand side G(t).Now the exponential estimate of the solution x(t) 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 β < cF * , k > 0, b > 0 for sufficiently small delay τ(t).Remark 5. Let us compare the difference of solution x(t) � col x 1 (t), . . ., x 5 (t)   of system (10) and solution y(t) � col y 1 (t), . . ., y 5 (t)   of system (31) with the same initial conditions x(0) � y(0).It is clear that Denote z(t) � y(t) − x(t).z(t) satisfies the following equation: e solution z(t) can be represented in the following form: It follows from the exponential estimates of C(t, s) and y(t) that z(t) tends to zero exponentially, when t ⟶ ∞.
In order to estimate z(t), we can use the results of the paper [25], where the elements of the Cauchy matrix C(t, s) were constructed.

Construction of the Cauchy Matrix and Stability of Model with Delay
In order to use eorem 2, we to obtain of max 1≤j≤5 sup t≥0  t 0 |C 2j (t, s)|ds.In this section, we present these estimates.