A discrete method for the identification of parameters of a deterministic epidemic model

A Discrete Method for the Identification of Parameters sf a Deterministic Epidemic Model*

GIOVANNI DI LENA Istituto di Anuiisi Mutemuku,

Puluzzo Ateneo, 70100 Ruri, Itu!y

AND

GABRIELLA SERlO Jstituto di Igiene - Policliriico Uniuei-sity of h-i, und Istituto di Anulisi Mutemuticu, Puluzzo A terreo, 70100 Buri, It@ Received 5 September I98O; revised 29 Jumtcy

1982

ABSTRACT We study a problem of identification of the parameters for a deterministic epidemic model of the Kermack-McKendlick type. Particular emphasis is put on the analysis of the conditions of numerical stability of the method of integratlljn used to calculate the solutions of the system of differential equations which describe the model. The nu~xxical me*.hodcan be regarded as a discrete model which reproduces the basic qualitative prop& ‘4 of the continuous model, which are positivity of the solutions, points of equilibrium, and the “ threshold theorem.” This allows us to identify the parameters with good reliability. by means of an iterative procedure to minimize the functional which is the measure nf discrepancy between the data observed and the data obtained from the discrete model. The initial estimate of the parameters is obtained by a direct method applied to the discretized system of equations.

1. INTRODUCTION In recent years a large amount of literature has developed concerning deterministic and stochastic mathematical models for the study of infectious diseases [ 1, 3, 8, 171. It is well known that an epidemic model can be described by a system of differential ecluations, dependent on parameters which are usually unknown, and related to the way in which the disease develops. Such parameters play a basic role in that, from their analysis, it is *Work performed under the auspices of the GNIM in the context of the “Progetto Finalizzato CNR di Medicina Preventiva.” MA THEMA TICA L BIOSCIENCES

60: 161- 175 ( 1982)

QElsevier Science Publishing Co., Inc., 1982 52 Vanderbilt Ave., New York, NY 10017

161 0025-5564/82/06 161+ 15$02.75

162

GIOVANNI

DI LENA AND GABRIELLA

SERIO

possible to suggest to the public health authorities the best sanitary policy with the aim of reducing and also controlling the spread of the disease. In this paper we propose to study, with optimization techniques, a problem of identification of parameters for a deterministic model, which belongs to the class which is known in the literature as the class of SIR models. The SIR models are those which, due to contagion, consider a decay from the class of susceptibles (S) to the infective class (I) and finally to the removed class (R), who have no possibility of returning to the state of susceptibility. It is known that these models are suitable for the description of diseases of viral origin, such as measles, chicken-pox, mumps, German measles etc. The models can be formulated in terms of a system of nonlinear ordinary differential equations, which describes the evolution of the vector x(t) = (S(t), I(t), R(t))‘, where S(t), I(t), R(t) respectively denote the proportion of susceptibles, infectives, and removed in the population at the time t. Such a system depends on the unknown parameters 6 = (A, a, y)‘, where A, (Y,y are specified in Section 2, in which an epidemic model is described that is a generalization of the well-known Kermack-MsKendrick model [4]. Our aim is to determine, on the basis of information gained from experimental data in an interval of time (0, T), the value of the unknown parameters 8 = (A, (x,y )‘. Such a problem is known as an inverse problem or a problem of identification of parameters. In order to resolve the problem of identification we construct a nonlinear functicmal J( Q, which is a measure of the discrepancy between the observed data and the corresponding values calculated from the model in the time interval (0, T). The problem of identification of parameters, by means of the minimization of the functional J(e), can present numerical instability, according to the algorithm used for the calculation of the functional itself. It is therefore indispensible that the numerical method for the integration of the system of differential equations should allow us to obtain numerically stable solutions. The method analysed by us in Section 3 not only ‘givesstable solutions but also reproduces the qualitative properties of the continuous model: positivity, the equililxium points, and the “threshold theorem,” which guarantees a “good” evaluation of the parameters to be identified. We have focused our attention on an implicit method of order 1 because, as proved in [3], numerical methods of a higher order do not guarantee the positivity of the solution and therefore the physical meaning of the problem. The minimization with respect to the parameters 8 of the functional J(0) is resolved by means of an iterative technique. Since in an iterative procedure the choiceof the initial point plays a fundamental role, especially in the case of the nonlinear functionals, in Section 4.3 we use a direct method to obtain a good initial estimate of the parameters. This is to guarantee that the minimum obtained by the minimizat$on procedure is that which is sought.

DETERMINISTIC

163

EPIDEMIC MODEL

By mea.us of numerical trials, we show that the initial estimate of the parameters is very close to the minimum point. This is because all the information relative to the observed data in the whole time interval (0. T) is used. Finally, the estimate of the parameters in the Kermack-McKendrick model is an obvious consequence of this analysis. 2.

THE MATHEMATICAL MODEL AND EPIDEMIOLOGICAL CONSIDERATIONS

Deterministic models for contagious diseases have been introduced in a systematic way by Kermack and McMendrick [ 12, 131. In the KermackMcKendrick (KMK) model, the population is divided into three disjoint classes of individuals: (S), the susceptible class, i.e. the class of those individuals who are capable of contracting the diseaslt: and becoming infectious; ( I), the infectious class, i.e. the class of those individuals who are capable of transmitting the disease to susceptibles; ( R ), the removed class, i.e. the class of those individuals who have had the disease and are dead, recovered! plermanently immune, or isolated. In addition, the total population, whicir is uniform and homogenously mixed, has a constant size M, which is sufficiently large. The generation of new infectious individuals occurs by contact between a susceptible and an infectious individual. On the basis of the previous classification and of epidemiological hypotheses, the epidemic process may be described by the following set of differential equations:

dS(t) = -g( -dt

I)S,

dW

yI,

-=g(r)s-

dt

dW -= dt

y[ ,

with initial conditions for t = 0:

s(0) = s() > 0, I(0) = I, > 0, R(0)

(2)

= 0.

The system (1) is a generalization of the original KMK form, as suggested in [4] (see also [7]).

GIOVANNI

164

DI LENA AND GABRIELLA

SERIO

By adding the equations in (1) and taking (2) into consideration we obtain $[S(t)+I(t)+R(t)]=O,

so that, for each moment in time t a0 belonging to the existence interval of the solutions, we have N=S(t)+I(t)+

R(t)=&+

IO.

Here, we have denoted by S(c), I(t), R(t) the number of individuals in the classes (S),(I),(R) respectively, at time t, and y is the removal rate: l/y is the average time of infectiousness of an infectious individual. According to [4], in this model (1) the function g(I) may be a linear increasing function OFthe number of the infectious individuals of the form

as in the KMK model (where A represents the rate of infection), or it may be a continuous nonlinear bounded function. We will refer to the case in which it eventually tends to a “saturation level” in the form:

as suggested in [4]. The introduction of an interaction term of the type g(I) in (4) makes it possible to take into account other i9a.rticula.reffects and has been suggested by the analysis of real epidemics. As usual, we shall limit our analysis to the first two equations in the system (l), which are independent of the third. For the KMK model and for its generalization (1) a “threshold theorem” has been found which contains most of the information concerning the evolution of the epidemic. According to it, we can state that in the phase plane (S, I) al3 the poi-sts on the S-axis are equilibrium poLrts; the points to the right of p* = y/g’(O) are uns,$ablepoints, while the points to the left of p* : are stable points. In the phase plane (S, 1; the system evolves in the triangle S + I G N; the maximum value of I(t) is reached on the “ threshold curve” S = yil/g( I), which becomes S = y/A if g(I) = AI (KMK). Moreover, if the initial point (S,, I, j’ is to the left of this “ threshold curve,” then I(t) is a decreasing function of t; if (So, IO)’is to the right of this curve, then I(t) is an increasing functbn of t up to a maximum value located on the curve itself, and then it decreases. Furthermore, the solution of the system of the differential equations {l), (2) always turns out nonnegative, and S(t) is a decreasing function like S(t)+ l(t).

DETERMINISTIC

165

EPIDEMIC MODEL

It is clear from this that we are able to control an epidemic by imposing suitable conditions on the parameters A, y,p*, and on the size of the population N, as well as on the initial number of the susceptibles So. On the basis of these considerations we conclude that it is very important to know the values of the parameters for a real system, on the basis of observed data, since through them we are capable of defining whether the System is in an endemic state or in an epidemic state, or if the disease is tending to extinction. Furthermore, if we know the values of the parameters, then via the system (1) we may describe the future evolution of the process being studied. 3.

QUALITATIVE ANALYSIS FOR THE DISCRETE MODEL

As observed in the introduction, the method for the numerical integration of the system of differential equations must reflect all the qualitative properties of the continuous model. Furthermore, it is impsrtant tc, choose h (the step of integration) in such a way that the numerical agproximation has the same stability character as the exact solution, independently ?f the vahes of the parameters [ 161. Then we consider the first two equations in ( 1)

4X= dt

dI -=-dt

XIS

-iqr XIS 1+pr

(5)

yr9

which are independent of the third one, subject to the kitial conditions

s(0) =

s, > 0,

(6)

I(0) = 10> 0,

where X > 0, y > 0, /3 = 1/a 2 0. If we apply the implicit Euler method, the system (5) becomes

s:+, =sn-h

AS,,+Jn+, I+/?&,+,

’

(7.1)

where S,? and In are the approximations of S( t,,), I( t,,) with t,, = nh. If we know the values of S,, and I,,, t;ilen by Equations (7) we can obtain approximations to Sn+ , and I,,+, . Furthermore, one may choose h in such a

GIOVANNI

166

DI LENA AND GABRIELLA

SERIO

way that for all nonnegative values of the parameters it gives the results

(8) The inequalities (8) give the stability conditions for the numerical method. Now, in order to prove, in our case, the stability character of Euler’s implicit method, we shall show that the following theorem holds. Eqwtions (7) satisfy the conditions(8) for ali values of the parameters. ProoJ

We observe that if So a 0, I, 3 0, Sn 2 0, and I,, 3 0, then the inequalities (8) are true for n + 1. We determine I,,+, by (7): al,:, , + bI,,+, + c= 0,

(9)

whertu=(l+hy)(/3+Ah), b=(l+hy)-(P+hX)I,,-h?& -c’-4 ’

If we look at the signs of a, 6, c, we can see that b’ -4ac 2 0; then (9) admits real point solution:;. Furthermore, if the solutions are nonzero, they must have opposite signs.’ Then it follows that

Th‘anks to this and to the conditions imposect on the parameters, by using Equation (7.1) we obtain immediately

Then by adding (7.1) and (7.2), we have

xl+,

+r,+,=s,+Q-W,,,,,

and by looking back at ( 10) we can prove the last inequality in (8). This concludes the proof of the theorem. ‘On the basis of this we choose the greatest solution for &+ *.

DETERMINISTIC

EPIDEMIC MODEL

167

Remark 1. Thanks to Theorem 1, the numerical solutions found by Euler’s implicit method are stable for every h. Moreover, since the implicit quantities Sn+ , and I,,+ 1 are calculated most accurately by means of (9) and (1 l), it is possible to combine Euler’s method with Romberg’s extrapolation to obtain a better approximation of the solutions of (5) and therefore of the functional J( 6). Furthermore, analogous with the continuous case and also the discrete case, we can determine the equilibrium points of the system, which give us information on the asymptotic behavior of the solutions. In fact the equilibrium points (S,*, In) of the system (7) are obtained by setting

Sn=Sn+,9

(12.1)

In =J,1+1.

(12.2)

If s”=s”+, = 0, it follows from (7) that 1, = 0, since h, ); ) 0. If S,, + 0, then taking into consideration (12.1), it follows directly from (11) thaL f_:= 0. Consequently the equilibrium points will all be of the ( Sn, 0) type. Also, if It, =O, then In+, is given by the maximum root of the equation (9). Such a root is I,,+, = 0 only when

Therefore the equilibrium points depend on the integration step h, arad all tend to those of the continuous model for infinitesimal h. In any case, whatever the value of h, the: stable equilibrium points for a continuous model are still the equilibrium points for a discrete model. The unstable equilibrium points of the continuous model are not easily identified in a real system, and therefore, the fact that they coot be extracted fro-m the discrete model does not mean any loss of information. On the basis of previous observations and taking into consideration what was proved in Theorem 1, we can state the following “threshold theorem”: THEOREM

2

The necessary and sufficient condition for the epidemic to tend to extinction (that is, 1”2 I,,+,) is that

S,,+(l+8r,)+yhl,,=jj.

(14

GIOVANNI DI LENA AND GABRIELLA SERIO

168

Proof. Considering (7), the values of &I and I,1corresponding to I,, a I,,+, can be obtained by calculating (9) in &. If the result is

then I,, 2 I,, + 1. Otherwise 1” < I,, + , . But (15) is equivalent to (14). The inequality (14) represents the “ threshold value’” in the discrete case. Contrary to the continuous case, this value depends on the integration step h, and it tends to that of the continuous case for h -+ 0. The terminal h ?I,, isn’t considered a perturbed terminal, because it gives information on what. happens after a time interval h = At. The function for the infectious individuals, I(t), is increasing for all t 3 0 if S,, > p and decreasing if S,, ~3. After a time interval At>& I(t+At) will tend to the same value ?asI(t). The “ threshold theorem” in the discrete case depends upon this condition. By means of (14) it is possible therefore to determine the length of the time interval after which the number of kfectious individuals decreases. By giving the values of S, and In, from (14) we obtain Ml, - Y( l+wl,)

At=h=

hY4

(16)

*

Then it is possible to choose values of the parameters X,,y, /? such that the epidemic tends to extinction by imposing the condition

which is related to the threshold value in the continuous case. Observe that all previous considerations are true also if we put in the equations a nonconstant h, i.e., h,,. 4. 4.1.

IDENTIFICATION MODEL

OF THE PARAMETERS

IN A DISCRETE

GENERAL FORMULA TION

We shall devote our attention to the problem of identification of the parameters for a discrete model which has the general form x n+l=X,,+~(XnYXn+,;h,‘e),

with initial conditions x0 = (xh.

n=O,l,..., l

l

xt )‘, and

Y,

(17)

where

8 is the vector of the unknown parameters, t, is the time interval which elapses between the points x, and h.=t,t+,X n+ 1,

DETERMINISTIC

EPIDEMIC MODEL

169

From now on let +n = +(xn,xn+ ,; h,,, 8). If we know the initial condition x0 and the values of the parameters 8, by (17) WC3a.n obtain the solution of the model. The problem of identification of the parameters for the model (17) is solved by the construction of a functional, which is a measure of the distance between the observed data W, for each time t, E (0, T) (n = 0, 1, . . , r) and the corresponding value obtained through Equation (17). In particular, we must minimize the following functional:

J(h,)

= i

w

(x,, -,a)’ 0

tl=

the gradient of which, with respect to the initial conditions and the parameters. is given by ilJ

(19)

ax- = 2qo9 0

where q. and qn are the sclutions in the time ir system: % =rlN+, +Qn+, r),=rl1+

.erval(0,

+ B,rln +e,,

WO

axo%+eo~

n=l,...,r

(21 .l)

n =O,

(21.2)

with initial conditions %Z+,=Qy and where 9, =x, -3,

for

n=O, l,...,r.

IFurthermore A,, B,,, C,l are the following matrices:

4.2.

APPLICATION

TO THE DISCRETE

T) for the following

EPIDE_MIC

MODEL

In the case in which x, = (&, In)‘, 0 = (A, Y,

B)‘, and

.

IA-’

(W+ 1)

,(%I+ 1)

(2’92)

\

I.

“SU

“su

-

“rd+r

W+1 i

’

“ru *(“Id+ I) ;I "SU 6

I

I-“y

=“a

I

,(W+1)

-

i

,“l”sU

i -14 = 4 “V

0

“I -

“Id+1

“rd+1 “r “s

II

T?--

‘Oa+‘&=%

‘()=u $4c*-.rZ‘I=U

(1.9’;‘)

‘ “‘3 + ,‘&Ug + I +q& = “&

uralsk ayl 30 suoqqos ayl an ‘LGpue Ok alay,ti I=24

w

z +

‘IkUy

1

‘OLZ = -Oxe

bz)

re

UIJOJxfs Sl?Y(81j @uopaun3 agi mm vmnww9

aw vNm Ia

6LI

INNVAOI~

DETERMINISTIC

EPIDEMIC MQDEL

171

tabulated at times t, which are different from the time where the solutions Sn and I,, of the system (7) are evaluated. For example, the observed data may be known for regular intervals of time such that h = 1, while, for better accuracy of the numerical integration method, the solutions of the systems (7) are evaluated for h = 0.1. On the basis of the previous remarks it is sufficient to define:
e,,=co,I,-in) if we know only the value of [, at times

t,;

if we know only the value of $ at times

t,, ;

if we know both $, and in at times the form

t,!.

Finally we can rewrite Equation (23) in

(e,,)2.

J(Qo> = i 0

=

0

The gradient of J is evaluated using the method described in this section. Remark 3. If we want to solve the problem of identification of the parameters for the continuous model (5), the form of the functional that we must minimize is [6] J’(~,So,Io)=S’([S(t)-~(tj]2+[I(t)-i(t)]2}dt, 0

where s’(t ) and f( t ) are the continuously differentiable solutions of the differential equation system (5) obtained by a regularization technique or by the interpolation of the observed data Si and Ii at times t, E (0, T) (i = 192,-.*,p).

In the case of the discrete model the regularization of the data is not necessary. Furthermore, in the discrete case, once having obtained the values of the parameters via a minimization technique applied to the functional J’, it is possible to evaluate the unique solution of the system (7); this solution can be considered itself the correct regulsrization of the experimental data.

GIOVANNI

172 4.3.

THE INITIAL,

ESTIMATION

DI LENA AND GABRIELLA

SERIO

OF THE PARAMETERS

“Z’heinitkd estimate of the parameters, which is necessary for the iterative procedure of the minimization of (23), is obtained in the following way. Let us suppose that we know a solution of the difference system (7): ($, I;,), which is evaluated at time t,, , where h ,, = t,, l - t, is not constant, while A, y , /3 are the unknown parameters. To solve the problem of initial estimates we obtain from Equations (7.1), (7.2)

s,.,+r,,,==s,,-+~,,-hnYL,*

(27)

By summing the equations (27) over n = 0, 1, . . . , r - I we obtain

c_

y-

(s,+w-(s,+r,) r-l

.

Furthermore, we rewrite the equation (7. I) in the form

(Srr+ 1 -s,)l,+,B+h,,&

,JI1+*+sn+l-S,~=0

(28)

and in the form

(S,,+,-S,)/?+h,JS,,,+s;A=O.

(29)

n+l

By summing both equations (28) and (29) over n = 0, I,. . . , r - I we obtain r-l

2

I1= 0

wn+,Wn+,P+

r-l

2 hnAS,+,~,,+,+S,-S,=O,

(30.1)

h,AS,,+,+

(30.2)

n=O

r -- 1

(S,-S,)p+

x

II= 0

‘jj’ s,,;,-” II = 0

=O.

ntl

The linear system (30) gives a solution for the initial estimate of the parameters X and 8. The identification of the parameters for the KMK model is obtained by applying the method described in Section 4, taking into account that f’cr the KMK model /3 I= 0. Rtymark 4.

DETERMINISTIC

5.

EPIDEMIC MODEL

173

NUMERICAL TRIALS

Now we show some numerical trials.” First of all we normalize (dividing by N, the dimension of the total population) the differential equation system (I) where g(I) has been given in the form (4). Then the experimental data are obtained for fixed values of the parameters via a numerical method of integration (we have used the RungeKutta method) applied to the system (1). The Runge-Kutta method was used with a time step h such that the numerical values all have a global error of 0.5 X 10V3; these values are displayed in Tables 1 and 2. For the data in Table I the initial estimates of the parameters obtained via the method described in Sections 4.3 are f=l.l27,

(where

?=0.30;’

h = I).

We observe that the value of the functional J( 0) corresponding to the true value of 8 =(1,0.3) is J&0.3) = 2X 10w4, wh ‘e if we put h = 0.05 in the minimization procedure we obtain (1.007, 0.299) as the minimum point for the functional J( 6) and correspondingly J( 1.007, 0.299) = 6 X ! C- ’ A comparison between these two values of the functional J( 8) clearly shows that a better estimate of 0 is not possible, with the assigned accuracy olr the data. Observe that if the experimental data are given with ten significant digits and h = 1, the initial estimates are X=1.122,

f = 0.306.

Furthermore, if we use a tabulation for the observed data in which h = 0.1 (i.e. 100 data with ten significant digits), the initial estimates are ii = 1.010,

ii = 0.3005,

which already are in a good agreement with the true values. For the data in Table 2 we have fixed h = I ; the initial estimate of the parameters are then A = 1.22,

_3= 0.35.,

,i?- 2.42.

If we use a tabulation with 100 data (i.e. h = 0.1 and ten significant digits for the solution), we have

x= 1.OOs,

f = 0.300,

,8 =

1.698.

2The numeric al routines have been implemented on the IBM the APL language.

370/l

58

in

GIOVANNI DI LENA AND GABRIELLA SERIO

174

TABLE I Numerical Simulation of the Experimental Data for the KMK modela T

__~

0 1

2 3 4 5 6 7 8 9 10

__

S

0.8

0.622 0.441 b.296 0.189 0.139 0.112 0.079 0.065 0.055 0.049

I

0.2 0.302 0.380 0.406 0.384 0.336 0.280 0.227 0.181 0.142 0.111

“T = time, S = susceptibles, 1= infectives. Values of parameters: X = 1, y = 0.3, /3 = 0. TABLE 2 Numerical Simulation of the Experimental Data for the Generalized Epidemic Modela T

S

I

0

0.8 0.679 0.564 0.464 0.381 0.314 0.262 0.222 0.191 0.167 0.148

0.2 0.253 0.286 0.299 0.292 0.274 0.247 0.218 0.188 0.160 0.134

I

2 3 4 5 6 7 8 9 10

“T = time, S = susceptibles, I = infectives. Values of parameters: A = 1, y = 0.3, p = 1.7.

The value of the functional corresponding to the true value of the parameters is J(1,0.3,1.7) = 2X lo- *. The functional at the minimum point, which is obtained via the minimization technique in which h = 0.05, is J( 1.004,0.300, 1.699) = 1 x 10-f We conclude by observing that if we lclavea large number of experimental data, the initial estimate of the parameter;; is very close to the minimum point and also to the true value of the parameters.. T.hen it is not strictly necessary

DETERMINISTIC

EPIDEMIC

hIODEL

175

to apply the minimization technique. Vice versa, in the case in which we do not have a large number of experimental data, or we have incomplete observed data, we need to apply the minimization of the functional .I( 8) using the initial estimate as thle starting point of this procedure of minimization. Finally we minimize the functional J( 8,x,) via a gradient-Newton method combined with the procedure of Goldstein [2]. REFERENCES I. 2. 3.

4. 5. 6. 7. 8. 9. IO. I I. 12. 13. 14. 15. 16. 17. 18.

N. T. J. Bailey, The Muthemuticuf Theory of Infections Diseases und Its Applicutions, Charles Griffin, London and High Wycombe, 1975. E. K. Blum. Numerical Anu.&sis and Computation: Theory und Pructice. AddisonWesley, London. 1972. C. Balky and M. Crouzeix, Conservution de Iu positicite Iors de lu JiscrPtisation des prohlt!wres d ‘evolution pwuhaliqws, R.A.I.R.O. Analyse Num&ique. 12. No. 3, 1978. pp. 237-245. V. Capasso and G. Serio, A generalization of the Kcrmack-McKendrick deterministic epidemic model, Muth. Biosci. 42:43-6 1 ( 1978). C. L. Chiang, An Introduction trr Stochustic Processes ucd Tlleir Applicutrons. Robert Krieger, Huntington, N.Y., 19813,Chapter 9. pp. ?f, I-268. G. Di Lena et al., Identification of parameters for a deterministic epidemic model of the Kermack-McKendrick type, Quaderni Istituto di Anal&i Ma;c-matica. Bari, 198 I. H. W. Hethcote. H. W. Stech, and P. van den Driessche, Stability analysis ~,CII. models of diseases without immunity, J. Muth. Biol. 1?:185 (1981). F. Hoppensteadt, lbfuthemuticul Theories of Popidutions: Dtwogruphit’s. Genc*tr~:vund Epidemics, SIAM. ‘I975. S. C. S. Jacoby., J. S. Kowalik, and J. T. Pizzo. Iterutioe Methods for Norrlirwr Optimixtion ProWems, Prentice-Hall, E972. ,S,stem Idnrrtifiwtion Methods und Appiicutions. Addison-Wesley, H. H. Kagiwad;, London, 1974. R. E. Kalaba and K. Spingarn. Optimal input system identification for non-linear dynamic system:;, J. Qptinz. Theop?*Appl. 2 I:9 I- 102 ( 1977). N. 0. Kermack and A. G. McKendrick, Contributions to the mathematical theory of epidemics, Part I, !‘roc. &I*. Sot. London Ser. A I 15:700-721 (1927). W. 0. Kermack and A. G. McKendrick, Contributions to the mathematical theoi?: of epidemics, Part II, Proc. Rcy. Sot. Londw Ser. A 138:55-83 (1932). A. P. Sage and 9. I.. Melsa, Swtem ldentificution, Academic, London. 197 I. Y. Sawaragi, T. Soeda, and S. Okato, Modeling Estinrutions und their Applitutions fol Distributed Purumeter .Qstems, Springer, New York. 1978. H. J. Stetter, Artu!\*sis of discretizution nwthods for ordinuy differentiul t qwtiora. Springer, New York, 1973. P. Waltman, Deterministic Threshold Models in the Theo?* of Epidemics. Lecture Notes in Biomathematics. Springer, New York, 1974. G. L. Yang and (2. L. Chiang. A time dependent Simple stochastic epidemic. in Proceedings of tha, Sixth Berkelgbv Symposium on Muthenwticctl Stutistits und Prohubifi~~~ (J. Neyman. Ed.). \Jniv. of California Press, Berkeley, 1971. Vol. 4, pp. 147- 1%.