key: cord-0004525-5ra3lgvw authors: Greer, Meredith; Saha, Raj; Gogliettino, Alex; Yu, Chialin; Zollo-Venecek, Kyle title: Emergence of oscillations in a simple epidemic model with demographic data date: 2020-01-29 journal: R Soc Open Sci DOI: 10.1098/rsos.191187 sha: 0e41575d99d67cf66f4611f3f1d44dca610db5d1 doc_id: 4525 cord_uid: 5ra3lgvw A simple susceptible–infectious–removed epidemic model for smallpox, with birth and death rates based on historical data, produces oscillatory dynamics with remarkably accurate periodicity. Stochastic population data cause oscillations to be sustained rather than damped, and data analysis regarding the oscillations provides insights into the same set of population data. Notably, oscillations arise naturally from the model, instead of from a periodic forcing term or other exogenous mechanism that guarantees oscillation: the model has no such mechanism. These emergent natural oscillations display appropriate periodicity for smallpox, even when the model is applied to different locations and populations. The model and datasets, in turn, offer new observations about disease dynamics and solution trajectories. These results call for renewed attention to relatively simple models, in combination with datasets from real outbreaks. Mathematical models for disease reappearance across time in a specified geographical area typically incorporate delay terms, age structure or periodic forcing. These modeller-selected features permit fine-tuned adjustment of oscillations for the purpose of matching a dataset, and periodic forcing terms, in particular, can permit the teasing apart of different contributors to disease incidence. Still, a chicken-and-egg question must be asked: should observed oscillations in the data cause us to construct a model in which periodic behaviour is guaranteed, or should a model written without the expectation of periodicity be permitted to display oscillation, or not, depending on its parameter values and underlying demographics? More questions closely follow: what can we learn from sustained oscillations in a very simple model? What is the role of historical data in determining model parameters? What might we discover about disease transmission or population dynamics by parametrizing a model with demographic data? And how do results from a mathematical model help to highlight unusual or especially interesting data points from the historical record? We address these questions using historical smallpox and demographic data from three geographically separated regions: the Hida district of Japan, British India and Sweden. We begin with an overview of modelling history and results, then employ a series of computational approaches that constrain parameter values for each of the three regions by comparing model, data, and observed outbreak periods. We next circle back to the data, demonstrating the ability of a calibrated model to illuminate unusual historical occurrences. Throughout, we show that oscillation is an emergent property of simple models involving demographics. Periodic outbreaks have piqued the interest of researchers for over two hundred years. In 1929, statistician Soper [1] commented on the periodicity of measles outbreaks. He applied a second-order differential equation model, of the kind regularly used for mass-spring systems, to represent oscillations and connect with data. Soper refers [1] to work approximately twenty years earlier by epidemiologist Sir William Hamer. Hamer later writes about Noah Webster's 1799 volume A Brief History of Epidemic and Pestilential Diseases [2] as saying that influenza is a 'transitory scourge only afflicting the peoples of the world at times separated by long intervals of freedom' [3] . Current research more typically draws from the classic 1927 paper by Kermack & McKendrick [4] that introduced compartmental models, now often known by initials such as SIR for 'susceptible-infectiousremoved'. Sustained oscillations in differential equations-based SIR and related models are frequently described using delay differential equations, periodic forcing terms involving sine or cosine functions, and/or age structure. The marvellous 2000 review paper by Hethcote [5] describes seven articles that model sustained oscillation in epidemics, and all seven use one or more of these approaches. It is also frequently the case that modellers keep total population size constant. The sustained oscillation discussion in Hethcote [5] does not refer to models with non-constant population size. The 2008 text Mathematical Epidemiology [6] features two chapters that thoroughly describe the topic of oscillation in compartmental models. In Chapter 1, Earn provides references to the analysis showing that a fixed-population-size SIR model with constant parameter values displays damped oscillations toward a stable endemic equilibrium. As a follow-up, he uses [7] to show that stochasticity in model demographicsa staple in any real-world dataset-can sustain oscillations indefinitely. In Chapter 11, Bauch describes models having sustained oscillations as involving either exogenous or endogenous mechanisms. An exogenous mechanism, typically periodic forcing of the transmission term, explicitly introduces oscillations into the model. Endogenous mechanisms instead destabilize the model's endemic equilibrium. The most frequent examples of endogenous mechanisms are age structure or delay terms. Other examples include stochastic approaches, or a non-constant contact rate of infectious and susceptible individuals. Before and since Mathematical Epidemiology summarized these findings, both exogenous and endogenous mechanisms have continued to be incorporated into models. Networks [8, 9] and cellular automata [10] [11] [12] provide additional ways to understand oscillations in SIR models. Each new approach extends the concept of the compartmental model in a way that yields new insight. By contrast, we step back, invoking Occam's Razor and asking what can be learned from a simpler model. Such a model keeps the transmission parameter constant, along with all other parameters that directly represent aspects of smallpox: the basic reproduction number R 0 , length of infectivity, percentage of deaths due to smallpox. Only the parameters for overall births and non-smallpox deaths in the population vary, their values determined by annual computations based on historical demographic data, with annual values interpolated using straight lines to prevent discontinuous inputs to the system of differential equations. Periodic or quasi-periodic oscillations from such a reduced and simple model can then point to the fundamental ingredients of oscillatory dynamics in disease outbreaks. In addition to this model in which the parameters for births and non-smallpox deaths vary with time-hereafter referred to as the demographically forced model-we present a corresponding autonomous model. A formula for oscillation period in the autonomous model is analytically determined, and parameter values for the autonomous model are obtained through a combination of wavelet analysis and simulation-based experiments. Ultimately, the analytically determined formula royalsocietypublishing.org/journal/rsos R. Soc. open sci. 7: 191187 for oscillation period, when computed with the obtained parameter values, is shown to be consistent with demographic data. Overall this study demonstrates that (i) a simple dynamical model with timevarying births and deaths is sufficient to explain sustained periodicity of outbreaks, and (ii) demographic and disease outbreak data can be used to constrain model parameters such as the infection parameter, β. We work with the autonomous model and its demographically forced counterpart The estimated per capita death rate for infectious individuals is given by e, the average rate at which infected individuals recover is γ, and μ = γ + e is the overall rate of loss from I(t). As a note, some models write e ¼ d þ (deaths due to smallpox). We agree that e > δ yet compute values of δ and δ(t) from historical datasets and compute e from government-provided smallpox information, as will be shown in § §4.1 and 4.2. Due to these distinct parameter estimation approaches, e is best described separately from δ or δ(t) in models (3.1) and (3.2) . In this model, all newborns are susceptible to smallpox, and disease is transmitted via mass action incidence. The smallpox incubation period, typically 10-14 days [13] , is neglected, regarding these individuals as remaining in the S compartment: people move to I only when they become infectious. Vaccination is not considered in this article because the large oscillations in our datasets occurred before extensive vaccine campaigns. This is one of the simplest possible compartmental models with demographic effects, yet its pairing with historical data shows fascinating insight into oscillation and determination of parameter values. When parameters are held constant, the system has two equilibrium points: a disease-free equilibrium ((α/δ), 0, 0) and an endemic equilibrium Stability analysis shows that solutions approach the disease-free equilibrium in the case that αβ < δμ. In the case that αβ > δμ, the endemic equilibrium is asymptotically stable, as shown by the eigenvalues provided below and in equation (3.4) . Note that on the left-hand side of the inequalities involving αβ and δμ, α indicates replenishment of the susceptible population and β determines disease incidence. On the right-hand sides of the inequalities, all terms show removal of people from the infectious compartment (μ) or from the population as a whole (δ), in both cases decreasing opportunity for infection to spread. Therefore, the interpretation of the model matches the mathematical analysis royalsocietypublishing.org/journal/rsos R. Soc. open sci. 7: 191187 regarding which equilibrium is favoured, which depends on relative parameter values. Furthermore, R 0 for this system is computed [14] to be and R 0 . 1 corresponds exactly to the case when the endemic equilibrium is stable. In some cases, approach to the stable endemic equilibrium is via damped oscillations. Note that the relevant eigenvalues of the linearized system are − δ and Àab + ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where δ > 0, α > 0, β > 0 and μ > 0, indicating that all eigenvalues have negative real part. Damped oscillations occur when with the period T of damped oscillations, that is, the interepidemic interval, given by While β is not easily computed from epidemiological data, the quantity R 0 is regularly estimated for outbreaks. Substituting (3.3) into the square root in (3.5) , the interepidemic interval is Parameter values for smallpox shed light on when to expect oscillatory versus asymptotic approach to the endemic equilibrium. In §4, three historical datasets, along with biological smallpox information, are used to determine realistic parameter value ranges. In §5, multiple analyses of periodicity in the historical demographic data further narrow down the possible combinations of parameter values for the autonomous model. Finally, in §6, we summarize results and discuss insights into historical data based on model and data analysis. Parameter ranges for R 0 , μ, γ and e are computed using studies from multiple smallpox outbreaks. Values of R 0 have been estimated to be 3-5 in developing countries before the global eradication campaign [15] , 3.5-6 if a contemporary outbreak were to occur [16] , and 5-7 in files from the Centers for Disease Control and Prevention (CDC) [17] . Combining ranges, R 0 % 3 À 7. Computing μ, γ and e starts with the average infectious period. The CDC lists several phases of smallpox infection, along with how contagious an infected person is during each phase [13] . Taking these sequentially, smallpox is most contagious for approximately 14 days, and the total contagious period could be as much as 24 days. This means individuals remain in the I compartment 14 to 24 days; on an annual time scale, individuals remain in I for 14/365 to 24/365 years. Therefore, μ ∈ [(365/24), (365/14)]. About 30% of infected individuals historically have died of smallpox [18] , signifying that about 70% of the individuals who leave I subsequently move to the R compartment. Therefore, set e = 0.3μ and γ = 0.7μ, or set e ≈ 0.3μ and γ = μ − e. A word of caution to the modeller: permitting CDC values for infection length to be used for γ alone, that is letting γ ∈ [(365/24), (365/ 14)], and simultaneously permitting e to equal an additional 30% of the value of γ, results in a combined removal rate from I that is higher than indicated by demographic data. Three datasets are analysed in this paper: Hida, Japan (1771-1851), British India (1870-1920) and Sweden (1774-1800). Hida data come from [19] : population, deaths overall, and deaths from smallpox, each royalsocietypublishing.org/journal/rsos R. Soc. open sci. 7: 191187 computed annually. Annual numbers of smallpox deaths for British India and Sweden are reported in [20, 21] , respectively. Total annual deaths and total annual population values for both British India and Sweden come from [22] . The regions' smallpox death data appear visually in figures 1a, 2a and 3a. Each dataset for smallpox deaths is accompanied by its wavelet power spectrum in (b) and its global wavelet spectrum and Fourier spectrum in (c). The graph in (c) includes the dominant period for the data: 6.6 years for Hida, 5.5 years for British India and 5.2 years for Sweden. To compute δ for each year t, smallpox deaths are subtracted from total deaths for that year, and the resulting number is divided by the population size P t that year: To compute α for each year, population change P t+1 − P t is set equal to the difference between births and deaths. Rearranging gives a t ¼ P tþ1 À P t þ (total deaths) t : (4:2) Note that α t and δ t are the sole demographically determined parameters in this model. Therefore, population irregularities such as high death counts due to famine are encapsulated in the values of α t and δ t . By contrast, the value of e draws from on-average global smallpox estimates, as shown in §4.1. This means that any periodicity of smallpox deaths observed in model simulations results from the model and the demographic data; periodicity is not predetermined by a periodically forced parameter such as β or e. To distinguish notation for demographic parameter values: α t and δ t are the discrete year-by-year values computed as shown above. When using model (3.2) for simulations, α t and δ t are linearly interpolated to produce α(t) and δ(t). This prevents discontinuous inputs to the differential equations of model (3.2) , which improves the stability of numerical results. Full datasets and programs used to convert data to α t and δ t appear in this article's electronic supplementary material. The value of β does not draw as directly from historical or epidemiological data. Instead, R 0 is instrumental in determining the range of viable values for β. Equation Combined parameter values appear in table 1. The phase portraits in figure 4 encapsulate the differences between the autonomous and demographically forced models. The model (3.1) mathematical analysis in §3 details the damped oscillation and periodicity expected for any autonomous set of parameters, resulting in the consistent and smooth phase space behaviour in the left image of figure 4 . The demographic forcing of model (3.2) permits sustained and irregular oscillation, as in the right image of figure 4 ; such trajectories provide the opportunity to calibrate parameter values using a breadth of historical data. The following three sections detail our approach to computing model parameters using historical data. In §5.1, a wide range of R 0 and μ values are tested against demographic data to determine Full biologically viable ranges of R 0 and μ are tested; these values appear in table 1. For each historical dataset, the ranges of R 0 and μ form a two-dimensional region, as shown on the axes in figure 5 , and this region is divided into a 25 × 25 grid of fixed pairs (R 0 , m). For each pair (R 0 , m) in the grid, a model simulation is run in which α(t) and δ(t) take on their historical values as determined in §4.2. Their interpolation for smoother integration is completed using the default Interpolation command in Mathematica version 11. Each model simulation produces a mean period for I(t) oscillations; these mean periods are depicted using colour in the plots of figure 5 . Computation of the mean periods shown in figure 5 uses the argrelextrema package from Python's SciPy library to determine local maxima in model output. Within each simulation, times are recorded from each local maximum to the next, and the mean of these times is plotted as the overall period length for the simulation. Overlaid on each graph in figure 5 is a black dashed line showing the dominant outbreak period for that historical dataset, as computed in the spectral analysis of data shown in figures 1-3. Pairs (R 0 , m) along this dashed line are good candidates for inclusion in an autonomous model. For any such selected pair (R à 0 , m à ) the corresponding constant value of β à is computed using equation (4.3) and median values of α and δ. To summarize, demographic data are used to determine likely parameter values with priority on matching the mean periodicity of the demographically forced model. This approach produces pairs (R à 0 , m à ) and their corresponding β à . Section 5.2 considers a wider range of corresponding β values for pairs (R à 0 , m à ). For our selected pair (R à 0 , m à ) and each of 50 fixed and equally spaced values β ∈ [β low , β high ], the demographically forced model is run using α(t) and δ(t). The power spectrum is then computed for the output of each of these 50 model simulations. Figure 6 shows the collection of power spectra for the three geographical regions of study. Each region's graph consists of 50 vertical strips, each a power spectrum, one for each β ∈ [β low , β high ]. A black curve encloses outputs in the 95% confidence interval; these are coloured more boldly than outputs outside the 95% confidence interval. Overlaid on the graphs in figure 6 , and computed directly from the datasets shown in figures 1-3, are solid horizontal lines showing the 95% confidence interval band. The dashed horizontal line shows the peak of the global wavelet spectrum. The vertical line overlaid on each figure 6 graph is the value of β à that corresponds to the pair (R à 0 , m à ) used to create the graph, with β à determined via the method in §5.1. With best values of β thus tied to values of R à 0 and μ à , §5.3 turns toward determining fixed values of α and δ for use in the autonomous model (3.1). In §5.1, we obtained pairs of parameter values (R à 0 , μ à ) and corresponding β à based on simulations of the demographically forced model. In §5.2, we generated a range of β values corresponding to each (R à 0 , μ à ) pair, again using the demographically forced model. In this section, we investigate how well the autonomous model reproduces outbreak periods with tuned parameters μ and β. This approach centres on equation ( The analyses in § §5.1, 5.2 and 5.3 show that the simple, autonomous model (3.1) can be calibrated to produce oscillations with periodicity matching that of a historical outbreak. Though such oscillations are damped in an autonomous model, introducing stochasticity via realistic data or even via noise [6] permits the oscillations to extend indefinitely. When that stochasticity arises from real-life data, or from data with values similar to historical data, and when other parameters are selected accordingly, the oscillations remain in a biologically realistic range of periods. Our results from § §5.1, 5.2 and 5.3 show consistent periodicity between historical data and both the demographically forced and autonomous models. This consistency shows that biological disease information and historical demographic data can produce autonomous parameter values for model (3.1) . Resulting simulations of model (3.1) match outbreak periodicity in the short term, with oscillations damped because the model is autonomous. The corresponding demographically forced model (3.2) can match outbreak periodicity in the long term, including sustained oscillations, without the need of a tuned periodic forcing or delay term to impose the oscillation period. Autonomous values of α and δ are determined via an analysis of data. Median values of each, computed from demographic data, are a good first approximation for use in the autonomous model. The results of §5.3 provide more precisely tuned combinations of α and δ that best match the periodicity of the historical outbreak. Values of R 0 and μ are biologically determined from information on smallpox. Web information is readily available from the CDC [13, 17] , for instance. Because R 0 is unitless and μ is a rate, their values change minimally across different outbreaks. Medical care differences and public health approaches have some effect on R 0 and μ, but their values are well represented by ranges provided by the CDC and other national or international health centres. Compared with all other parameters in models (3.1) and (3.2) , the infection parameter β is the least easy to estimate from datasets or biological disease information. The value of β depends on population size and density, as well as various dynamical aspects of disease transmission such as types and frequencies of interactions. Indeed, the units of β show dependence on the number of people in the population. Therefore, the computational work in §5 is required to systematically narrow down a range of values of β based on model outputs and agreement with outbreak and demographic data. One outcome of our analysis is an overlapping set of values relating directly to the infection parameter. While the numbers found for β itself have very different orders of magnitude, rescaling by population size N = S + I + R produces values β 0 = Nβ that can be compared; here, N is the median population across all years of data. Figure 8 shows the range of values of β 0 ∈ [β low , β high ] for all three countries. It is worthwhile to note the correspondence between the geographical scale of data and the range of values for β 0 . A second outcome of our analysis is identification of parameter combinations (α, δ, μ, β) for which the autonomous model best matches the periodicity of historical data. These parameter combinations appear as bold dashed (α, δ) lines in figure 7 for which μ and β are fixed, with most historical pairs (α t , δ t ) lying close to these lines of maximum global power, and almost all historical pairs lying within the 95% confidence interval. This comparison of analysis with data shows that an autonomous model can match periodicity well for the short term. Longer term, historical birth and death rates introduce stochasticity, which sustains oscillations, yet almost always producing oscillation periods that match well with the signal from historical data. Of further interest are the exceptions to historical (α t , δ t ) pairs within the spectral analysis 95% confidence interval. British India has no exceptions, and Sweden's exceptions are few. Hida's exceptions include two years of famine, marked in red on figure 7. Hida is also the only case in which historical (α t , δ t ) combinations produce a non-negative square root in equation (3.4) . A nonnegative square root is a rare occurrence, to the extent that textbook analyses (such as [23] ) typically neglect the possibility so that formula analysis is more tractable: for example, the oscillation period then has a more tidy solution. A third outcome of our analysis is that its success is enhanced by longer datasets covering more homogeneous populations and geographical regions. In particular, the dataset for Hida is 81 years, the dataset for British India is 51 years, and the dataset for Sweden spans only 26 years. Further advantageous for Hida is that it is a relatively small and homogeneous population in a small geographical region. For these reasons, the Hida dataset appears to best lend itself toward this analysis. Hida's peaks of smallpox incidence are clearest, rising the most in comparison with nonpeak years and spanning just one to two calendar years. For contrast, while oscillation is noticeable in the British India and Sweden datasets, it is yet true that smallpox incidence in those regions encompasses a much larger geographical area, so that outbreaks spread spatially and therefore take up multiple years apiece. This smooths and spreads out some of the effects of oscillation. Length of dataset affects the oscillation periods modelled using demographic data, such as in figure 5 , where the pattern is quite similar in period lengths for Hida and British India, but the smaller Sweden dataset leads to more variation in period lengths as R 0 and μ change. These second and third outcomes, combined, support the interest in outlier points in the Hida dataset. Hida's data are well suited to these analyses due to both number of data points and homogeneity of population. As further support for this model, research [19] shows that, aside from famine in two separate years, smallpox was the main contributor to untimely deaths during the years of our dataset. Therefore, the (α t , δ t ) pairs lying outside the 95% confidence interval, or even outside the region of parameters for which model (3.1) oscillates for otherwise-appropriate (μ, β) combinations, may hold important information about the presence of smallpox in Hida. By using simple models to emphasize the connections between data and mathematical properties such as oscillation, we can thus highlight historical data points of interest for further study. This cycle, from data to model calibration to model simulation to re-emphasis on data, has the potential to reflect back to us historical data points of special interest. The fourth and final outcome of our analysis is that periodic forcing, delays and other exogenous mechanisms are not necessary ingredients for generating periodic behaviour in a model. Instead, the application of demographic forcing alone, via rates of live birth and of death, can produce sustained oscillations with similar periods as observed infection rates. This further suggests that the underlying dynamical mechanisms in disease outbreaks can be captured well by simple SIR models. Meanwhile, datasets from historical outbreaks worldwide have become widely available for study. We, therefore, call for a renewed interest in relatively simple mathematical models, in conjunction with data from real outbreaks. This combination of models with data will bring new insights into models, especially models involving realistic demographics. In turn, these models can shed light on special or unusual data points, for improved understanding of the relationship between stochastic real-world dynamics and the varying incidence of endemic disease. Data accessibility. The code and the datasets supporting this article have been uploaded as electronic supplementary material. Authors' contributions. M.G. estimated parameter values, co-developed and mathematically analysed the model and drafted the manuscript; R.S. estimated parameter values, designed and performed computational analyses and helped draft the manuscript; A.G. collected epidemiological and demographic data, curated information from many sources, carried out a literature review of prior work on periodic disease outbreaks, and provided critical direction and content for the manuscript; C.Y. characterized theoretical conditions in which model oscillations occur, estimated parameter values and initial conditions, and contributed critical intellectual content to the manuscript; K.Z.-V. developed initial model formulation and analytical solutions, collected epidemiological and demographic data, carried out a literature review of prior work on periodic disease outbreaks, and provided critical direction and content for the manuscript. All authors gave final approval for publication and agree to be held accountable for the work performed therein. The interpretation of periodicity in disease prevalence 1799 A brief history of epidemic and pestilential diseases The crux of epidemiology A contribution to the mathematical theory of epidemics The mathematics of infectious diseases Mathematical epidemiology Stochastic population models in ecology and epidemiology Modeling epidemics on adaptively evolving networks: a data-mining perspective Complex dynamics of epidemic models on adaptive networks Oscillations in an epidemiological model based on asynchronous probabilistic cellular automaton Big cities: shelters for contagious diseases Self-sustained oscillations in epidemic models with infective immigrants Centers for Disease Control and Prevention. Smallpox: Signs and Symptoms Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission Parasitic infections as regulators of animal populations: the dynamic relationship between parasites and their host populations offers clues to the etiology and control of infectious disease Transmission potential of smallpox in contemporary populations Smallpox: History of Smallpox Epidemics and mortality in early modern Japan A concise history of small-pox and vaccination in Europe Smallpox and vaccination in British India during the last seventy years Mathematical models in population biology and epidemiology Competing interests. We declare we have no competing interest. Funding. No funding has been received for this article.