key: cord-0755851-a9fds45e
authors: Bombardt, John N.
title: Congruent epidemic models for unstructured and structured populations: Analytical reconstruction of a 2003 SARS outbreak
date: 2006-06-09
journal: Math Biosci
DOI: 10.1016/j.mbs.2006.05.004
sha: cc1473a772ad2b656313caf604e9e2db41c4a90a
doc_id: 755851
cord_uid: a9fds45e

Both the threat of bioterrorism and the natural emergence of contagious diseases underscore the importance of quantitatively understanding disease transmission in structured human populations. Over the last few years, researchers have advanced the mathematical theory of scale-free networks and used such theoretical advancements in pilot epidemic models. Scale-free contact networks are particularly interesting in the realm of mathematical epidemiology, primarily because these networks may allow meaningfully structured populations to be incorporated in epidemic models at moderate or intermediate levels of complexity. Moreover, a scale-free contact network with node degree correlation is in accord with the well-known preferred mixing concept. The present author describes a semi-empirical and deterministic epidemic modeling approach that (a) focuses on time-varying rates of disease transmission in both unstructured and structured populations and (b) employs probability density functions to characterize disease progression and outbreak controls. Given an epidemic curve for a historical outbreak, this modeling approach calls for Monte Carlo calculations (that define the average new infection rate) and solutions to integro-differential equations (that describe outbreak dynamics in an aggregate population or across all network connectivity classes). Numerical results are obtained for the 2003 SARS outbreak in Taiwan and the dynamical implications of time-varying transmission rates and scale-free contact networks are discussed in some detail.

The presence and likely impact of contagion are critical issues that surround a bioterrorist attack, the natural emergence of a new disease, or the accidental release of a dangerous pathogen into the natural environment. Even if effective countermeasures are available, a failure to completely contain a contagious disease may give rise to a far-reaching epidemic with considerable social disruption and economic impact. Both civilian and military homeland biodefense planners are thus faced with a broad spectrum of contagious disease threats and a sound conceptual framework is necessary to shape and plan threat responses. Arguably, well-designed epidemic modeling studies can contribute significantly to such a conceptual framework.

Severe acute respiratory syndrome (SARS) is the first new contagious disease to emerge in the 21st century. Initial SARS cases occurred in southern China during November of 2002 and international air travel quickly spread the disease through Asia, North America and Europe. On the 5th of July 2003, the World Health Organization announced the global containment of SARS cases. This international outbreak was responsible for about 8100 cases and 775 deaths. Moreover, Asian economies (for example) suffered losses in the neighborhood of 15 billion dollars (plus additional billions of dollars for economic stimulus packages) [1] .

The SARS outbreaks of 2002 and 2003 are instructive for several reasons. First, modern pharmaceuticals could neither prevent SARS infections nor abort the natural course of SARS cases. Each affected country controlled the spread of SARS in traditional ways: (a) case identification and contact tracing, (b) case-patient isolation (hospital infection control), and (c) exposure management (quarantine). Second, disease surveillance shortcomings and public health capacity limitations constrained SARS responses of Asian governments. Asian public health systems were clearly under stress and international assistance may have been critical in avoiding much larger outbreaks. Third and last, a great deal of SARS epidemiological information has already been published and more SARS data is likely to appear in the scientific literature.

The available SARS outbreak data has prompted researchers to reexamine the question of how best to quantitatively describe outbreak dynamics when traditional outbreak controls are the only viable countermeasures [2] [3] [4] [5] . Good answers to this question have been (and will be) essential to the development of sound biodefense requirements and plans.

Epidemic models have long been valuable tools for studying the dynamics of contagious diseases in human populations. Assuming an unstructured population and the standard incidence, disease transmission occurs by means of homogeneous mixing, where each contagious individual is free to contact and infect any susceptible individual. But if the population is structured according to cultural, socio-economic, demographic or geographic factors, there is a mixing matrix that constrains opportunities for disease-causing contacts. Proportionate and preferred mixing matrices have enhanced the explanatory and predictive powers of epidemic models since the late 1980s [6] . It is also true that sufficient data is seldom accessible to quantitatively define key mixing matrix elements for a large and highly structured population.

Viewing a structured population as a contact network, nodes or vertices in this network represent individuals and links or edges indicate potential disease-causing contacts. Mounting evidence shows that various real-world networks are scale-free; i.e., the probability distribution for the number of links per node (or the node 'degree') follows a power law [7, 8] . This power law is associable with network evolution or growth. More specifically, in contrast to classical random graphs, a scale-free network is a growing open system that accommodates a new node by attaching it to certain preexisting network nodes. The growth process for a scale-free network gives rise to another distinguishing feature, correlated degrees of connected nodes [9] . In the absence of node degree correlation, network-based epidemic models often entail proportionate (random) mixing assumptions.

May and Lloyd seem to have been the first researchers to incorporate scale-free network characteristics in a basic SIR (susceptible-infectious-removed) model for structured populations [10] . These investigators relied upon a random mixing matrix and treated the discrete node degree as a continuous variable. They derived mathematical expressions and obtained numerical results describing (a) the infected fraction of the entire population as a function of a basic reproductive parameter and (b) the infected fraction by node degree for selected parameter values. These researchers focused on epidemic threshold issues and they did not discuss how a scale-free network could affect outbreak dynamics.

Barthelemy, Barrat, Boguna, Pastor-Satorras and Vespignani have published a series of interesting and informative papers that explore network effects on outbreak dynamics [11] [12] [13] [14] [15] [16] . They have mathematically defined important properties of scale-free networks and examined epidemiological implications of these properties using susceptible-infectious (SI), susceptible-infectioussusceptible (SIS) and SIR models. In their network-based SI, SIS and SIR analyses, Barthelemy and coauthors routinely assumed an unbounded total population and employed a uniform initial condition (forcing all node degree or connectivity classes to be initially infected at the same level); additionally, they concentrated on an epidemic's exponential growth phase. The principal finding of these analyses is that, once high connectivity classes become infected, a cascade of infection occurs and affects lower and lower connectivity classes over time. Lastly, with regard to a non-uniform initial condition, the investigators found that confining all initially infected individuals to a single high-connectivity class would lead to an outbreak involving a faster rise and a faster invasion of the network.

Two published studies describe stochastic SARS outbreak simulations based upon lattices or grids and 'small world' linkages. Masuda, Konno and Aihara [17] argue that their simulations of the 2003 SARS outbreak in Singapore (as well as the general nature of SARS epidemics) are inconsistent with clustering properties and epidemic thresholds of scale-free networks. On the other hand, findings of the present author indicate that a finite scale-free network with node degree correlation is consistent with the 2003 SARS outbreak in Taiwan. Small and Tse [18] discuss their simulation of the 2003 SARS outbreak in Hong Kong and they make the following observations: (a) a small world network with long distance links can exhibit scale-free behavior and (b) further work must be done to ascertain whether their simulation results reveal that type of behavior. Although clear evidence of scale-free contact networks is at hand for sexually transmitted diseases [19] , similar published evidence of SARS-related contact networks is currently unavailable.

The main objectives of the present study are to: (a) develop congruent epidemic models for unstructured and structured populations that can better exploit currently available information on SARS progression and outbreak controls; (b) demonstrate the utility of these models by carefully reconstructing the dynamics of a historical SARS outbreak; and (c) employ historical time-varying rates of SARS transmission to assess effects of control measures in possible future outbreaks. Zoonotic infections and ordinary human activities in Asia may again be responsible for major SARS outbreaks, but it is also conceivable that a bioterrorist attack (for example, human vectors and intentional person-to-person SARS transmission) would cause even more human suffering and deaths. Epidemic models are clearly needed to quantitatively understand both natural and unnatural outbreaks of contagious diseases and to help improve or delimit outbreak response capabilities of health care systems.

This paper puts forward two complementary epidemic models for unstructured and structured populations. The SEIR-like compartmental framework 1 for both models is deterministic and encompasses the standard incidence. An 'X' compartment (with X S , X E and X I sub-compartments) is also a feature of this framework and it facilitates quantitative characterizations of outbreak control measures (quarantine, contact tracing and patient isolation). Most importantly, in order to enhance epidemiological realism, the two models accommodate a time-varying transmission rate and probability density functions (PDFs) for several important time intervals: namely, latent (incubational), contagious, onset-to-diagnosis, admission-to-discharge (from a hospital) and onset-to-death.

When exponential PDFs govern the latent and contagious periods, only mean values of these random variables enter a typical SEIR model. But alternative SEIR formulations make the assumption of exponentially distributed time intervals unnecessary. This is to say, by replacing differential equations (DEs) with integro-differential equations (IDEs), previous investigators have utilized non-exponential PDFs in fully characterizing latent and contagious periods [20, 21] . The present author adopts a similar approach and obtains IDEs containing general PDFs for all key time intervals. Although numerical algorithms for these IDEs can require substantially more computer time than those for exponential PDFs and associated DEs, the higher level of epidemiological realism helps to confidently resolve critical dynamical aspects of network-based models.

The finite and discrete scale-free network of primary interest is the product of a rigorous mathematical analysis and the derived node degree PDF is in good agreement with results of numerical simulations [22] . Network characteristics enter our epidemic model for a structured population through (a) the node degree PDF and (b) the Pearson correlation coefficient for degrees of connected nodes. The size of the total population and the node degree PDF determine sizes of all connectivity classes, which are the selected sub-populations for the epidemic model. Similarly, the aforementioned correlation coefficient can be related to the preferred mixing parameter for our structured population. This correlation coefficient has been evaluated for some scale-free social networks [23] , but it is essentially an unknown quantity in the realm of SARS outbreaks and 1 Typical SEIR (susceptible-exposed-infectious-removed) models allocate the fixed total population among compartments in accordance with a constant disease transmission rate and mean incubation and contagious periods. In these models, 'exposed' individuals are infected and non-contagious, 'infectious' people shed the pathogen and can infect others, while those in the 'removed' category are immune and no longer sources of infection. Removed individuals can be alive or dead and, if alive, they may still need medical care. contact networks. Numerical results in the main body of this paper thus cover a range of correlation coefficient values as well as different locations (connectivity classes) of the index case.

Utilizing a network-based model to analytically reconstruct a historical outbreak indicates the model's explanatory value and also ties its dynamical behavior to real epidemiological circumstances. For instance, if the analytical reconstruction of a 2003 SARS outbreak demonstrates a cascade of infection (from higher to lower connectivity classes) within the selected contact network, this resultant chain of disease transmission is at least connected with (a) the observed progression of one specific disease, (b) specific social (including health care) settings and (c) modern implementations of traditional control measures.

Since it has received relatively little attention from the epidemic modeling community at large, the 2003 SARS outbreak in Taiwan is analyzed herein. Ying-Hen Hsieh at the National Chung Hsing University (along with colleagues at other universities and government organizations in Taiwan) co-authored some interesting descriptive and analytical papers dealing with the Taiwan SARS outbreak. Hsieh and colleagues used SARS hospitalization and fatality data spanning one month (May 5 through June 4, 2003) to obtain parameter values for an epidemic model (i.e., a system of ordinary DEs) without quarantine or hospital isolation [24] . Subsequently, these researchers included a two-level quarantine and hospital isolation in another system of ordinary DEs and they performed a stability analysis that identified necessary quarantine rates for outbreak containment [25] .

A number of investigators have developed epidemic models (systems of ordinary DEs) that match certain temporal data from SARS outbreaks. In their study of SARS transmission and control in China, Zhou, Ma and Brauer assumed an inexhaustible supply of susceptible individuals, formulated a system of linear difference equations, and followed a process of trial and error in obtaining a time-varying transmission rate that enabled the model to closely match data on diagnosed cases per day [26] . A more complex model of the SARS outbreak in China was later put forward by Zhang, Lou, Ma and Wu [27] . These authors first quantified their 'basic adequate contact rate' using about three weeks of data on diagnosed SARS cases and then they obtained model results for the daily number of SARS patients (in good agreement with epidemiological data). Lastly, Gumel and co-authors studied a system of ordinary DEs with constant parameters and the chosen values of the two free parameters (transmission rates for contacts involving contagious individuals within and without hospital isolation) provided the 'best' agreement between calculated and reported cumulative SARS fatalities over time. After obtaining values of the two free parameters for four SARS outbreaks (Toronto, Hong Kong, Singapore and Beijing), calculations of cumulative probable cases over time generally conformed to reported data for each outbreak [28] .

The three previously-summarized studies are just a few examples of how simpler SARS outbreak models (systems of ordinary DEs) have been structured and parameterized to produce results that fully or partially agree with selected epidemiological data. Perhaps it goes without saying that any epidemic modeler must work with the data that is accessible. In other words, even though an epidemic model is capable of emulating a particular temporal feature of an outbreak, this does not mean that the underlying modeling approach is necessarily a systematic or suitable vehicle for analyzing outbreak dynamics on the whole. To examine the relative value of the proposed epidemic models of intermediate complexity, the present author also formulates and numerically evaluates a parallel or companion system of ordinary DEs.

The dashed curve in the lower graph of Fig. 1 identifies daily onsets of illness during the Taiwan outbreak [29] , but this curve excludes the index case to draw attention to the new infections. 2 These 670 onsets of illness (without the index case) were originally classified as probable new SARS cases. But the Center for Disease Control in Taiwan later published separate epidemic curves for both probable and laboratory-confirmed SARS cases [30] . Only 346 probable SARS cases were confirmed through laboratory testing. Arguably, the true number of SARS cases in Taiwan could have been larger (even substantially larger) than 346. As a consequence, an epidemic curve representing a total of 671 probable onsets in Taiwan is deemed to be adequate for our modeling purposes. One way of inferring new infections from symptom onsets is to uniformly shift these onsets backward in time according to the mean latent period. On the other hand, the Monte Carlo method is a more rigorous way to obtain average new infections per day. A useful Monte Carlo algorithm is straightforward and a single trial involves a tractable number of computations. First, obtain random latent periods for all onsets of illness (excluding the index case). Second, backtrack in time to determine when all infections began. And third, compile the total score for each time step. Averaging scores per day for a large number of Monte Carlo trials yields the desired infections per day.

The solid curve (vice points) in the upper graph of Fig. 1 comes from 50 000 Monte Carlo trials. In these trials, the assumed PDF for the latent period was a Gamma distribution, G 1 (x). The mean value (l 1 ) and standard deviation (r 1 ) of the latent period are in Table 1 with a supporting note and reference. This table likewise displays mean values and standard deviations of several additional random time intervals that, by assumption, follow Gamma distributions [31] [32] [33] [34] [35] [36] .

Solid curves in the upper and lower graphs of Fig. 1 are related by means of a mathematical convolution. Treating the time t as a continuous independent variable, the desired convolution integral can be expressed in the following manner:

The function p(u) is the new infection rate at time u, whereas G 1 (t À u) du is essentially the probability of developing SARS at time t after becoming infected at time u. Thusly, Eq. (1) defines the expected onsets of illness per unit time [37] . Regarding the numerical evaluation of this convolution integral, series of discrete p and G 1 values and fast Fourier transforms (FFTs) led to the calculated onsets (solid curve) in the lower graph of Fig. 1 . 'Nested' convolution integrals appear in many of the equations that follow and the notation of Eq. (1) can be extended to cover integrals of this type. Consider four continuous functions of time:

. Two examples using these functions demonstrate the extended notation for nested convolution integrals, namely, The new infection rate, previously defined PDFs, the size of the fixed total population, N 0 , and the number of initial or primary infections, E(0), are key elements of the basic epidemic model for an unstructured population. A general mathematical description of this model begins with a set of IDEs:

Remaining initial conditions become

Laplace transforms simplify the integration of Eqs. (5)-(7) and the results can be stated in terms of the time-dependent cumulative number of follow-on (new) infections, y(t):

Let z be the independent variable in the Laplace domain, let an underlined letter (e.g., E) denote a function in the Laplace domain and let L À1 signify an inverse Laplace transformation. Eq. (5) in the Laplace domain becomes zEðzÞ À Eð0Þ ¼ ÀEð0ÞG 1 ðzÞ þ pðzÞ À pðzÞG 1 ðzÞ À zX E ðzÞ þ X E ð0Þ;

or

Several inverse Laplace transforms are requisites for writing the last equation in the time domain: specifically, L À1 ð1:0=zÞ ¼ 1:0ðfor t > 0Þ;

Allocations of individuals among the E, I and R compartments are now defined by

EðtÞ ¼ Eð0Þ

IðtÞ ¼ Eð0Þ

and RðtÞ ¼ Eð0Þ

Integrals in Eqs. (11)-(13) are readily interpretable and the parallel progression of initial and follow-on infections is apparent. For instance, the first term on the right-hand side of Eq. (11) is the number of initially infected individuals who still remain in the E compartment at time t, while the difference between the second and third terms is the number of individuals whose follow-on infections began before time t and who have not yet moved into the I compartment.

Eqs.

(2)-(4) and (8)-(13) establish the general form of the basic epidemic model. Although the X E sub-compartment is not part of our analytical reconstruction of the Taiwan SARS outbreak, X S and X I sub-compartmental dynamics are indeed required to emulate, respectively, the quarantine of susceptible individuals and the isolation of SARS patients in hospitals. And because R(t) is just the time-dependent cumulative number of individuals who are no longer contagious, extensions of the basic epidemic model are necessary to evaluate time-dependent cumulative deaths and recoveries (or recovered SARS case-patients who are discharged from a hospital). At the end of an outbreak, the total number of individuals in the R compartment must be the same as the total number of SARS-related deaths plus the total number of SARS recoveries.

Six parameters characterize the sub-compartment for quarantined susceptible individuals, X S . The parameters u 1 and u 2 denote the respective proportions of initial and follow-on infections that result in hospital isolation upon symptom onset. 3 Half of the quarantine parameters are time constants: (a) delay time, t D , extending from day number 0 up to the start of contact tracing, (b) average time, t T , to trace and quarantine the contacts of a case patient entering hospital isolation and (c) quarantine period, t Q . Also, the average number of quarantined susceptible people per case-patient is k. Values of all quarantine parameters for the SARS outbreak in Taiwan and supporting references are in Table 2 .

If the diagnosis-to-isolation time is (on average) small, the time derivative of X S basically depends on diagnoses per unit time. Letting q 1 and q 2 be time-dependent numbers of diagnosed initial and follow-on cases per unit time, respectively, the delay differential equation for X S is

and

Eqs. (8) and (14)- (16) , along with the average new infection rate in Fig. 1 and parameter values in Tables 1 and 2 , are sufficient to evaluate the quarantine rate, which is input for Eq. (4). An early detection of symptoms or a later clinical diagnosis can trigger hospital isolation in the basic epidemic model. In other words, the hospital isolation sub-compartment, X I , accommodates initial and follow-on case-patients who are identified by either tracing activities or clinicians. Generally, a case-patient stays in hospital isolation until he or she is no longer contagious; the patient then continues to receive medical care for some period of time (ending in hospital discharge or death). During the 2003 SARS outbreaks, efforts to isolate contagious people were not always successful. Nevertheless, in the model, everyone with a case of SARS is eventually admitted to a hospital and each case-patient within the X I sub-compartment is assumed to be adequately isolated during the contagious period. The model does allow people with undiagnosed cases of SARS (who are outside the X I sub-compartment) to infect other individuals (who may be in hospitals or elsewhere).

Since a SARS diagnosis and the associated hospital isolation are assumed to take place almost simultaneously, case-patients in the X I sub-compartment fall into four distinct categories: (a) initial infection and diagnosis upon symptom onset, (b) initial infection and diagnosis after symptom onset, (c) follow-on infection and diagnosis upon symptom onset and (d) follow-on infection and diagnosis after symptom onset. As previously indicated, when case-patients in these categories cease to be contagious, they leave the X I sub-compartment. The time-dependent number of case-patients in the X I sub-compartment is below and it encompasses all four of the previously mentioned categories: 

Cumulative time-dependent numbers of case-patient recoveries, H R (t), and deaths, H F (t), are the final elements of the basic epidemic model. To mathematically define these cumulative quantities, proportions of initial and follow-on cases resulting in death are useful parameters: respectively, h 1 and h 2 . (See Table 2 for values of h 1 and h 2 .) Cumulative case fatalities are easily defined, i.e.,

Since the discharge of a case-patient from a hospital signifies a recovery, H R depends on the above four case-patient categories and the desired mathematical definition is

The complete model for an unstructured population spans Eqs. (2)- (4) and (8)- (19) . A known new infection rate and finite-difference replacements for Eqs. (4) and (14) enable numerical evaluations of X S and S at successive (full-day) time steps. FFTs make it possible to efficiently and accurately convolve two or more Gamma PDFs, but convolution integrals with p (or its definite integral, y) in the integrands call for a standard integration algorithm like the trapezoidal rule. Numerical results for compartments and sub-compartments in the basic epidemic model are discussed in a subsequent section of this paper.

Because average new infections per unit time are derivable from a given epidemic curve, the quantitative analysis of populated compartments and sub-compartments in the basic epidemic model does not invoke the standard incidence. But once results of this analysis are in hand, assuming the standard incidence allows us to characterize a historical outbreak by means of a time-varying rate of disease transmission, b(t), where

Under certain circumstances, historically derived time-varying rates of disease transmission may have predictive value [38] . The rate of disease transmission is essentially the product of (a) the probability of infection, given an 'adequate' contact, and (b) the rate of adequate contacts. That is, an adequate contact is a necessary (but not sufficient) condition for disease transmission. The conditional probability of infection could vary during an epidemic or outbreak if, for example, the disease-causing microorganism mutates and becomes more or less able to overcome the host's defensive mechanisms. Even in the absence of mutations, changes in environmental conditions might alter a microorganism's infectivity. Similarly, an actual rate of adequate contacts is likely to fluctuate with day-to-day human activities, imposed restrictions on these activities (e.g., quarantine) and a growing public awareness of an ongoing outbreak. An epidemic curve and an appropriately-derived time-varying rate of disease transmission may both be stepping-stones in resolving the fine structure of outbreak dynamics.

If a safe and effective SARS vaccine were to become widely available, public health officials might well rely on targeted vaccination programs (instead of quarantines) in controlling future SARS outbreaks. Minor modifications of Eq. (14) and the supporting parameters yield the following mathematical description of a populated sub-compartment that retains successfully vaccinated people:

where Eqs. (15) and (16) Tables 1 and 2 ; namely, l 2 , r 2 , l 4 , r 4 , l 5 , r 5 , h 1 and h 2 . Interestingly, l 2 and r 2 are the only chemotherapy-related parameters affecting disease transmission in the epidemic models at hand and the efficacy of chemotherapy need not be specified explicitly. But until there is accessible data on a specific therapeutic drug, establishing values of chemotherapy-related epidemic modeling parameters is largely a matter of guesswork.

Embedding a scale-free contact network in our basic epidemic model is a reasonably efficient way to explore disease transmission patterns in a structured population. That is, major functions of the desired network-based epidemic model include reconstructing historical outbreak dynamics with some fidelity and capturing disease transmission patterns at the sub-population level. These network-based modeling capabilities ought to facilitate the development of effective contact tracing, quarantine and/or vaccination plans.

Krapivsky and Redner [39, 40] have carefully analyzed a fundamental model of network growth. In that model, the growth process for a network includes (a) an initial condition, (b) the addition of one node or vertex at each time step and (c) the attachment of a new node to a single pre-existing (or ancestral) network node. The initial condition of interest is a 'dimer,' involving a primordial node that is linked to the first network node. 4 Most importantly, when the probability of a new node attaching to a pre-existing network node is proportional to the degree of the latter node, a scale-free network emerges in the manner of Barabasi and Albert [41] .

As the number of nodes in a scale-free network becomes infinite, the discrete PDF for nodes of degree j is

wherein j runs from 1 to infinity. With regard to correlating the degrees of linked nodes in an infinite scale-free network, researchers have obtained the following relative number of nodes of degree j that are linked (attached) to ancestral nodes of degree i:

with values of j as before and i running from 2 to infinity. The relative number of nodes of degree j with links to ancestral nodes of degree 1 is undefined.

Of particular interest in this paper is a finite scale-free network with N 0 nodes and N 0 links. The Krapivsky-Redner PDF for a finite scale-free network deviates from Eq. (21), especially at higher node degrees:

where erfc stands for the complementary error function and

If N 0 = 10 6 nodes, Eqs. (23) and (24) tell us that virtually no nodes are in degree classes above j = 1482. 5 In passing, this author is unaware of a published counterpart to Eq. (22) that delineates node degree correlation in a finite scale-free network. By and large, the basic epidemic model in the previous section is easy to disaggregate. Eqs. (2)-(4) and (8)-(10) can be restated for nodes (or individuals) in the jth node degree class (or jth connectivity class): viz., N j ¼ S j ðtÞ þ E j ðtÞ þ I j ðtÞ þ R j ðtÞ þ X j ðtÞ ¼ N 0 P ðjÞ; ð25Þ X j ðtÞ ¼ X S;j ðtÞ þ X E;j ðtÞ þ X I;j ðtÞ; ð26Þ _ S j ðtÞ ¼ Àp j ðtÞ À _ X S;j ðtÞ; ð27Þ

and y j ðtÞ ¼ Z t 0 p j ðuÞ du ¼ S j ð0Þ À S j ðtÞ À X S;j ðtÞ:

Invoking the standard incidence and employing the preferred mixing concept, the new infection rate for the jth connectivity class, p j (t), and corresponding transmission matrix elements, b jk (t)/N k , can be written as

Summation indices in Eqs. (31) and (32) range from 1 up to 1482, e is the combinative mixing parameter, a j (t) is the average number of disease-causing contacts per unit time per individual in the jth connectivity class, and d jk is the Kronecker delta [42] . When e = 0, members of different connectivity classes mix proportionately (randomly) and, when e = 1, there is no mixing between members of different connectivity classes. Furthermore, since each individual in the jth connectivity class is linked to j other individuals, a j (t) becomes

where x(t) is the average time-varying rate of disease transmission. The transmission matrix elements can now be restated as Given an individual in the jth connectivity class, P(kjj) is the conditional probability that he or she is linked to an individual in the kth connectivity class [43] . Utilizing Eqs. (32) and (33) 

Simply indexing Eqs. (14)-(16) is problematic, essentially because this action would quarantine the same number of susceptible individuals from any connectivity class containing a newly identified case-patient, regardless of that group's size. 6 Over 150 000 susceptible people were quarantined during the 2003 SARS outbreak in Taiwan and equating connectivity class PDFs for the quarantined susceptible population and the overall susceptible population is a reasonable approximation. Thus, in the jth connectivity class, the approximate number of quarantined susceptible individuals per unit time is

Eqs. (15) and (16) still define the respective aggregate numbers of diagnosed initial and follow-on cases per unit time, q 1 and q 2 . Of course, summing over the index j in Eq. (37) yields Eq. (14) .

Indexed versions of Eqs. (11)- (13) and (17)- (19) entail nothing more than substituting E j (0) for E(0) and y j (u) for y(u). Such straightforward substitutions obviate the need for a full listing of every indexed equation in the network-based epidemic model. In summary, the model for a structured population encompasses Eqs. 13) and (17)- (19) .

The last item of business with respect to the network-based model is a relationship between the combinative mixing parameter, e, and the Pearson correlation coefficient, r. Consider two linked individuals in a heterogeneous network, one of which is in the jth connectivity class. Pastor-Satorras and Vespignani [44] have shown that the average connectivity class (K nn ) for the remaining individual (nearest neighbor) is

and, in this paper,

In addition, Ramasco, Dorogovtsev and Pastor-Satorras [45] defined r as

so that Eq. (38) then allows us to write

The preferred mixing concept is related to the assortivity [46] of a contact network in that the combinative mixing parameter depends linearly on the Pearson correlation coefficient. A positive value of r (an assortative network) implies that individuals in higher connectivity classes tend to be linked to other individuals in high connectivity classes, while a negative value of r (a disassortative network) indicates individuals in lower connectivity classes are often linked to highly connected individuals. For our scale-free network with N 0 individuals and N 0 links, the second term on the right-hand side of Eq. (41) is 0.047. A null value of the combinative mixing parameter means that r is negative and, therefore, purely proportionate mixing is slightly disassortative. Values of r for certain real-world collaboration networks range from about 0.1 up to 0.4.

Wallinga and Teunis developed a likelihood-based algorithm that relies upon symptom onset data to determine the mean effective reproduction number, R e , over time [47] . Denoting L ij as the relative likelihood that case j was the source of infection for case i, the Wallinga-Teunis algorithm for SARS outbreaks is comprised of

In Eq. The mean basic reproduction number, R 0 , for the 2003 Taiwan outbreak is related to (a) the total number (670) of follow-on infections, (b) reported daily onsets of illness in Fig. 1 and (c) mean effective reproduction numbers in Fig. 2 . If O(n) symbolizes the reported number of illness onsets on day number n, the defining relationship for R 0 becomes

where n = 4, 14 and 114 are the respective day numbers for the onset of the index case, the immediately following onsets and the final onset. The resulting value of R 0 is 2.23, which conforms to comparable findings of Lipsitch and co-authors.

Since DEs (versus IDEs) represent a more conventional approach to the reconstruction of outbreak dynamics, a system of parallel or analogous DEs can serve as a baseline for assessing the merits of the basic and network models put forward here. The desired DEs encompass many of the parameters in Tables 1 and 2 and the principal free parameter is still the SARS transmission rate. This rate is assumed to be a fixed non-zero constant within a time interval corresponding to the duration of SARS transmission (in Fig. 2 ) and its numerical value assures a total of 670 follow-on infections. Employing lower-case labels for populated compartments and sub-compartments, the initial conditions for the DEs are

eð0Þ P 1; ð45Þ

Note that s, x s , e, c, c iso , r, h r and h f are, respectively, analogs of S, X S , E, I, X I , R, H R and H F . A system of eight DEs is mathematically defined below: _ cðtÞ ¼ ð1 À u 2 ÞeðtÞ=l 1 À cðtÞ=l 3 ; ð51Þ _ c iso ðtÞ ¼ u 2 eðtÞ=l 1 þ cðtÞ=l 3 À ð1=l 2 þ 1=ðl 2 À l 3 ÞÞc iso ðtÞ; ð52Þ

Although Eq. (20) still defines the functional form of the new infection rate in the above equations, an epidemic curve no longer constrains the time history of p(t) or, for that matter, b(t).

Of particular interest is the following specification of the disease transmission rate:

where the adjusted value of b FIX is 0.4899, U stands for the Heaviside step function and t off equals 110 days. Lastly, since definitions of all parameters in Eqs. (46)-(55) are already at hand, the reasoning behind compartmental entry and exit rates (as well as delays) is evident. Fig. 3 contains two sets of numerical results from the basic epidemic model (Eqs. (2)-(4) and (9)- (19) ) and one set of numerical results from the system of parallel DEs (Eqs. (44)-(55) ).

Regarding the basic epidemic model, red curves signify that both the backtracking Monte Carlo algorithm and the IDEs relied upon Gamma PDFs, while green curves are connected with an input of Exponential PDFs. Differences between red and green curves for the same compartment or sub-compartment are generally quite small; the exceptions are results for X I (n) and H R (n). But even small differences between red and green curves for either p(n) or I(n) can lead to more noticeable differences between red and green curves for b(n). This is to say that b(n) tends to follow the ratio of p(n) over I(n), assuming a small outbreak in a large population.

Black curves in Fig. 3 describe numerical results from the system of parallel DEs. The three curves for p(n) yield the same total number of follow-on infections, but the black curve is in stark contrast to the companion red and green curves. Of course, the slow rise and sharp decline of the black curve for p(n) are due to the assumed rectangular shape of the transmission rate. Slow rise times are likewise distinguishing aspects of the other black curves in Fig. 3 .

Since the basic epidemic model and Gamma PDFs may have yielded the most realistic results in Fig. 3 , some detailed aspects of the red curves are worth mentioning. As many as 3000 people were quarantined in a single day ($D + 95) and the maximum number of people in quarantine (on any given day) was found to be 50 000 ($D + 115). Similarly, the calculations show that nearly 175 (50) SARS case-patients were (were not) in hospital isolation on or about D + 85. And after D + 24 or so (20 days after onset of the index case), the calculated waveform of the time-varying transmission rate is close to the waveform in Fig. 2 . In passing, note that red curves in Fig. 3 and Eq. (20) gave rise to the points (vice solid curve) in the upper graph of Fig. 1 .

Focusing on the network-based model, the total number of new infections in the jth connectivity class is the subject of Fig. 4 . The contents of this figure originate in Eq. (36) , where the new infection rate, p j (t), is summed over all time steps (t = nDt and Dt = 1 day). With regard to the upper graph of Fig. 4 , the combinative mixing parameter is fixed at a value of 0.147 and the initial infection is either distributed across all connectivity classes (j SEED = PDF 7 ) or localized within a single class (j SEED = 3, 10, 30 or 100). In our network-based reconstruction of the Taiwan SARS outbreak, the disposition of the initial infection seldom influenced allocations of new (follow-on) infections to connectivity classes. But placing the initial infection in a higher connectivity class (specifically, j SEED = 100) did significantly increase the number of new infections in that class 8 (leading to the small blue spike in the upper graph of Fig. 4) .

The location of the initial infection is fixed (j SEED = 3) and the combinative mixing parameter takes on different values (e = 0.147, 0.247, 0.347 and 0.447) in the lower graph of Fig. 4 . As the value of e increases in this graph, 'shoulder' regions (between j ffi 30 and j ffi 400) of the curves become more pronounced; also, as j decreases from about 10 down to 1, smaller values of e give rise to slightly larger numbers of new infections. To better illustrate how the combinative mixing parameter affects the distribution of new infections within our contact network, consider the following definitions:

and g jk ¼ X n p jk ðnÞ: ð58Þ

The function p jk (n) is the time-dependent new infection rate in the jth connectivity class due to infectors in the kth connectivity class, while g jk is the resultant number of new infections. Density plots in Fig. 5 are associated with j SEED = 3 and selected values of e, and each of these plots shows variations of g jk for j and k ranging from 1 to 50. The lower-right density plot contains the color legend for the entire figure; e.g., a dark blue pixel indicates g jk 6 1 and a pure red pixel means g jk P 3. A 'pool' of red pixels (k 6 2 and j 6 12) and a few more red pixels along the connectivity class is relatively small, then the magnitude of the early spike is relatively large. Regarding the lower four graphs and a fixed location (j SEED = 3) of the initial infection, previously-considered values of e engender four time histories of the new infection rate in each connectivity class. In the 3rd (100th) connectivity class, lower (higher) values of the combinative mixing parameter produce higher (lower) peak values of the new infection rate. These influences of the combinative mixing parameter are consistent with curves in the lower graph of Fig. 4 . Fig. 6 does not reveal a cascade of infection affecting lower and lower connectivity classes over time. However, Fig. 7 displays new infection rates in very high connectivity classes (j = 300, 500 and 700) that are compatible with the concept of an infection cascade. More specifically, the maximum new infection rate in the 700th (500th) connectivity class precedes the peak rate in the 500th (300th) connectivity class. Perhaps it's worth noting again that sizes of very high connectivity classes are minuscule in our network-based model.

Numerical results for the time-varying rate of SARS transmission, x(n), are in Fig. 8 . Transmission rates for different values of j SEED span two time frames (n 6 25 and 25 < n 6 120) in the two upper graphs, and additional transmission rates for different values of e are displayed similarly in the two lower graphs. The fixed value of e for the upper graphs is 0.147 and the fixed value of j SEED for the lower graphs is 3. The seeding assumption shapes the early transmission rate, whereas the combinative mixing parameter is most influential at later times. And because all complete time histories in Fig. 8 produce the same number (670) of new infections, the lower graphs show that smaller e values entail higher late-time transmission rates. As e approaches a value of 0.047, the Pearson correlation coefficient vanishes and the connectivity classes of linked individuals become uncorrelated.

A network-based SARS modeling investigation by Meyers, Pourbohloul, Newman, Skowronski and Burnham [48] warrants some discussion before considering further variations of the present author's modeling parameters. Meyers and co-authors constructed a detailed contact network by focusing on households, schools, workplaces and hospitals and by drawing upon available data for 1000 households ($2600 people) in the city of Vancouver, British Columbia. Percolation theory and the node degree distribution for this contact network (with node degree correlation) enabled the investigators to theoretically and numerically predict outcomes of outbreaks for various fixed values of the transmissibility (transmission rate). Similarly, they obtained contrasting results for Poisson and scale-free contact networks (without node degree correlation). It is also worth noting that the modeling approach of these investigators yielded no information about outbreak dynamics.

Revisiting the upper-left graph in Fig. 8 , it's apparent that the connectivity class of the seed can strongly influence the early (6 day number 13 or so) transmission rate for the selected historical outbreak; i.e., to attain the same basic reproduction number, the early transmission rate for a seed in the 3rd connectivity class generally dominates the accompanying early rate for a seed in the 100th connectivity class. The clear implication is that the basic reproduction number for an unconstrained outbreak in a contact network depends on the connectivity class of the initial infection. In this vein, Fig. 8 complements and partially supports the impact of initial conditions as described by Meyers and co-authors.

Summing over connectivity classes in the network-based model generates aggregate populated compartments and sub-compartments that are identical to counterparts in the basic epidemic model. This confluence of epidemic models for structured and unstructured populations is due to a common modeling framework, common parameters and linked time-varying transmission rates (deriving from the same epidemic curve). Most importantly, under certain constraints, dynamical impacts of parameter variations in the basic epidemic model will also be aggregate dynamical impacts of the same parameter variations in the network-based model. The main constraint calls for maintaining the correspondence of transmission rates in the two models. A variation of parameters in Table 1 , for example, may significantly alter the dynamics of the Taiwan outbreak through changes in the disease transmission rate, which could be derived anew using the basic or network-based model. Alternatively, variations of parameters in Table 2 could be explored by utilizing either model with the appropriate previously-derived rate of disease transmission. In bringing this section of the paper to a close, quantitative results from the basic epidemic model elucidate changes in outbreak dynamics due to variations of selected parameters in Tables 1 and 2 .

The onset-to-diagnosis PDF (G 3 ) governs the random time interval extending from the onset of illness up to a clinical diagnosis. In the models under consideration, this time interval is fundamentally the amount of time that a contagious individual is free to infect other people before he or she enters hospital isolation. It's noteworthy that onset-to-diagnosis periods (l 3 for either the Taiwan or Singapore (l 3 = 3.09 and r 3 = 2.50) outbreak. This information leads to the following question: How would longer (Hong Kong or Singapore) onset-to-diagnosis periods have affected the dynamics of the Taiwan outbreak?

The upper graph in Fig. 9 shows onset-to-diagnosis PDFs and CDFs for Taiwan, Hong Kong and Singapore and the lower graph displays the resultant SARS transmission rates for the Taiwan been a substantially lower rate of SARS transmission and a concomitant increase (decrease) in the maximum number of contagious individuals who were not (were) in hospital isolation. One parameter that accounts for the effectiveness of both contact tracing activities and medical surveillance is the proportion (u 2 ) of follow-on infections entering hospital isolation upon symptom onset. A reasonable estimate of u 2 for the Taiwan outbreak is 45/670, which appears in Table 2 and comes from published information on probable cases and quarantine. To assess how the Taiwan outbreak would have been altered by larger values of u 2 , the basic epidemic model is implemented with a single 'underlying' transmission rate. (See the red curve in the lower graph of Fig. 9 .) Calculated time histories of new infection rates and compartmental populations appear in Fig. 11 for four values of u 2 (45/670, 90/670, 135/670 and 180/670). Unsurprisingly, more effective contact tracing and medical surveillance could have substantially reduced the total number of infections and the total number of deaths. The proportion (u 1 ) of multiple initial infections entering hospital isolation upon symptom onset is another parameter dealing with contact tracing and medical surveillance. In the event of a single initial infection and a new emerging infectious disease, the prompt isolation of the index case-patient is unlikely and a null value of u 1 is appropriate. But multiple individuals who are simultaneously infected with a familiar pathogen would have different incubation periods and sufficient information may be at hand to promptly identify and isolate some of the index case-patients. The basic epidemic model with one underlying transmission rate is readily implemented for multiple initial infections and non-zero values of u 1 . Fig. 12 contains time histories of new infection rates and compartmental populations for 10 initial infections and three non-zero u 1 values (0.2, 0.4 and 0.6). In terms of the total number of infections and the total number of deaths, the potential benefits of quickly finding initially infected individuals are also large.

The 2003 SARS outbreaks have yielded much epidemiological data covering both the course of the disease and modern implementations of traditional outbreak controls. The amount and quality of this data is unusual, spawning many simple and complex SARS outbreak models. Since few of those models were designed to exploit currently available data in a systematic manner, an epidemic modeling approach of moderate complexity has been described here that accommodates either an unstructured or a structured population.

A key thrust of this study is to describe the dynamics of an outbreak in new informative ways, especially when traditional outbreak controls are the only viable countermeasures. The epidemic models incorporate PDFs for a number of random time intervals that define the progression of initial and follow-on SARS infections as well as the implementation of traditional outbreak controls. Another significant dynamical aspect of these models is the semi-empirical time-varying rate of disease transmission. Using probabilistic reasoning to account for important random time intervals and deriving time-varying disease transmission rates from epidemic curves augment the realism and explanatory power of outbreak analyses.

Simpler SARS outbreak models in the literature often rely on a system of ordinary DEs with constant coefficients and models of this type have produced many important results concerning basic reproduction numbers and other dominant features of SARS outbreaks. To highlight the utility or value of the basic epidemic model, the present author formulated and evaluated a parallel simple system of DEs. The basic epidemic model and the parallel system of DEs generated (a) disparate new infection rates for the 2003 Taiwan outbreak (even though the total numbers of new infections were identical) and (b) very different time histories for each compartmental population.

When disease progression and outbreak controls are well-characterized, the basic epidemic model seems better suited to the task of reconstructing the overall dynamics of a historical outbreak.

The mathematical theory of scale-free networks is a work in progress and the full epidemiological utility of these networks remains to be seen. To date, ideal scale-free networks (infinite number of nodes and/or continuous node degree distribution) have been incorporated into relatively simple epidemic models (SI, SIS and SIR). The present author chose a finite and discrete scale-free network (with a rigorous mathematical foundation), developed a network-based deterministic model (with a realistic progression of infections) and implemented the preferred mixing concept to account for correlated connectivity classes of linked individuals.

The matter of cascading infections within a scale-free contact network warrants further discussion. Even though the initial infection is in the 3rd connectivity class, Fig. 7 displays new infection rates peaking at times that decrease (from $D + 75 down to $D + 53) as the number of links per infected individual increases (from 300 up to 700). Fig. 7 thus implies that individuals in very high connectivity classes can become infected well before those individuals in somewhat lower connectivity classes. But the upper four graphs in Fig. 6 suggest that the first generation of follow-on infections emerges by D + 25. (As the degree of the initial infection gets smaller and drops below 30 or so, the total number of infections in that connectivity class gets large and the first generation of follow-on infections becomes harder to identify.) After the first generation of follow-on infections, the upper half of Fig. 6 reveals remarkably similar exponential phases and general wave shapes for all four locations of the initial infection. And lastly, keeping the initial infection in the 3rd connectivity class, the lower half of Fig. 6 shows that different values of the combinative mixing parameter primarily affect magnitudes (rather than exponential phases) of new infection rates. In summary, our network-based analytical reconstruction of an actual outbreak exhibits no signs of cascading infections in connectivity classes containing more than about one individual.

Basic and network-based epidemic models under consideration here allow outbreak analyses at two levels of resolution, respectively, aggregate compartmental populations and connectivity classes. The common analytical framework for these two models assures congruence at the level of aggregate compartmental populations. Whereas the basic epidemic model is efficient in assessing the general dynamical impacts of disease progression and outbreak control parameter variations, the network-based model provides insights into disease transmission within and between connectivity classes. Both models were designed to support the development and/or evaluation of biodefense scenarios, requirements and investment alternatives. It's also worth emphasizing that biodefense investment decisions must often be made years before the formulation and implementation of new outbreak control strategies.

Coherent analyses of historical outbreak dynamics at multiple levels of resolution can broaden our knowledge of both disease transmission and the effectiveness of outbreak control measures. Perhaps more importantly, in establishing biodefense requirements and making biodefense investment decisions, retrospective questions ('what would have happened if. . .') deserve as much attention as questions about possible future events ('what will happen if. . .'). Historical outbreaks are obviously tied to particular settings and epidemiological circumstances that may or may not appear again in future outbreaks. Nevertheless, sound analyses and interpretations of historical outbreak dynamics could play a stronger supporting role in the biodefense planning process.

The epidemic modeling framework in this paper is somewhat flexible and adaptations of the SARS models for other diseases and/or other control measures (including prophylaxis and chemotherapy) are feasible. In reconstructing a primary pneumonic plague outbreak, for instance, log normal PDFs should be good choices for characterizing the progression of primary pneumonic plague infections and the insertion of these PDFs in the above models is an easy first step. On the other hand, determining a suitable scale-free contact network for primary pneumonic plague and getting enough data to evaluate the key parameters (Tables 1 and 2 ) could be more challenging tasks. It is interesting to note that the projection of respiratory droplets within a radius of 2 m or so is the main mechanism for transmitting both SARS and primary pneumonic plague from person to person. If the scale-free network in this paper is a reasonable vehicle for analyzing SARS transmission, then this same network might also be helpful in understanding primary pneumonic plague transmission.

The Corporate Research Program at the Institute for Defense Analyses funded this work and the author thanks Dr. Victor Utgoff (Institute for Defense Analyses) for his interest and encouragement. Opinions and findings of the author are not necessarily endorsed by the Institute for Defense Analyses or the US Department of Defense.

Emerging Infectious Diseases: Asian SARS Outbreak Challenged International and National Responses

Curtailing transmission of severe acute respiratory syndrome within a community and its hospital

Transmission dynamics of the etiological agent of SARS in Hong Kong: impact of public health interventions

Transmission dynamics and control of severe acute respiratory syndrome

Factors that make an infectious disease outbreak controllable

Epidemiological models for heterogeneous populations: proportionate mixing, parameter estimation, and immunization programs

Statistical mechanics of complex networks

Evolution of random networks

Organization of growing random networks

Infection dynamics on scale-free networks

Epidemic spreading in scale-free networks

Epidemic outbreaks in complex heterogeneous networks

Epidemic dynamics in finite size scale-free networks

Absence of epidemic threshold in scale-free networks with degree correlations

Velocity and hierarchical spread of epidemic outbreaks in scale-free networks

Dynamical patterns of epidemic outbreaks in complex heterogeneous networks

Transmission of severe acute respiratory syndrome in dynamical small-world networks

Clustering model for transmission of the SARS virus: application to epidemic control and risk assessment

Scale-free networks and sexually transmitted diseases: a description of observed patterns of sexual contacts in Britain and Zimbabwe

Disease extinction and community size: modeling the persistence of measles

Destabilization of epidemic models with the inclusion of realistic distributions of infectious periods

Finiteness and fluctuations in growing networks

The structure and function of complex networks

Modeling intervention measures and severity-dependent public response during severe acute respiratory syndrome outbreak

A discrete epidemic model for SARS transmission and control in China

A compartmental model for the analysis of SARS transmission patterns and outbreak control measures in China

Modelling strategies for controlling SARS outbreaks

Use of Quarantine to Prevent Transmission of Severe Acute Respiratory Syndrome

Memoir of Severe Acute Respiratory Syndrome Control in Taiwan

Epidemiology, transmission dynamics and control of SARS: the 2002-2003 epidemic

Viral shedding patterns of coronavirus in patients with probable severe acute respiratory syndrome

Quarantine for SARS

Model parameters and outbreak control for SARS

Epidemiological determinants of spread of causal agent of severe acute respiratory syndrome in Hong Kong

The transmission dynamics of human immunodeficiency virus (HIV)

Potential Influenza Effects on Military Populations, Institute for Defense Analyses Paper P-3786

Connectivity of growing random networks

Finiteness and fluctuations in growing networks

Emergence of scaling in random networks

Modeling heterogeneous mixing in infectious disease dynamics

Epidemic spreading in correlated complex networks

Dynamical and correlation properties of the internet

Self-organization of collaboration networks

Assortative mixing in networks

Different epidemic curves for severe acute respiratory syndrome reveal similar impacts of control measures

Network theory and SARS: predicting outbreak diversity