lg This is to certify that the thesis entitled SYSTEMS MODELLING, ANALYSIS AND SIMULATION OF TYPE A INFLUENZA EPIDEMICS presented by Roy Gardner has been accepted towards fulfillment of the requirements for Ph. D. degreein Systems Science $10 [Jag/Q Major professor Date February 19, 1980 0-7639 __ __ , _ 7 . OVERDUE FINES ARE 25¢ PER DAY PER ITEM Return to book drop to remove this checkout from your record. SYSTEMS MODELLING, ANALYSIS AND SIMULATION OF TYPE A INFLUENZA EPIDEMICS By Roy Gardner A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Electrical Engineering and Systems Science 1980 ABSTRACT SYSTEMS MODELLING, ANALYSIS AND SIMULATION OF TYPE A INFLUENZA EPIDEMICS By Roy Gardner A pandemic results from the emergence of a new subtype of type A influenza virus. Following the pandemic, epidemics occur in two to three year cycles due to the structural change of the same subtype virus. The graphic representa- tion of disease incidence from one pandemic to the next resembles underdamped oscillation. A review of the liter- ture reveals no existing dynamic model of influenza epi- demics which depicts the behavior of the virus and its interaction with the pOpulation. This paper develops a model for such a system using systems theory and application. As a basis for describing the dynamic behavior, three epidemic models are analyzed in the context of systems analysis and simulated. All simulations in this research are performed using DYNAMO. Variables which are of interest to epidemiologists, such as the peaks and the durations of epidemics, are defined, and the formulas for computing these variables are derived. The numerical values obtained from analysis are compared with the simulation results. The influenza epidemic model develOped using the basic epidemic models is a nonlinear system with four state variables, one of which is the structural state of the virus. _The virus was modelled using the idea of a logistic curve, the asymptote of which is the immunity level of the pOpulation. VThe line- arized model was used to analyze system behavior. With a proper choice of the parameter associated with the-virus, the system provides an epidemic curve which resembles under- damped oscillation. A range of the parameter required for this underdamped oscillation was computed using a sufficient condition, expressed in terms of the coefficients of the characteristic equation, for the third order linear system to be stable, oscillatory, and underdamped. This sufficient condition was derived based on the relationship between the roots location and the coefficients of a cubic equation. Simulation results of the system for four test cases are presented and compared. With this model three epidemics can be shown. The attempt to fit the model outputs with the historical data on the mortality rate of three selected epidemics of recent times was carried out as model vali- dation. The model was modified to include the population growth and the seasonal variation of influenza outbreaks. Simulation results show a reasonable fit, considering the complexity of the real world system. To minimize the severity of the epidemic, an immunization model was developed as a subsystem of a control system. This model consists of the vaccinated pOpulation as the state variable. The state equation was constructed using the idea of a goal seeking curve. A comparison of simulation outputs of the system without immunization and with immunization is given. Simulation shows that, with this control system, in- stead of the expected epidemic, an endemic state results. This control system, however, is an idealized system where only individuals in the susceptible pOpulation are immunized, and it does not consider the efficacy of the vaccine. To interpret the simulation results more realistically, a simple formula is given to compute the percentage of the total pOpulation that requires immunization to prevent the epidemic. This percentage varies depending on the efficacy of the vaccine. To Martha and Kenzo. ii ACKNOWLEDGMENTS This work is a result of the culmination of knowledge I have acquired from each member of my guidance committee: I wOuld like to thank Dr. Robert Barr, my chairman and thesis advisor, for his direction and for his excellent teaching. A special thanks is due to Dr. Robert Moon, Department‘ of Microbiology and Public Health, for patiently teaching me about the fascinating subjects of viruses and epidemics. I wish to thank Dr. Gerald Park for guiding me toward a more interdisciplinary approach to this dissertation. I also wish to express my gratitude to Dr. Connie Shapiro, Department of Statistics, for her constant concern and encouragement. Finally, I want to thank my wife Martha for editing my manuscript and coming up with suitable sentences at times when nothing seemed to come out of my head. iii TABLE OF CONTENTS LIST OF TABLES 0000.0000.000.000.000000000000000000.00. LIST OF FIGURES 000.00.00.000000000000.0..000.00.000.00 CHAPTER 1.1 1.2 1.3 CHAPTER 2.1 2.2 2.3 CHAPTER 3.1 3.2 3.3 CHAPTER h.l h.2 l-INTRODUCTION 0000000000000.0.0.0.0000.0..00 Statement of the Problem ..................... Description of the Influenza Epidemic SYStem 00000.00000000000000.000000000000.00... Purpose, Background and Methodology Of ResearCh 000000.0..000.000.000..000.0000... 2 - BASIC EPIDEMIC MODELS ..................... Simple Epidemics Model ....................... General Epidemics Model ...................... Endemic Model ................................ 3 - INFLUENZA EPIDEMIC MODEL .................. Formulation of the State Model ............... Causal Loop Analysis ......................... System Analysis .............................. 3.3.1 Equilibrium Points Analysis ........... 3.3.2 Linearized Model ...................... 3.3.3 Computation for Range of Drift Rate 00000000.00.0000...00.000.000.000. A - SYSTEM SIMULATION AND MODEL VALIDATION .... AnalYSis Of Drift 000.00.000.00000000.0000000. Threshold of Epidemics ....................... iv Page vi vii 12 13 21 3h L6 #6 50 51 56 58 68 71 76 79 Page A.3 Causal L00p Diagram of Parameters ........... 82 h.h Model Validation ............................ 82 CHAPTER 5 ‘ IMMUNIZATION MODEL AND CONTROL SYSTEM .... 91 5.1 Causal Analysis ............................. 92 5.2 State MOdel of the Immunization System ...... 93 5.3 Simulation Results .......................... 97 CHAPTER 6 - SUMMARY AND CONCLUSIONS .................. 102 APPENDIX A - DERIVATION OF THE FORMULA FOR THE FUNDAMENTAL MATRIX OF THE 2X2 MATRIX A IN WHICH THE EIGENVALUES ARE COMPLEX .... 106 APPENDIX B - DYNAMO LISTINGS.......................... 108 BIBLIOGRAPHY .0000.00.000.000.000000000.0000000...0000 11L Table 1.1 01.2 4.1 L..2 z..3 5.1 B.1 3.2 3.3 B.h 8.5 3.6 LIST OF TABLES Page Major Epidemic and Pandemic Years and Their Intervals 00.0000.00.000.00.0000000000... 2 Subtypes of Human Influenza A Viruses and the Year ISOlated 0000.00.000.000000000000.000. ’0 Relationship Between Virus Drifts and Epidemic Peaks .00.0000.000.0.0000000000000000. 78 Annual Birth and Death Rates in the United States (per 1,000) ocoooooooooooooooooooooooooo 87 Monthly Mortality Rates in the Three Epidemics Used in MOdel Validation ............ 88 Percentage of Total POpulation Required to Be Immunized 00.000.000.0000000000.00.00.00...0 101 Simple Epidemics Model ........................ 108 General Epidemics Model ....................... 109 Endemic Model ................................. 110 Influenza Epidemic Model ...................... 111 Model Validation .............................. 112 Immunization MOdel 00000000000000.00000000.000. 113 vi Figure 1.1 1.2 2.1 2.2 2.3 2.h 2.5 2.6 2.7 2.8(a) 2.8(b) 2.9 2.10 2.11 2.12 2.13 2.1h 3.1 3.2 3.3 3.h LIST OF FIGURES Page Diagram of Influenza Virus Section ............ 5 The Pattern of Influenza Pandemics ............ 10 Causal Loop Diagram of Simple Epidemics ....... 15 Block Diagram of Simple Epidemics ............. 16 Flow Diagram of Simple Epidemics .............. 16 Simple Epidemics MOdel ........................ 17 Causal Loop Diagram of General Epidemics ...... 23 Block Diagram of General Epidemics ............ 25 ,Flow Diagram of General Epidemics ............. 25 General Epidemic Model: p LOO, cl = 950 .... 26 400, c1 = 350 .... 27 General Epidemic Model: p Phase Portrait of General Epidemics Model ..... 28 Causal Loop Diagram of Endemic Medel .......... 36 Block Diagram of Endemic Model ................ 38 Flow Diagram of the Endemic MOdel ............. 39 Endemic Mbdel: p1 = 800, p2 = 62.5 ........... LO Endemic Model: p1 = #00, p2 = 12.5 ........... Ll Logistic Curve of the Virus Strain ............ L8 Causal L00p Diagram of Influenza Epidemic MOdel 00000000000000.000000...0.000000000000000 52 Simplified Causal L00p Diagram of the Influenza Epidemic MOdel 00000000000000.0000... 53 Block Diagram of Influenza Epidemic Model ..... 5h vii Figure 3.5 3.6 3.7 h.1 h.2 h.3 h.h h.5 h.6 h.7 h.8 5.1 5.2 5.3 5.h 5.5 5.6 5.7 viii Flow Diagram of Influenza Epidemic MOdel ...... System Response of a Third Order Linear system 0.0000000000000000000.00000000000000000. Curves of Cubic Polynomials Having One Real Root and Two Complex Roots, (a) With No Relative Maximum or Relative Minimum, (b) With a Relative Maximum and a Relative Minimum 000000000000000000000000000.0000.00.00. Influenza Epidemic Model With the Drift Rated=1 .000.00.0000000000000.00000000.00.00 Influenza Epidemic Model With d 05 0000000000 03 0000...... Influenza Epidemic Model With d Influenza Epidemic MOdel With d .1 .......... Drift and Epidemic Curve ...................... Non Occurrence of Initial Epidemic ............ Weekly Mortality Rate of Influenza- Pnemonia 00.00000000000000.00000000000000.0000 Model Validation: Graphical Comparison of mortality Rates 0000000000000...0.0.0.0000.0.00 Block Diagram of the Control System ........... Causal L00p Diagram of the Control System ..... The Goal Seeking Curve of the Vaccinated Population 0.0.00.0000...00.00000000000.0000.0. Immunization Period ........................... Simulation Outputs Before Immunization ........ Vaccination Rate and Total Vaccinated Population 0..00...0.000.000.0000....000000.0.. Simulation Outputs After Immunization ......... Page 55 65 67 72 73 7h 75 77 81 8h 90 91 93 9h 96 98 99 100 CHAPTER 1 INTRODUCTION 1.1 Statement of the Problem In 1979 smallpox was declared eradicated from the face of the earth. Once—devastating infectious diseases such as cholera, typhoid fever, and yellow fever are now under con- trol through improved sanitation, vaccination, and quaran- tine. Yet, control of influenza epidemics and pandemics (world wide epidemics) is not in sight. This is because the infectious agent, a virus, is constantly changing with the resultant consequence of invading a pOpulation which has little or no immunity for the changed agent. This disser- tation concerns the deve10pment, analysis, and simulation of models of influenza epidemics. Influenza can be just as devastating as cholera. It is estimated that the pandemic of 1918 killed 20 million peopb world wide and over 500,000 in the United States in a period of just a few weeks. The occurrence of pandemics and major epidemics in recent times (Kilbourne, 1975, p.h9h) is as given in Table 1.1. Table 1.1. Major epidemic and pandemic years and their intervals. Pandemic Major Epidemic Pandemic Interval Interval 1889 + 1900 11 1918 + 29 18 1929 11 19h6 + 28 17 1957 + 11 11 1968 + 11 11 + indicates occurrence of an unquestioned pandemic. It is interesting to note the recurrence of the 11 year intervals in the above table. The exact reasons for the cyclical prOperties of the epidemics are not known, but the relationship between the immunity level of the population and the occurrence of an epidemic is an important factor to be considered. Another interesting theory of the cyclical nature of pandemics is the return of the same virus subtype every sixty years. The pandemic of 1918 is believed to have been caused by the "swine flu" virus. Isolation of the swine flu virus from soldiers at Fort Dix, New Jersey, in 1976 prompted a national immunization program in anticipation of another pandemic. The fact that the pandemic did not materialize points out the need for a model to more accu- rately forecast future influenza epidemics. Such a model requires an extensive knowledge of the basic three elements in the system -- the virus, the host, and the environment -- and the complex interactions among them. In many socio- economic and ecological systems, the degree of 3 uncertainty, in essence, determines the validity and hence the usefulness of the model. In the case of the influenza epidemic system, this uncertainty is compounded by the limited knowledge of the virus itself. As the secrets of the virus are unraveled, prediction and con- trol of epidemics may become reality. 1.2 Description of the Influenza Epidemic System To build a realistic system dynamics model of any phenomena, it is fundamental to have a good understanding of the real world system the model is to represent. The following is a brief but essential description of the various elements and interactions involved in influenza epidemics as they are understood today. Influenza is an infectious disease caused by a class of viruses called myxoviruses, or interchageably known as influenza viruses. There are three types of influenza viruses, types A, B, C, which are classified according to the immunity each produces by infection. Type A has sub- types, and each subtype has variants. Type B has no sub- type, but has variants. Type C virus is uncommon and does not cause epidemics. Type A is the only type found in animals and birds as well as humans. Epidemics of type A occur every two to three years, and epidemics of type B occur every four to six years. Appearance of a new subtype of type A can result in a pandemic. In fact, one author defines an influenza l, pandemic as an emergence of a new subtype (Beveridge, 1977). A remarkable fact about the subtypes of type A virus is that when a new subtype appears, the old sub- type rapidly disappears. Hence, only one subtype of type A predominates at any given time. However, the type B virus can coexist with type A virus. The table below (Volk, 1978, p. 5h2) lists the subtypes of human viruses isolated up to the present time: Table 1.2. Subtypes of human influenza A viruses and the year isolated. Year Isolated Subtype 193a HONl 1947 HlN1 1957 H2N2 1968 HBNZ In immunological terms, the change in subtype is called antigenic shift. In simple terms, antigenic shift is a major change in the type A virus which may lead to a pandemic, while antigenic drift reflects minor changes in the type A or B virus which may cause more localized epi- demics. Understanding the structure of the virus makes it possible to understand shifts and drifts: The virus is 80-100nm in diameter and is roughly spherical in shape. It consists of a helical symmetric core which contains the genetic material of the virus called ribonucleic acid (RNA) and two types of spikes, hemagglutinin and neuraminidase; The RNA is in eight separate pieces, each being a single genetic unit. The Spikes are attached to the membrane which encloses the core (Beveridge, 1977, p. 69). N = Neuraminidase H = Hemagglutinin C = Core M = Membrane Figure 1.1. Diagram of influenza virus section. When the virus attacks a person, it first must attach to host cells. This attachment is facilitated by the H and N spikes and is necessary for the virus to pene- ‘trate the host cells. Once inside the cell, multiplica- ‘tion.of the virus takes place through the process of translation and transcription of the genetic codes. An antigen is defined to be any substance which causes the body to produce antibodies. There are four antigens in the influenza virus: the two spikes H and N, the core, and the membrane. Antibodies produced in reSponse to the core antigen and the membrane antigen are not pro- tective and do not prevent infection. The core antigen is significant, however, because the types A, B, and C are classified according to the differences in the antibodies induced by the core antigen. Antibodies induced by the N spike are not protective, either, but they may play some role in reducing the spread of infection, since they inter- fere with the release of the virus from the infected cell. Antibodies induced by the H spike are protective and neutralize infectivity of the virus by coating the surface of the cell so that the H spikes cannot be attached (Volk, 1978, p. 5&0). Basically, the coating of the surface of the cell by the antibodies results in immunity to the specific virus. However, if the chemical composition, and hence the anti- genicity, of the H or N spikes change, the body does not recognize the change, and the person can be infected again. Slight changes may occur by mutation during the multiplication process, and these changes in the spikes are called antigenic drift. If there is a major change so that there is negligible immunity in the pOpulation, then the change is called antigenic shift. Antigenic shifts are believed to result from the creation of a hybrid between human and animal influenza viruses. This is the recombination theory which states that a human influenza virus may infect lower animals or 7 birds and then be recombined with an animal virus present in the animal host to produce a new subtype of the human influenza (Laver, 1979). The recombination theory plays a significant role in the control of influenza epidemics. First, it implies that future pandemics will be difficult to control since vaccines cannot be prepared by simply manipulating the human influenza viruses. Second, it widens the whole dimension of the influenza epidemic system by involving the ecology of animals and birds. Early identification of a new variant or subtype is crucial for the preparation of a prOper vaccine to immunize against a possible epidemic. The World Health Organization maintains a world-wide network for surveillance of in- fluenza viruses. At present, early detection and vaccin- ation seem to be the only effective tools against influenza epidemics. 1.3 Purpose, Background and Methodology of Research The purpose of this research is to build a large- scale or aggregated,model of an influenza epidemic system using systems theoretic concepts which shows the relation- ship between the dynamics of the epidemic and the changing structural state of the influenza virus. In the systems approach, the identification of the virus and vaccination are controllable inputs to the influenza epidemic, while shifts and drifts are uncontrollable inputs. Exogeneous inputs include seasons, susceptibility due to genetic fac- tors and standard of living, and general health and ages which affect morbidity and mortality rates. POpulations with different immunity levels and the number of suscep- tibles and the number of immunes may be considered design parameters. Desirable outputs are a high level of immunity in populations and low morbidity and mortality rates. Complications such as Guillain-Barre syndrome and high economic costs are undesirable outputs. In reality many of the variables identified above are probabilistic or stochastic in nature. For instance, not all individuals who have contact with an infective person come down with the disease. The complexity of the human body and its interaction with the invading virus can never be adequately described by mathematical expressions. In a stochastic approach, disease incidence at any given time period is associated with a probability. Stochastic models are more apprOpriate if one is investigating the system behavior in which the statistical fluctuation is significant, such as in infectious outbreaks involving a few individuals. In a large-scale model this statistical fluctuation averages out and becomes less significant. Deterministic models are adequate for describing epidemics involving a large number of pe0ple and will be used in this research. Epidemic models using “system dynamics" techniques have been found only as examples and illustrations in the texts of Goodman (197A, pp. 85-88, 365-375) and Sage (1977, pp. 212-213). No influenza epidemic models were found in any of the systems and simulations literature. Bailey (1975) presents an extensive literature review of mathe- matical models of epidemics from the beginnings of research in that area up to and including recent develOpments. Elveback, Fox, Ackerman (1975) use a stochastic discrete time model to simulate an influenza epidemic in a community of one thousand peOple and deduce an Optimal strategy for vaccination. . Several Russian researchers (Baroyan, 1971) constructed a computer simulation model to forecast the outbreak of influenza epidemics in forty-three cities in the U.S.S.R. This model is based on a migration model and uses trans- portation data as one of the parameters. Using the parameter estimation technique developed by the above Russian researchers, Spicer (1979) shows that a stochastic discrete time model conforms to the historical data on influenza deaths in the United Kingdom from the years 1952 to 1973. Kilbourne (1975, p. A96) gives a graphic representa- tion of the pattern of pandemics in recent years, depicting the relationship between the immunity level of the popu- lation and the drifts and the shifts of the virus as shown in Figure 1.2. This graph resembles underdamped oscillation and gives some insight into the structure of the influenza epidemic model which is the subject of this research. Three basic models of epidemics described in Bailey 10 mammaammwn N< maca> m4 masunsm .moflsmvcma mucosamcfi mo chopped 6:9 .N.H shaman maths N< waspnam Hmowumpoth Hmofiponuoahs mo :oauosnoppsH mo soapospouch , mama» cw mafia NH Ha DH 0 m h m? m 4 m N H r Annanmv ApMHucv Apmfiucv Aamwnvv Annapvv nowanmv ma \\ :N< =N< =~< .m< m< \ x g \ \ x \ \ \ \ \ \,\ \xl/llxx / \ / \\\||ll//III\\\ // xx chameagm / \\ / \ /I| \\ fikfiq vowumm oasmucmahoucH u. owsmpcmm oHEmncmm maha> m< m> apopwpcm soapmasaon Ho Ho>oH cmos.ll:lil maha> N< m> hponwucm :oapmasaoa mo Hm>oa cams nnnnnnn Disease Incidence 11 (1975) will be discussed in Chapter 2. These models lead to an influenza epidemic model develOped in Chapter 3 which includes the structural state of the virus as one of the state variables. The model will be simulated, and the results compared to historical data provided by the Center for Disease Control in Chapter A. Finally, in Chapter 5, an immunization model will be developed as a subsystem of a control system, and the percentage of the pOpulation required to be immunized to prevent epidemics will be determined. All models will be discussed in the context of systems analysis, focussing on such t0pics as stability, lineariza- tion, phase portrait, and transient response. A qualitative modeling approach, namely causal 100p diagrams, will show the relationship among the various variables. Normally, the causal 100p diagram does not show the mathematical processes relating variables. However, to make the descrip- tion more precise, in the simple models such definitions taken from the state model will be included alongside each variable defined. Block diagrams as well as DYNAMO flow diagrams will be used to show the parameters, the flow of quantities, and the rate and level variables. All simulations will be done by DYNAMO, a well-accepted dynamic systems simulation program. DYNAMO was chosen because of its versatile capa- bilities of plotting and its ease of programming. CHAPTER 2 BASIC EPIDEMIC MODELS Before one attempts to build a large-scale socio- economic or biological model, it is useful to learn the in- depth behavior of the simple model which may be used as'a basis of a more complex one. The dynamic model of the basic system should contain a minimum number of the state variables, or levels, which provide outputs that represent a reasonable similarity to the real world situation. In the case of epidemic models, the system usually consists of three state variables: (1) the size of the susceptible pOpulation (population capable of contracting the disease); (2) the size of the infective pepulation (population. capable of transmitting the disease); and (3) the size of the removed p0pulation (either immune, recovered, or dead, and hence no longer infectious). In this chapter, three basic epidemic models which were presented by Bailey (1975) have been arranged from simplistic to more realistic and will be discussed in the context of systems analysis and the system dynamics approach. Discussion of each model will be organized as follows: First, the state model (a set of differential equations) and the assumptions made in formulating the model will be 12 13 presented. From these assumptions the causal 100p diagram will be constructed. A brief interpretation of feedback loops will be given. Next, the block diagrams, DYNAMO flow diagrams and typical DYNAMO simulation outputs will be shown. Analysis will follow. The variables which are of particular interest to epidemiologists will be examined. These variables are defined as follows: (1) Epidemic curve E(t) = the curve representing disease incidence versus time. (2) The peak of the epidemic Emax = the maximum disease incidence. (3) The peak time tmax = the time at which the disease incidence is maximum. (A) The duration of the epidemic td = the time at which the epidemic curve falls below the initial value in the first two models discussed; in the third model, which is oscillatory, td is taken as the two percent settling time. We will make one assumption common to all three models. This assumption is that the total population, the sum of the susceptible, the infective, and the removed pOpulations, is constant. 2.1 Simple Epidemics Model The first epidemic model, called the simple epidemics model, consists only of the size of the susceptible pOpu- lation and the size of the infective pOpulation as the state 1h variables. Yet, the model provides interesting results about the epidemic curve which are similar to the real world behavior of epidemics in general. Basically, we assume that the disease is transmitted through contact between an infective and a susceptible individual. No allowance for incubation period is made in the model. The rate of infection is then assumed to be proportional to the product of the susceptible and the in- fective pOpulations. Once infected, the infected indivi- dual remains infective until all the susceptible individuals are exhausted. We further assume that the susceptible pOpulation is homogeneous in the sense that every one has an equal chance of contracting the disease from the infec- tive population. We let x1 = size of susceptible pOpulation, x2 = size of infective population, b = infective (or contact) rate, N = total pOpulation. Then, these assumptions may be formulated into the state model: x.- H II -b x x 1 2 (2.1) x2 = b x1 x2 with the initial conditions xl(0) = Cl’ x2(O) = c2, and the boundary conditions x1 + x2 = N, O 5 x1, 0 5 x2. The causal 100p diagram is shown in Figure 2.1. 15 + _ _ Susceptible - -._______Total (x1 I N x2) POpulation POpulation (N) (-) + (£2 = b x1 x2) Rate of Infection 4. (+) + . (x ) Infective 2 POpulation Figure 2.1. Causal 100p diagram of simple epidemics. An epidemic curve is a graphic representation of disease incidence (the rate of disease occurrence per unit of time). In this case we will assume disease incidence; that is, the epidemic curve, to be the rate of infection. E(t) = 12(t) In the causal loop diagram, the rise of the epidemic curve is represented as the positive feedback 100p imbedded in the) negative feedback 100p, which is the fall of the epidemic curve. The block diagram and the flow diagram are shown in Figure 2.2 and Figure 2.3. DYNAMO simulation output for the following test case is given in Figure 2.h: N = 1000 C1 = 950 b = .002 c2 = 50 Goodman (197h. pp.‘85-88), in the discussion of an epidemic growth model, gives DYNAMO simulation output similar to 16 4;“ SPOP Tr r'jdt’ IPOP Jir SPOP = Susceptible POpulation (xi) IPOP = Infective POpulation (x?) b = Infective Rate (Contact Rate) Figure 2.2. Block diagram of simple epidemics. SPOP : '- IPOP \ f \\ I \ / IPOPR / I” I _nf. BETA SPOP = Susceptible Population IPOP = Infective Population IPOPR = Disease Incidence BETA = Infective Rate (Contact Rate) Figure 2.3. Flow diagram of simple epidemics. l7 . 0 O - B A - ' Susceptible xpPopplation . O 0 ' i - ' O - II 0 - | U ' 0 .50“ to :v a... .u. .._s, _ I w :m. .— ...v. . 0 0 0 0 0 o . 0 0 0 0 0 0 0 . 0 0 0 0 0 0 0 . - 0 0 0 0 0 0 run ho...- ' I I I I ........ m00~ !!!!!!! H. | I I I I no. thN. #000. Zedxofiu Hi «3.: Mimosa Time Simple epidemics model. Figure 2.h. 18 Figure 2.L. DYNAMO uses the Euler integration method to approximate numerical solutions of differential equations. The choice of the size of DT, the time increment, is an important step in simulation. The wrong choice of DT may cause instability of the system arising from the simulation, rather than the model. Also the size of DT must be small enough to maintain the bounds of each variable; in this case 0 s xl(t) s.N, O s x2(t) s N. There is no precise formula for determining DT other than the rule of thumb given in the DYNAMO User's Manual (Pugh, 1976, p. LA), which is to take one half the smallest first order delay in the model. In the test case cited above, when DT was set equal to .l and again set equal to .01, the graphs produced by simulation were close to the one plotted from the analytical solution; however, the tabulated results showed inaccuracies between the peak and the end of the epidemic curve. When DT = .001, the tabulated result was almost exactly the analytical solution. Since the simple epidemics model is a conservative system; that is, no quantities are lost through dissipation, and the sum of the state variables is constant, the dimension of the system can be reduced by one by writing one variable as a function of the other. Hence, consider the single state model reduced from the state model (2.1) by eliminating x1. x2 = b(N - x2)x2 19 with the initial condition x2(0) = c2. Separating variables and integrating, we have Jr dx2 _ jf dt b(N - x57x2 - ° The left hand side of the equation can be integrated using the partial fraction decomposition. l l (1. 1 T———)— = - -—+—---—-). Integration yields t = 6N [log x2 - log (N - x2)] + constant. Applying the initial condition and noting N - c2 = c1 we obtain the solution x2(t) = c N . 1 +‘—l e"bNt c2 From the relation xl N - x2 we have xl(t) = c N . 1 + 3% ebNt 1 The epidemic curve E(t) which we defined to be the rate of infection is given by E(t) = b xl(t) x2(t) = bNZ c1 -bNt c2 bNt ° (1 +-E— e ) (1 + E_ e ) 2 1 This symmetric prOperty of the epidemic curve can be shown by verifying the equality E(tmax + t) = E(tmax - t). 20 We now calculate the peak of the epidemic Emax’ the peak time tmax’ and the duration of the epidemic td. The initial value of the epidemic curve is given by 2 E(O) = CEN c2 (1 4- 3'2") (1 + F1.) = b c1 c2 (since C1 + c2 = N). Computation of E(t) is easier if we use the substitution C x = 3g ebNt 1 then bN2 E(t) = 1 (1 +12) (1 + x) Differentiating the above expression with respect to t using the chain rule, we obtain E(t) = W bNZ, (1 + x)3 Setting E(t) = 0 gives the solution x = 1, and hence c log —; °2 t = —— max bN ’ and 2 _ - EN. Emax “ E(tmax) I h ° Using the symmetric prOperty of the epidemic curve, the duration of the epidemic is given by td = 2 tmax c 2 log El 2 - bN 0 21 From the equations for th and Emax’ we note that the ax larger the infective rate, the shorter the peak time and the higher the peak. Therefore, we can say the epidemic is more intense and of shorter duration if the infective (or contact) rate is large. We also note that the peak of the epidemic is independent of the initial conditions c1 and c2, which is somewhat unrealistic. However, the peak time depends on the initial conditions as well as the infective rate and the total pOpulation. If the ratio of the initial susceptible p0pulation to the initial infective population is large, the peak time is longer and hence the duration of the epidemic is longer. The comparison of the analytical calculation and the DYNAMO output for the peak of the epidemic, the peak time, and the duration of the epidemic for the test case is as follows: Analytic DYNAMO Emax 500.00 L99.6L tmax 1.h722 1.5 td 2.9hhh 3.0 2.2 General Epidemics Model The second model takes the more realistic view that after a certain period of time, the infective pOpulation are removed from circulation and become inactive throughout the rest of the epidemic period. In the case of influenza, 22 most of the infected peOple will recover and become immune hence the term immune will be used in place of removed as usually denoted in the literature. Let x3 be the third state variable, namely, the size of the immune population. Then we have the following state model: b x1 x2 - rx2 (2.2) N 00 II I'X2 with the initial conditions x1 (0) = c1, x2 (0) = c2, x3 (0) = c3, and the boundary conditions x1 + x2 + x3 =IV, 0 5 x1, 0 5 x2, 0 5 x3. An additional assumption made in this model is that the infective pOpulation decreases at the rate of r, which is called the removal rate. We note that the infective population increases only when x2 > 0; that is, x2 = (bx1 - r)x2 > O equivalently, x1 > B = p. p is called the relative removal rate. This means that the epidemic occurs only when the initial size of the susceptible papulation is greater than p. This phenomenon is referred to as the threshold phenomenon. Thus, if we know the threshold of a particular epidemic disease to which everyone is susceptible, we only need to vaccinate up to the threshold level to prevent an outbreak of the .23 epidemic. In the simple epidemic model we considered x2 to be the epidemic curve. However, in reality, it is difficult to know the instant when an individual becomes infected. It is more_convenient to equate the rate of immunity *3‘With the disease incidence, because the disease incidence is normally associated with the number of individuals entering a hospital or seeing a doctor in a given time period, and that quantity is more readily observable. Hence, we define the epidemic curve to be E(t) = x3 = rxz. The causal loop diagram is shown in Figure 2.5. ‘ + (x = N - x - x ) Susceptible " - “~Total l 2 3 POpulation - POpulation (N) 4. (bx x ) Rate of Infection 1 2 + (+) (-) + (x2) Infective Pepulation 4. (£3 = rxz) Rate of Immunity 4. (X3) Immune Population Figure 2.5. Causal 100p diagram of general epidemics. 2b In this model the positive feedback of the infectious cycle is now imbedded in two negative feedback loops; hence, the epidemic is slowed down by the removal rate r. The block diagram and the flow diagram are shown in Figure 2.6 and Figure 2.7. The output from the simula- tion for the following test case is given in Figure 2.8(a): N = 1000 b = .002 C1 = 950 r = .8 c2 = 50 p = LOO c3 - O In the case where the initial susceptible pOpulation is below the threshold level p, the initial infective p0pula- tion steadily decreases to zero; hence, the epidemic does not occur. The output from the simulation of this case is shown in Figure 2.3(b). We now consider the reduced state model by eliminating x3 from the state model (2.2). x1 = -b x1 x2 22 = b x1 x2 - rx2 with the initial conditions xl(0) = Cl, x2(0) = c2. We will analyze this system by the isocline and the phase portrait technique (Olinick, 1978). The isocline is x2 - bxl x2 - rx2 s=§;- 'bxixz l is 25 RPOP IPOP SPOP IPOP RPOP b I' Figure 2.6. Susceptible Population Infective Pepulation Immune POpulation Infective Rate (Contact Rate) Removal Rate Block diagram of general epidemics. SPOP ‘37 IPOP ‘37 RPOP \ ///f /* \~ ///fi\\\ \\\ // \\\~ "" IPOPR "" RPOPR ,( I / / / _0'_ _n’__ BETA GAMMA SPOP = Susceptible POpulation IPOP = Infective POpulation IPOPR = Rate of Infection RPOP = Immune POpulation RPOPR = Rate of Immune (Disease Incidence) BETA = Infective Rate (Contact Rate) GAMMA = Removal Rate Figure 2.7. Flow diagram of general epidemics. 26 .on...-nuuun.uuouoaocuou ..,_|....... o .Omm u H0 .004 u a 0 I- o o 1.. .... .1..2..:..:.w$ma . oases“ m n u - . .qawpwmfifinom...:................ ocseeH - "HmuoE oHEouwao Hmpocmu 9 r clan-.o-oo-enso mane .Amvm.m ondmfim a.) .7 c. I I. 000 o r w r ... "or“. .ouuoacuoo-acuo-n-us-nn-o-apuo-nnoouuouuan on... ..... 0° 00 . . . . l . . . . . . .1. . . . . .l. . . .M . ,. .0 . . . la . . . . . . . . . . . .2 . . .o“¢ a. . . . . .J . . . _ . . . . . . . . . . o .t . . . .N l . . . . .OL cones-non...- oo-new-a-ao-cuooo.o.uuoul .cu.-oo.u¢.--“m . . _ . . H u ” genumasaom. _ . . . papapamomum . . . . . . i." all . . . . . .... use.uunto-auo-ce-o-oono'uo-nncoon-.co-o-oc-oo-o- .4 a o .000..-no-oao-ocnucnun-o-o-onuca.0000. . .9 .c ... 2.1. . A. hun— N 315 .I. "-3494: 124-4: I-aoor s=aoas 27 ll CL .Omm u H0 .004 “Hence owEocHQm Hmuocoo .Anvm.m mpzwfim w u. A. u. w. a w .... s. u. w...1......l£.._..°.1....,.4.....3.0.0....001“.3 .0000.0.009.00000000000000000000oooooooa...on000°..ocoosuccess-00.000.000.00... . . .. .x.u;.£...:.;.;.::.:;.:::.Iu::..;.::i.;.;e.:.;.. . . . o . v0. . u - o . vL 2...“..172 .. 2A. .91. .. 1.5.... s A... u o o . ... . o o a o o u. A o a 3.1. g 1.1 . . . .. . . . GOHQMHZQOQ . n o o n o n o .1. 1.» t3!“ 1 1.. (4., . mi. 3.. . :5. .. 5 3s... . fldnflpgwomjm . u o o a o c. t. c 4 t J . .0... c o 0.. I 0 0 O o 0 o 0 0.00.. 0 0 .0 00.000000 0000.000. .00 0 00.0 0 000 00000.. 0 0 0. 000—0 0 0 .0. 000 0.000 0 000 . C . 8.000.000 0.000 0 0 0.. u o u o a a . p o . c c o c . o u o o o J g o o n - m>hao . n - o o . . . . owamnflam. . . . . o o c c o c o o o o . o u c o o c o o a s o o o u o o o o a u . o o p . o c s u o . . . o a o o o o . o o c a o o o c n n o . c o o u q o a c o c o o o o o u o o . u c c n o o c u o . 3 .. o n o. O O 0 0 0 0 0 0 0 .— 0 0 0 0 0 0 0 0 0 0. 0 0 0 0 0 0 0 0 . 0. 0 0 0 0 0 . 0 . 0 0. 0 0 0 0 0 0 0 0 0 0. 0 0 0 0 I . 0 0 .. 0. .v .. :0 0 0 0.. 0 0 0 0 0 0 0 0 0. 0 0 0 0 n ..0 0 0 0. 0 . . 0 0 0 O 0 VA LL .nA . J: . ~ . . 4 n .14.. u. . o . . J o o o o . . . . COfiQmflzgom . . . . ... o o o . o c o a c u o o o . 0::EEH . - o o o I o a c n u u o a o o c o c o o o o o o . 000000000000000000000000.000000000000coo.oooooooooooosocoo-coo...coco0000.00000.000000000000000... a... n I 000°C? 1051' N311“! Isl“! 1:491! Studs LOEF' ' nuv+ 000‘09 1111' 090:0} ' ' ' ° our: 34 1'15 L 28 In this case, by separating variables, we can integrate and solve for x2 in terms of xi: .. ._ .1. fdxz— f(1+pxl)dx1 which gives x2 = -x1 + p log x1 + constant. From the initial conditions, we obtain x2 = el + c2 - x1 + p log :i. (2.3) The trajectory of the test case, cl = 950, c2 = 50, c3 = 0, p = #00, is shown in Figure 2.9. Note that the trajectory moves in the plane bounded by x1 + x2 = N, x1 = O and x2 = 0. When x1 = p, x2 attains maximum since the isocline at that point is horizontal. 1000+ 750 X2 500 0"254) (950,50) + 250 P500 750 1000 x 250 - I I (100,0)\0 ! f 1 Figure 2.9. Phase portrait of general epidemics model. The steady state for x1, denoted by S, can be found by setting x2 = 0 in Equation (2.3) and solving for x1: 29 x1 0 = c1 + c2 - x1 + p log 3; . Equivalently, S.is the root of the equation 1 - - (c + c - x) cl e p 1 2 - x = O. (2.h) The root can be approximated by Newton's method. From the relation x1 + x2 + x3 = N, the steady state for x3, denoted by R, is then ‘ R = N - S. A non-zero steady state for xl means that there is always a portion of the susceptible p0pulation which escapes being infected during the epidemic. We define the immunity level of the population A(t) to be the ratio of the immune pepulation to the total pepu- lation: A(t) =3??- . The steady state of the immunity level, denoted by I, is then -3 I - N . This ratio can also be considered the attack rate of the epidemic referred to in the literature. Hence, if we are given the attack rate, we can calculate the steady state for x3 and consequently the steady state for x1. From this information, together with the initial conditions, the relative removal rate p can be computed from Equation (2.4) and the general epidemics model which satisfies the steady 30 state can be constructed. The explicit time solution for the epidemic curve E(t) cannot be obtained, although the parametric form of the exact solution was worked out by Kendall (1956). An approx- imate form of the epidemic curve, in the case c3 = O, is given by Bailey (1975). To show this, we need to express x1 in terms of x3. This expression can be obtained from Equation (2.3). We rewrite Equation (2.3) as x1 x1 + x2 = C1 + c2 + p log-E— . Since x3 = N - (x1 + x2) and c3 = N - (c1 + c2), we have x1 x3 = c3 - p log 3; . Therefore, we can write X1 in terms of X}: 1 X1 = C1 9 Ip(XB CB). If we assume the initial condition c3 = 0, then 1 -‘- x _ p 3 x1 - cl e . The epidemic curve, then, can be written as E(t) = 13 = rx2 r(N - xl - x3) 1 P x3) r(N - x3 - c1 e We now approximate the exponential term by Taylor expansion to the third term: 1 2 e px3=1-3§2+fz_, p 2 Hence, we obtain the state equation for x3 in terms of the quadratic expression of x3 given by °1 c1 2 x3 = r[N - c1 + (if - l)x3 - 5;: x3 ]. Separating variables and integrating, we have 1‘ cl 2 x3 ] C 1 p 3 2p? The left hand side of the equation can be integrated using a table of integrals.. Because the algebraic manipulation is tedious, we give the final expression without stating the intermediate steps. 2 c x3(e)=L {4-1. a tanh(% art-«p )} c1 p where c 2c c . _ 1 _ 2 1 2 a C ¢ = tranh-l l'-(-p:-L-----1). O. Differentiating x3(t) with respect to time, we obtain the epidemic curve 2 2 E(t) a £3qu sechz (ears -¢ ). (2.5) tmax’ Emax’ td can be calculated from Equation (2.5). The curve y = sechzx is symmetric with respect to the Y 32 axis and attains maximum at the origin, and y = y(0) = 1. max Hence, setting icxrt - ¢ = O in Equation (2.5) and solving for t, we obtain N e- tmax = ' p *3 And, therefore, rel2 2 Emax = E(r'max) = 2cl ' I and td = 2 tmax“ However, the exact value of Emax can be obtained from Equation (2.3). _ _ - IL Emax - rxZmax - r(c1 + c2 p + p log (:1). Another approximate time solution, t in terms of xi, can be derived by the use of the quadratic approximation of Equation (2.3). We take three points on the curve in the phase portrait, (c1, c2), (p, ), (S, O), and determine x2max . 2 the coefficients of the parabola x2 = a2 x1 + a1 x1 + a0 by solving three simultaneous linear equations. We,.then, substitute the parabolic equation into the state equation and obtain _ _ 2 x1 - b x1(a2 x1 + a1 x1 + a0). Separating variables and integrating, we have 1] dx1 j.d con-b- 2 = to xl(a2 x1 + a1 x1 + a0) 33 The left hand side of the equation can be integrated using a table of integrals. Integration yields -2 x a 2a x + a t = - % [—l— log —l- + ————— tanh 1( 2 l 1)] + C, (2.6) 1 2 80 x2 aO ’q "q 2 where q = a1 - haz a0. From the initial condition t = 0, x1 = c1 we determine the constant C. c a 2a c + a C = % [—1— log -l- + .tanh 1( 2 l 1)] 1 ‘CZ aoJ—a ‘ V3 At the peak time, t max’ the infective population x2 18 maximum and the susceptible pOpulation x1 is p, and at the duration time td, by definition the infective population is equal to its initial value. To compute tmax' we substitute x2max and p into x2 and x1, respectively, in Equation (2.6). To compute td, we substitute c2 and h, which is the sus- ceptible pOpulation corresponding to c2, into x2 and x1, respectively. Then, 2 a 2a p + a _ _ 1 1 p 1 -1 2 1 tmax - b [-3- log + tanh ( )] + C 2 x O 2max aO / q \I q 2 1 h -——— log —— + 2aO c2 a 2a h + a l tanh-l( 2 r 1)] + C - - 1 b 30,]?!- «C,— td- where h is the smaller root of the equation 2 - azx + alx + a0 — c2. 3h The comparison of the two analytical approximations and the results obtained from the simulation for the test case shown in Figure 2.8 with respect to the peak of the epidemic, the peak time, and the duration of the epidemic is as follows: Taylor Expansion Quadratic Simulation Emax 167.37 203.20 203.18 tmax 2.1299 2.0872 2.7 td A.2578 6.3177 6.3 The quadradic approximation and the simulation reflect the skewness of the epidemic curve with reSpect to peak time, which is descriptive of many actual epidemic curves. 2.3 Endemic Model In the first two models discussed, the epidemic curve asymptotically approaches zero. The third model, which is more reflective of influenza epidemics, results in an endemic state (the state at which the disease incidence is constant) after an epidemic. Suppose we consider that a certain portion of the immune population is put back into the susceptible pool at a con? stant rate e. The cause of such a transfer is not specified in this model. One example for this transfer to take place is the loss of immunity in infants. A newborn baby has acquired from her mother immunity against many diseases (this type of immunity is referred to as passive immunity) 35 and loses it within six months to a year. Another example which gives the same effect as the transfer is the migration of the susceptible individuals into the system and at the same time the migration of immune individuals out of the system in such a way that the total pOpulation remains con- stant. A specific cause of the transfer; that is, once immune individuals becoming susceptible because of the changes in the virus, will be discussed in the next chapter. Using the above assumption of the constant transfer from the immune population to the susceptible pOpulation, we formulate the following state model: x1 = -b x1 x2 + e 22 = b x1 x2 - rx2 (2.7) 13 = rx2 - e with the initial conditions X1(O) = Cl: X2(O) = C2: X3(O) = C3, and the boundary conditions xl+x2+x3=N, 05x1, 05x2, 05x3. 36 + (x1 = N - x2 - x3) Susceptible “ “Total (N) P0pulation - Population ”\ (-) + (bx1 x2) Rate of Infection 4. (+) + b Infective (x2) POpulation (-> + b Disease (rxz) Incidence Susceptible Rate (e) + /\ Immune POpulation (x3) Figure 2.10. Causal loop diagram of endemic model. The causal loop diagram, Figure 2.10, is the same as that of the general epidemics model with the addition of the susceptible rate e. The increase in the susceptible rate causes the susceptible pOpulation to build up through the decrease of the immune population. When the susceptible population crosses the threshold, we can expect another epidemic to occur. Hence, this causal 100p diagram suggests that oscillation may occur. The block diagram and the flow diagram are shown in Figure 2.11 and Figure 2.12. The DYNAMO output for the test 37 case: N = 1000 b = .001 pl = 800 C1 - 900 r = .8 p2 = 62.5 c2 = 50 e = 50 c3 = 50 is given in Figure 2.13, and for the test case: N = 1000 b = .002 p1 = #00 c1 = 950 r = .8 p2 = 12.5 c2 = 50 e = 10 c3 = O is given in Figure 2.1L. We now consider the nonlinear autonomous state model X = F(X) reduced from the state model (2.7) by eliminating x3: x1 = -b x1 x2 + e 22 = b x1 x2 - rx2 with the initial conditions x1(0) 2 c1, x2(O) = c2. This model has a unique equilibrium point given by '1 I‘D 01"! We will analyze this system using the linearization tech- nique. If the system is structurally stable, then the behavior of the nonlinear system in the vicinity of the equilibrium point is similar to that of the linearized system (Aggarway, 1972) . 38 IPOP " . ~ RPOP 390- + + , TI' 2 dt ~— 6 Sdt SPOP IPOP Susceptible Population Infective POpulation RPOP = Immune Population b Infective Rate (Contact Rate) r = Removal Rate e = Susceptible Rate Figure 2.11. Block diagram of endemic model. 39 / \ / \ / ...\.Q_ EPSLN SPOPR 22> SPOP : — IPOP : RPOP T 7 \ \ I \ \\ // \ “~ IPOPR / ‘~. RPOPR / x l / / / __n;_ __aL_ BETA GAMMA - SPOP = Susceptible Population IPOP = Infective POpulation RPOP = Immune POpulation SPOPR = Rate of Transfer IPOPR = Rate of Infection RPOPR = Disease Incidence BETA = Infective Rate (Contact Rate) GAMMA = Removal Rate EPSLN Susceptible Rate Figure 2.12. Flow diagram of the endemic model. #0 Y. C C 0 C .-.-.. . .. .21. . . . o . a... 0- . - — . . . . - ’0'06 chemo " '09 . A....A. 0 1.1.1.]. 2 t .103 A. 31.1. . ..V.. .vA .. oust-.ocuuuc soc-nunuouun . u - - c u . ucuo-cuvuon- - . a . a - can-.00... - . _ o . n . soot-on... . . . . e o _ sou-cc...- . . - - - . — me .oom u an "Haves oHEoccm .mH.N enemas 9 a: h I 7— r. C. Wu 0 .v o o a 05: .0. “h H. my up. ".5 ”EMU? IIOQIl-IIIIOQ .IIIIIIIII0-t.lt0nl0l .0 I. . Inl-Illil........'...luIIIOIIMHW . — . - - u . L 031.1!5ddHQNHH—Qbiz 00113.01? 0...... 00.13: .. 0: . 0. .\ .mcsssH . .. 0., . . ,;n... . a . . . _ l . I . . .. ... . . . 0. 0 a n o c n o . . -.3 u - g u n - .HL“‘ .6- QIIOI.l000.|1 In IOOOIIAOQI‘IIIII unlctotltl .IIIOOCIIOIOIQI '0‘) . . u c c - .530 . . . . . . . 4. u o - n u — u u u _ u . ~ -. o a - o x n a . . . . . o .5 F3 - a c - o - .o - . Id.’ .0i0l0ll.0000.00....l00l00l 'IO0DII'00....IIOOIIO000. 00......Il. 00laurv . Fuo - - . u u u c l. o g u o u u o . . . . m>h50 . . . . . . oflEmUHQm . . u o - o n u a n. . . . . u - _.3 - g - u o . . . .oc.‘ I..I.l.....lbllll.0.0uu...3:00....I‘00IIIIOIOQOII 1'... .0... 0.00.6.) . . . 1 . .. . .0,” . . - . . . . c - n a I a- u.— o - :owpmasaom. . . . . . wanwpdmomdm . . . . — n - o - - av; SI - - o o . — .Mverv 0.00-0.00...-ensue-o...on.no...0.-uu.-no:00.00.00...sausages-noououaaco er- I ... A ,0 ‘U v: N313"! ”can“ Isdod‘r 3:096. hl - o a u o a c .0.‘ u o a u o o 0.0. a c... . o g - . .00. o o c n - a c 0.0. ... N. O I O O O O O O O . .ob to... 0' .m.ma n mg .000 . o p c . c . - . . . . . OOOOQQOOOOOOOOOOOOOO . . c o J - a 4 q o . o . c u 0 O 0 O O O O O O O O O o O o O o O o O a _ o . c c . ¢ . n c . . . . OOIOOOOOOOOO 0...... - . o c u . c . o . o o u . O D O O 0 v D O O O O O O O O O O I O O .3 0...; . CIA“ IWAWIW' 'H\.ll C u an "Hmuos afiEmvcm .JH.N mgsmwh u m 9 a z o o o a o 3 vi m. n. 1 M. . c o u . o c o a . mSpSQ . . . o a u u c . owswufiam . . . . . . . . .1. 2.....emSmSaQFE: . . , manwuamomsm .. u o o o u . . ... o o u n - o . , . g o o . o . o . . .. . o p c . o c u c a a . . - c o c . . . o .u . . . .a. Q. I. O .0. O... O fiowwmdamom . . . . madeeH : c c . . g o o o . . g o . o . . o . . . o . o o . o o . . a M nugmmmmumugnuu unnummmmmaummm.mmunm.mum “an M M ace 0 O O ~7de n§nv Nana radon 12d“: Sand: #2 Matrix A of the linearized system X = AX is the Jacobian matrix evaluated at Xe: F293. _r A = ‘31— - r 0 3x1 9.9. 0 Xe L r The characteristic equation and the eigenvalues are 2+2; r‘ s + be = O S -be be 2 —- + (—) - Abe S1,2 Since parameters b, r, e are positive, the real part of the eigenvalues are negative in all cases, and, hence, the linearized system is stable in the sense of Liapunov. Oscil- lation occurs when the eigenvalues are complex; that is, when the discriminant of the characteristic equation is negative: (pf-)2 - l+be < O, or equivalently, e z. - + I?) where H = u2 + v2 , and y! = tan-1 % . Therefore, t _ - (‘¢ + V’) max ' v ’ and, Emax = E(t'max) ut . = Ke max Sln (-‘? ). td is the 2 percent settling time in this model, as defined earlier, and is given by L td = IuI . To apply the linearized system for computation of E #5 max’ t and td, one must first transfer the coordinate of xi, max x2 plane so that the origin (the equilibrium point of the linearized system) corresponds to the equilibrium point of the nonlinear system; in this case (%, new coordinate x1- In particular, 2). Hence, for the 1" system 1" e x1 ‘ B ' x2 ‘ x2 ‘ F ° 2 ' _ _ 2 C1 ‘ b ' C2 ‘ c2 r e = 1'(x2 ‘ F) Comparison of the computation from the analytical re- sults of the linearized system and the tabulated results from the simulation in the first test case (Figure 2.13) is as follows: max max Linearized Model Simulation 67.8013 70.3h 8-h5h2 8.00 128 129 CHAPTER 3 INFLUENZA EPIDEMIC MODEL The underdamped oscillation exhibited by the endemic model in Chapter 2 can be used to illustrate Kilbourne's concept of the cyclical nature of influenza epidemics (Figure 1.2) when that model includes the agent; that is, the virus, and its interactions with the population dynamics. Using the endemic model as a basis, we will now develop an influenza epidemic model and analyze it as in Chapter 2. 3.1 Formulation of the State Model In the endemic model discussed in Chapter 2, the oscilla- tion and the endemic state result from the fact that some individuals are transferred from the immune population into the pool of the susceptible population. The influenza epidemic model described in this chapter specifies the cause of that transfer to be the structural change of the virus. Indivi- duals who have been infected by the virus have gained immunity against that specific strain. As the virus changes due to mutation (drift), a certain percentage of the immune pOpula- tion becomes susceptible to that changed virus. To construct a state equation for the virus strain, we consider the relationship between the immunity level of the A6 L7 pOpulation and the virus strain. Kilbourne (1979, p. L96) and more recently Pereira (1979, p. 10) state that antigenic drift is the result of the selection process in which the virus mutates in order to overcome a high immunity level of the population built up through a series of epidemics. Pereira further reports that the occurrence of antigenic drift has been demonstrated in the laboratory by growing the virus in the presence of a specific antibody. In addition to the three state variables in the endemic model we now introduce the fourth state variable and let xh = the virus strain. The virus strain x," is defined to be a qualitative measurement of the structural state of the virus. A numerical value of xh may be a measurement of a chemical composition of the virus structure which indicates a degree of mutation; that is, the higher the number, the greater the degree of mutation resulting in a higher susceptibility level, or equivalently, a lower immunity level of the pOpulation. The above description of drift suggests that the immunity level is the constraint (or the boundary) toward which the virus structure changes to overcome that constraint. We therefore assume that the numerical value of the structural state of the virus xh changes in such a way as to approach asymptotically to the immunity level. Such a process may be described by a state equation which produces a logistic curve. In Chapter 2 we defined the immunity level of the pOpulation: "3 A(t) = T . A8 If the immunity level is constant, A(t) = A, then the state equation in = d(A - Xh)xh (3.1) with the initial condition xh(0) = Ch gives a solution A x = A A - c _ 1 + h e Adt C1. The graph of the solution is a logistic curve where XL is asymptotically approaching A as shown in Figure 3.1. Based on equation (3.1), we construct the state equation for the virus strain to be in d(A(t) - xh)xh x d(j? - xh)xh, where the constant d is called the drift rate. ct. t Figure 3.1. Logistic curve of the virus strain. We also assume that the rate of the transfer of A9 individuals from the immune to the susceptible pOpulation is prOportional to the product of the immune pOpulation and the virus strain. That is, infection resulting from contact between an individual in the immune pOpulation and the new strain of virus indicates susceptibility to that virus as well as the occurrence of drift. Hence, the state equation for the immune pOpulation is given by x3 = rx2 - ex3 xh. The constant e is called the susceptible rate. The transfer of a certain portion of the immune pOpula- tion to the susceptible pOpulation now takes place with the introduction of the drift rate d. Individuals who are now susceptible to the changed virus are still immune to the original strain. Therefore, the susceptible p0pulation is now a mixture of those who are susceptible to a new strain of the virus as well as those who are susceptible to both new and old strains. Thus, we can speak of the susceptible population collectively as having different levels of susceptibility at different times. » Let S(t) be the susceptible level of the population and define S(t) = 1 — A(t). We then assume that the infective rate x2 is proportional to the product of the susceptible and the infective populations and the susceptible level. In effect we have replaced the constant infective rate b in the endemic model with the time varying infective rate bS(t). Then the state equations for 50 the susceptible pOpulation and the infective population are given by x1 = -bS(t) x1 x2 + ex3 xh x = -b(l - 3%) x1 x2 + ex3 xh - _ 1‘}. - x2 — b(l N) x1 x2 rx2. Now we write the state model of the influenza epidemic system: x x1 - -b(l - 3%) x1 x2 + ex3 x1+ x x2 = b(l - 3%) x1 x2 - rx2 x3 = rx2 - ex3 XL x in = d(1% - xh)xh (3.2) with the initial conditions X1(O) = cl’ x2(0) = C2: x3(0) = 039 xh(0) = CL, and the boundary conditions x1 + x2 + x3 = N, O 5 x1, 0 5 x2, 0 5 x3. 3.2 Causal Loop Analysis The detailed causal loop diagram is shown in Figure 3.2. To show more clearly the relationship between the virus strain and the pOpulation, a simplified causal loop diagram, Figure 3.3, is obtained by leaving only the level variables and the two auxiliary variables, the immunity level A(t) and the 51 susceptibility level S(t). This diagram shows that the virus strain is the direct cause of the changes in the immune population. The appearance of a drift reduces the size of the immune pOpulation which in turn increases the size of the susceptible pOpulation. When the size of the susceptible pOpulation exceeds the threshold, an epidemic occurs. The size of the epidemic is now dependent not only on the con- tact rate b, but also on the susceptibility level S(t) which is a function of the virus strain as seen in the causal 100p diagram. The increase in the size of the immune pOpulation, as a consequence of the epidemic, causes the higher immunity level of the population. In order for the virus to survive, the structure of the virus changes in response to the high level of immunity in the pOpulation, which causes another drift, and the epidemic cycle begins again. 3.3 System Analysis The influenza epidemic system described by the state model (3.2) is classified as a nonlinear autonomous system i = F(X) with the dimension A. This system can be analyzed by linearization about its equilibrium point provided the linearized system is stable in the sense of Liapunov. With the prOper choice of the drift rate d, the linearized system of the influenza epidemic model satisfies the stability condition; therefore, the linearized technique will be used for this analysis. 52 .Hmooe owsonwdo mucosamcw Ho Empwmwv QOOH Hmmsmo .N.m ouzwwm + Lommcmha mo 09mm Aaxmxmv + a :owumasdom m + ocsssH A xv 4 a 4 mwcmno A xA xlfipv mo mama . + I AmocmcwocH ommmmwav Alv + ocsssH mo mama Amxhv + Asxv , camtsm mata> A+V cowpmasaom A+v + m>wpommcH Amxv E . z Ho>og N H ANluAuvoq :owpmazdom + H Axyvauauxpvmv spaaanapawomsm . ospauamomsm A xv + sz :ofipmHSQom Hmpoa 53 _ Susceptible Susceptibility Population Level (-) (-) + Infective + Virus Population Strain + (-) + (-) Immune - + Immunity POpulation \ 4/ Level Figure 3.3. Simplified causal loop diagram of the influenza epidemic model. 51+ .Hmcos owsmvwam mucosamca mo Emhmmfic xooam .4.m ohzwfim chLum msuw> u H> Hm>ma spacsseH u a>azH comm phage u n Hm>oq manwpdmomsm u q>qm mumm manwpdmomsm u m cowumasaom mGSEEH u momm mumm Hm>osmm u L cofipmasaom m>wuommcH u momH Ampwm pompcoov mama o>wpomMCH u n :owpmasdom manflpdmomsm u moan ‘I moeomm mDmH> _ mOBOmm MZDEEH MOBomm OHEWDHQM 55 ’-‘ v EPSLN ‘ ,I’ \ O ~—-————-_ A v K I; / \ IPOPR / RPOPR ,’ // I, / / / // // / / / / / _.C{_ // / ’ BETA GAMMA // //’ / z / / \ / I N \\_d SPOP = Susceptible POpulation SLVL Susceptible Level of POpulation IPOP = Infective POpulation IMLVL = Immunity Level of RPOP = Immune POpulation Population SPOPR = Rate of Transfer BETA = Infective Rate (Contact Rate) IPOPR = Rate of Infection GAMMA = Removal Rate RPOPR = Disease Incidence EPSLN = Susceptible Rate VIR = Rate of Virus Strain DELTA = Drift Rate Figure 3.5. Flow diagram of influenza epidemic model. 56 3.3.1 Equilibrium Points Analysis The equilibrium points for the system are solved by setting each state equation in the state model (3.2) equal to zero; that is, F(X) = 0. They are given by F r 7 r b l-c N-03 % ch 0 X: = . X: = (3.3) No c k c .J L O , where c is some constant. Since we are concerned with a system which produces an endemic state; that is, the steady state of the infective pOpulation is non zero, this analysis considers the linearization about the first equilibrium point Xi. Moreover, the second equilibrium point x: can be shown to be structurally unstable for all values of c. The constant c must be such that the boundary conditions are met; that is: xe1 + xe2 + xe3 = N and O s < N, 0 s s N. xe1 - s N, O s Xe2 xe3 Hence we have the relations: m+chz+Nc=N and O < c < 1. ”Sim As before, let p1 = % be the relative removal rate and p2 = be the relative immune rate. Then we have P1 2 FE-i-pZNC +NC=N. The above expression can be written as - 1 N - p 3 _ E£____ 2 _ 2L _____l = o c p2 c p2 c + p2N . (3.4) Therefore, c is the real root of the cubic equation which lies between 0 and 1. The existence of such c is shown by the intermediate value theorem of calculus. We let p - 1 N - p f(x)=x3--£-—-x2--2-x+-———l. 92 P2 p2N Then, N ‘ P1 . f(0) = _E;N__ > 0 (Since N > pl), and 2N + p l f(l)=--—--—- <0. P2N Since the sign of f(x) at x = 0 and x = l are opposite, and f(x) is a continuous function, the graph of f(x) must cross the X axis between x = 0 and x = 1. From the above analysis, given p1, p2, there exists an equilibrium point given by (3.3) where the constant c is the real root of the cubic equation (3.h). In most cases this equilibrium point is unique because a chance of having three real roots being in the interval x = 0 to x = 1 is relatively small. Since the virus strain is bounded by O and l, the struc- tural change of the virus (drift) resulting from one epidemic 58 to the next may be measured in terms of percentage differences. Suppose that the state of the virus strain causing the ini- tial epidemic were .1 and at the beginning of the next epidemic had changed to .3. Then we would say that the drift was 20 percent, the quantitative measurement of the struc- tural change of the virus. (A more precise definition of drift will be given in Chapter A.) As a series of epidemics occurs, the virus strain changes in such a way as to approach the immunity level of the pOpulation A(t) provided the system is stable. Let I be the steady state of the immunity level of the pOpulation as defined before. Then I = lim A(t) = lim 51%‘1 t9” t-yoo = Nfic- = C. Hence, the steady state of the virus strain is equal to the steady state of the immunity level, and the reciprocal of the drift rate d is analogous to a delay constant in an exponential growth curve. 3.3.2 .Linearized Model The relations x1 + x2 + x3 = N enables one to reduce the dimension of the state model (3.2) by one because any one of the variables can be expressed as a function of the other. We let x3 = f(xl x2). Then, we can eliminate the state equation for x3 from the state model (3.2) and obtain the following reduced state model X = F(X): 59 X x X in = d (-Nl -- x1.) *1. (3.5) with the initial conditions xl(0) = cl, x2(0) = c2, xh(0) = Ch' Note x3 is now a function of x1 and x2. The Jacobian matrix Eg— for the state model (3.5) is given by axi r-b[(1-%l>x2+x1Nx2]—exh -b[(l-%}')xl+:]§:gJ-exh ext: b[(l—;%)x2+f%;g] b[(l-§%)xl+f%;g]-r 0 k “K xh "% x4 ’2dxu) The A matrix of the linearized model i = AX is obtained by evaluating the Jacobian matrix at the equilibrium point IP13 l-c p2Nc2 and has the form FU - ec V - ec Nec A = -U -(V+r) o (3.6) where 60 P1P2 U = -bC2 [(1 -C) p2N + 7:;- J p p 02 v = -b (p1 + J-i-EF—L (3.7) We now derive a formula which directly gives the coef- ficients of the characteristic equation of the A matrix without having to compute det(sI - A) = 0 each time as is the normal procedure for obtaining the characteristic equation. The computational advantage of using such a formula becomes quite clear when one is comparing the behavior of the linearized system for a number of different sets of para- meters. The inverse of the matrix A (3.6) is calculated as A": Ad j A det A F 2dc(V+r) dc(2V-3ec) Nec(V+r) 3 = de% A ~2ch -dc(2U-3ec) -Nch d d Afich—(V+r)] fic(U-V) U(V-ec)-(V+r)(U-ecy where det A = -dc[U(2V - 3ec) - (V + r)(2U - 3ec)]. In the case of a three by three matrix A, it can be shown that the characteristic equation is given by (Wilkinson, 1965) s3 - (tr A)s2 + (tr Adj A)s - det A = o (3.8) where tr A = trace of matrix A tr Adj A = trace of the adjoint matrix A. From the equation (3.8) we obtain the coefficients of the characteristic equation of the matrix A (3.6), 6l 83 + a2 52 + a1 5 + a0 = 0, given by a2 = 2dc - (U - ec) + (V + r) a1 = dc [2(V + r) - (2U - 3ec)] + U(V - ec) - (V+r)(U-ec) a0 = dc [U(2V - 3ec) - (v + r)(2U - 3ec)]. (3.9) Note that if we are given p1, p2, then the coefficients of the characteristic equation are expressed as functions of d: s3+f(d)32+f(d)s+f(d)-0 1 2 3 ‘ ° In a linear system, the location of the roots of the characteristic equation determines the system response. Kilbourne's observation of the influenza epidemic cycles (see Figure 1.2.) suggests that the system response for the type A virus epidemic system is stable and underdamped oscil- lation. The drift rate constant d can be chosen in such a way as to determine the desired response described above. We will consider first the choice of d required for stability, then for oscillation, and finally, for underdamped oscillation. A linear system is stable in the sense of Liapunov if all the roots of the characteristic equation lie in the left half plane. The Routh stability criterion gives the condition for which the linear system is stable (Saucedo, 1968). The Routh array for the characteristic equation is formed as follows: 62 3 s 1 a1 2 5 a2 a0 1 30 1 a2 0 3 a0 0 The Routh criterion states that the number of roots with positive real parts is equal to the number of sign changes in the first column of the Routh array. In this case each element of the first column must be positive for all three roots to lie in the left half plane. Also, from the theory of equations, all coefficients must be positive. Hence, the system is stable if the coefficients of the charac- teristic equation satisfy the inequalities a a. > 0, i = 0,1,2; and a1> '55 In an oscillatory system one of the roots of the characteristic equation has to be real and the other two complex. From the theory of equations, the necessary and sufficient condition for the cubic equation 3 2 x + 32 x + a1 x + 30 = O to have one real root and two complex roots is that the discriminant 63 _ _ 3 2 2 _ 3 _ 2 A - 18 azalaO A a2 a0 + a2 al A al 27 a0 is negative. Here we have the characteristic equation whose coefficients are functions of the drift rate d. To find the range of d which gives a set of coefficients satisfying the above condition is computationally impractical. Therefore, we now seek a simpler sufficient condition for the system to be oscillatory. Consider the characteristic polynomial f(s) = 53 + azs2 + als + a0. The graph of the cubic polynomial crosses the real axis only once if there is no relative minimum or relative maximum, or if the critical point is the point of inflection. The critical points of f(s) are solved by taking the derivative and setting f’(s) = 0. Hence, the critical points are the roots of the quadratic equation 352 + 2a2s + a1 = 0 and are given by 2 _ -2a2.: V[ha2 - 12al (3 10) 51,2 "’ 6 ° ° If the discriminant is non positive; i.e., La 2 - 12a < 0 2 1 “ or equivalently a1 2 -§- , (3-11) then the roots are either complex or multiple roots. If the roots are complex, the cubic polynomial has no critical point; and if the roots are multiple, the critical point is the 6h point of inflection. Therefore, a sufficient condition for the system to be oscillatory is that the coefficients of the characteristic polynomial satisfy inequality (3.11). The system is both stable and oscillatory if the coef- ficients of the characteristic equation satisfy the inequal- ity 2 a1 a» max 2% , 'i§- . (3.12) Although such a choice of coefficients guarantees the oscillatory, stable system, the real root of the character— istic equation may lie to the right of the real part of the complex roots. In order to have underdamped oscillation, we require that the real root be on the left side of the real part of the complex roots, preferably as far to the left as possible. (The influence of the real root is negligible if it lies six times as far to the left as the real part of the complex roots.) Figure 3.1 shows the location of the roots and the system response for the third order linear system. We now derive a sufficient condition for the system to be underdamped. This condition will be such that the widest range of d can be easily obtained. The coefficients of the characteristic equation 3 2 .. s + azs + als + a0 — O with the roots $1, 52 3 = u i jv have the following relation- 9 ships: 65 c(t) Im 31 Re 7‘ I >< t =7 c(t) Im X ;< Re X’ t r: c(t) I . m .X x Re ‘X t : Figure 3.6. System response of a third order linear system. 66 - -(sl + $2 + 83) SD N I a1 - 8182 + $183 + 8283 m 0 I Assume all roots lie in the left half plane; i.e., 51 < O, u <:O. Then all coefficients are positive. Since a2 is equal to the negative of the sum of the three roots, the magnitude of the real root is greater than one third of a2 whenever the real root lies to the left of the real part of the complex roots. That is, $1.< - :5 if and only if sl‘< u. Note that the real part of the critical points given in Equation (3.10) is - é; . Hence, if the coefficients of the characteristic polynomial f(s) satisfy Inequality (3.11) and if Sl‘< u, then we must have f(- 2%) > 0. (See Figure 3.7 (a).) This inequality can be simplified to yield the relation a 2 2 O \OIN a < 3 -— + 1 a2 Combining the above inequality with inequality (3.12), we obtain a sufficient condition for the system to be stable, oscillatory and underdamped. This condition is given by 2 a a a max{ -9- , £- < a1 < 3 ~- + 5' 322. (3-13) 67 f(s f(s) f(s)=sB+a $2+a s+a 2 l O ":10 a O 2 2 3 -— +‘— a - l I I l | I s 51 :2. 3 (a) (b) Figure 3. 7. Curves of cubic polynomials having one real root and two complex roots, (a) with no relative maximum or relative minimum, (b) with a relative max— imum and a relative minimum. So far we have dealt only with the case where the characteristic polynomial has no relative maximum or relative minimum. We now relax this condition and allow it to have a relative maximum and a relative minimum (Figure 3.7 (b)). We note from Equation (3.10) that the point of inflection sp a is equal to - —g Hence, if 31 < u, then the relative 3 0 minimum must be positive. To obtain a simple relationship among the coefficients of the characteristic polynomial which satisfy the above condition, we draw two lines, one tangent at the point 68 a - a (- 7&3 f(- ?§)) and the other tangent at the point (0, a0). The point of intersection of the two tangent lines is given a a by (- ig’ 9 3% - a1). Thus, a sufficient condition for the relative minimum to be positive is given by a0 9-55-a1>0 or equivalently, Since the existence of the relative maximum and the relative 2 minimum implies that al < a2 , the sufficient condition for the system to be stable, oscillatory and _underdamped is given by 2 a a a O . 2 O 33 o. The graph of the above parabola does not cross the real axis; i.e., the roots of the quadratic equation are complex 70 (d1,2 = -.1775 i J .LBL). Hence, the solution is d >00 (3’16) - a 0 2 2 . a1 <-3 -— +.§ a2 yields a 2 .0525 d2 - .2388 d + .0252 < 0. The parabola crosses at d1 = .1080 and d2 h.hh20. Hence, the solution is , .1081 < d <1h.h420. (3.17) a0 a1 < 9-—— yields a2 .1505 d2 - .8106 d + .01 < 0. 5-3357. -Hence, The parabola crosses at d1 = .0h98 and d2 the solution is .0198 < d < 5.3357. (3.18) Taking the union of the solutions (3.17) and (3.18) and then taking the intersection of the union and the solution (3.16) we obtain the range of d: 1.081 < d < 5.3357. CHAPTER h SYSTEM SIMULATION AND MODEL VALIDATION A number of simulation outputs for the influenza epidemic model developed in Chapter 3 will be presented here, and the relationship between each parameter and its corresponding system behavior will be analyzed. In particular, the peaks and the thresholds of the epidemics, the intervals between epidemics, and the drift of the virus will be discussed. Then the system output will be compared with historical data on influenza epidemics. We have computed the equilibrium point and the range of drift rate d for underdamped oscillation for the test case in Chapter 3. The range of d is between .1080 and 5.3357. Figures 1.1, h-Z, 4.3, and h.h show the DYNAMO out- puts for d = 1, d = .5, d = .3, d = . 1 respectively, with the initial conditions: 01 = 900 (susceptible p0pu1ation) c2 = 50 (infective p0pu1ation) c3 = 50 (immune p0pu1ation) cl+ = .001 (virus strain) 71 72 .H u c oumh champ one saw: Hooos oHSocHQo wucoaamcH .H.¢ M 5 9 1 9 s o r... 0 O O 0 r .....: . m. .. m. m. m. . o n o . - u . u c . . c - o c . . . p . . u o c c c c o o o n o . o o o o o c c o O. a)” O I 0 O O O 0).- 05010 I O o O O O 0.. O O O I O O O). 0.. O O O o O O C 0.0 I O I O O o O O O o 0 O C O o o u c . . . . o>mso oaamnaam . o c o g g c ”I .p I I . r .r . I .f o o p o o . . . Hmsmg snaggesH . . . . . . . o u a . o o o o o c o u . n o o u c . .. . . . . o c . . a o . . . . . . .00 0.0...O‘DOODOOOOOOCOOOOOO.0.0.00.0.00.0..OOOOIOOIOOOOO.v. . o o c o c . o . c-Hmhpm msufl> . o o o o a c u c o . c c . o . u u . mzucoe :a mafia 0.0.... .009... C. 35 C .U o O a» a 0.... ....... OI . c . q . o . c o 0.0.0.0.003 . o . u o g g g n c . a \ o 0 so. 0 ‘00 o a - o . u o u o c n no 90.000000. . c . . c c o o . . o o n - whamHm .0.GI A-ov.ufi.w O o o O .v.-.vf~ .4. 93 u o UAV nit I“-U Go; . 6.“ AU YD AU 0 o 0 739° 09.0 0.“ HIV .530 0540 O . C O - E - O . O - . 9 O 90 0 a co. ”Pg-V 00......00......OOOOOOOOOOOOOOOOOC.U-b at A c-Ido AV» ,4.” hum $19M: 73 .mo nu we 6 8 I. 3 o s... A... o o o o 3 D. .U a.» 00' ...... 0.00.0000000-00000000000n o O . . . . . . . . . . o o o . . . . . . .o . o o a . . u c O... C0.0.00.0IOOQOI‘IOOOOOOIII. a. . L up). I ..> . o . . o . . . . a. . . o>hso owsouwam. . . o .11 1 -1. 11111111. I1 I... . o . . . .00....e0.0...OOOOOOOIOIOOIOOOOOOO.I... . . . u - . . o . . . . o . . c . . . . . . . . . . - . O. 00. OOOOODOOOOOOOOCOOO0.0.0....O... . o n n . o . . (II. IIIII .11 A JIIII. III. - 1 . . o . . . a . o . o u . . . . 9 p! O 9 o o 0 0 . ..... o. co... . - . . o . . . . . o o . . 00.... .0... o . . o . . . . c - . . ooooooooooooo . o u . o . . . O. .0000... u .r\ c o . . o . ooooooonooooo o Sufi: Hocos oHEmcHQo mucosamcH mnucoe :H oEHB P. .0 n0 .~.4 «human 0°52 0°02 dflmhpm . u “.8... n H H H Howe; . . suaczeeH 0.0.. II An. OOOOOOOOOOOOOOOOOOOOO.0000... ..... h\.uoyu c wrung I .5 O o I .3,—,3 . o o 1' v '- Cr any)» km '8 M‘MI NSI'JNI 7h 0 fl\ 0 II M“ .3 8 A... AV .0 . . . . . . . . . . . . oypso oasflvwam O. O O O O O Q Q .- 0 .0. C .0 .O C O C o O . . o . . . . . . . . . . . . - u c . . . . . . I. C I 0 O ICOfiOmOLOOQm. O . 0 O O .1. I . C O O O O O . .mSLHS .nn. . .. . . "A. . . .4 . . . . . . . .. :C'Cg ‘0’09 msuco " E LC'OS . O . - - - - - . ..- - - - - - - . . -:- - - - - - - . - ‘ I :« mafia. p :uflz Hence owsmcfino mucosaucH .m.4 opswwm C. 3 T. u a (v 0 O I v #9 o I O 006...... dm>oq hpacsesH .HAH . .Lccark 0 AU ufs . 03L. . . o c 7.. g .L a v o. O O ”1.5.. O O O o O rv~,v...w nL-.. AV" “=13!!! AdA Vs‘MM 75 .H. u c npflz H0605 owEmcwdo mucozamcH 3.4 magma '0 .U '0 a L P, g f C. 2 cs 5» C c 5 c o as O 5 0 .5C.... 0 o o o o o o o o o o... as r.» nu A]. IV 0 av a 0 wk...- ) I d1“) 55 I).4L)) L I I11 I 4 III]! I 1 4 4L 4 )2 I .IOOQOucI 02.3.51 III .v~... - . o c . a g o u o. . . . . . . . . . o . n c o g o o . . . o u g o o c u o c c o . u . . o o c u o c 7. p) o o . n c n u . a c .3 J 0.. o . o . o u . o . c u at...» 0.. coo. cocooooooooooooooOoooo..pooooonooo. u...oooo..oo.ono.po.ocoooooooonoooooocoooooooooo o .v a . ~I O o O K. O n . . O -.\.u\ d0 . . _ . :.mhpm mama. . . . . . . .V c . n c o - c o o . . . ... . . . . . o>h=o case an . o u n o o a o o 0 ”Com . l v. . . c . o u . o o a .2. . L ..- . o o o a - u o o o .00. O r5”..s coco-coococoooooooooooooooooo onnoose-00.000.000.00.coo-coco...ooooooo.00000000000000000000000 ...4£.t KuUH‘ . . o c. . o . . o c .v. o o o c o o a u - a o o o . o . . o o c . o . . u o - o . c c o . . . . . .5 . . . . . . I . a 2. . . o . o no. u . . g . up, rs 0.9 o u . o u o o o c o a .T...)... OOOOOOOOOOIOOOO COCIOOIODOODOOOOOOOOOO OOOOOIOOCOOOO~I. OIIOOO0.00I000......OOOCOOOOOOOOOOOOO COCO-Lune . . o . . c . o o o .33: v o c o u . o g .. on. o n a u . Ur. ' s c o o c o u c u . . . . . . Ho>mg apHGSEEH. . . . .1 . . . . . . . . . . I70 . u . c . u o o u . ..v..- .\D p c o o o o c . c . .9.- pv.vu.$ nono.0.0000000000000000000000coco...ocoo-0.00.0.0no...a...00.00000...coo.0000.00000000000000.9000no.0...rnv ‘ A!“ A31}. Vafl‘hfl' N113 NI 76 h.1 Analysis of Drift In the case of oscillation (Figures L.1, h.2, and h.3) there are three distinct epidemics before the system reaches its steady state. Also, the virus strain peaks twice, each peak occuring between the epidemics. We note that the time intervals between epidemics (measured from one peak time to the next peak time) is longer as the reciprocal of d is larger. Hence, we can consider T =.% as a delay constant similar to the time constant of the exponential growth as mentioned in the discussion of the equilibrium point analysis in Chapter 3. Using the basic unit of time as one month, we find that the time intervals between the epidemics for T = 1 month are 2A and 16 months; and for T = 2 months, 38 and 17 months; and for T = 3 1/3 months, 56 and 19 months. Likewise, the virus strain reaches its peak in longer time when T is larger. In this model, drift is not only the measurement of the difference in the virus strain, but also a function of the time interval between the two virus strains. As the virus strain grows, a continuous transfer of indivi- duals from the immune p0pu1ation to the susceptible popula- tion takes place. Hence, even though the difference of the virus strain from time tl to time t2 may be small, if the interval t2 - t1 is large, the total transfer of individuals from the immune to the susceptible p0pu1ation is large, thus causing a severe epidemic. As the size of the immune p0pu1ation becomes smaller through this transfer, the immunity level of the p0pu1ation decreases which causes the growth of 77 the virus strain to slow down or even decrease. The occurrence of another epidemic builds up the immune p0pu1a- tion which in turn causes the virus strain to grow again. From this causal description of the virus and the epi- demics, we define the drift D to be D = 21-1— (vmax — v ). (1.1) min where Vma is a peak of the virus strain and Vmi is the X n minimum value of the virus strain preceding V (I.e., max' the difference of the virus strain is measured from one valley to the next peak. See Figure 4.5.) The time of the occur- rence of drift TD is defined to be the time at which the virus strain attains its peak. Virus q ? I Epidemic Curve 0-3 ___._......___........._.___._ min D Figure h.5. Drift and epidemic curve. 78 guano u a :Hmppm wzpfl> u > :fimhpw 02.6, mo 0nd; Essfidwso 0?“:0 owE0uwao u m msgw}\owsocfim0 mo x000 x Anacosv 05w» n 0 NmN. 40* 55 00m. 40* mm own. 00* N4 544. num.x 0m on Hm. omm.x mm 0m Noa. mmm.x as am 04H. mHHx mm NBH. 50H: H4 JON. ow: hm qu.o an em moa.o am on Hoa.o 4m mm mo. noa.x o 08 004. mam.x 0 mm mmm. 0mm.x a pa HOG. mmax N HOG. Omar m Oma. Omax m H006 08 0 H00; 04 0 H00. 0 04 o a > m p a > m a a > m p m. u m. H u .mxm0n owEocHQ0 0:0 mumwpn mspfi> 2003909 afismcowpwa0m .H.4 0Hnms 79 Table 4.1 shows the numerical values of the epidemic peaks and the drifts for the oscillatory cases (Figures h.1, 4.2, and 4.3). The relationship between the drift calculated using the above definition and the following peak of the epidemic curve confirms the reasonableness of the definition; that is, the higher the drift, the more severe the subsequent epidemic. A.2 Threshold of Epidemics In the general epidemics model the threshold level of the epidemic is equal to the relative removal rate p. When- ever the size of the susceptible population exceeds that level, the size of the infective p0pu1ation increases, thus causing the epidemic. The difference between the initial susceptible p0pu1ation and the threshold is the measurement for the severity of the epidemic. If the initial susceptible population is very large compared to the threshold, the peak of the epidemic is high. On the other hand, if the initial susceptible population is smaller than the threshold, the infective p0pu1ation decreases asymptotically to zero as shown in Figure 2.8 (b). In the case of the influenza epidemic model, the threshold is a time varying one. Setting x2 = 0 in the state model (3.2) we obtain the threshold H(t) as 80 “(‘5’ = 8m N _ x - (4.2) The threshold is inversely proportional to the susceptible level. Even if the composition of the initial p0pu1ation is such that the susceptible p0pu1ation is below the threshold, it will not necessarily prevent epidemics because the sus- ceptible p0pu1ation may cross the threshold and stay above it. If we consider the infective p0pu1ation negligible, the size of the susceptible population required to be equal to the threshold level can be calculated from Equation (0.2) and is given by x1 = J Npl . Figure 4.6 shows the simulation output for the test case with the initial conditions: cl = 600 (susceptible p0pu1ation) c2 = 50 (infective population) c3 = 350 (immune p0pu1ation) CL = .001 (virus strain). The initial threshold is H(0) = 615. An Optimal immunization strategy to minimize the peak of the epidemic will be treated in the next chapter. 81 .0fisonwa0 HmflprH mo 0oc0upsooo :oz a. 6 3 .L a o o o W... 0 0° OOOOOWUOOOOOOOOOHO 00.... o o a o c o . o o o a a c o c c . c c . . o c . to...00.0.000000000.0.0....0.00.00.00.000000 . o n u . . . o o . o u c o I o I . q a I I IlIbaIII I no". .0 . . o o. o . a o a . I -I--.- --I--I-HIII... .- . .0>nzo owaocwam . . c o a o o c o o o o . . o o o 0.0 00.0.0 0..000.00.00.00000000.0.0.0000... . o o c I III II IUHI IIUIIHII...I. . . . Cfimhbm mfihw>. - . o o . u o c a c o . eo.cooooooocooooooooooooooooooooooooooooooooooooooooococo-coo.ooooooooooooooooooooo 9 o O a 0000...... o o n o o u . 90.0.00... 0 o o a o o 1 - - - . - - u‘ 0 ‘ -‘C - C - C - - . . C . nausea :H 0EH9 g o o 0 000000.... 0 c u c a c 0 00.0000... . o . I. 00.4 0usmfim 3'32 $2.3 . hpacseeH IO... 9.. o 0 av 000. o \ a c o u c o 0000000 no . o g o . g o 000...... o o c o . .I . c . 00000000.. . . a . c . . 090000000. a"! 90 o o o c o a o o o 0.620 0 can cavco Q o o o O o o c o O o O O C . tau-v O .Uafio . Ooh...- - o a o o O coc 3.0..» LVN Axum Va1flflfl1 usxant 82 h.3 Causal Loop-Diagram of Parameters Simulations are performed more efficiently in terms of the number of runs if the causal effects of the parameters on the system outputs are known. A desired response may be obtained fairly quickly by adjusting one appropriate para- meter rather than making adjustments on all of them. Such effects are observed from the results of the simulation of a test case obtained by varying one parameter and fixing all others. The causal effects on the peaks and the intervals of epidemics by each parameter is demonstrated in the causal loop diagram Figure h.6. This causal relationship was used to generate the system outputs for the model validation. h.h Model Validation The Center for Disease Control (CDC) in Atlanta, Georgia maintains the weekly statistics on the number of deaths resulting from influenza pneumonia. These data are collected from 121 cities with a p0pu1ation in excess of 100,000. This mortality chart is also broken down into nine geographi- cal regions. Figure h.7 shows the mortality data from September 1968 to September 1977. The expected curve is taken from the least square method of fitting the curve to the data for the previous five years excluding epidemic periods (periods in which the number of deaths exceed the threshold). 83 .mh090smpmq mo Empwmflo mooa Hmmsmo .©.¢ 0uswflm A0000 0czeswv 0 A0000 agency 0 1' Amofisoowa0 mo mam>popcwv H + A0000 H0>oE0pv 0 III Ill! A0pmh 0::EEa 0>Huma0hv m." Na Amowsonfimo mo mxm0dv xmsm + A0000 0>fluo0mcwv n.IIIIlIl\‘.Aopmn H0>oE0h 0>H90H0hv n I am Hi 8h PNEUMCNIA-IHFLUENZA - React-tad Dam: in 121 Sdocztd Cities. Unind Sta-m. Septumbu 1563~30mmbn 1977 I. In! in: awe v . an: Dunn on. a... 0" on .kfl: - ... 0 q A "' r . 3. . nut . an. I. ”“0 N I. bad” x -. an: I .1 on. .0”. \\N. . .u. . .1 .80 .N.» Q... I o 1 an. /~ . .. ‘0. '.. ‘~ I .1. ‘.. Ill-1 ‘LLA‘J 'I <- .‘I "a. u...- 53 H u” r m It). __ O 1111'_._L..I ' I O 1 n... A. If f. l.-H.' -a At . .a.- J L‘.. 11114 -JIII‘ 1‘14 .4 WJ_1L ‘ 22:32:22.:u Ger-awe I?! . hda .. nu. . ”no AHO- . «a. (an. I \x ,f". 'V 4! ll. In. at I” . .. MLAV . .w an; if... .dummm . e: “IP‘dO 1 Weekly mortality rate of influenza-pneumonia. 1 From the Annual Summary 1977 of Morbidity and Mortality Weekly Report. Figure 4.7. 85 From numerous simulation results of the model, it be- came apparent that the Oscillation dies down quickly after the third epidemic. Hence, with this model we can generate three distinct epidemics and two drifts. There were six epidemics from the onset of the 1968 Hong Kong influenza (an occurrence of the virus shift) up to 1977 excluding the small epidemic in July 1972. Of these six epidemics, we will attempt to produce the system outputs of the model to fit the data for three epidemics and compare the results with the actual figures. The three epidemics we will consider are: the outbreak of November 1968 - March 1969, that of January 1972 - March 1973. and that of January - March 1975. The difficulty in validating the model by fitting the his- torical curve is that data for the initial values and the parameters are not available. Nevertheless, we will estimate these values in order to generate the historical curve. The model will be modified to accomodate the population growth and the seasonal variation of influenza. We let POP = total p0pu1ation SPOP = susceptible p0pu1ation IPOP = infective population BRTH = birth rate (pe0ple/month) DTH a death rate (peOpIe/month). Assuming the exponential growth of the p0pu1ation, we have 3%- POP = (BRTH - DTH) * P0? and POP = SPOP + IPOP + RPOP. 86 The newborns are put into the susceptible p0pu1ation,and the death rate is equally applied to the three groups of population. Table h.2 gives the annual birth and death rates of the United States from 1968 to 1976. The birth rate BRTH and death rate DTH for the simulation is calculated to be the average of the actual death and birth rates of the eight year period, and this figure is converted into the monthly rate. These rates are BRTH 1.36 x 10‘3 DTH .775 x 10‘3. The seasonal variation of influenza is considered by making the infective rate sinosoidal, a sine curve which starts at the initial time equal to September. The sine curve represents the fact that influenza occurs most fre- quently during the fall and the winter and least frequently during the summer. We let INF infective rate (fraction/peOple month) CR = contact rate (fraction/people month) PRD = period = 12 months H = height of the sine curve (dimensionless). Then we can form the infective rate as INF = CR * (1+H + SIN((6.25/PRD) * TIME)). Integrating the above changes to the influenza epidemic model and simulating by varying parameters, we obtain a reasonable fit with the historical data. The comparison of simulation output and the actual figures for the three epidemics are given in Table h.3 and 87 .ptonom neapneunpm Hmom> manage: sons H m.« o.« n.m m.m «.0 0.5 o.m N.m m.u 000m gamma 5.8H m.4H 0.4a 0.4H o.ma m.sH 4.ma n.5a «.5H 000m spasm onH name name mesa mama Hema owma mooa mommy .HAooo.H p0dv m0pmpm v0pflca 0:» ca 00000 20000 020 zones Hmscc< .m.4 0Hnme 88 Table 4.3. Monthly mortality rates in the three epidemics (used in model validation. Epiggmic 1 Actuali Simulation Nov 68 159 ' 919 DEC 1068 1184 JAN 69 1359 _ 1304 FEB 677 1236 MAR 673 ' 1037 APR 695 759 MAY 545 574 Epidemic 2 Actual1 Simulation DEC 72 508 966 JAN 73 830 1029 FEB 904 1006 MAR 701 917 APR 398 796 Epidemic 3, Actugll Simulation DEC 75 479 691 JAN 76 819 710 FEB 824 714 MAR 606 702 APR 467 678 1 From Monthly Vital Statistics Report. 89 plotted in Figure 4.8. The values of the parameters and initial conditions used to obtain these results are as, follows: I Parameters b = .07 x 10-6(1 + .03 sin .52t) (infective rate) r = .8 (removal rate) e = .5 (susceptible rate) d = .6 (drift rate) Initial conditions 01 = 19.19 x 106 (susceptible population) c2 = .5 x 106 (infective p0pu1ation) c3 = .5 x 106 (immune population) ch = .001 (virus strain). The composition of the initial p0pu1ations is 80 percent susceptible, 10 percent infective, and 10 percent immune. Since the monthly mortality data is from the 10 percent sample of the total U.S. p0pu1ation, the simulations were performed using that figure as the initial total p0pu1ation. Figure 4.8 shows that the simulation results in terms of the epidemic peaks and the peak times are reasonably close to the actual data. This demonstrates a possible use of the model for prediction provided that the initial conditions and the values of the parameters can be measured. mopmu hpwampuos mo somwhmdsoo Hmoflsamhw “coeooomam> Homo: .m.e onomem nso— eso— mun. «No— eco— coo— ¢m< 52¢ cm“. 25. U8 «.2 «<2 mm“. 25. Ho >2 ¢m< :52 my“. 24.. U3 >02 fcom float rcoo room 1 TOGO— 100N— .ooe. 3503 ..<3pU< 20:.<.=..<§m I I I I I CHAPTER 5 IMMUNIZATION MODEL AND CONTROL SYSTEM In the influenza epidemic model deve10ped in Chapter 3, the cycle of epidemics begins as a result of a virus shift. The size of the susceptible population in relation to the threshold determines the severity of the initial epidemic; that is, if the difference between the initial susceptible p0pu1ation and the threshold is large, the peak of the epidemic curve is high. Subsequent epidemics are caused by drifts. The peaks of these epidemics are also determined by the difference between the susceptible p0pu1ation and the threshold. We now construct an immuniZation system to pro- vide the control input to the influenza epidemic system. Figure 5.1 shows the block diagram of this control system. INFLUENZA ‘2 EPIDEMIC :- SYSTEM . IMMUNIZATION SYSTEM Figure 5.1. Block diagram of the control system. 91 92 The goal of the control system is to minimize the severity of the second epidemic caused by drift. One method of pre- venting the occurrence of the second epidemic is to remove a sufficient number of individuals from the susceptible p0pu- lation so that the size of the susceptible p0pu1ation is always below the threshold level. This method could be costly in terms of a percentage of the total p0pu1ation required to be immunized. We, therefore, deve10ped an immunization system which accomplishes the goal of the control system by immuniz- ing as few individuals as possible. First we will examine the relationships among the variables which determine the size of the epidemic peak; then construct the state model for the immunization system; and finally, present the results of simulation. 5.1 Causal Analysis After the initial epidemic, the susceptible population starts building up through the transfer of individuals from the immune population due to the structural change of the virus. The decrease of the immune population causes the immunity level to fall, thus increasing the susceptible level. As the susceptible level increases the threshold decreases and the difference between the susceptible population and the threshold becomes larger causing the higher peak of the next epidemic. To remedy this situation, we add the vaccinated p0pu1ation in the causal loop diagram (Figure 5.2). We now decrease the susceptible population and follow through the 93 dotted lines. The causality indicates that the difference between the susceptible p0pu1ation and the threshold becomes smaller, causing a smaller epidemic peak. Epidemic (E ) Peak max 4. (x1) Susceptible Difference (xl-H) P0pulation _ ._ I I I I I - 7 (V) p1 (x3) Immune Vaccinated Threshold (H: 3—) POpulation Population Level I, - I / / x +v - // (A: —%—d Immunity’ + Susceptible (3:1-1) Level ~\\\\\‘-—_—fl’gsvel Figure 5.2. Causal loop diagram of the control system. 5.2 State Model of the Immunization System We construct the state model in such a way that the vaccinated population behaves like a goal seeking curve. This goal is the number equal to the difference between the susceptible p0pu1ation and the threshold level. We let v vaccinated p0pu1ation G = goal. Then 94 c(t) = x1(t) - H(t) where x1(t) H(t) = threshold level. susceptible p0pu1ation If the goal is constant; i.e. C(t) = G, then the state equa- tion for the goal seeking curve is given by I = k (0 - v) (5.1) with the initial conditidn v(0) = 0. The constant k is the vaccination rate. The solution for the state equation (5.1) is v = G (1 ~ e-kt) The graph of the vaccinated p0pu1ation is shown in Figure 5.3. t Figure 5.3. The goal seeking curve of the vaccinated p0pu- lation. Using the time varying goal C(t), the state model for the immunization system is then obtained as v = k (S(t) - V) (5.2) 95 with the initial condition v(0) = 0. The rate of immunization must be positive; hence, we consider applying the model only when the susceptible p0pu- 1ation exceeds the threshold. After the initial epidemic the susceptible p0pu1ation crOSses over the threshold just before the drift time TD.and the next epidemic begins just after the virus strain reaches the relative minimum. There- fore, the immunization system will be initiated at the drift time TD and continue up to the time the virus strain attains relative minimum Tmin (see Figure 5.4). This procedure can be accomplished by using a rectangular function (the difference of two step functions) to simulate the model. The state model (5.2) becomes v = k (S(t) - v)Tr where TT is the rectangular function with the height equal to 1 and the length equal to the period of the immunization. The immunity level A(t) is now the sum of the immune population and the vaccinated population divided by the total pOpulation: X + V A(t) = “l—N—n As before, the susceptible level S(t) is defined to be S(t) = 1 - A(t). Since the susceptible p0pu1ation decreases at the rate of immunization, the state equation for the susceptible popu- lation in the state model (3.2) is changed to x1 = -b S(t) x1 x2 + ex3 x4 - v. HW Hmuou cam mpmu :Ofipmcwoom> .o.m muswwm msuzos :0 mafia tIIOOOOOI. 0.0 I n. 00 L I: 5 W a a o a o o 0 0 0 0 0 0 0 5 0 0 O 0 a . 0 0 0 0 0 0 0 0 0 0 0 0 0 0 o o 0 . 0 0 0 0 0 0 0 o 0 0 0 0 0 0 0 0 0 . . . 0 . . . . . . . . . . 00.000.000.000...000.000.0000000000000000000000000.00.0000... . 0 0 0 0 o o 0 0 o 0 0 o 0 c 0 0 . o 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 000000000000000000000000000000000000000000000000.000000000000 a4. 4 JIM d 0‘44...“ 0‘10““44 J‘J‘JI“‘II‘JII“.‘““‘1{‘4D‘J N . 0 0 0 0 0 0 . . .cowpmasnom umumcwoom> . . . 0 0 0 0 0 o 0 0 0 0 0 0 0 . . . . 0 . . . . . . o o . 0 o 0 . 0 0 0 0000000000000000000000000000000000000000000000000000000000000 . 0 0 0 0 0 0 o 0 c . 0 0 0 o 0 0 0 0 n o 0 0 0 0 0 0 0 0 0 0 0 0 0 0 o 0 0 . 0 0 0 o 0 0 0 0 0 o 0 . '0°O§ 0 i 4114 d- 44 fl 0 . 0 o 0 0 0 0 0 0 0 0 0 o 0 0 00.0 0 0 0 0° 0 0 0 0“ O 0r...............0.0.0.0.... 0 0 .0 0 0 0 . .mnflmm . . .0 0 . cowmewoom> 0 0 0 0 0 0 0 0m 0 0 .0 .0 o 000 000000000000000000000000000O o 0 0 0 o 0 0 0 0 0 o 0 u 0 0 0 0 0 0 0 0‘ o 0 0 0 0° 0 0 0 0 0. 00000000000000000000000000000000 . 0 0 .0 0 0 . o 0 0 0 0 0 0 0 0 . . o 0 . .7 0 0 0 0° 0 0 . o .0 . 000000000000000000000. o UtOAhflfihfi-cflh0lfinofinnflnn0) “AAAAAAAAAAAAAfiAAAAAAh-AAAAI. linden» Ib‘uou .cowpmuwzssefi 909mm wusapso cowpmHSEHm .m.m mhswfim mnucos :0 mafia 100 I u u u m. u u n u u n .8 r .o no r r r r w r m. rm... I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I v I I I I I I I I I I I .0 . 00 . . . . . . . . . . L . . . . . . . . . . . . . 0 . . - . . . .0 . . . . 059:0 oflvaflam . . . . . s 0 0 0 0 0 0 0 0 o 0 0 0 0 0 0 0 0 0000: I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I . Ive . . . . . . , 0“ < I; J _ . . . .. . . mapwpawomzm. . . . . . . . . 0 0 0 0 0 0 0 0 . . . . . . Uflormwhgh. . . .N . . . . . . . . 0 ...0u.. I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I r I u" 0' . . . . . . . . . . '- 0 0 0 0 0 - 0 0 0 0 0 . . . . . . . . . . . 0 0 0 0 0 - 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0'0 . . . . . . . . . . .Wn. . . . . . . . . . 0 .00. I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I a, . . . . . . . . . . .00 1 . . . . 0 . . . . . . . . . . . . . . . . . . . 0 . 0 . . . . . . 0 0 0 0 0 0 0 0 0 0 0 9. 0 0 0 0 0 0 0 0 0 0 0o 09.. 0 . . a 0 . 0 0 . - .0 0 an I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I ”m 0 3 .0 5 . I. H h n n. ”s H NIIJKI H=HL 534 MS 101 Therefore, to compute the percentage of the total p0pu1ation required to be immunized, taking into consideration the effi- cacy of the vaccine and assuming that the immune p0pu1ation and the susceptible population are vaccinated proportionally, we use the formula: 1 VPOP % POP = EFF x SPOP x 100 where % POP = percentage of total population required to be immunized EFF = efficacy of the vaccine (expressed in decimal) VPOP = vaccinated p0pu1ation from the simulation SPOP = the size of the susceptible p0pu1ation at the beginning of the immunization period. The percentage of the total p0pu1ation required to be immunized are calculated for each test case using the efficacy of the vaccine to be #0 and 80 percent. The results are given in Table 5.1. Table 5.1. Percentage of total population required to be immunized. d D VPOP SPOP EFF % POP EFF % POP 1.0 .338 178 682 .4 65 .8 33 .5 .L96 205 763 .4 67 .8 34 .3 .650 215 795 .h 68 .8 3a d drift rate D ; drift CHAPTER 6 SUMMARY AND CONCLUSIONS The influenza epidemic model was systematically built from the basic three epidemic models described by Bailey (1975). In the basic epidemic models, the variables which are of particular interest to epidemiologists, such as the peaks and durations of epidemics, were defined and analyzed. The numerical values of these variables computed from analy- sis were compared with those obtained from simulation. The results show that the Euler integration method used in DYNAMO was reasonably accurate for simulating these nonlinear models. The influenza epidemic model is an aggregated system with four state variables which includes the structural state of the virus as one of the state variables. The state equation for the virus strain was constructed using the idea of a logistic curve, the asymptote of which is the immunity level of the p0pu1ation. With a prOper choice of the drift rate, the system provides an epidemic curve which resembles under- damped oscillation. Linearization was used to analyze the system behavior and to compute the range of the drift rate required for this underdamped oscillation. The computational procedure for finding this range is based on a sufficient condition, expressed in terms of the coefficients of the 102 103 characteristic equation, for the third order linear system to be stable oscillatory and underdamped. With this model three epidemic cycles can be shown. The initial epidemic results from virus shift which is represented by the presence of low immunity level in the p0pu1ation. The following two epidemics result from drifts which cause the susceptible p0pu1ation to build up above the threshold level. For the first three epidemics, the dynamic behavior of the, virus, the immunity level, and the epidemic curve are similar to those graphically represented by Kilbourne. Model validation was performed using the data taken from the Vital Statistics Report of the United States. The model was modified to accomodate p0pu1ation growth and the seasonal variation of the infective rate. The attempt to fit the historical data on the mortality rate of the three selected epidemics in recent time to the simulation results was reasonably successful, considering the complexity of the real world system. The model analysis and simulation of the aggregated sys- tem provide a good insight into the behavior of the essential variables. In particular the threshold is one such variable whose behavior may not have been so obvious. The threshold is a time varying quantity which is inversely prOportional to the susceptible level of the p0pu1ation. The size of the differencebetween the susceptible population and the threshold prior to an epidemic determines the severity of that epidemic. A formula used to calculate the numerical value of 10h drifts was developed. It was shown from simulation results that the value of the drift is an indication of the severity of the following epidemic, measured in terms of the peak of the epidemic curve. The drift time, the time at which a drift occurs, was also defined. Knowing the drift time, the beginning of the next epidemic could be estimated. In this model the second epidemic occurs if the sus- ceptible population exceeds the threshold. To control this epidemic, an immunization system was deve10ped as a sub- system of a control system. The immunization model consists of the vaccinated population as a state variable. The state equation was constructed using the idea of a goal seeking curve, the goal of which is the difference between the sus- ceptible p0pu1ation and the threshold. The immunization period begins with the drift time and ends at the time the virus reaches relative minimum. The control system is an idealized model. It takes no consideration of the efficacy of the vaccine, and only the susceptible p0pu1ation is immunized. To translate the total vaccinated population obtained from simulation into more realistic figures, a simple formula was used to estimate the percentage of the total p0pu1ation required to be immunized. The following are suggestions for building a more complex model from this aggregated model: (1) The susceptible level of the population may be made more precise by defining it to be a function of general health and age as well as antibody level. 105 (2) The infective rate may be a function of season, mobility, the density of the p0pu1ation and the effective contact rate. (3) Geographical Spread of infection may be modelled using partial differential equations. (A) Population may be divided into three or four age brackets with each group having different infective removal and immune rates. (5) Subclinical cases may be considered by having a separate immune p0pu1ation for those cases. (6) Surveillance systems may be established to identify drifts and shifts of the virus. (7) Immunization systems may include the efficacy of the vaccine which is a function of the virus and the effect- iveness of surveillance. They may also include economic factors such as production, distribution and manpower costs. (8) A more precise mathematical representation of the virus may be developed as research in this area progresses. APPENDICES APPENDIX A DERIVATION OF THE FORMULA FOR THE FUNDAMENTAL MATRIX OF THE 2X2 MATRIX A IN WHICH THE EIGENVALUES ARE COMPLEX Let s1 = u + jv, 32 = u - jv be the eigenvalues. Then, from the theory of the function of matrices, s t s t eAt = Z1 e l + 22 e 2 , where 21’ 22 are the constituent matrices given by _ A - s2 I Zl ’ s - s l 2 Writing 21’ 22 in terms of the real and the complex parts, we obtain zl=§[I-j%(A-u1)] (A.1) and Z2 = El.(conjugate of Zl)' Hence, 2 Re 2 1 Now, using the Euler Formula we can write 106 107 = e(u + jv)t eut ert = eut (cos vt + j sin vt). (A.2) Multiplying Equations (A.l) and (A.2) and taking twice the real part, we obtain the formula eAt = eut [(cos vt) I + (% sin vt) (A - uI)]. APPENDIX B DYNAMO LISTINGS Table 8.1. Simple epidemics model. * SIMPLE EFIQZHICS NSDEL NOTE ' NOTE 2‘ WE! (7 L SPOP. V‘"“"’ . -.."‘-'D *3? QP'7. \JH. SUSCEPTIB L IPOP. 3"."- -‘z 3):: g." QT: IEDO'J' JK INFECTIVE TIB IVE P ..:.=-fi:¥}23PoP KPIPDP.K SUSCEP LE R IP0PP.P;.us:: ...:0P KPIPDP.K INFECT ' NOTé H INITlRLIZdTION LE POPULATION OPULATION C S=950 C 1300 NOTE . NOTE P 33‘ :-‘.' 1‘.- 3'7"" C BETA=.OOP INFECTIVE RATE PRINT SPOP ..'.,F’OHH PLOT SPO.— “.{%h”~mylku.fl N SPEC DT=.C01/:L Ln../LENGTH=5/PRFPE? -.l RUN BASE *EOF 108 109 Table B.2. General epidemics model. * GENERAL EPIDEMICS MODE NOTE NOTE LEVE1.? L SPOP K=SP P U+OIPEPOPR.UR SUSCEPTIBLE POPULATION. L IPOP VPIPOP.J+DTPIPQPR.JK INFECTIVE POPULATION hOREOP K=RPOP.J+DT%RPDPR.JK IMMUNE POPULATION NOTE R AT. - :5 R SPOP..RL:-OETA -EPOP. KRIPOP K R IPOPR RL=P ETAIEPOP RaszOP. R- —GAMNA*IPOP K R RPOPR.RL=OANNAPIPOP K NOTE NOFE INITIALIZATION N SPCIPI=S N IPOPPI N RPOP=R c 5:950 G 1:50 c R=o NOTE NOTE PARAMETERS c BETA=.002 INFECTIVE RATE c GAMMA=.8 REMOVAL RATE PRINT SPOP.IPOP,RPOP,RPOPR PLOT SPOP=E.IPOP=I.IMPOP -R/RPOPR=N SPEC DT=”001/PLTPE.=.1/LENGTH=10/PRTPER=.1 BET‘I‘FP. ("f 7.1 110 Table 8.3. Endemic model. * ENDEHIC MO} EL NOTE NOTE L ’ WEI. 13 L SPOP K: EFOF 5+DT’"~PPQ UK SUSCEPTIBLE POPULATION L IPO..nnIPC+ U+DT¢IPC.: .OH INFECTIVE POPULATION L RPOP.M=RPOP.J+DT H'PP? um IMMUNE POPULATION NOTE NOTE Rfifffi R SPOPP.RL=—BETP.IZPOP.KPI POP. HiEPSLfi R IPOPR.EL?EETA% SPOP KPEPOP. KP'L AAMAPIP .K R RPOPR.KL=( A'MIIIPb‘.&-EP5LN NOTE NOTE DISEASE INCIDENCE RPPANNAPIPOP.R I\3IFIDLIZATIDN F’ARAI”IE C BETA=.OQI C GAMMA=.S C EPSLNPSO TERS PRIN-I Sr CF 1.95:3; RP “3’ INCI PLOTSPOP-5.1POP=I,RPOP=§/1 SPEC DT-= 1/DLTPEH21/LENGTH2 RUN BASE . *EUR 0 8:950 C R=O C BETA=.PC: C EPSLT‘E‘ l '5‘: _ - SPEC HLTPERPE/LENEIHPEDO RUN CH1 *EOF INFECTIVE RATE REMOVAL RATE SUSCEPTIBLE RATE 7PRTPER=I 111 Table 3.4. Influenza epidemic model. * INFLUENZA EP1DEHIC EE DUEL NOTE NOTE LL: “SLR ' L SPOP. “.3 SFE'S. JJUT’SFOPR.JK SUSCEPTIBLE POPULATION L IPCIU. K32PQ:) J‘DT* IPLPR JK INFECTIVE POPULATION L RP UP. {QnPPEZ‘EP' J" I“: r5 FEVDPR. 4'34. IMMUNE POPULATION L V]. K=VL ‘ HEV‘IW -JE\ VIRUS STRAIN NUTF. NOTE. PAT" I‘: R SPUPR.KL=-PEE R IPOPR E ‘SETI R RPUPP.EE~< EE H LVL. K.*SP3P. h*IPCP. K+EPSLN*VI. K*RPOP.K TSE SVL. KKSPOP. K%IPDP. K“ GAMMA*IPDP. K «IFUP. P”EPSLN*VI. K*RPOP. K §D¥&R.PL ~uFIE JMEEVL K VI. K)§VI. K NOTE AUX/SUPPL VARIABLES A IMLVL.K*RPCP H/IOOO IMMUNITY LEVEL A SLVL. ' '=-E.-'13‘-EL./E. . SUSCEPTIBLE LEVEL S TH. K=EIUOU*(bwhnA/B TA))/(IOOO-RPDP K) THRESHOLD S INCI.h~”AVHHKIPU“ K DISEASE INCIDENCE NOTE INITIALIZATIDN C BETA=.OCE INFECTIVE RATE C GAMMA? r': RENO mm PM" a E C EPSLN”.3 SUB C‘ZPTI3'.E RATE C DELTA£=. DPIFK RniE PRINT SP”: 1FSP-PPOP,TH.IHCI VI IMLVL PLOT ]'( rM/ WLJi= n/VI=V SPEC D72 LJFL.‘LH~1/LENUTH loo/PRTPER- =1 RUN BfSE *EDR C DELTF? {'7 RUN CH1 C DELTA$.3 RUN CH2 C DELTéa.1 RUN CH3 C =6OO C 19:93.53; \ 112' Table 8.5. Model validation. * MODEL VALIDATTDH NOTE NOTE . LC ’i'fl {'3' L SPOP.K=SP§P.J+DTPSPOPR.JK SUSCEPTIBLE POPULATION L IPOP.K=IPOP.\'+DT*IPO R.JK INFECTIVE POPULATION L RPOP K=PPOP J+D'%RPCPR.JK IMMUNE POPULATION L VI.K=“1.Ur~l“V R. uK VIRUS STRAIN NOTE NOTE udiLS R SPOPR.AL=-IPF.H%SLVL.K%SP OP.KfilPOP.K+EPSLN*VI.K*RPOP.K X +BR‘H-4-é-P23-i-D. E'x-"Q'e‘i-éi-T .=’.-S.-‘{2P. K R IPOPP.AL=IAP h~fLVL.HhEisOP.K*IPOP. K- GAMMA*IPOP. K- DTHI*IPOP.K R RPU~~.PLnCAm:,»:Ku. K- EPSLN¥VI. K*RPOP. -DTHR*RPOP. K R VIR.KLrDELTA a : .K- VI. K)*VI. K NOTE NOTE AU>;/313”PL VARIABLES A INP.K%CPP(1+H*SIN€(5. QS/PRD)#TIME. K)) INFECTIVITY S MORT.~—¢7HZ&IP" !& MORTALITY RATE A POP.K=SP3P. Hfi IP OP."+RPOP. V TOTAL POPULATION S INCI.P~~AH.A:1POP h. DISEASE INCIDENCE A IMLVL. H: TF -DP. K/POP V IMMUNITY LEVEL NOTEVL' K=L*INLVL.K SUSCEPTIBLE LEVEL NOTE SEASONAL VARIATION C CT=.O?E~A CONTACT RATE C =.03' HIGHT STEP FUNC C PROP 2 PERIOD SINE NOTE - NOTE INITIALIZATION N SPUP=S N IPOP=I N RPOP=P N VI= C 831?.1“L C I=.5Eb C R=.5Eé C V=.OOJ NOTE . NOTE TAP~W’¢L}S C GAMNA=.B REMOVAL RATE C EPSLN2.5 SUSCEPTIBLE RATE C DELTA= é DRIFT RATE C BRTH=1.3u-—3 BIRTH RATE C DTHS=7.731 A DEATH RATE C DTHI=7.7E "4 C DTHR==7. '7‘.‘w.‘---~’!. PRINT POP.BPOP.IPOP’RPOP.INF,VI,IMLVL;MORT.INCI PLOT PCP=P,EPOPnS,RPOP=R.IPOP=I/VI=J/IMLVL=A/MORT=M/INCI=N SPEC DT:.1ffiLTFERxlfLEMGTH=120/PRTPER=1 RUN .BAEL 113 Table 8.6. Immunization model. * IMMUNIZATION MODEL AND CONTROL SYSTEM NOTE NOTE LEVELS L SPOP.K=SPOP. J+DT*SPOPR.JK SUSCEPTIBLE POPULATION L IPOP.V.= IPOP. J+DTAIPOPR.JK INFECTIVE POPULATION L RPOP.K=RPOP.J+DT£RPOPR.JK IMMUNE POPULATION L V1. K=VI. J+DT*V1R.UK VIRUS STRAIN NOTEbP' K: VPOP.xJ+DT*VPOPR.JK VACCINATED POPULATION NOTE RATES SPOPR.‘ =A. V-wW‘R" J34 IPOPR. #3: =3- TA '53! ...VL KKSPOP. K*IPOP. K-GAMMA‘H'IPOP. K RPOPR. 31.3 AG A. “’33’3. -‘ " IPOP. KmEPSLNfi-‘I'I. Kit-RPOP. K VPO'JR 3.4 _‘SKEAI‘KWDIFF. K~VPOP. KHi'VAC. K VIR. KL. =1C~EZI-"2".-‘\3':( INLVL. KNVI. K)*VI. K :TJX/SUPPL. ‘PVQIABLE S A. K=-BET‘ A;: VL. 3i*SPOP K*IPOP. K+EPSLN*VI. K*RPOP.K DD 4.; mm XSET?=:T:P(HGHT1 STTM1)-STEP(HGHT1 STTNE) VACCIN PERIOD 1 =.- STTM3=27 STTN2=LR IMLVL. K=€RPOP. K+VPOP.K)/1000 IMMUNITY LEVEL SLVL. V.=1-IPHLVL. SUSCEPTIBLE LEVEL" TH. K=RHO3/E :3”!_ THRESHOLD DIFF V=SPOP .n-TH K SUSCEP-THRESH INCI.H=GANMA%IPOP.K DISEASE INCIDENCE IN i '3' I Al- I ZAT I OH OTE P AR A3-33'STE P. S RHOI=40§ ELATIVE REMOVAL RATE cocoonzzooooozzzzzggbbbbbooopbzzmmmmm ‘44 mm KSAIO=I'. “ " VACCIN RATE PLOT SP‘sz. * ---- W/LMCI—i! SPEC DT-= 1/PLTPER=I/LEi"GTH=IOO/PRTPER=1 RUN BASE *EOR PLOT VP: — 3:: ‘13"43’ 3 RUN CHI” BIBLIOGRAPHY Aggarway, J. K. 1972. Notes on nonlinear systems. Van Hostrand Reinhold Company, New York.' . Bailey, Norman T. J. 1957. The mathematical theory of epidemics. Hefner Publishing Company, New York. . 196A. The elements of stochastic process. John Wiley and Sons, Inc., New York. . 1975. The mathematical theory of infectious disease. Hafner Press, New York. Barker, D. J. P., and G. Rose. 1976. Epidemiology in medical practice. Churchill Livingston, New York. Baroyan, O. V., gt. g1. 1971. Computer modeling of influenza epidemics for the whole country (USSR). Adv. Appl. Bartlett, M. S. 1960. Stochastic population models. Methuen and Co., LTD., London. Beveridge, W. I. B. 1977. Influenza. Prodis, New York. Center for Disease Control: Influenza surveillance report no. 91, issued July 1977. Center for Disease Control: Morbidity and mortality weekly report, annual summary 1977, issued September 1978. g9 (53): 61-63 . Dayananda, P. w. A., and W. L. Hogarth. 1978. Optimal health programs of immunization and isolation for some approximations to chain-bronomial epidemic models. Mathematical Biosciences. 51: 2h-251. D'Azzo, John J., and Constantine H. Houpis. 1975. Linear control system analysis and design. McGraw-Hill, New York. Dickson, L. E. 1922. First course in the theory of equations. John Wiley and Sons, Inc., New York. 11h 115 Elveback, Lila R., John P. Fox, and Eugene Ackerman. Stochastic simulation models for two immunization problems, p. 90-103 in Ludwig, D., and K. L. Cooke, eds. 1975. Epidemiology: Proceedings of a SIAM Conference on Epidemiology. Alta, Utah, July 8-12, 197A. Society for Industrial and Applied Math, Philadelphia, Pennsylvania. Fenner, F. J., and D. 0. White. 1970. Medical virology. Academic Press, New York. Fenner, F. J., gt, a1. 197A. The biology of animal viruses. 2nd ed. Academic Press, New York. Forrester, J. w. 1968. Principles of systems. wright- Allan Press, Cambridge, Massachusetts. Fox, J. P., gt. a1. 1968. Epidemiology. Macmillan, London. Goodman, Michael R. 1974. Study notes in system dynamics. Wright-Allan Press, Cambridge, Massachusetts. Hinshaw, V. S., gt. g1. 1978. The prevalence of influenza viruses in swine and the antigenic and genetic related- gess of influenza viruses from man and swine. Virology. _g: 51- 2. Holland, W. W., ed. DATA handling in epidemiology. 1970. Oxford University Press, London. Kendall, D. B. 1956. Deterministic and stochastic epidemics in closed populations. Proc. Third Berkeley Symp. Math. Statist. and Prob. 5: 1h9-65. Univ. California Press, Berkeley. Kilbourne, E. D. 1975. Epidemiology of influenza, p. 483- 538 in Kilbourne, E. D., ed., 1975. The influenza viruses and influenza. Academic Press, New York. Kojima, K., ed. 1970. Mathematical topics in population genetics. Springer-Verlag, New York. Laver, W. 0., and R. G. Webster. 1979. Ecology of influenza viruses in lower mammals and birds. British Medical Bulletin 12 (1): 29-33. May, R. M. 197h. Model eco-systems. 2nd ed. Princeton University Press, Princeton, New Jersey. Olinick, Michael. 1978. An introduction to mathematical models in the social and life sciences. Addison-Wesley, Reading, Massachusetts. Oxford, J. 8., and J. D. Williams. 1976. Chemotherapy and control of influenza. Academic Press, New York. 116 Pereira, M. S. 1978. 'Flue: the mutant virus. World Health, July 1978, WHO, Geneva, Switzerland. . 1979. Global surveillance of influenza. British Medical Bulletin 12 (l): 9-lh. Perkins, F. T., and R. H. Regamey, eds. 1977. International symposium on influenza immunization. International Association of Microbiological Societies, New York. Pugh, Alexander. 1976. DYNAMO user's manual. The MIT Press, Cambridge, Massachusetts. Sage, A. P. 1977. Methodology for large scale systems. McGraw-Hill, New York. Sanders, M., and M. Schaeffer. 1971. Viruses affecting man and animals. Warren H. Green, Inc., St. Louis, Missouri. Saucedo, Robert, and Earl E. Schiring. 1968. Introduction to continuous and digital control systems. Macmillan, New York. Savageau, M. A. 1976. Biochemical systems analysis. Addison- Wesley, Reading, Massachusetts. Selby, P., ed. 1976. Influenza: Virus, vaccines, and strategy. Sandoz Institute for Health and Socio- Economic Studies, Academic Press. Smith, Gordon C. E. 1976. Epidemiology and infections. Meadowfield Press, Busbey, England. Spicer, C. C. 1979. The mathematical modeling of influenza epidemics. British Medical Bulletin 22 (1): 23-28. Stuart-Harris, D. H., and G. C. Schild. 1976. Influenza. Publishing Sciences Group, Inc., Littleton, Massachusetts. Supplement. 1977. The Journal of Infectious Diseases. 136 (December 1977). Susser, Mervyn. 1973. Causal thinking in the health sciences. Oxford University Press, New York. Tyrrell, D. A. J., and J. W. G. Smith. 1979. Vaccination against influenza. British Medical Bulletin 12 (1): 77-85 0 Volk, W. A. 1978. Essentials of medical microbiology. J. B. Lippincott Company, New York. 117 Waltman, P. 197a. Deterministic threshold models in the theory of epidemics. Springer-Verlag, New York. Wilkinson, James H. 1965. The algebraic eigenvalue problem. Clarendon Press, Oxford.