- Email: [email protected]
A mathematical model for the global spread of influenza
A Mathematical
Model for the Global Spread of Influenza
LEONID A. RVACHEV Gomulq~u Institute of Epidemiology and Microbiology, Acudemy of Medical Sciences, Moscow, U.S..% R.
AND IRA
M. LONGINI, JR.*+ Depurtment of Statistics und Biometty, Emory Universit.v, Uppergute House, Atlanta, Georgiu 30322 Received
7 January
1985; revised IO April I985
ABSTRACT A mathematical model is presented for forecasting the global spread of influenza based on information from the initial city in the transportation network to experience the disease. This model represents the natural extension of nearly 20 years of Soviet work on modeling the geographic spread of influenza in the U.S.S.R. and Bulgaria. However, the work presented in this paper is the first attempt at applying the methods on a global scale. A description of the model formulation is given along with a method for estimating critical parameters. The model is then applied to forecasting the global spread of the “Hong Kong” pandemic of 1968-1969 based on estimated parameters from Hong Kong, the initial city to report the appearance of the then new strain of influenza A. The forecast is shown to reproduce Health
the general
Organization
time-space
spread
of the actual
epidemic
as documented
by World
sources.
INTRODUCTION Influenza is the last of the classic “great plagues” of the past which has yet to be brought under control. Major influenza pandemics (see Glossary for an explanation of epidemiological terms) occur every lo-12 years, and more localized “epidemics” occur almost amurally during the interpandemic periods [l]. The use of killed virus vaccines has had little effect on controlling
*Part of this paper was written while I. Longini and Epidemiology at the University of Michigan. + Please address reprint requests to I. Longini. MATHEMATICAL
BIOSCIENCES
was in the Departments
of Biostatistics
3
75~3-22 (1985)
QElsevier Science Publishing Co., Inc., 1985 52 Vanderbilt Ave., New York, NY 10017
0025-5564/85/$03.30
4
LEONID
A. RVACHEV
AND IRA M. LONGINI,
JR.
influenza spread. Effective control strategies may ultimately depend on greater understanding of the geographic transmission dynamics of influenza virus, which is caused by the flow of infected individuals traveling through the global transportation network. This paper reports on progress made by Soviet and Western researchers on the mathematical modeling of the geographic spread of influenza among the world’s major cities [2]. Although the model was originally constructed for influenza spread, it could be used to model the geographic spread of any agents that are transmitted from person to person. Furthermore, this research represents the first attempt at modeling the global spread of an infectious agent. It is the direct outgrowth of modeling studies of the geographic spread of influenza in the territory of the U.S.S.R. which were carried out as early as the 1960s [3, 41. During the 1970s numerous monographs and articles were published describing applications of such models to the territory of the U.S.S.R. (see [5] for a summary of this work) and to the territory of Bulgaria (e.g., see [6]). Reports in English on this work have been written by Baroyan et al. [7-lo], Bailey [ll] (see Chapter 19), and Fine [12] (see Section 2.4.2). Most recently, a refined model has been applied to the territory of the U.S.S.R. [13]. The objective of the present work is to model the pandemic spread of a single strain of influenza among 52 of the world’s major cities located on all the continents (see the list of cities in Table 1). The original global model is given in Baroyan et al. [lo], but the refinement and application of this work is covered in detail by the first set of algorithms in the previously cited manuscript [2]. This first set of algorithms describes the construction of the forecast on a computer. Rather than give the algorithms here, we instead give the formulation of the model, in detail, along with a description of the methods used for parameter estimation. This is followed by an example of how the model was used to forecast the 1968-1969 “Hong Kong” influenza pandemic based on initial epidemic data from Hong Kong, i.e., the first major city to report the then new variant of influenza A. This is followed by a general discussion of the method as well as of future applications and refinements of the model.
THE MODEL The model is formulated as a system of difference equations in a continuous state space and discrete time domain, where the time unit is one day. Intracity epidemics can occur in each of the n cities in the model. Individuals traveling through the transportation network can spread infection from city to city, creating the intercity pandemic; the cities are located on several continents.
.. -
_,
POP.
_
___
Mean Daily Number
_
of Airline Passengers
-
m
_
;
;
-
_
,
i
TABLE 1 between Modeled Cities during the 1968-1969
Pandemic
-_ .--
Influenza
-
-_,.
6
LEONID
DYNAMIC
A. RVACHEV
AND IRA M. LONGINI,
JR.
VARIABLES
For modeling purposes, the population is partitioned into four explicit, disjoint disease states within each city, i = 1,2,. . , n, which are as follows: x, (t)= number of susceptible individuals on day t, u,( r, t)= number of latent (incubating) individuals on day t who were infected on day t-r, y, ( T, t) = number of infectious (ill) individuals on day t who were infected on day t-r, z, (t) = number of recovered (immune) individuals on day t. The population size, p,, of each city is assumed to be time invariant purposes of modeling influenza. Therefore, we have x,(t)+
g 7=0
u,(7,t)+
5 Y,(~,f)+zr(t)=P,
for every t,
7=0
for the
(1)
where ri is the maximum length of the latent (incubation) period and r2 is the maximum length of the infected period, i.e., the latent plus infectious period, taken from the initial day of infection. The variable z,(t) can be considered redundant and, therefore, a total of n [l + (TV+ l)+ ( r2 + l)] = n (TV + TV + 3) equations are required to model the pandemic process. Insofar as the length of the latent and incubation periods can be considered equal for influenza [l], the latent and incubating states are identical. For epidemiological purposes, two additional dynamic variables are used. They are as follows: w,( t)=number of individuals who became infectious (ill) on day t, i.e., daily morbidity or daily incidence of illness, b,(t)= computed (or modeled) number of new ill individuals reported to the health registry on day t. DISTRIBUTION
OF INFECTION
IN INDIVIDUALS
The infection distribution for individuals is defined for a set of empirically derived probability distributions. For an infected individual, let L, be a random variable for the length of that individual’s latent period, and L, be a random variable for the length of that individual’s infectious period. Also, let .& be a set which is the statistical ensemble of individuals who were all initially infected at the same time. Then the infection states (i.e., latent, infected, removed) of the individuals in S? form a partition with the following discrete probability distributions: f(r)=P(&>r),
7’0,l
g(7)=P(L1<7<&),
r=O,l,...,
T2,
(3)
h(T)=P(Lz<7),
r=O,l,...,
72+1,
(4)
7.. . ,711
(4
MODEL
FOR GLOBAL
7
SPREAD
and f(7) + g( 7) + h(7) = 1, T 2 0; pi and T* are the maximum lengths of the latent and infectious periods, respectively. The distributions f(7), g(T), and h(7) can be thought of as the fractions of individuals in S? who are in the latent, infectious, and removed states, respectively. It is assumed that f(0) = 1 and therefore g(0) = h(0) = 0, since the latent period should be at least one day. Using Equations (2)-(4) and standard probability arguments, we derive the following two transition probabilities:
f(T)-ff(T+l) f(T) ’
y(T)=P(L1T)=
f(T)
‘0,
T=o,l,...,
T1, (5)
h(~+l)--h(7)
S(T)=P(L*~T+~~L~~T
g(T)
’ T=O,I
g(T)‘0 ,..., T2. (6)
Then Y(T) is the probability that a latent individual becomes infectious (ill) on day T + 1, given that the individual was still latent on day T; while S(T) is the probability that an infectious (ill) individual recovers on day T + 1, given that the individual was still infectious (ill) on day 7. INFECTION
PROCESS
Within each city, susceptibles and infectives are assumed to mix homogeneously. The daily infectious contact rate h is defined as the average number of individuals with whom an infectious individual will make sufficient contact (to pass infection) in a day. Therefore, the average number of susceptibles infected by an infectious individual in city i on day 1 is Xx, (t)/p,.The number of infectious individuals on day t is given by
Therefore, the number in city i on day t is
of new latent individuals
u,(0, t) = y TRA NSPOR
TA TION
g ’
due to infectious
individuals
U,(O,t-T)g(T).
7=1
OPERA TOR
Individuals move through the transportation network with daily intensity a,,. We define q, as the daily passenger flow from city i to city j, i.e., the average number of individuals that travel from city i to city j in a 24-hour
8
LEONID
A. RVACHEV
AND IRA M. LONGINI,
JR.
period. The n x n matrix u with elements a,, is assumed to be symmetrical, i.e., u,, = a,,, i,j=l , . . . , n. It is further assumed that the number of individuals in each disease state traveling from a particular city on day t is proportional to the fraction of individuals in each disease state in the city at that time and that the population size of each city p, remains constant over time. Ill individuals are assumed not to travel. Based on the above assumptions, a transportation operator Q is constructed to operate on any dynamic variable A,(t): i-@,(t)]
=A,(t)+
2 /=I
A modification of D has been constructed to prevent premature epidemics in cities due to small (fractional) movements of latent individuals. This modification is discussed in the section on forecasting the 196881969 influenza pandemic. STATE
EQUATIONS
for the global system (i = 1,2,. . , n) are as follows:
The state equations
x,(t+l)=D[x,(t)]-u,(O,t),
(9) r=O,l,...,
u,(7+1,t+l)=[l-Y(~)]52[u,(T,t)], y,(7+1,t+1)=
W,(t+l)=
rl-1,
Y(T>Q[u,(7,t)]+[I--(T)]y,(T,t), [I-6(T)]y,(T,t),T=T1+1,T1+2
2 Y(T)~[U,(T,t)]t
(IO)
r=o,l,...,Tir ,...,
T2-1,
(11)
(12)
T=O
h,(t+l)=Pw,(t+l),
(13)
where p is defined as the fraction of newly ill individuals that are reported to the health registry, i.e., reported illnesses incidence. Note that X and p are the same for all cities. The boundary conditions are given by Equation (7) and y, (0, t) = 0. LOCAL.
INFLUENZA
AND
THE INITIAL
CONDITIONS
The city which experiences the earliest epidemic peak is taken as the index city i,. The epidemic process is first modeled locally, i.e., without the transportation operator, in order to estimate critical parameters from data and to index the starting time of the forecast. The usual data are the observed registered number of ill individuals reported { a,(,( t)} in the initial city. Other
MODEL
FOR GLOBAL
SPREAD
9
data can also be used, such as the total number of reported illnesses up to some point in time, or the reported influenza mortality [14]. The equations for local influenza are
x,& i + 1) = X,“(i> - U,“(O,i)
(14)
1
2 u,o(o,i-7)g(T),
q(O, i) = y 10
(15)
r=l
b,(,(i+l)=P; [f(~)-f(T+l)lu,“(o,t-T).
(16)
r=O
where i is a shifted time domain, starting at i = 0, that is lined up with the actual calendar time domain f. Since the early stages of the epidemic are not observable, the most reliable point to make the alignment of t^and t is that as that time time, t,,, , when {a,,(r)} reaches its peak. If we define i,, when {b,&i)} reaches its peak, then we align times t,,, and i,,, and simultaneously count back in time until we reach that time where i = 0; the corresponding calendar time will signify the initial date of the epidemic in city i,. We initialize the local influenza epidemic by specifying time i = 0 as that time when u,&O, i) = 10m5p,0, since the model only includes cities with populations of at least 100,000. In order to get the full set of initial conditions on u,&O, +), we assume that the first appearance of new latent individuals in city i, follows a geometric progression with parameter p > 1, i.e., x,&t) = x,&O), during the early stages of the epidemic. Then it follows that
u&x i + 1) p= for i not too large. Iterating
(17)
U,“(O,i)
on (17) yields
u,O(o, i - T) = UJO, i) p-T. Using the fact that u,~(O,O) = 10-5p,0 and substituting conditions
u,()(O,- T) =loVp,“p-‘, The initial number
of susceptibles
7=0,1,_..,
(18) into (18) yields the
T2.
(19)
is given by
where 0 < OLI 1 is the initial fraction of city i,‘s population
that is suscepti-
LEONID
10
A. RVACHEV
AND IRA M. LONGINI,
JR.
ble to the agent. (Note: (Yis assumed to be the same for all the cities, and when (Y< 1 some of the population may be immune from a previous epidemic.) In order to initialize the pandemic process, we will define t, as the first day at which at least one latent individual (at any stage of infection) travels from city i, to another directly connected city. Since we are still dealing with the early stages of the epidemic in city i,, we substitute Equation (18) into
which yields
Then it follows from the above definition
of f, that
The forecast is made using the state equations (7), (9)-(13) over the time domain t = t,, t, + 1,. . . , to + T, where T is the forecasting horizon. For the forecast. the initial conditions are set at
(24)
and
X,(b) =ap,. ESTIMATOR
OF FREE
(26)
PARAMETERS
The free parameters h and a are estimated from the ascending epidemic curve in initial city i,. The estimation problem is to find estimates ^x and & which minimize the mean squared deviation between the daily reported morbidity { a,$ t)} and the computed morbidity { b,O(f)]~,~} over the time horizon of the available data. The estimated procedure involves minimizing the mean squared deviation A(k), which is a function of p [given in (17)] in the sense that the least squares estimates ^x and & are found at different values of h.
MODEL
FOR GLOBAL
SPREAD
11
The two-dimensional estimation procedure is reduced to a one-dimensional procedure (to a variation of CL)by addition of the constraint ah = S, where s is a constant. The value of s, at fixed values of CL,is found by substituting (18) into (15) yielding
u(o,t)=* :
u(0,t)p-‘g(7).
(27)
7=1
(Note: The subscript i, has been dropped in this section, since it is clear that we are referring to the initial city.) Simplifying and using the initial condition x(0) = ap in (27) yields .S=
1 72 c P-k(T)
(28)
r=l
Now, in order to estimate X and cx, we let A =X/r and a! = r-i?, where x and 5 are arbitrary numbers (with 15 = s) and r is to be estimated. Then x=x/i and &=?C; also ^x&=s. Generally, we set A=s and Z=l. The estimate F is that value which minimizes the sum of squares
where the summation is taken over all values of t where {a(t)} is defined, and b(t) varies within the class of functions generated by Equations (14)-(16) at various pairs of A =x/r and a = r5 subject to Xa!= s. Thus, we compute b(t) using b(t)=rb(t)
for every t,
where b(t) is found by stepping through Equations and a = 5. The least squares estimator
(30) (14)-(16)
when J =h
(31)
is found by solving the equation
gp-(r)b(r)l’l I
r=i
=o.
(32)
12
LEONID A. RVACHEV AND IRA M. LONGINI. JR.
Since s is a function of p, then so are different “best” estimates of A and a when k[‘), k=1,2 ,..., I, e.g., pikJ E (1.05, 1.06,.. and OLare associated with the fixed value squared deviation
x and &, and we get a class of varying the fixed values of CL,say .,1.50}. Within this class, the X of p which minimizes the mean
(33)
where b(t) is b(r) evaluated at 1 and &, and ? is the total number of days over which {a(t)} is defined. There may be gaps in {a(f)}, since morbidity data may not be registered over weekends and holidays. The algorithm is summarized as follows: 1. 2. 3. 4. 5. 6. 7. 8.
Set up the values ).L[‘], k =1,2 ,..., 1. For the present pL’] compute s from (28) and set Cu= 1, x = s. Substitute x and Z into the state equations (14)-(16) and step through them, generating the sequence {b(t)}. Substitute b(f) into (31) and obtain i. Find 1 =I/?, 8= Gi, and {g(r)} = i{%(f)}. Substitute {b(t)} into (33) and obtain A(p). Return to 2 for next value of pLkl until k = 1. Choose ^x and & associated with the smallest A(k) as the best estimates.
It is useful to plot A(p) against p and to approximate the plot with a smooth curve; the minimum of this curve should be close to the same p, as is the minimum A(p). In those rare cases when this does not hold, one should consult Reference [5]. When other data than {a(t)} are available, the algorithm can be modified accordingly. THL
HEMISPHERIC
SEASONAL
SWING
P&gents causing upper respiratory illness tend to occur during the colder months, when conditions are best for person-to-person spread. Accordingly, the timing of epidemics in cities in the northern and southern hemispheres is affected by the beginning and ending of the colder months. In order to model this :easonal swing in the most parsimonious fashion, the cities modeled are distrbuted into three zones: northern, southern, and tropical, with the variable Sn, coded 1, - 1, and 0, respectively, indicating in which zone city i is locited. For cities in the northern and southern zones, seasonal alterations which affect spread of the agents are assumed to be pronounced; for cities in the tnpical zone they are assumed to be absent. Fo the northern and southern zones, the dates Ti and T2 correspond to the btginning and end of the epidemic season in the northern zone. Con-
MODEL
FOR GLOBAL
13
SPREAD
versely, T, and Tl correspond to the epidemic season in the southern zone. The daily infectious contact rate is reduced by a factor of 0.1 during the nonepidemic period, producing the following seasonal pattern in each zone: Northern
zone (Sn, = 1):
A+
Southern
Tropical
x
for
TV
0.1X
for
TV [T,,T,).
0.0
for
f~ [T,,T,),
x
for
TV [T,,T,).
zone (Sn, = - 1):
zone (Sn, = 0): A+A
For influenza pandemics, the calendar October and 1 April, respectively. ALGORITHMS
[T,,T,),
for all t.
dates for Tl and
T2 are set at 1
FOR COMPUTATION
The forecast is accomplished by stepping through the first set of algorithms given in [2]. We give a brief description of these algorithms here. Algorithm 1 (infection process): Sets up the distributions (S), (6), given the distributions (2)-(4). Algorithm 2 (local influenza): Steps through Equations (14)-(16), (19)-(20) in the initial city i,. Estimates the free parameters LY,X from Algorithm 3 (free parameters): data for the initial city i, using the algorithm given after Equation (331. Computes t, from (23) and sets mitial Algorithm 4 (initial conditions): conditions for algorithm 6 using Equations (24)-(26). Algorithm 5 (transport operator): Applies the operators E and D as given by Equations (37) and (8) respectively. Algorithm 6 (global influenza): Steps through Equations (7) (C)-(13), i=1,2 ,...,n. A seventh coordinating algorithm also is given to direct the timng and movement of flow as the computer passes through the above six algorithms. Given this set of algorithms, it is fairly easy to write down the compllter code in any standard computer language, e.g., FORTRAN or PL/I.
14
LEONID
FORECASTING
OF THE 1968-1969
A. RVACHEV
AND IRA M. LONGINI,
INFLUENZA
JR.
PANDEMIC
In July 1968, a new strain of influenza A, i.e., influenza A(H3N2), caused a major epidemic in Hong Kong. This influenza virus, which probably originated in southern China, subsequently swept the globe, terminating its first wave late in 1969 [15]. It was subsequently given the name Hong Kong influenza. The model described in the previous section is used to forecast the influenza pandemic of 1968-1969 in 52 cities based only on background data and the first half of the epidemic curve in Hong Kong. THE
DATA
The population sizes for each city, p,, in 1968, were taken from demographic reference books and are given in Table 1. Also given in Table 1 are the estimates for the transportation matrix taken from international air transport statistics [16-221. The upper triangle of this symmetric matrix has 440 nonzero entries; e.g., an average of 68 individuals traveled between London and Sydney during each 24hour period in 1968. The distribution of the length of the latent period, f(7), was the same as that used by Elveback et al. [23] for simulating influenza A epidemics. The distribution of the length of the infectious period, g(7), was taken from past Soviet work [5, 241. Then h(7) was computed by subtraction, and Y(T) and &(T) were computed using Equations (5), (6), respectively. The distributions are given in Table 2, and we see that T, = 2 and TV = 8. ESTIMATION
OF
FREE
PARAMETERS
Hong Kong was taken as the initial city, and therefore only source of morbidity data for Hong Kong is that of the weekly morbidity from nine state clinics is given. sequence {a(t)} from the weekly totals does not produce
we set i, = 10. The Chang [25], where Extrapolating the a stable sequence
TABLE 2” Infection
Distribution
f(T) g(7) h(T) Y(T) S(7)
Process Distributions
0
1
2
3
4
5
6
7
8
9
1.00 0.00 0.00 0.30 0.00
0.70 0.30 0.00 0.70 0.10
0.20 0.77 0.03 1.00 0.20
0.00 0.82 0.18 0.00 0.34
0.00 0.54 0.46 0.00 0.44
0.00 0.30 0.70 0.00 0.50
0.00 0.15 0.85 0.00 0.60
0.00 0.06 0.94 0.00 0.83
0.00 0.01 0.99 0.00 1.00
0.00 0.00 1.00 0.00 0.00
“From Belova et al. [2]. ‘Tlrne since initial day of infection.
MODEL
FOR GLOBAL
SPREAD
15
for estimating h and a, although such a sequence is useful for approximating p. Chang [25] reported 250 notifications during the first week of the epidemic and 1150 during the second week. Using this information, we extrapolate fi = (1150/250)“’ = 1.2436. Thi s value of p compares well with the value estimated for the influenza A(H3N2) epidemic during January-February of 1969 in the Soviet city of Kuibyshev. From the Kuibyshev data, F = 1.24 was estimated from a proper sequence {a(t)} (see, [5, p. 2041 and, [24, p. 531). This estimate was subsequently used to produce a very good forecast of the daily influenza morbidity in 40 Soviet cities (see, [5, p. 4771 and [24, p. 561). The closeness of these two estimates is remarkable considering that different methods were used to estimate p for the two cities. We further note that Kuibyshev is in the northern zone while Hong Kong is in the tropical zone, and that the estimates were made during different seasons. The estimate F = 1.24 was thus taken as the best estimate and used for the estimation of h and cy. Substituting p =1.24 into Equation (28) yields s = 0.641 and therefore Cu= 0.641. The population of Hong Kong in 1968 was about 3.9 million people. Thus, from Equation (19) we set St0 (0,O) = urO(O,O)h, = 39, and we also computed the initial conditions ~(0, - r), r = 1,2,. . . ,8, from (19). Then was computed via Algorithm 2. From (13) the sequence {W(i)} = { w(i)h6} it is easy to see that the sequences W(i) and b(i) peak at the same time, I,,,, . This time was found to be i,, = 44. The calendar time, I,, , when (a(t)} reached its peak was extrapolated to be 27 July 1968. By aligning t,, with and counting back, the calendar time corresponding to i = 0 was found Lx to be 13 June 1968. In order to estimate X and a, the computed (from Algorithm 2) number of reported illness in Hong Kong up to the peak of the epidemic, Sb( i,,) = C&b(i), was set equal to the observed number Sa( t,,). The latter number was estimated to be 292,000 individuals from Chang’s data [25]. From Equations (13), (30) we have
b( t) = B4 t>
(34)
w(t)=fi(f),
(35)
and
where w(t) is computed using Algorithm 2 with X =x = 1 and OL= 5 = 0.641. From Equations (34) and (35) we have
,,=b(r)= w w(t)
c,>
~~(bulX)’
(36)
LEONID A. RVACHEV AND IRA M. LONGINI, JR.
16 and we estimate
r by
w cmx) i = jL%i?(f,,) =
(0.3;:~~000) z 0.95,
and ^x=x/i=l.OS, &=?&=0.61. The value of t, was computed from Equation (23) and was found to correspond to 12 July 1968. Accordingly, the initial conditions (24)-(26) were set. The forecasting time horizon was set at T = 440 days, which yields a forecast up to 26 September 1969. However, since data were used from Hong Kong up to 27 July 1968, the actual forecast was over a 425 day period. During the first attempt at modeling the 1968-1969 influenza pandemic it was observed that epidemics in cities began too early because of small fractions of latent individuals, e.g., 0.001 latent individuals, passing through the transportation network and initiating epidemics. The problem did not arise previously when the model was applied within the USSR territory and the duration of the epidemic process usually did not exceed three months. This problem was circumvented by constructing an operator, E, that allows latent individuals at a particular stage to be present in a city i only if the approximate expected sum of these individuals in city i and those entering from the other connected cities is at least one. The form of this operator is as follows:
I
!2[[u,(~,t)] if
E{fhht)l)
=
Q[u,(T,t)]
Cf(T)21, 7=0
0
if
Q[u,(~,t)]
i f(~) 7=0
(37) cl.
The physical interpretation of E is most apparent if we assume that the epidemic is peaking at time t = t* in city i. This implies that U,(O,1*)=~,(O,t*-I)= since
71
...
is small. Then, the total number
=U,(0,t*-T71),
(38)
of latents in city i at time t = t* is
2 U.(T,f*)=,~o~(T)u,(O.l*-T). (39)
T=O
Substituting
(38) into (39) yields ; 7-o
U,(T,t*)
=
U,(0,f*)7~oj(T).
(40)
MODEL
FOR GLOBAL
SPREAD
17
Thus, the factor C:,,f( T) is the expected number of latents in all classes for each new latent when the epidemic is peaking. In order to forecast the 1968-1969 influenza pandemic, Equations (lo)-(12) were modified by applying the operator E to all 3[ u, (7, t)]. THE FORECAST The results of the forecast in terms of predicted daily reported morbidity, { b,(t)}, versus observed daily reported morbidity, { ai( t)}, are given in Figure 1. The observed values were taken from World Health Organization sources [26-301. For some cities the data were quite detailed, but for many cities only the duration and peak period of influenza activity were available. The most conspicuous disadvantage of the data is the absence of empirical information on influenza activity in many of the cities modeled. Probably, in the majority of cases, this means that there was no epidemic in these cities. It is hard to imagine that no information would have reached the World Health Organization had epidemics occurred at least in the territories around the cities in question. Only Havana was predicted not to experience influenza activity over the time period of the simulation. Nevertheless, the occurrence of cities not experiencing significant influenza activity during the pandemic is consistent with the basic premises used in building the model, since the infectious contact rate in cities is decreased during the “nonepidemic” period in the northern and southern zones. In Figure 1, we see that on the space scale of the globe and the time scale of a year the predicted epidemic course generally corresponds to the outline of the real course of the pandemic.
DISCUSSION The original global model was formulated as a system of integrodifferential equations with partial derivatives [lo]. In order to develop the first step of computer algorithms briefly discussed in this paper, the system of difference equations (7)-(13) was formulated as the discrete analog of the original system of equations. A second set of algorithms also has been developed to forecast the global spread of an infectious agent which spreads like the influenza virus but exhibits a higher case fatality rate like the agent of Legionnaire’s disease. In this case, the rapid isolation of exposed and infective individuals is important. Accordingly, the model is expanded to incorporate isolated states, and vital dynamics are added to account for high case fatality rates, i.e., p, becomes a dynamic variable. More advanced methods are used to estimate the free parameters A and (Y where p is not treated as a constant but rather as a continuous function. Further details are available in [2]. With regard to the forecast of the 1968-1969 influenza pandemic, the most time-consuming task was the construction of the transportation matrix.
0
26
I
,...
Kinshasa
I * Empirical dataamjointwit/lIhesurro”“di”g repion.
Cairo
..r
.,
..-
,
..
. . . . . . --++++,+--...
.. ....-- .++f+r+
..,.
... .. .
1
When a symbol falls directly
on next puge)
indicated.
(Continued
on a border
between
months.
this indicates
occurred
during
the four days two days in each month.
b,(t) < 10, - for 10 I b,(r) I 100, + for 100 > b,(r), and i_ when the peak in morbidity
FIG. 1. A schematic plot of the forecasted h(f) and actual a(r) course of the 1968-1969 influenza pandemic. Each of the following symbols represents the daily forecasted morbidity incidence per lo5 over four calendar days: 0 for
1
26
z _-
3
B k s
.
peaks were not reported.
(Redrawn
.
from Belova et al. 121.)
. _
solid line indicates that significant influenza activity was reported from the city or region indicated; parallel lines indicate that a peak in morbidity was reported. If only the peak period was reported, then the parallel lines are surrounded by short solid lines. The arrows indicate the direction, in time. towards the period of peak activity when
If the symbol touches the border from the left or from the right, this indicates that three of the four days were in the month to the right or left, respectively. The prevalence of the actual pandemic is given by the lines above the symbols. A
FIGURE1 (continued)
20
LEONID
A. RVACHEV
AND IRA M. LONGINI,
JR.
For this first application, only the daily intercity airline passenger turnover was used for the 52 cities modeled. Numerous technical difficulties had to be overcome. For example, many passengers going from city i to j may have a short intermediate stop in some other city, k. However, the airline statistics would record the flow from cities i to k and k to j separately, thus inflating the estimated number of individuals who are actually staying in the cities involved. This difficulty, as well as others, had to be resolved by considering that the transportation matrix is not a separate empirical source of data but rather a set of parameters that have to work in concert with the state equations and the estimated population sizes of cities of departures. In order to make the forecast more realistic, it is necessary to expand the number of cities modeled from 52 to at least 152. It is also desirable to refine the transportation matrix to include other sources of transportation such as ground transportation and the movement of freight personnel. There also should be a separate transportation matrix for each day of the week in order to simulate the weekly periodicity of the world timetable of airflights. Another complication involves the modeling of infected individuals. It was assumed that all infectious individuals were ill. However, a certain fraction of infected individuals are asymptomatic for influenza [23]. Since it would be hard to separate this factor (i.e., pathogenicity) from the estimation of /3, asymptomatic individuals were not modeled separately. The model with asymptomatic individuals is given elsewhere [13]. In conclusion, the model was successful in forecasting the broad patterns of the geographic spread of Hong Kong influenza. This lends credence to the theoretical correctness of the concepts used to build the model. The length of the forecasting interval, which proceeded from 27 July 1968 to 26 September 1969, may indeed be a milestone since hardly anybody has ever predicted, even roughly, any daily course of events in the world society for such a period-about 425 days. GLOSSARY Agent.
OF EPIDEMIOLOGICAL
TERMS
An entity which causes disease (e.g., an influenza virus). Unusually frequent occurrence of disease compared to background level (e.g., outbreak of disease in a single city). Incubation period. The period of time beginning when an individual first harbors an agent and ending when that individual begins to manifest symptoms of disease. Infected. Harboring an agent of disease. Infectious period. The period of time when an infected individual can spread the agent to other individuals. This period begins at the end of the latent period and usually ends with recovery. Lutent period. The period of time beginning when an individual first harbors an agent and ending when that individual becomes infectious. Epidemic.
MODEL
FOR GLOBAL
21
SPREAD
Pandemic. An epidemic which includes significant outbreaks in many regions of the world. Puthogenicity. A measure of the ability of an agent to cause illness once the agent has infected an individual. L. A. Rvachev would like to express his gratitude to V. E. Kvitka, Deputy Director and M. L. Petryanova, Senior Researcher, of the Research Institute of the Ministty of Civil Aviation, U.S.S.R.; E. N. Chernvochkin, Department Manager of Automated Control Systems, of the Sklifosovsky Institute; V. A. Becketova, G. P. Busuyok, Y V. Gudkova, Research Workers at the Gamaleya Institute; and V. Z. Sotnikov at the Institute of Applied Mathematics of the Academy of Sciences of the U.S.S.R. I. Longini would like to thank P. E. M. Fine for introducing him to the work described here. REFERENCES 1 2
E. D. Kilbourne, The lnfluenzu L. A. Belova. J. W. Donovan, Rvacbev, modelling
V. A. Shashkov, (Part
I-influenza),
and
Viruses und Influenzu, Academic, P. E. M. Fine, D. W. Fraser, V. I. Vasilyeva, Unpublished
Experiment
manuscript
New York, 1975. M. B. Gregg, L. A. on pandemic
(in Russian
and
process English),
4
1983. 0. V. Baroyan and L. A. Rvachev, Deterministic epidemic models for a territory with a transport network (in Russian), Kiherneriku, No. 3, pp. 67-74 (1967). L. A. Rvachev, Computer modelling experiment on large-scale epidemic (in Russian),
5
Dokl. Akud. Ncluk SSSR 180(2):294-296 (1968). 0. V. Baroyan, L. A. Rvachev, and Y. G. Yvannikov,
3
6
7
Modeling und Forecusting
of
Influenzu Epidemics for Terriroty of the U.S.S.R. (in Russian), Gamaleya Inst. of Epidemiology and Microbiology, Moscow, 1977. 0. V. Baroyan, L. A. Rvachev, K. Frank, et al., Model of influenza epidemics spread through the territory of Bulgaria and its significance (in Bulgarian), J. Eprdemiol. Microhiol. Jnfecr. Dis. (Soficr) 15(3):168-173 (1978). 0. V. Baroyan, V. V. Basilevsky, V. V. Ermakov, et al., Computer modeling of influenza epidemics for large-scale systems of cities and territories (in Russian and English), Working Paper for W.H.O. Symposium on Quantitative Epidemiology, Moscow, 23-27 Nov. 1970.
8
0. V. Baroyan, L. A. Genchikov. L. A. Rvachev, and V. A. Shashkov, large-scale influenza epidemic modelling by means of a computer, Epidemiol. Assoc. 18:22-31 (1969).
9
0. V. Baroyan, L. A. Rvachev. V. V. Basilevsky, et al., Computer modelling of influenza epidemics for the whole country (U.S.S.R.), Ado. Appl. Pro/x&. 3~224-226 (1971). 0. V. Baroyan. G. A. Mironov, and L. A. Rvachev, An algorithm modeling global epidemics of mutant origin, Progrumming and Comput. Software 6(5):272-277 (1981) English transl. from Programmirovunie, No. 5, pp. 73-79 (1980) (in Russian). N. T. J. Bailey, The Muthemutical Theoy of Infectious Diseases and its Applrcutrons, Charles Griffin, London, 1975.
10
11
An attempt at Bull. Inrernur.
22
LEONID
A. RVACHEV
AND IRA M. LONGINI,
JR.
12
P. E. M. Fine, Background paper, in Influenza Models - Prospects for Development and Use (P. Selby, Ed.), MTP Press, Boston, 1982, pp. 18-85.
13
L. A. Belova et al., New influenza epidemic model for the U.S.S.R. territory Russian), presented at Seventh International Symposium on Influenza, Lenningrad, Dec. 1982.
14
C. C. Spicer, 35:23-28
15 16
The mathematical
modelling
of influenza
epidemics,
hit.
(in 10
Med.
Bull.
(1979).
W. C. Cockbum, P. J. Delan, and W. Ferreira, Origin and progress of the 196881969 Hong Kong influenza epidemic, Bull. M/id. Hlth. Org. 41:345-348 (1969). Air Trunsport Stntistics,
Australia,
international
1980. Dept.
Air Transport,
of Trans.,
Canberra,
June 1981.
17
J. W. Donovan, ing Paper.
Australian
18 19 20
ARC World Airyvs Guide: (a) Part One, (b) Part Two, ABC Travel Guides, Ltd., Dunstable. Beds, U.K., Aug. 1981. Official Airline Guide, North American ed., OAG Inc., Oak Brook, Ill., Oct. 1981. Officrul Airline Guide, Worldwide ed. OAG Inc., Oak Brook, III.. Dec. 1981.
21
Data
22
G. P. Busuyok. V. A. Becketova, and V. P. Charlamov, The methods of epidemiological modelling of the U.S.S.R. passenger traffic (in Russian), presented at Seventh International Symposium on Influenza, Leningrad, Dec. 1982.
23
L. R. Elveback, J. P. Fox, E. Ackerman, et al.. An influenza immunization studies, Amer. J. Epid. 103:152-165 (1976).
from the International
V. Baroyan
International
Airport
Civil Aviation
and L. A. Rvachev.
24
0.
25
Moscow, 1977, pp. l-64. W. K. Chang, National influenza
Traffic,
Organization
Mathematics
Dynamics
1968-1980.
Work-
for the years 1968-1982.
simulation
and epidemiology
model
for
(in Russian),
in
Znanie,
Org. 41:349-351 26
International
experience
in Hong
Kong,
1968,
Bull.
Wld.
Hlth.
(1969).
Conference
on Hong Kong Influenza,
Bull. W/d. Hlth.
Org. 41:345-414
(1969). 27
Week!v Epidemiological
Record. No. 34-52. Wld. Hlth.
28 29
Week!v Epidemrologicul
Record, No. l-40,
Org.. Geneva, 1968. Wld. Hlth. Org.. Geneva, 1969.
W/d.
Hlth.
Stutist.
Rep. (Geneva)
22. No. 7 (1969).
30
Wld.
Hlth.
Statist.
Rep. (Geneva)
23, No. 7 (1970).
Model for the Global Spread of Influenza
LEONID A. RVACHEV Gomulq~u Institute of Epidemiology and Microbiology, Acudemy of Medical Sciences, Moscow, U.S..% R.
AND IRA
M. LONGINI, JR.*+ Depurtment of Statistics und Biometty, Emory Universit.v, Uppergute House, Atlanta, Georgiu 30322 Received
7 January
1985; revised IO April I985
ABSTRACT A mathematical model is presented for forecasting the global spread of influenza based on information from the initial city in the transportation network to experience the disease. This model represents the natural extension of nearly 20 years of Soviet work on modeling the geographic spread of influenza in the U.S.S.R. and Bulgaria. However, the work presented in this paper is the first attempt at applying the methods on a global scale. A description of the model formulation is given along with a method for estimating critical parameters. The model is then applied to forecasting the global spread of the “Hong Kong” pandemic of 1968-1969 based on estimated parameters from Hong Kong, the initial city to report the appearance of the then new strain of influenza A. The forecast is shown to reproduce Health
the general
Organization
time-space
spread
of the actual
epidemic
as documented
by World
sources.
INTRODUCTION Influenza is the last of the classic “great plagues” of the past which has yet to be brought under control. Major influenza pandemics (see Glossary for an explanation of epidemiological terms) occur every lo-12 years, and more localized “epidemics” occur almost amurally during the interpandemic periods [l]. The use of killed virus vaccines has had little effect on controlling
*Part of this paper was written while I. Longini and Epidemiology at the University of Michigan. + Please address reprint requests to I. Longini. MATHEMATICAL
BIOSCIENCES
was in the Departments
of Biostatistics
3
75~3-22 (1985)
QElsevier Science Publishing Co., Inc., 1985 52 Vanderbilt Ave., New York, NY 10017
0025-5564/85/$03.30
4
LEONID
A. RVACHEV
AND IRA M. LONGINI,
JR.
influenza spread. Effective control strategies may ultimately depend on greater understanding of the geographic transmission dynamics of influenza virus, which is caused by the flow of infected individuals traveling through the global transportation network. This paper reports on progress made by Soviet and Western researchers on the mathematical modeling of the geographic spread of influenza among the world’s major cities [2]. Although the model was originally constructed for influenza spread, it could be used to model the geographic spread of any agents that are transmitted from person to person. Furthermore, this research represents the first attempt at modeling the global spread of an infectious agent. It is the direct outgrowth of modeling studies of the geographic spread of influenza in the territory of the U.S.S.R. which were carried out as early as the 1960s [3, 41. During the 1970s numerous monographs and articles were published describing applications of such models to the territory of the U.S.S.R. (see [5] for a summary of this work) and to the territory of Bulgaria (e.g., see [6]). Reports in English on this work have been written by Baroyan et al. [7-lo], Bailey [ll] (see Chapter 19), and Fine [12] (see Section 2.4.2). Most recently, a refined model has been applied to the territory of the U.S.S.R. [13]. The objective of the present work is to model the pandemic spread of a single strain of influenza among 52 of the world’s major cities located on all the continents (see the list of cities in Table 1). The original global model is given in Baroyan et al. [lo], but the refinement and application of this work is covered in detail by the first set of algorithms in the previously cited manuscript [2]. This first set of algorithms describes the construction of the forecast on a computer. Rather than give the algorithms here, we instead give the formulation of the model, in detail, along with a description of the methods used for parameter estimation. This is followed by an example of how the model was used to forecast the 1968-1969 “Hong Kong” influenza pandemic based on initial epidemic data from Hong Kong, i.e., the first major city to report the then new variant of influenza A. This is followed by a general discussion of the method as well as of future applications and refinements of the model.
THE MODEL The model is formulated as a system of difference equations in a continuous state space and discrete time domain, where the time unit is one day. Intracity epidemics can occur in each of the n cities in the model. Individuals traveling through the transportation network can spread infection from city to city, creating the intercity pandemic; the cities are located on several continents.
.. -
_,
POP.
_
___
Mean Daily Number
_
of Airline Passengers
-
m
_
;
;
-
_
,
i
TABLE 1 between Modeled Cities during the 1968-1969
Pandemic
-_ .--
Influenza
-
-_,.
6
LEONID
DYNAMIC
A. RVACHEV
AND IRA M. LONGINI,
JR.
VARIABLES
For modeling purposes, the population is partitioned into four explicit, disjoint disease states within each city, i = 1,2,. . , n, which are as follows: x, (t)= number of susceptible individuals on day t, u,( r, t)= number of latent (incubating) individuals on day t who were infected on day t-r, y, ( T, t) = number of infectious (ill) individuals on day t who were infected on day t-r, z, (t) = number of recovered (immune) individuals on day t. The population size, p,, of each city is assumed to be time invariant purposes of modeling influenza. Therefore, we have x,(t)+
g 7=0
u,(7,t)+
5 Y,(~,f)+zr(t)=P,
for every t,
7=0
for the
(1)
where ri is the maximum length of the latent (incubation) period and r2 is the maximum length of the infected period, i.e., the latent plus infectious period, taken from the initial day of infection. The variable z,(t) can be considered redundant and, therefore, a total of n [l + (TV+ l)+ ( r2 + l)] = n (TV + TV + 3) equations are required to model the pandemic process. Insofar as the length of the latent and incubation periods can be considered equal for influenza [l], the latent and incubating states are identical. For epidemiological purposes, two additional dynamic variables are used. They are as follows: w,( t)=number of individuals who became infectious (ill) on day t, i.e., daily morbidity or daily incidence of illness, b,(t)= computed (or modeled) number of new ill individuals reported to the health registry on day t. DISTRIBUTION
OF INFECTION
IN INDIVIDUALS
The infection distribution for individuals is defined for a set of empirically derived probability distributions. For an infected individual, let L, be a random variable for the length of that individual’s latent period, and L, be a random variable for the length of that individual’s infectious period. Also, let .& be a set which is the statistical ensemble of individuals who were all initially infected at the same time. Then the infection states (i.e., latent, infected, removed) of the individuals in S? form a partition with the following discrete probability distributions: f(r)=P(&>r),
7’0,l
g(7)=P(L1<7<&),
r=O,l,...,
T2,
(3)
h(T)=P(Lz<7),
r=O,l,...,
72+1,
(4)
7.. . ,711
(4
MODEL
FOR GLOBAL
7
SPREAD
and f(7) + g( 7) + h(7) = 1, T 2 0; pi and T* are the maximum lengths of the latent and infectious periods, respectively. The distributions f(7), g(T), and h(7) can be thought of as the fractions of individuals in S? who are in the latent, infectious, and removed states, respectively. It is assumed that f(0) = 1 and therefore g(0) = h(0) = 0, since the latent period should be at least one day. Using Equations (2)-(4) and standard probability arguments, we derive the following two transition probabilities:
f(T)-ff(T+l) f(T) ’
y(T)=P(L1
f(T)
‘0,
T=o,l,...,
T1, (5)
h(~+l)--h(7)
S(T)=P(L*~T+~~L~~T
g(T)
’ T=O,I
g(T)‘0 ,..., T2. (6)
Then Y(T) is the probability that a latent individual becomes infectious (ill) on day T + 1, given that the individual was still latent on day T; while S(T) is the probability that an infectious (ill) individual recovers on day T + 1, given that the individual was still infectious (ill) on day 7. INFECTION
PROCESS
Within each city, susceptibles and infectives are assumed to mix homogeneously. The daily infectious contact rate h is defined as the average number of individuals with whom an infectious individual will make sufficient contact (to pass infection) in a day. Therefore, the average number of susceptibles infected by an infectious individual in city i on day 1 is Xx, (t)/p,.The number of infectious individuals on day t is given by
Therefore, the number in city i on day t is
of new latent individuals
u,(0, t) = y TRA NSPOR
TA TION
g ’
due to infectious
individuals
U,(O,t-T)g(T).
7=1
OPERA TOR
Individuals move through the transportation network with daily intensity a,,. We define q, as the daily passenger flow from city i to city j, i.e., the average number of individuals that travel from city i to city j in a 24-hour
8
LEONID
A. RVACHEV
AND IRA M. LONGINI,
JR.
period. The n x n matrix u with elements a,, is assumed to be symmetrical, i.e., u,, = a,,, i,j=l , . . . , n. It is further assumed that the number of individuals in each disease state traveling from a particular city on day t is proportional to the fraction of individuals in each disease state in the city at that time and that the population size of each city p, remains constant over time. Ill individuals are assumed not to travel. Based on the above assumptions, a transportation operator Q is constructed to operate on any dynamic variable A,(t): i-@,(t)]
=A,(t)+
2 /=I
A modification of D has been constructed to prevent premature epidemics in cities due to small (fractional) movements of latent individuals. This modification is discussed in the section on forecasting the 196881969 influenza pandemic. STATE
EQUATIONS
for the global system (i = 1,2,. . , n) are as follows:
The state equations
x,(t+l)=D[x,(t)]-u,(O,t),
(9) r=O,l,...,
u,(7+1,t+l)=[l-Y(~)]52[u,(T,t)], y,(7+1,t+1)=
W,(t+l)=
rl-1,
Y(T>Q[u,(7,t)]+[I--(T)]y,(T,t), [I-6(T)]y,(T,t),T=T1+1,T1+2
2 Y(T)~[U,(T,t)]t
(IO)
r=o,l,...,Tir ,...,
T2-1,
(11)
(12)
T=O
h,(t+l)=Pw,(t+l),
(13)
where p is defined as the fraction of newly ill individuals that are reported to the health registry, i.e., reported illnesses incidence. Note that X and p are the same for all cities. The boundary conditions are given by Equation (7) and y, (0, t) = 0. LOCAL.
INFLUENZA
AND
THE INITIAL
CONDITIONS
The city which experiences the earliest epidemic peak is taken as the index city i,. The epidemic process is first modeled locally, i.e., without the transportation operator, in order to estimate critical parameters from data and to index the starting time of the forecast. The usual data are the observed registered number of ill individuals reported { a,(,( t)} in the initial city. Other
MODEL
FOR GLOBAL
SPREAD
9
data can also be used, such as the total number of reported illnesses up to some point in time, or the reported influenza mortality [14]. The equations for local influenza are
x,& i + 1) = X,“(i> - U,“(O,i)
(14)
1
2 u,o(o,i-7)g(T),
q(O, i) = y 10
(15)
r=l
b,(,(i+l)=P; [f(~)-f(T+l)lu,“(o,t-T).
(16)
r=O
where i is a shifted time domain, starting at i = 0, that is lined up with the actual calendar time domain f. Since the early stages of the epidemic are not observable, the most reliable point to make the alignment of t^and t is that as that time time, t,,, , when {a,,(r)} reaches its peak. If we define i,, when {b,&i)} reaches its peak, then we align times t,,, and i,,, and simultaneously count back in time until we reach that time where i = 0; the corresponding calendar time will signify the initial date of the epidemic in city i,. We initialize the local influenza epidemic by specifying time i = 0 as that time when u,&O, i) = 10m5p,0, since the model only includes cities with populations of at least 100,000. In order to get the full set of initial conditions on u,&O, +), we assume that the first appearance of new latent individuals in city i, follows a geometric progression with parameter p > 1, i.e., x,&t) = x,&O), during the early stages of the epidemic. Then it follows that
u&x i + 1) p= for i not too large. Iterating
(17)
U,“(O,i)
on (17) yields
u,O(o, i - T) = UJO, i) p-T. Using the fact that u,~(O,O) = 10-5p,0 and substituting conditions
u,()(O,- T) =loVp,“p-‘, The initial number
of susceptibles
7=0,1,_..,
(18) into (18) yields the
T2.
(19)
is given by
where 0 < OLI 1 is the initial fraction of city i,‘s population
that is suscepti-
LEONID
10
A. RVACHEV
AND IRA M. LONGINI,
JR.
ble to the agent. (Note: (Yis assumed to be the same for all the cities, and when (Y< 1 some of the population may be immune from a previous epidemic.) In order to initialize the pandemic process, we will define t, as the first day at which at least one latent individual (at any stage of infection) travels from city i, to another directly connected city. Since we are still dealing with the early stages of the epidemic in city i,, we substitute Equation (18) into
which yields
Then it follows from the above definition
of f, that
The forecast is made using the state equations (7), (9)-(13) over the time domain t = t,, t, + 1,. . . , to + T, where T is the forecasting horizon. For the forecast. the initial conditions are set at
(24)
and
X,(b) =ap,. ESTIMATOR
OF FREE
(26)
PARAMETERS
The free parameters h and a are estimated from the ascending epidemic curve in initial city i,. The estimation problem is to find estimates ^x and & which minimize the mean squared deviation between the daily reported morbidity { a,$ t)} and the computed morbidity { b,O(f)]~,~} over the time horizon of the available data. The estimated procedure involves minimizing the mean squared deviation A(k), which is a function of p [given in (17)] in the sense that the least squares estimates ^x and & are found at different values of h.
MODEL
FOR GLOBAL
SPREAD
11
The two-dimensional estimation procedure is reduced to a one-dimensional procedure (to a variation of CL)by addition of the constraint ah = S, where s is a constant. The value of s, at fixed values of CL,is found by substituting (18) into (15) yielding
u(o,t)=* :
u(0,t)p-‘g(7).
(27)
7=1
(Note: The subscript i, has been dropped in this section, since it is clear that we are referring to the initial city.) Simplifying and using the initial condition x(0) = ap in (27) yields .S=
1 72 c P-k(T)
(28)
r=l
Now, in order to estimate X and cx, we let A =X/r and a! = r-i?, where x and 5 are arbitrary numbers (with 15 = s) and r is to be estimated. Then x=x/i and &=?C; also ^x&=s. Generally, we set A=s and Z=l. The estimate F is that value which minimizes the sum of squares
where the summation is taken over all values of t where {a(t)} is defined, and b(t) varies within the class of functions generated by Equations (14)-(16) at various pairs of A =x/r and a = r5 subject to Xa!= s. Thus, we compute b(t) using b(t)=rb(t)
for every t,
where b(t) is found by stepping through Equations and a = 5. The least squares estimator
(30) (14)-(16)
when J =h
(31)
is found by solving the equation
gp-(r)b(r)l’l I
r=i
=o.
(32)
12
LEONID A. RVACHEV AND IRA M. LONGINI. JR.
Since s is a function of p, then so are different “best” estimates of A and a when k[‘), k=1,2 ,..., I, e.g., pikJ E (1.05, 1.06,.. and OLare associated with the fixed value squared deviation
x and &, and we get a class of varying the fixed values of CL,say .,1.50}. Within this class, the X of p which minimizes the mean
(33)
where b(t) is b(r) evaluated at 1 and &, and ? is the total number of days over which {a(t)} is defined. There may be gaps in {a(f)}, since morbidity data may not be registered over weekends and holidays. The algorithm is summarized as follows: 1. 2. 3. 4. 5. 6. 7. 8.
Set up the values ).L[‘], k =1,2 ,..., 1. For the present pL’] compute s from (28) and set Cu= 1, x = s. Substitute x and Z into the state equations (14)-(16) and step through them, generating the sequence {b(t)}. Substitute b(f) into (31) and obtain i. Find 1 =I/?, 8= Gi, and {g(r)} = i{%(f)}. Substitute {b(t)} into (33) and obtain A(p). Return to 2 for next value of pLkl until k = 1. Choose ^x and & associated with the smallest A(k) as the best estimates.
It is useful to plot A(p) against p and to approximate the plot with a smooth curve; the minimum of this curve should be close to the same p, as is the minimum A(p). In those rare cases when this does not hold, one should consult Reference [5]. When other data than {a(t)} are available, the algorithm can be modified accordingly. THL
HEMISPHERIC
SEASONAL
SWING
P&gents causing upper respiratory illness tend to occur during the colder months, when conditions are best for person-to-person spread. Accordingly, the timing of epidemics in cities in the northern and southern hemispheres is affected by the beginning and ending of the colder months. In order to model this :easonal swing in the most parsimonious fashion, the cities modeled are distrbuted into three zones: northern, southern, and tropical, with the variable Sn, coded 1, - 1, and 0, respectively, indicating in which zone city i is locited. For cities in the northern and southern zones, seasonal alterations which affect spread of the agents are assumed to be pronounced; for cities in the tnpical zone they are assumed to be absent. Fo the northern and southern zones, the dates Ti and T2 correspond to the btginning and end of the epidemic season in the northern zone. Con-
MODEL
FOR GLOBAL
13
SPREAD
versely, T, and Tl correspond to the epidemic season in the southern zone. The daily infectious contact rate is reduced by a factor of 0.1 during the nonepidemic period, producing the following seasonal pattern in each zone: Northern
zone (Sn, = 1):
A+
Southern
Tropical
x
for
TV
0.1X
for
TV [T,,T,).
0.0
for
f~ [T,,T,),
x
for
TV [T,,T,).
zone (Sn, = - 1):
zone (Sn, = 0): A+A
For influenza pandemics, the calendar October and 1 April, respectively. ALGORITHMS
[T,,T,),
for all t.
dates for Tl and
T2 are set at 1
FOR COMPUTATION
The forecast is accomplished by stepping through the first set of algorithms given in [2]. We give a brief description of these algorithms here. Algorithm 1 (infection process): Sets up the distributions (S), (6), given the distributions (2)-(4). Algorithm 2 (local influenza): Steps through Equations (14)-(16), (19)-(20) in the initial city i,. Estimates the free parameters LY,X from Algorithm 3 (free parameters): data for the initial city i, using the algorithm given after Equation (331. Computes t, from (23) and sets mitial Algorithm 4 (initial conditions): conditions for algorithm 6 using Equations (24)-(26). Algorithm 5 (transport operator): Applies the operators E and D as given by Equations (37) and (8) respectively. Algorithm 6 (global influenza): Steps through Equations (7) (C)-(13), i=1,2 ,...,n. A seventh coordinating algorithm also is given to direct the timng and movement of flow as the computer passes through the above six algorithms. Given this set of algorithms, it is fairly easy to write down the compllter code in any standard computer language, e.g., FORTRAN or PL/I.
14
LEONID
FORECASTING
OF THE 1968-1969
A. RVACHEV
AND IRA M. LONGINI,
INFLUENZA
JR.
PANDEMIC
In July 1968, a new strain of influenza A, i.e., influenza A(H3N2), caused a major epidemic in Hong Kong. This influenza virus, which probably originated in southern China, subsequently swept the globe, terminating its first wave late in 1969 [15]. It was subsequently given the name Hong Kong influenza. The model described in the previous section is used to forecast the influenza pandemic of 1968-1969 in 52 cities based only on background data and the first half of the epidemic curve in Hong Kong. THE
DATA
The population sizes for each city, p,, in 1968, were taken from demographic reference books and are given in Table 1. Also given in Table 1 are the estimates for the transportation matrix taken from international air transport statistics [16-221. The upper triangle of this symmetric matrix has 440 nonzero entries; e.g., an average of 68 individuals traveled between London and Sydney during each 24hour period in 1968. The distribution of the length of the latent period, f(7), was the same as that used by Elveback et al. [23] for simulating influenza A epidemics. The distribution of the length of the infectious period, g(7), was taken from past Soviet work [5, 241. Then h(7) was computed by subtraction, and Y(T) and &(T) were computed using Equations (5), (6), respectively. The distributions are given in Table 2, and we see that T, = 2 and TV = 8. ESTIMATION
OF
FREE
PARAMETERS
Hong Kong was taken as the initial city, and therefore only source of morbidity data for Hong Kong is that of the weekly morbidity from nine state clinics is given. sequence {a(t)} from the weekly totals does not produce
we set i, = 10. The Chang [25], where Extrapolating the a stable sequence
TABLE 2” Infection
Distribution
f(T) g(7) h(T) Y(T) S(7)
Process Distributions
0
1
2
3
4
5
6
7
8
9
1.00 0.00 0.00 0.30 0.00
0.70 0.30 0.00 0.70 0.10
0.20 0.77 0.03 1.00 0.20
0.00 0.82 0.18 0.00 0.34
0.00 0.54 0.46 0.00 0.44
0.00 0.30 0.70 0.00 0.50
0.00 0.15 0.85 0.00 0.60
0.00 0.06 0.94 0.00 0.83
0.00 0.01 0.99 0.00 1.00
0.00 0.00 1.00 0.00 0.00
“From Belova et al. [2]. ‘Tlrne since initial day of infection.
MODEL
FOR GLOBAL
SPREAD
15
for estimating h and a, although such a sequence is useful for approximating p. Chang [25] reported 250 notifications during the first week of the epidemic and 1150 during the second week. Using this information, we extrapolate fi = (1150/250)“’ = 1.2436. Thi s value of p compares well with the value estimated for the influenza A(H3N2) epidemic during January-February of 1969 in the Soviet city of Kuibyshev. From the Kuibyshev data, F = 1.24 was estimated from a proper sequence {a(t)} (see, [5, p. 2041 and, [24, p. 531). This estimate was subsequently used to produce a very good forecast of the daily influenza morbidity in 40 Soviet cities (see, [5, p. 4771 and [24, p. 561). The closeness of these two estimates is remarkable considering that different methods were used to estimate p for the two cities. We further note that Kuibyshev is in the northern zone while Hong Kong is in the tropical zone, and that the estimates were made during different seasons. The estimate F = 1.24 was thus taken as the best estimate and used for the estimation of h and cy. Substituting p =1.24 into Equation (28) yields s = 0.641 and therefore Cu= 0.641. The population of Hong Kong in 1968 was about 3.9 million people. Thus, from Equation (19) we set St0 (0,O) = urO(O,O)h, = 39, and we also computed the initial conditions ~(0, - r), r = 1,2,. . . ,8, from (19). Then was computed via Algorithm 2. From (13) the sequence {W(i)} = { w(i)h6} it is easy to see that the sequences W(i) and b(i) peak at the same time, I,,,, . This time was found to be i,, = 44. The calendar time, I,, , when (a(t)} reached its peak was extrapolated to be 27 July 1968. By aligning t,, with and counting back, the calendar time corresponding to i = 0 was found Lx to be 13 June 1968. In order to estimate X and a, the computed (from Algorithm 2) number of reported illness in Hong Kong up to the peak of the epidemic, Sb( i,,) = C&b(i), was set equal to the observed number Sa( t,,). The latter number was estimated to be 292,000 individuals from Chang’s data [25]. From Equations (13), (30) we have
b( t) = B4 t>
(34)
w(t)=fi(f),
(35)
and
where w(t) is computed using Algorithm 2 with X =x = 1 and OL= 5 = 0.641. From Equations (34) and (35) we have
,,=b(r)= w w(t)
c,>
~~(bulX)’
(36)
LEONID A. RVACHEV AND IRA M. LONGINI, JR.
16 and we estimate
r by
w cmx) i = jL%i?(f,,) =
(0.3;:~~000) z 0.95,
and ^x=x/i=l.OS, &=?&=0.61. The value of t, was computed from Equation (23) and was found to correspond to 12 July 1968. Accordingly, the initial conditions (24)-(26) were set. The forecasting time horizon was set at T = 440 days, which yields a forecast up to 26 September 1969. However, since data were used from Hong Kong up to 27 July 1968, the actual forecast was over a 425 day period. During the first attempt at modeling the 1968-1969 influenza pandemic it was observed that epidemics in cities began too early because of small fractions of latent individuals, e.g., 0.001 latent individuals, passing through the transportation network and initiating epidemics. The problem did not arise previously when the model was applied within the USSR territory and the duration of the epidemic process usually did not exceed three months. This problem was circumvented by constructing an operator, E, that allows latent individuals at a particular stage to be present in a city i only if the approximate expected sum of these individuals in city i and those entering from the other connected cities is at least one. The form of this operator is as follows:
I
!2[[u,(~,t)] if
E{fhht)l)
=
Q[u,(T,t)]
Cf(T)21, 7=0
0
if
Q[u,(~,t)]
i f(~) 7=0
(37) cl.
The physical interpretation of E is most apparent if we assume that the epidemic is peaking at time t = t* in city i. This implies that U,(O,1*)=~,(O,t*-I)= since
71
...
is small. Then, the total number
=U,(0,t*-T71),
(38)
of latents in city i at time t = t* is
2 U.(T,f*)=,~o~(T)u,(O.l*-T). (39)
T=O
Substituting
(38) into (39) yields ; 7-o
U,(T,t*)
=
U,(0,f*)7~oj(T).
(40)
MODEL
FOR GLOBAL
SPREAD
17
Thus, the factor C:,,f( T) is the expected number of latents in all classes for each new latent when the epidemic is peaking. In order to forecast the 1968-1969 influenza pandemic, Equations (lo)-(12) were modified by applying the operator E to all 3[ u, (7, t)]. THE FORECAST The results of the forecast in terms of predicted daily reported morbidity, { b,(t)}, versus observed daily reported morbidity, { ai( t)}, are given in Figure 1. The observed values were taken from World Health Organization sources [26-301. For some cities the data were quite detailed, but for many cities only the duration and peak period of influenza activity were available. The most conspicuous disadvantage of the data is the absence of empirical information on influenza activity in many of the cities modeled. Probably, in the majority of cases, this means that there was no epidemic in these cities. It is hard to imagine that no information would have reached the World Health Organization had epidemics occurred at least in the territories around the cities in question. Only Havana was predicted not to experience influenza activity over the time period of the simulation. Nevertheless, the occurrence of cities not experiencing significant influenza activity during the pandemic is consistent with the basic premises used in building the model, since the infectious contact rate in cities is decreased during the “nonepidemic” period in the northern and southern zones. In Figure 1, we see that on the space scale of the globe and the time scale of a year the predicted epidemic course generally corresponds to the outline of the real course of the pandemic.
DISCUSSION The original global model was formulated as a system of integrodifferential equations with partial derivatives [lo]. In order to develop the first step of computer algorithms briefly discussed in this paper, the system of difference equations (7)-(13) was formulated as the discrete analog of the original system of equations. A second set of algorithms also has been developed to forecast the global spread of an infectious agent which spreads like the influenza virus but exhibits a higher case fatality rate like the agent of Legionnaire’s disease. In this case, the rapid isolation of exposed and infective individuals is important. Accordingly, the model is expanded to incorporate isolated states, and vital dynamics are added to account for high case fatality rates, i.e., p, becomes a dynamic variable. More advanced methods are used to estimate the free parameters A and (Y where p is not treated as a constant but rather as a continuous function. Further details are available in [2]. With regard to the forecast of the 1968-1969 influenza pandemic, the most time-consuming task was the construction of the transportation matrix.
0
26
I
,...
Kinshasa
I * Empirical dataamjointwit/lIhesurro”“di”g repion.
Cairo
..r
.,
..-
,
..
. . . . . . --++++,+--...
.. ....-- .++f+r+
..,.
... .. .
1
When a symbol falls directly
on next puge)
indicated.
(Continued
on a border
between
months.
this indicates
occurred
during
the four days two days in each month.
b,(t) < 10, - for 10 I b,(r) I 100, + for 100 > b,(r), and i_ when the peak in morbidity
FIG. 1. A schematic plot of the forecasted h(f) and actual a(r) course of the 1968-1969 influenza pandemic. Each of the following symbols represents the daily forecasted morbidity incidence per lo5 over four calendar days: 0 for
1
26
z _-
3
B k s
.
peaks were not reported.
(Redrawn
.
from Belova et al. 121.)
. _
solid line indicates that significant influenza activity was reported from the city or region indicated; parallel lines indicate that a peak in morbidity was reported. If only the peak period was reported, then the parallel lines are surrounded by short solid lines. The arrows indicate the direction, in time. towards the period of peak activity when
If the symbol touches the border from the left or from the right, this indicates that three of the four days were in the month to the right or left, respectively. The prevalence of the actual pandemic is given by the lines above the symbols. A
FIGURE1 (continued)
20
LEONID
A. RVACHEV
AND IRA M. LONGINI,
JR.
For this first application, only the daily intercity airline passenger turnover was used for the 52 cities modeled. Numerous technical difficulties had to be overcome. For example, many passengers going from city i to j may have a short intermediate stop in some other city, k. However, the airline statistics would record the flow from cities i to k and k to j separately, thus inflating the estimated number of individuals who are actually staying in the cities involved. This difficulty, as well as others, had to be resolved by considering that the transportation matrix is not a separate empirical source of data but rather a set of parameters that have to work in concert with the state equations and the estimated population sizes of cities of departures. In order to make the forecast more realistic, it is necessary to expand the number of cities modeled from 52 to at least 152. It is also desirable to refine the transportation matrix to include other sources of transportation such as ground transportation and the movement of freight personnel. There also should be a separate transportation matrix for each day of the week in order to simulate the weekly periodicity of the world timetable of airflights. Another complication involves the modeling of infected individuals. It was assumed that all infectious individuals were ill. However, a certain fraction of infected individuals are asymptomatic for influenza [23]. Since it would be hard to separate this factor (i.e., pathogenicity) from the estimation of /3, asymptomatic individuals were not modeled separately. The model with asymptomatic individuals is given elsewhere [13]. In conclusion, the model was successful in forecasting the broad patterns of the geographic spread of Hong Kong influenza. This lends credence to the theoretical correctness of the concepts used to build the model. The length of the forecasting interval, which proceeded from 27 July 1968 to 26 September 1969, may indeed be a milestone since hardly anybody has ever predicted, even roughly, any daily course of events in the world society for such a period-about 425 days. GLOSSARY Agent.
OF EPIDEMIOLOGICAL
TERMS
An entity which causes disease (e.g., an influenza virus). Unusually frequent occurrence of disease compared to background level (e.g., outbreak of disease in a single city). Incubation period. The period of time beginning when an individual first harbors an agent and ending when that individual begins to manifest symptoms of disease. Infected. Harboring an agent of disease. Infectious period. The period of time when an infected individual can spread the agent to other individuals. This period begins at the end of the latent period and usually ends with recovery. Lutent period. The period of time beginning when an individual first harbors an agent and ending when that individual becomes infectious. Epidemic.
MODEL
FOR GLOBAL
21
SPREAD
Pandemic. An epidemic which includes significant outbreaks in many regions of the world. Puthogenicity. A measure of the ability of an agent to cause illness once the agent has infected an individual. L. A. Rvachev would like to express his gratitude to V. E. Kvitka, Deputy Director and M. L. Petryanova, Senior Researcher, of the Research Institute of the Ministty of Civil Aviation, U.S.S.R.; E. N. Chernvochkin, Department Manager of Automated Control Systems, of the Sklifosovsky Institute; V. A. Becketova, G. P. Busuyok, Y V. Gudkova, Research Workers at the Gamaleya Institute; and V. Z. Sotnikov at the Institute of Applied Mathematics of the Academy of Sciences of the U.S.S.R. I. Longini would like to thank P. E. M. Fine for introducing him to the work described here. REFERENCES 1 2
E. D. Kilbourne, The lnfluenzu L. A. Belova. J. W. Donovan, Rvacbev, modelling
V. A. Shashkov, (Part
I-influenza),
and
Viruses und Influenzu, Academic, P. E. M. Fine, D. W. Fraser, V. I. Vasilyeva, Unpublished
Experiment
manuscript
New York, 1975. M. B. Gregg, L. A. on pandemic
(in Russian
and
process English),
4
1983. 0. V. Baroyan and L. A. Rvachev, Deterministic epidemic models for a territory with a transport network (in Russian), Kiherneriku, No. 3, pp. 67-74 (1967). L. A. Rvachev, Computer modelling experiment on large-scale epidemic (in Russian),
5
Dokl. Akud. Ncluk SSSR 180(2):294-296 (1968). 0. V. Baroyan, L. A. Rvachev, and Y. G. Yvannikov,
3
6
7
Modeling und Forecusting
of
Influenzu Epidemics for Terriroty of the U.S.S.R. (in Russian), Gamaleya Inst. of Epidemiology and Microbiology, Moscow, 1977. 0. V. Baroyan, L. A. Rvachev, K. Frank, et al., Model of influenza epidemics spread through the territory of Bulgaria and its significance (in Bulgarian), J. Eprdemiol. Microhiol. Jnfecr. Dis. (Soficr) 15(3):168-173 (1978). 0. V. Baroyan, V. V. Basilevsky, V. V. Ermakov, et al., Computer modeling of influenza epidemics for large-scale systems of cities and territories (in Russian and English), Working Paper for W.H.O. Symposium on Quantitative Epidemiology, Moscow, 23-27 Nov. 1970.
8
0. V. Baroyan, L. A. Genchikov. L. A. Rvachev, and V. A. Shashkov, large-scale influenza epidemic modelling by means of a computer, Epidemiol. Assoc. 18:22-31 (1969).
9
0. V. Baroyan, L. A. Rvachev. V. V. Basilevsky, et al., Computer modelling of influenza epidemics for the whole country (U.S.S.R.), Ado. Appl. Pro/x&. 3~224-226 (1971). 0. V. Baroyan. G. A. Mironov, and L. A. Rvachev, An algorithm modeling global epidemics of mutant origin, Progrumming and Comput. Software 6(5):272-277 (1981) English transl. from Programmirovunie, No. 5, pp. 73-79 (1980) (in Russian). N. T. J. Bailey, The Muthemutical Theoy of Infectious Diseases and its Applrcutrons, Charles Griffin, London, 1975.
10
11
An attempt at Bull. Inrernur.
22
LEONID
A. RVACHEV
AND IRA M. LONGINI,
JR.
12
P. E. M. Fine, Background paper, in Influenza Models - Prospects for Development and Use (P. Selby, Ed.), MTP Press, Boston, 1982, pp. 18-85.
13
L. A. Belova et al., New influenza epidemic model for the U.S.S.R. territory Russian), presented at Seventh International Symposium on Influenza, Lenningrad, Dec. 1982.
14
C. C. Spicer, 35:23-28
15 16
The mathematical
modelling
of influenza
epidemics,
hit.
(in 10
Med.
Bull.
(1979).
W. C. Cockbum, P. J. Delan, and W. Ferreira, Origin and progress of the 196881969 Hong Kong influenza epidemic, Bull. M/id. Hlth. Org. 41:345-348 (1969). Air Trunsport Stntistics,
Australia,
international
1980. Dept.
Air Transport,
of Trans.,
Canberra,
June 1981.
17
J. W. Donovan, ing Paper.
Australian
18 19 20
ARC World Airyvs Guide: (a) Part One, (b) Part Two, ABC Travel Guides, Ltd., Dunstable. Beds, U.K., Aug. 1981. Official Airline Guide, North American ed., OAG Inc., Oak Brook, Ill., Oct. 1981. Officrul Airline Guide, Worldwide ed. OAG Inc., Oak Brook, III.. Dec. 1981.
21
Data
22
G. P. Busuyok. V. A. Becketova, and V. P. Charlamov, The methods of epidemiological modelling of the U.S.S.R. passenger traffic (in Russian), presented at Seventh International Symposium on Influenza, Leningrad, Dec. 1982.
23
L. R. Elveback, J. P. Fox, E. Ackerman, et al.. An influenza immunization studies, Amer. J. Epid. 103:152-165 (1976).
from the International
V. Baroyan
International
Airport
Civil Aviation
and L. A. Rvachev.
24
0.
25
Moscow, 1977, pp. l-64. W. K. Chang, National influenza
Traffic,
Organization
Mathematics
Dynamics
1968-1980.
Work-
for the years 1968-1982.
simulation
and epidemiology
model
for
(in Russian),
in
Znanie,
Org. 41:349-351 26
International
experience
in Hong
Kong,
1968,
Bull.
Wld.
Hlth.
(1969).
Conference
on Hong Kong Influenza,
Bull. W/d. Hlth.
Org. 41:345-414
(1969). 27
Week!v Epidemiological
Record. No. 34-52. Wld. Hlth.
28 29
Week!v Epidemrologicul
Record, No. l-40,
Org.. Geneva, 1968. Wld. Hlth. Org.. Geneva, 1969.
W/d.
Hlth.
Stutist.
Rep. (Geneva)
22. No. 7 (1969).
30
Wld.
Hlth.
Statist.
Rep. (Geneva)
23, No. 7 (1970).