Assume that the model consists of a system of ordinary differential equations
where
zt E R n is a vector of state variables, a E
R '
is a vector of coefficients.Furthermore, assume a linear parameter structure such that
where F ( z ) is an appropriately smooth n x .t function.
Let z t ( a ) denote the solution of equation ( 3 . 1 ) for t E [0
, TI.
It was assumed previously that the experimental data have the following form:However the real t r a j e c t o r i e s of the state variables presumably have a stochastic character and c a n n o t be described within the framework of a deterministic model.
The stochastic character of the trajectories depends not only on errors of measurements but also various internal and external factors which influence process dynamics.
4 . 1 . M a z i m u m Likelihood Estimation
The stochastic character of trajectories can be described by the introduction of a small random perturbation for the model parameters. For each trajectory zi E
Xm
assume that there exists a piecewise continuous function a! such that % ( a { ) = fi,
t E 8. A set of these functions can be considered a s a set of realizations of some stochastic processMoreover
E a t =
I
a t ( w ) d P ( w ) = 6 , V t E [ 0 ,T ] .
n
In this case a set X , can be considered as a contraction on 8 of the set of realizations of the stochastic process
which satisfies the stochastic equation
Then the solution t o the problem of model coefficient estimation reduces t o maximizing the likelihood function
max
@(2,,
a ).
a (4.4)
This problem is discussed in detail in [38], [64].
4.2. Adjoint Estimation of Model Parameters
In the previous section the problem of stochastic estimation of model parameters is discussed. It is a difficult problem, because the likelihood function depends on parameters of the model in the implicit form.
Fortunately an effective numerical algorithm can be constructed which uses the ad- joint equations. As an example, consider the simple deterministic task [50]-[52].
Let the model be represented by (4.1). For the sake of simplicity, we assume that the initial values z(0
,
a) = q E R 3,
R 3 = {z E R nI
z 2 0) are known. Denote z0 G z ( t,
ao) the solution of equation (4.1) satisfying the initial condition z(0,
ao) = q.This solution is said t o be a non-perturbed or a reference one.
Assume that
(a) statistical errors in the measurements are eliminated by appropriate preprocessing of the data,
(b) within the given accuracy f(0) w z(0
,
a0+
e 6 a ),
0 E 8 where z ( t,
a.+
c 6 a ) isthe true or perturbed solution of equation (4.1) satisfying the initial condition z(0
,
a.+
c 6 a ) , a. is a known vector, e>
0 is a small parameter.The problem of evaluating the coefficients of the model using the available data reduces t o that of determining the variation of the coefficients of the model 6 a I
.,
j = 1,..
.,t, which have been chosen, for example, from the condition112 (0) - z ( 0 , a.
+
c 6 a )(I2 -
minAlternatively, a sequence
has t o be defined such that
1)
2 ( 8 ) - ~ ( 8,
a"+
6 a U)I2
-+ min.
as u -+ oo where a. E R L is a given vector.
Let us write the perturbed solutions of equation (4.1) z ( t
,
a.+
r 6 a ) as a series in powers of a small parameter r such that 0<
r<
ro (ro > 0 is a fixed number)Substitute (4.5) into (4.1); expand the right-hand side of the above equality in powers of the small parameter r > 0 up t o terms of the order of O ( r N ) , and equate the terms with the same powers of r > 0 t o obtain a recursive system
Neglecting terms of the order of r2 or higher, for 6zt x zt(ao
+
6 a ) - zt ( a o ) , we have - d 62, = A (t)6zt+
B ( t ) 6 a,
dt (4.7)
6 z O = O , t E [ 0 , TI
,
where
For system (4.7) write the adjoint system
- dt d Y t = - A T ( t ) y t
+
~ ( t ),
Y ( T ) = 0 , t E [ O , T I ,
where p ( t ) is an appropriate function which will be defined below.
Taking the scalar product of (4.7) by yt and (4.8) by 6zt, integrating from 0 t o T, adding together and using the relation
< A ( t ) 62
,
yt > - < A T(t) yt,
6zt > = 0,
we obtain
=
/
< B ( t ) 6 a , yt > dt+ /
< p ( t ) , 6 z t > d t ,0 0
If we choose the function p ( t ) as follows:
where 6(t - 8) is the Dirac delta function, 1 5 k 5 n
,
then (4.10) can be rewritten in the formwhere
and k(t) satisfies equation (4.8) for k = 1
,
. . .,
n.Let 6 a. be the solution of system (4.12) for t E 8. a* = a
+
6a is an unknown vec- tor. As a result of our computation we have 6 ao. Therefore al = a.+
6 a. is a first a p proximation for the unknown vector a*. Now we can write an iterative process for es- timating a*.Actually this is the Newton process with the convergence rate [38]
I
a" - a*I
5 c-l(cI
a. - a* ( )2s,
c = const .4.3. Simple Ezample
Let us consider the following zero system [52]
dz
- dt =
f ( ~ , ,
t E [ O , tk],
~ ( 0 , a) = z O , where
This system of equations in a simplified form describes the change in the number z2 of T-lymphocytes (z2) which occur during the immune response t o non-reproductive an- tigens, e.g. sheep red blood cells in CBA mice. The d a t a on the number of T-cells (helpers) were obtained by R.N. Stepanenko.
Here, 22 (Of) is an average number of T-cells at the instant of time t = Of
,
t = 1,...
5 The non-perturbed solutions of equation (4.15) has the form0 5 2 -is1 t
z2(t
,
5) = z2 exp {T (1 - e ) -63 t}Let us write the equation for the linear part of the variation,
and the corresponding adjoint equation
where
yi(tk) = 0
,
i = 1,2, t E [0,
tk],
8 E (0,
tk).
The solution of problem (4.18) has the form
~ l ( t ) =
The estimate of the variation 6a of the coefficients of the model can be found from the condition
11
J - ~ 6 ~- m i n , ~ 1 1 ~ (4.20) y2(t) =where
&2 0
- z2 exp El
exp
0 , f o r t L 8 , (4.19)
6 2 - - -
- 1 - e 0 1
'1
-z3
81 (1 - e-Ol (B-t))61
- 4
-[
e -%t - e -El Ba 1
,
for t < 40 , f o r t 2 8
and D is a 5 x 3 matrix whose entries are
4.4
Stochastic CaseIn this case a ( t ) = a.
+
ec(t),{ti,
~ E [ O , T ] ) is a stochastic process with E c t = 0 and6 > 0 is a small parameter. A vector of deviations 6zi(a) has random character so t h a t T
6zjk ( a ) = - <6a, $ i ( a ) >
+ $
<BY;, dw,> (4.21)0
which has approximatly a Gaussian probability density function. If the perturbations are independent, t h e mathematical expectation and dispersion have forms of
where
r
is a vector of intensity of perturbation. Estimation of the coefficients of this model can be obtained from t h e likelihood functionIn [38] it is proven t h a t t h e iterative process (6au, r Y ) = arg min 4 ( a , 60, I?)
a
,r
is a quasi-Newton process with first-order convergence rate. T h e estimators, computed by this method, converge t o t h e true values a
*, r *,
with probability one.References
[I] Hoffmann, G. W. and Hraba, T . (Eds.), "Immunology and Epidemiology", Lecture Notes in Biomathematics, Spri ~ger-Verlag, New York, 1986.
[2] Sigmund, K. (Ed.), "Evolution and Control in Biological Systems", J. Acta Applic. petent Cell and AntigennFolia Microbiol. 1 5 , 294-302,1970.
[14] Jilek, M. and Sterzl, J., "Modeling of the Immune Processes", in Morphological a n d F u n c t i o n a l A s p e c t s of I m m u n i t y , Adv. Exp. Med. Biol. 12, K. Lindahl- Kaessling, G. Alm, and M.G. Hanna, Jr., eds., Plenum Press, New York, 1971.
[15] Sercarz, E. and Coons, A.H., "The Exhaustion of Specific Antibody Producing Capa- city During a Secondary Response", in M e c h a n i s m s of Immunological Toler- Optimiz. Conf., Springer-Verlag, New York, 1977.
[I91 Belykh, L.N., Analysis of M a t h e m a t i c a l Models in I m m u n o l o g y a n d Medi-
Mohler, R.R., Barton, C.F., and Hsu, C.S., "T and B Cells in the Immune System",
Adam, G. and Weiler, E., "Lymphocyte Population Dynamics During Ontogenetic Generation of Diversity", in The G e n e r a t i o n of A n t i b o d y D i v e r s i t y , A.J. Cun- ningham, ed., Academic Press, London, 1976.
Hiernaux, J., 'Some Remarks on the Stability of the Idiotypic Network", I m m u n o - chem. 14,733-739, 1977.
Jerne, N.K.: Clonal Selection in a Lymphocyte Network, in Cellular Selection and Regulation in the Immune Response, Edelman, ed., Raven, New York, 1974.
Bell, G.I., "Model for the Binding of Multivalent Antigen to Cells", N a t u r e 248, 1980 (Russian). English version, Optimization Software, Publication Division, distri- buted by Springer-Verlag, New York, 1984.
Mohler, R.R., Nonlinear D y n a m i c s a n d C o n t r o l , Prentice-Hall, Englewood Walsh Functions", IEEE Trans. Auto. Cont. AC-23, 704-713, 1976.
Mohler, R.R., Barton, G.F. and Karanan, V.R., "BLS Identification by Orthosonal
Mathematical Modeling in Immunology and Medicine
-
North-Holland Publishing Co. Amsterdam, 1983, 396 p.Yashin, A.I., Merton, K.G., Stallard, E., E v a l u a t i n g the Effects of Observed a n d U n o b s e r v e d Diffusion Processes in S u r v i v a l Analysis of Longitudinal D a t a , Mathematical Modelling, Vol. 7, 1986, pp. 1353-1363.
(631 Yashin, A.I., Dynamics in Survival Analysis, IIASA, 1984, 30p.
[64] Zuev, S.M., Statistical Estimation of the Coefficients of Ordinary Differential Equations using Observational Data, Sov. J . Numc-. and Math.
Modelling, Vol. 1, 1986, pp. 235-244.
[65] Romanycha, A.A., Mathematical Model of V i a l Hepatitis B. Data Analysis Constructing a Compartmental Model. Preprint
,
Dept.
Numer. Math. Acad.Sci.,USSR, Moscow, 1985,68p.
[66] Romanycha, A.A., Bocharov, G.A., Numerical Identification of Coefficients of the Mathematical Model of Anti-viral Immune Response. Preprint, Dept.
Numer. Math. and Sc., Moscow, 1987, 74p.