=Paper=
{{Paper
|id=Vol-1662/mod4
|storemode=property
|title=Adaptive Wegstein method for a coefficient inverse problem for one model of HIV infection
|pdfUrl=https://ceur-ws.org/Vol-1662/mod4.pdf
|volume=Vol-1662
|authors=Irina F. Iumanova,Svyatoslav I. Solodushkin
}}
==Adaptive Wegstein method for a coefficient inverse problem for one model of HIV infection==
https://ceur-ws.org/Vol-1662/mod4.pdf
Adaptive Wegstein method for a coefficient inverse
problem for one model of HIV infection
Irina F. Iumanova1 Svyatoslav I. Solodushkin1,2
irina.iumanova@urfu.ru s.i.solodushkin@urfu.ru
1 – Ural Federal University (Yekaterinburg, Russia)
2 – Krasovskii Institute of Mathematics and Mechanics (Yekaterinburg, Russia)
Abstract
We consider a coefficient inverse problems for one model of HIV in-
fection. The problem is formulated as an minimization problem of a
quadratic residual functional. The last one is turned to the fixed point
problem. We study two approaches to obtaining approximations such
as the adaptive Wegstein method to solutions of the fixed point prob-
lem. One way is the vector approach, this method uses the ratio of the
norms of residuals for finding parameters of the method. The other way
is the componentwise approach that shift tracking of approximations
in each coordinate more subtle. Numerical results are presented.
1 Adaptive Wegstein method for fixed point problem
There are an impressive number of methods for the approximate solutions of nonlinear scalar equations and
systems. The search for new methods of solving such equations and ways to improve the efficiency of the well-
known and practically proven methods remains important. Newton’s method and simple iterations method,
being the simplest to build and structure and most capable of an exhaustive study, often serve as the starting
point for this.
One way to improve the convergence and extend the scope of applicability of linearly convergent iterative
processes is to build processes of the Fejer type [1] which are called Mann iterations (see [2, 3, etc.]). Suppose
that the problem of a fixed point is defined by the equality
x = Φ(x). (1)
Here Φ : Rn → Rn is a given continuous nonlinear mapping. These iterations are based on the formula
x(k+1) = λΦ x(k) + (1 − λ)x(k) (k = 0, 1, 2, . . . ). (2)
Here λ is a scalar parameter or a sequence of parameters λ = λ(k) . Different approaches to the method of
fixing this parameter lead to different specific processes of the family (2). The Wegstein method [4] for solving
one-dimensional equations (1) is one of sufficiently effective members of this family. In elaboration of this method,
the authors of this paper proposed in [5] to fix the coefficients of the linear combination of the points x(k) and
Copyright c by the paper’s authors. Copying permitted for private and academic purposes.
In: A.A. Makhnev, S.F. Pravdin (eds.): Proceedings of the 47th International Youth School-conference “Modern Problems in
Mathematics and its Applications”, Yekaterinburg, Russia, 02-Feb-2016, published at http://ceur-ws.org
261
Φ x(k) to be inversely proportional to the values of taking residuals of these values. This method is defined by
the following set of formulas:
λ(k) 1
x̃(k+1) = x̃ (k)
+ x (k+1)
, x (k+1)
= Φ x̃ (k)
, (3)
1 + λ(k) 1 + λ(k)
(k) x(k+1) − Φ x(k+1)
λ = , k = 0, 1, 2, . . . (4)
x(k+1) − x̃(k)
The resulting iterative method is identical to the Aitken ∆2 method, in which the formula for acceleration is
applied every other step of the method of simple iterations. The Aitken method in the form (3), (4) is convenient
for studying the conditions of its quadratic convergence and gave rise to an adaptive algorithm in which the
fulfillment of these conditions is checked and corrected in the process of its implementation [6].
When solving scalar equations, the parameters can be calculated using the formula
(k) x(k+1) − Φ x(k+1)
λ = , (5)
x(k+1) − x̃(k)
which allows extending this method to the case of systems of nonlinear equations.
The method (3), (5) has more limited conditions for applicability, but the method (3) with the parameters
(k) kx(k+1) − Φ x(k+1) k
λ = (6)
kx(k+1) − x̃(k) k
is formally suitable for solving systems of nonlinear equations in the form (1). A careful study of the one-
dimensional case [6] shows that, to ensure fast convergence of the method (3), when fixing the parameter λ(k) at
each iterative step, it is important to consider the nature of convergence/divergence of simple iterations delivering
the next intermediate value x(k+1) . In particular, it is important to know whether the points x(k+1) and x̃(k)
are on the same side or on the opposite sides with respect to the fixed point x∗ . In an n-dimensional case, the
behavior of the starting and intermediate points during each approximation can besimilarly taken into account
(k) k (k)
for each coordinate separately. The scalar parameter λ becomes a vector λ = λi for this purpose. The
Mann process is implemented by coordinates as follows.
Let the given nonlinear system (1) in an expanded form be
x1 = ϕ1 (x1 , x2 , . . . , xn ),
x2 = ϕ2 (x1 , x2 , . . . , xn ),
(7)
...
xn = ϕn (x1 , x2 , . . . , xn ).
(k) (k+1)
The transition from the value x̃(k) = x̃i to the value x̃(k+1) = x̃i is carried out based on formulas
(k+1) (k) (k)
xi = ϕi x̃1 , x̃2 , . . . , x̃(k)
n , (8)
(k+1) λ(k) (k) 1 (k+1)
x̃i = x̃ + x (k = 0, 1, 2, . . . ), (9)
1 + λ(k) i 1 + λ(k) i
The matrix-vector form of the method (8), (9) is
x̃(k+1) = Λ(k) Λ̄(k) x̃(k) + Λ̄(k) x(k+1) ,
(10)
x(k+1) = Φ x̃k
(k = 0, 1, 2, . . . ),
( )
(k)
n
(k) (k) (k)
o 1 1 1
where Λ = diag λ1 , λ2 , . . . , λn , Λ̄(k) = diag (k)
, (k)
, ..., (k)
.
1 + λ1 1 + λ2 1 + λn
(k)
We use two options for calculating values λi required in (10).
262
(k) (k+1) (k+1)
β (k) , if
sgn x̃ − x = sgn ϕi (x(k+1)
) − x , x(k+1) − Φ(x(k+1) )
(k)
1. λi = i i
i
where β (k) = .
(k) (k+1) (k+1)
−β (k) , if sgn x̃i − xi = sgn xi − ϕi (x(k+1) ) , x(k+1) − x̃(k)
(k+1) (k)
(k) xi − xi
2. λi = ; i = 1, . . . , n.
x̃i (k−1) − x̃i (k)
Theorem 1. Let a vector-valued function Φ be continuously differentiable in some domain M ⊆ Rn , that
contains its fixed point ξ. Let there exist nonnegative constants ak , bk , k = 0, 1, . . . , such that
n
1 (k)
X ∂ϕ(x)
max (k)
6 ak , max λi + 6 bk .
i 1 + λi i
j=1
∂xj
If there exist q, that the following inequality holds
a k bk 6 q < 1 (11)
(k)
and all members of the sequence x̃ and x(k) , defined by (10), do not leave the M, then the sequence x̃(k)
converges to ξ provided that the initial point x̃(0) ∈ M.
(k)
P r o o f. Regardless of the way in which parameters λi , i = 1, . . . , n is chosen, we deduce from equality (10)
that h i
ξ − x̃(k+1) = ξ − Λ(k) Λ̄(k) x̃(k) − Λ̄(k) x(k+1) = Λ̄(k) ξ − x(k+1) + Λ(k) ξ − x̃(k) ,
because
! ! !
(k)
(k) (k) (k) 1 1
(k)
1 λi
Λ̄ + Λ̄ Λ = diag (k)
+ diag (k)
· diag λi = diag (k)
+ (k)
=
1 + λi 1 + λi 1 + λi 1 + λi
!
(k)
1 + λi
= diag (k)
= E,
1 + λi
where E is the identity matrix.
Using the equality ξ = Φ(ξ) and the mean-value theorem for each row vector function ξ − x(k+1) , we obtain
h i
ξ − x̃(k+1) = Λ̄(k) Φ (ξ) − Φ x̃(k) + Λ(k) ξ − x̃(k) = Λ̄(k) Γk + Λ(k) ξ − x̃(k) , (12)
∂ϕ1 ∂ϕ1 ∂ϕ1
∂x1 ∂x2 . . .
∂xn
∂ϕ2 ∂ϕ2 ∂ϕ2
...
(k) (k)
where Γk = Γ θ1 , . . . , θn = ∂x1 ∂x2 ∂xn is the matrix, in which the partial derivatives are
... ... ... ...
∂ϕn ∂ϕn ∂ϕn
...
∂x1 ∂x2 ∂xn
(k) (k)
calculated with x = θi , i = 1, . . . , n. If the following inequality holds ξ − θi 6 ξ − x̃(k) , equality (12)
implies the inequality
ξ − x̃(k+1) 6 Λ̄(k) · Γk + Λ(k) · ξ − x̃(k) 6
6 ak bk ξ − x̃(k) 6 q ξ − x̃(k) .
Successive iterations of this inequality yield
ξ − x̃(k+1) 6 q 2 ξ − x̃(k−1) 6 · · · 6 q k ξ − x̃(0) .
We have (11). It is now obvious that the above estimate for the error of the approximation x̃(k+1) implies the
convergence of the sequence x̃(k) to ξ for each x̃(0) ∈ M.
263
Theorem 2. Let a vector-valued function Φ be continuously differentiable in some domain M ⊆ Rn , that
contains its fixed point ξ. Let there exist nonnegative constants a, b c, such that
n
1 (k)
X ∂ϕ(x)
max (k)
6 a, max λi 6 b, max 6 c.
i 1 + λi i i
j=1
∂xj
If the following inequality holds
a(b + c) < 1
(k)
and all members of the sequence x̃ and x(k) , defined by (10), do not leave the M, then the sequence x̃(k)
converges to ξ provided that the initial point x̃(0) ∈ M.
To illustrate good properties of the elaborated method we apply it for solving one inverse coefficient problem
for a model of HIV dynamics. An observer measures the concentration of viral particles in blood as well as
a serum level of immunocompetent cells. These data are used to estimate coefficients of the model. This is
a typical inverse dynamic problem and to deal with it a quadratic residual functional is introduced and the
identification problem is formulated as a minimization problem. As a necessary condition for extremum dictates
the last problem is reduced to a nonlinear system to be solved numerically. To accelerate a convergence of
iterative process we use method (10). The implementation of these ideas in a formal and precise way is given in
a paragraph 3 while a concise description of the HIV model is given in a paragraph 2.
2 Description of the model
Mathematical models provide a means to understand the human immunodeficiency virus (HIV)-infected immune
system as a dynamic process. Models formulated as differential equations for the dynamic interactions of CD4+
lymphocytes and virus populations are useful in identifying essential characteristics of HIV pathogenesis and
chemotherapy [7, 8]. The equations for the model with treatment are as follows:
dT (t) S2 VS (t) λ1
= S1 − − µT T (t) + T (t)V (t) − (kS VS (t) + kr Vr (t)) T (t),
dt B S + V S (t) C + VS (t)
dTS (t) λ2
= kS VS (t)T (t) − µTi TS (t) − TS (t)V (t),
dt C i + VS (t)
dT (t) λ2
r
= kr Vr (t)T (t) − µTi Tr (t) − Tr (t)V (t), (13)
dt C i + VS (t)
dVS (t) λ3 GS V S
= (1 − q) TS (t)V (t) − kV T (t)VS (t) + ,
dt Ci + VS (t) B + V (t)
dVr (t) λ3 Tr (t) λ3 TS (t) Vr (t)
= V (t) + q V (t) − kV T (t)Vr (t) + Gr (V (t)) ,
dt Ci + VS (t) Ci + VS (t) B + V (t)
where V (t) = VS (t) + Vr (t), Gr = (GS e(V −V0 )p )/(1 + e(V −V0 )p ), the initial conditions are T (0) = T0 , TS (0) = 0,
Tr (0) = 0, V (0) = V0 , the dependent variables are T is uninfected CD4+ T-cell population, Ts is CD4+ T-cell
population infected by virus sensitive to the treatment, Tr is CD4+ T-cell population infected by virus restrictive
to the treatment, Vs is infectious HIV population sensitive to the treatment, Vr is infectious HIV population
restrictive to the treatment.
In equation 1 the term S1 − S2 VS (t)/(BS + VS (t)) represents the external input of uninfected CD4+ T-cells
from the thymus, bone marrow, or other sources. It is assumed that there is a deterioration of this source as
the viral level increases during the course of HIV infection, µT is the death rate of uninfected CD4+ T-cells
whose average lifespan is 1/µT . The term λ1 T (t)V (t)/(C + VS (t)) represents CD4+ T-cell proliferation in the
plasma due to an immune response that incorporates both direct and indirect effects of antigen stimulation,
where C is a saturation constant. The form assumed here idealizes the growth mechanisms of CD4+ T-cells,
since subpopulations of antigen specific CD4+ T-cells are not modeled. In equation 1 kS is the infection rate of
CD4+ T cells by virus (it is assumed that the rate of infection is governed by the mass action term kS VS (t)T (t)).
In the absence of virus the CD4+ T-cell population converges to a steady state of S1 /µT .
In equation 2 there is a gain term kS VS (t)T (t) of CD4+ T-cells infected by drug-sensitive virus, a loss term
µTi TS (t) due to the death of these cells independent of the virus population, and a loss term λ2 TS (t)VS (t)/(Ci +
VS (t)) dependent on the virus population due to bursting or other causes. The dependence of the loss term
264
λ2 TS (t)VS (t)/(Ci +VS (t)) allows for an increased rate of bursting of infected cells as the immune system collapses
and fewer of these cells are removed by CD8+ T-cells. The structure of equation 3 is the same.
In equation 4 the virus population is increased by the term λ3 TS (t)V (t)/(Ci +VS (t)). This term corresponds to
the internal production of virus in the blood. The dependence of this term on TS (t) allows for a decreased rate of
viral production in the plasma when the infected CD4+ T-cell population in the plasma collapses. Since most of
the plasma virus is contributed by the external lymph source, the plasma virus population still increases steeply
at the end stage of the disease. In equation 4 the virus population is decreased by the loss term kV T (t)VS (t),
which represents viral clearance. In equation 4 there is a source of virus from the external lymphoid compartment,
which is represented by the term GS VS /(B + V (t)) (B is a saturation constant). This term accounts for most of
the virus present in the blood. the structure of equation 5 is the same.
A complete list of parameters and their estimated values for model (13) is given in [8]. Figure 1 shows a
schematic diagram of the entire immune response process.
Figure 1: Schematic diagram of the working immune system
Figure 2 represents the progression to AIDS, the total number of viruses increased rapidly after 3000 days,
while the population of CD4+ T-cells decline to zero.
3 A coefficient inverse problem and numerical experiments
Let tk , k = 1, . . . , m, t1 < t2 < · · · < tm are fixed time moments of measurements. We also defined state vector
x = x(t) = (T (t), TS (t), Tr (t), VS (t), Vr (t)) and the vector of the parameters a∗ = (a1 , . . . , ap ) to be estimated,
1 6 p 6 16. We denote by x(a) = x (t; a) the solution obtained with the given parameters a, xk (a) = x (tk ; a) .
Let us denote by yk the results of measurements at moment tk , and the accuracy of measurements satisfies
|yi, k − xi, k | 6 ξi , i = 1, . . . , 5, where ξi are given.
Let the initial approximation a0 = a01 , . . . , a0p is known, it is required to estimate the parameters a∗ . We
consider a quadratic residual functional
5 X
m
X 1 2
Φ(a) = (xi, k (a) − yi, k ) . (14)
i=1 k=1
ξi
A minimization problem for a functional (14) is turned to solving of nonlinear equation
∂Φ(a)
= 0, (15)
∂al
265
Figure 2: This is the numerical solution of model (13). Parameter values used to generate this figure can be
found in [8]
∂Φ(a) P5 P m 1 ∂xi, k (a)
it is easy to check that =2 (xi, k (a) − yi, k ) , l = 1, . . . , p.
∂al ξ
i=1 k=1 i ∂al
Method (10) is applied for solving the problem (15). We use a priory information to set the initial approxima-
tion of the vector a0 . Calculation of the function in (15) implies a numerical solving of the differential equations
(13). Because of stiffness we use implicit fourth order Runge-Kutta method each step of which implies solving
of an appropriate nonlinear system, and adaptive method (10) is used for this too.
Figure 3: This is the results of parameters estimation for model (13)
To illustrate good properties of the proposed method (10) and efficiency of the approach described above
we performed several numerical experiments. Let B and BS are unknown parameters to be estimated, so
a = (B, BS ). To simulate the system dynamics and obtain the observations a∗ was defined according to [8] as
a∗ = (2, 13.8). The initial approximation was taken as a0 = (3, 18).
The system dynamics for true parameters value is depicted by solid lines, for initial value a0 by green and blue
266
stars, and as it clearly seen in the picture it is rather far. Numerical algorithm (10) found a = (2.1207, 13.4800),
the corresponding dynamics is depicted by red and black stars and is extremely closed to the true one.
4 Conclusions
The advantage of the proposed method (10) for solving nonlinear systems is that only one function calculation
is required on each iterative step, while the calculation of the derivatives and matrix inversion is not required at
all. In one-dimensional case the proposed method has the quadratic convergence; in general case according to
the results of numerical estimations method (10) converges not slower than Newton’s method.
Acknowledgements
This research is supported by the program 02.A03.21.0006 on 27.08.2013, Russian Science Foundation (RSF)
14-35-00005, RFBR 14-01-00065.
References
[1] V. V. Vasin, I. I. Eremin. Operators and iterative processes of the Fejer type. Theory and Applications.
Yekaterinburg: Ural Branch of RAS, 2005.
[2] W. R.Mann. Mean value methods in iteration. Proc. Amer. Math. Soc., 4: 506–510, 1953.
[3] Ş. M.Şoltuz. The equivalence between Krasnoselskij, Mann, Ishikawa, Noor and multistep iterations. Math-
ematical Communications, 12: 53–61, 2007.
[4] J. H.Wegstein. Accelerating convergence of iterative processes. Communications of the ACM, 1(6): 9–13,
1958.
[5] V. M. Verzhbitskiǐ and I. F. Yumanova. An Analog of the Wegstein Method for Accelerating the Convergence
of Iterative Processes. Intell. Sist. Proizvodstv., 1(15): 18–28, 2010.
[6] V. M. Verzhbitskiǐ and I. F. Yumanova. On the Quadratic Convergence of the Aitken ∆2 Process. Compu-
tational Mathematics and Mathematical Physics, 51(10): 1659–1663, 2011.
[7] D. Kirschner. Using mathematics to understand HIV immune dynamics. Notices of the AMS, 43(2): 191–
202, 1996.
[8] D. E. Kirschner and G. F. Webb. Resistance, Remission, and Qualitive Differences in HIV Chemotherapy.
Emerging Infectious Disiases, 3(3): 273–283, 1997.
267