Complex multiannual cycles of Mycoplasma pneumoniae: Persistence and the role of stochasticity

TL;DR

This paper explains the complex outbreak cycles of Mycoplasma pneumoniae using mathematical models and shows how randomness in host behavior helps sustain these cycles.

Contribution

The study introduces a novel explanation for M. pneumoniae's multiannual cycles using quasicycles and stochasticity.

Findings

Mycoplasma pneumoniae's 5-year cycle is explained by quasicycles and loss of immunity.

Environmental stochasticity stabilizes multiannual cycles in smaller populations.

Cycles can be disrupted by large perturbations like NPIs from the COVID-19 pandemic.

Abstract

The bacterium Mycoplasma pneumoniae causes respiratory illness worldwide and is known for its complex and poorly understood multiannual outbreak cycles. By fitting a mechanistic mathematical model to data from Denmark spanning several decades, we find a parsimonious explanation for the persistent 5-y periodicity, in the form of quasicycles. A deterministic model explains the periodicity and shape of the cycles, while environmental randomness—such as variability in host contact patterns—is needed to sustain them. The work has practical implications for predicting the epidemiological dynamics of M. pneumoniae across diverse settings and suggests that noise-induced dynamics should be considered as a potential driver of complex cycles in other endemic pathogens.

The epidemiological dynamics of Mycoplasma pneumoniae is characterized by poorly understood complex multiannual cycles. The origins of these cycles have long been debated, and multiple explanations of varying complexity have been suggested. Using Bayesian methods, we fit a dynamical model to half a century of M. pneumoniae surveillance data from Denmark (1958 to 1995, 2010 to 2025) and uncover a parsimonious explanation for the persistent cycles, based on the theory of quasicycles. The period of the multiannual cycle (approx. 5 y in Denmark) is explained by susceptible replenishment due, primarily, to loss of immunity. While an excellent fit to shorter time series (a few decades), the deterministic model eventually settles into an annual cycle, unable to reproduce the persistent cycles. We find that environmental stochasticity (e.g., varying contact rates) stabilizes the multiannual cycles and so does demographic noise, at least in smaller or incompletely mixing populations. The temporary disappearance of cycles during 1979 to 1985 is explained as a consequence of stochastic mode-hopping. The circulation of M. pneumoniae was recently disrupted by COVID-19 nonpharmaceutical interventions (NPIs), providing a natural experiment on the effects of large perturbations. Consequently, the effects of NPIs are included in the model and medium-term predictions are explored. Our findings highlight the intrinsic sensitivity of M. pneumoniae dynamics to perturbations and interventions, underscoring the limitations for long-term prediction. More generally, our findings provide further evidence for the role of stochasticity as a driver of complex cycles across endemic and recurring pathogens.

Figures

Model fits to M. pneumoniae incidence data. A seasonal SIRS model is fitted to recent (2010 to 2025) and historical (1958 to 1995) M. pneumoniae testing data from Denmark using Bayesian inference. (A) The seasonal SIRS model explains observed M. pneumoniae dynamics 2010 to 2025. The data are PCR positivity rates. The Inset shows a full cycle (the interval highlighted in blue in the main plot) in the (S,I) configuration space. Markers are 1 wk apart and the color indicates time of year. (B) The same type of model captures key features of observed M. pneumoniae dynamics 1958 to 1995. The fitted data are variance-stabilized numbers of samples which are simultaneously CA and MpCF-positive. The Shallow Bottom panel shows a smoothed version of the time series, to emphasize the multiannual peaks (Data source: Statens Serum Institut and ref. 2).

Demographic stochasticity excites multiyear cycles in smaller populations. (A) Theoretical power spectral density (PSD) of multiannual fluctuations excited by demographic noise (without seasonality), and avg. power spectrum of 50 stochastic simulations (with seasonality). See SI Appendix, Fig. S3 for plots of the dependence of the theoretical PSD on individual parameters. (B) Stochastic simulation of a seasonally forced SIRS model parameterized with M. pneumoniae parameters as determined by Bayesian fit to 2010 to 2025 data. In small populations (initial population size N0=20,000), demographic stochasticity leads to noisy multiyear cycles. (C) In intermediate populations (N0=250,000), multiyear cycles are prominent and fairly regular. (D) In larger populations (N0=2×106), annual dynamics dominates as the impact of demographic stochasticity is smaller.

Moderate environmental noise stimulates multiannual cycles, especially approximately 5-y cycles. Average Fourier spectrum across 40,000 realizations (for each noise level indicated in the legend) of the prevalence I(t) of the seasonal SIRS model with parameters from the fit to 2010 to 2025 data from Denmark. The noise source is random fluctuations in the transmission rate β(t), implemented as Gamma-distributed multiplicative noise with mean 1 and variance σ2 (Materials and Methods and ref. 59). 80 y of each simulation time series were analyzed, and a 100 y “transient” was discarded.

One-off exogeneous perturbations (“kicks”) transiently excite multiannual cycles. (A) System trajectory depicted in (S,I) configuration space before and after a perturbation. The marker color denotes the time of year, and markers are 1 wk apart. The approximately vertical line starting from the red cross shows an exogenous “kick” (a discrete change in prevalence, in this case) which causes the system to transition from the annual (thick yellow line) to 5-y cycle (thin purple line). (B) After a discrete change ΔI in the prevalence, the system may temporarily settle into a 5-y cycle, as in A. The plot shows how the size of the perturbation and the time of year affect the strength of the resulting cycle. The color bar indicating the power of the 5-y cycle in the Fourier spectrum (relative to the highest peak) is shared between B, C, and D. Spectra were based on the first 20 y post-perturbation. The red cross indicates the perturbation implemented in A. (C) If the perturbation ΔI is proportional to the preperturbation prevalence I0, the prominence of the 5-y cycle ceases to depend on time of year. The red cross indicates the perturbation implemented in A. (D) Effects of perturbations in S(t), rather than the prevalence.

Prediction is hard. (A) Prediction using the deterministic seasonal SIRS model, parameters fitted to 2010 to 2025 data. (B) Prediction using the seasonal SIRS model with environmental stochasticity (noise in transmission, σ=0.12, fitted value). Noise can significantly affect long- and medium-term predictions, but moderate to high M. pneumoniae activity is predicted in the 2025/2026 season. The 95% CI is shown in light purple, and the 85, 75, and 50% intervals are shown in successively darker shades.

Tables

Table 1.. Estimated and assumed parameters from fit to 2010 to 2025 M. pneumoniae surveillance data from Denmark

Parameter	Description	Estimated value [95% HPD]
β0 (determ.)	Mean transmission rate	0.65/wk [0.58; 0.73]
a (determ.)	Amplitude of seasonality	0.18 [0.17; 0.19]
δ (determ.)	Immunity waning rate	1/(8.2y) [1/12.0; 1/6.4]
β0 (stoch.)	Mean transmission rate	0.71/wk [0.65; 0.77]
a (stoch.)	Amplitude of seasonality	0.18 [0.15; 0.20]
δ (stoch.)	Immunity waning rate	1/(8.7y) [1/12.0; 1/6.9]
σ (stoch.)	Level of stochasticity	0.12 [0.11; 0.13]
Parameter	Description	Assumed value
μ	Birth/death rate	1/(80y)
γ	Recovery rate	1/(2.5wk)

Full text

The bacterial pathogen Mycoplasma pneumoniae is simultaneously common and unusual. It lacks a cell wall and is among the smallest free-living bacteria (1). The infection it causes has been given names such as (primary) atypical pneumonia (2), cold pneumonia, and walking pneumonia (3) due to its tendency to cause less severe disease compared with typical bacterial pneumonia. It is a very common respiratory pathogen, especially in children and adolescents, accounting for 15 to 20% of all cases of pneumonia (4). Disease courses are typically mild, but severe outcomes do occur. About 5 to 10% of upper airway infections with M. pneumoniae progress to pneumonia, and there are likely many asymptomatic cases (5), suggesting a high degree of underascertainment. Extrapulmonary symptoms are increasingly recognized (6, 7) and include complications such as encephalitis and RIME/Stevens–Johnson syndrome (8, 9), with M. pneumoniae being the most common infectious cause of the latter in children (10, 11). As many as 25% of M. pneumoniae cases develop some extrapulmonary symptoms (12), most commonly dermatological or gastrointestinal. Sometimes cardiovascular (13) and neurological (14) complications are seen, the latter in as many as 1 to 10% of those M. pneumoniae cases which require hospitalization (12, 15, 16).

From an epidemiological dynamics perspective, a remarkable feature of M. pneumoniae are the observed complex outbreak cycles which generally exhibit a 3- to 7-y period (17??–20). In Denmark (2), as well as England and Wales (21), 4 to 5 y cycles have generally been observed. This is in contrast to some classical, strongly immunizing childhood diseases, such as measles, which tend to exhibit simpler annual or biennial dynamics (22). However, there are exceptions to that rule, for example rubella (23, 24). Despite being such a common pathogen, the mechanisms behind the complex epidemiological dynamics of M. pneumoniae remain poorly understood, complicating prediction and forecasting. In general, even short and medium term predictability are important for timely resource allocation and deployment of relevant screening and diagnostic testing schemes. However, when the year-to-year incidence exhibits wide variability, planning for outbreaks becomes simultaneously more important and more difficult. Discerning and understanding patterns in this variability is thus paramount for effective planning, not least for M. pneumoniae, where severe outcomes are typically preventable if diagnosed and treated in a timely manner.

Since early 2020, the multiannual cycles of M. pneumoniae have been largely disrupted (25, 26). Such polymicrobial effects of the COVID-19-related nonpharmaceutical interventions (NPIs) have been observed across many pathogens (27). The most striking example is perhaps the apparent extinction of the Yamagata strain of influenza B virus (28, 29).

Over 3 y, M. pneumoniae circulation was drastically curtailed in many locations, and Denmark is no exception. A similar multiyear delay in M. pneumoniae resurgence has recently been reported for Europe and Asia (30) as well as the United States (25, 31). Such large-scale perturbations effectively provide a “natural experiment,” demonstrating how the dynamics of M. pneumoniae responds to a prolonged disruption and the associated build-up of susceptible individuals. A recent study of respiratory syncytial virus (RSV) and human metapneumovirus (hMPV) found that including the NPIs and post-NPI rebound could aid model parameter estimation (32), since including data from outside the regular, cyclic regime tends to improve identifiability.

While the origins of the multiannual cycles of M. pneumoniae are generally regarded as unclear (21, 33), a number of different mechanisms of varying complexity have been proposed in the literature. Among the earlier ones was the suggestion by Lind et al. (2) that immunity following infection simply has approximately the same duration as the observed cycle. This explanation, however, has difficulty accounting for how cycle periods can differ somewhat between countries. Without accounting for susceptible dynamics in detail, it is also unclear why the immunity duration should have to match the cycle period. Zhang et al. (34) suggested a combination of strain dynamics, contact network structure, and waning immunity as the explanation, while Omori et al. (35) explored three different possible explanations: school-term forcing, a low-dispersion distribution of immunity duration, and interference between M. pneumoniae serotypes. They find that, of these three explanations, a less variable duration of immunity (compared with the typically assumed exponential waning) best captures the generally observed patterns in that it allows for incidence oscillations. Recently, the authors of ref. 26 used a TSIR model to describe the post-COVID rebound of M. pneumoniae in Denmark as well as four UN regions (Europe, Asia, the Americas, and Oceania) but found that they required at least a 90% reduction in the transmission rate due to nonpharmaceutical interventions in order to fit the data.

Here, we use a mathematical model that incorporates seasonality, vital dynamics (births/deaths), immunity loss, and the effects of NPIs to study recent (2010 to 2025) and historical (1958 to 1995) M. pneumoniae incidence data from Denmark. The simple deterministic model is found to capture the observed dynamics well over a few decades but eventually reverts to annual outbreaks, failing to reproduce the prominent multiannual cycles long-term. We show that the presence of environmental noise can indefinitely sustain the multiyear outbreak patterns. In doing so, we highlight the important role that fluctuations in contact patterns, environmental noise, and demographic processes play when forecasting pathogens with complex dynamics such as M. pneumoniae.

Lessons from Contemporary and Historical Data

We utilize two different datasets as correlates of M. pneumoniae incidence, both from Denmark (Statens Serum Institut, SSI). One has M. pneumoniae-positive PCR tests (and the raw number of tests performed) with a weekly reporting frequency and covers the period 2010 to 2025. This dataset forms the primary basis for parameter estimation, since the test data are highly M. pneumoniae-specific and temporally well-resolved. Secondarily, we use a dataset covering the period 1958 to 1995. Beginning in 1946, patient blood samples collected at hospitals and general practitioners for the diagnosis of Primary Atypical Pneumonia (later identified as M. pneumoniae) were sent to SSI to undergo the cold agglutinin (CA) reaction. This reaction, unfortunately, has only low to moderate sensitivity [ $[eqn]$ 20 to 80% (2, 36)] and moderate specificity [ $[eqn]$ 82 to 92% (2, 36)] for M. pneumoniae. However, beginning in 1958, CA-positive samples were also subjected to M. pneumoniae complement fixation (MpCF), which was found to have higher specificity (37, 38). We thus also treat those tests from the 1958 to 1995 dataset which are simultaneously CA- and MpCF-positive as an incidence proxy on which we base (secondary) parameter estimates.

The duration of immunity following M. pneumoniae infection is known not to be life-long, but the precise duration is unknown, with estimates ranging from a few years to more than a decade (2, 25, 39). For this reason, we opted for an SIRS (susceptible-infected-recovered-susceptible) compartmental model which allows previously infected individuals to become susceptible, and thus potentially reinfected. The model includes seasonal modulation of the transmission rate and equal birth/death rates (which are assumed constant within each data fit, for simplicity). The duration of immunity was treated as a parameter to be estimated during fitting to data.

The influence of NPIs is implemented as a slowly varying multiplicative effect on the transmission rate from March 2020 onward, representing a reduction in overall contact rates, the magnitude of which is an estimated parameter (SI Appendix, Fig. S7).

We estimate a basic reproductive number of $[eqn]$ (95%CI [1.45, 1.82]), in good agreement with ref. 40, a mean duration of immunity of 8 to 9 y and a relative amplitude of seasonality in transmission of $[eqn]$ . All estimated parameters and their uncertainties can be found in Table 1.

With a sinusoidal seasonality function,* the model fits the 2010 to 2025 data well (with the partial exception of the 2017 season in which the apparent incidence was somewhat higher than modeled), see Fig. 1A. The historical 1958 to 1995 data, which is overall of much lower fidelity, is captured more loosely by the model, see Fig. 1B. Interestingly, here we had to allow for a slow change in the mean transmission rate and ascertainment (allowed to change only slightly every 4 y), otherwise the model would tend to settle into an annual pattern.

Model fits to M. pneumoniae incidence data. A seasonal SIRS model is fitted to recent (2010 to 2025) and historical (1958 to 1995) M. pneumoniae testing data from Denmark using Bayesian inference. (A) The seasonal SIRS model explains observed M. pneumoniae dynamics 2010 to 2025. The data are PCR positivity rates. The Inset shows a full cycle (the interval highlighted in blue in the main plot) in the (S,I) configuration space. Markers are 1 wk apart and the color indicates time of year. (B) The same type of model captures key features of observed M. pneumoniae dynamics 1958 to 1995. The fitted data are variance-stabilized numbers of samples which are simultaneously CA and MpCF-positive. The Shallow Bottom panel shows a smoothed version of the time series, to emphasize the multiannual peaks (Data source: Statens Serum Institut and ref. 2).

This leads us to a more general observation: the SIRS model with M. pneumoniae parameters is able to capture the patterns reasonably well over a few decades, but in the absence of perturbations, the modeled dynamics always decays into an annual pattern eventually. This is of course at odds with the ubiquity of complex multiannual cycles in M. pneumoniae (2, 17???–21) and other SIRS-like epidemiological systems such as parvovirus B19 (41???–45). In the next section, we discuss how the time-scales of the M. pneumoniae system may give rise to an approximately 5-y cycle, at least transiently. We then go on to investigate the role that stochasticity plays in sustaining such cycles over decades and in reigniting them even after periods of approximately annual dynamics. In the final part of the paper, we fit a stochastic model to the data and study projections of future M. pneumoniae incidence with and without stochasticity.

The Origins of Intrinsic Cycles.

The seasonally forced SIRS system with vital dynamics contains two comparable oscillatory time-scales: the 1-y seasonal period and the time-scale set by the intrinsic epidemiological dynamics. To ascertain the latter time-scale, the corresponding autonomous system may be studied (that is, the SIRS system with vital dynamics but no seasonal forcing). This system will eventually reach an endemic equilibrium but will perform underdamped oscillations as it approaches this equilibrium. A linearization of the autonomous system around its endemic equilibrium yields an intrinsic oscillation frequency as well as an exponential decay rate (46), given by a combination of the transmission rate $[eqn]$ , immune waning rate $[eqn]$ , recovery rate $[eqn]$ and birth/death rate $[eqn]$ . See Eq. 5 in Methods and Materials for the analytic expression for the intrinsic period. Using the parameter values estimated from the Bayesian fit to Danish data (2010 to 2025), the intrinsic period is found to be $[eqn]$ (95%CI: [4.96, 5.10]). In SI Appendix, Fig. S6, we explore the pairwise parameter dependence of the intrinsic period. The parameter values estimated by fitting to the 2010 to 2025 data are indicated with black crosses. The intrinsic period depends strongly on the transmission rate ( $[eqn]$ ), infectious period ( $[eqn]$ ), and immunity duration ( $[eqn]$ ), while the birth/death rate ( $[eqn]$ ) plays only a minor role. Since the mean duration of immunity and infectious period are primarily a function of intrinsic host–pathogen interactions, we conclude that the primary determinant of periodicity in different countries is likely to be the transmission rate.

Grossman remarked (46, 47) that the SIR equilibrium (w. vital dynamics) tends to be only weakly stable, since the real part of the eigenvalue (the decay rate) is much smaller than the imaginary part (the frequency). We find the same phenomenon for the SIRS system with M. pneumoniae parameters, with the ratio $[eqn]$ (with $[eqn]$ being the relevant eigenvalue of the Jacobian matrix for the linearized system).

A simple Fourier transform of the 2010 to 2020 (pre-COVID) data yields a dominant period of 5.1 y, while a Lomb–Scargle (LS) periodogram (48), which has finer frequency resolution, places the peak at 4.9 y. For the MpCF-based 1958 to 1995 time series, the LS periodogram shows the annual signal as being the strongest, which makes sense given the observation of Lind et al. (2) that the multiannual outbreak pattern was somewhat overtaken by annual dynamics during the late 1970s and early 1980s. However, the strongest multiannual peak for the whole period (1958 to 1995) is at 4.1 y (LS)/4.2 y (Fourier). If the 1958 to 1975 period is viewed in isolation, the strongest peak (eclipsing the annual signal as well) is at 4.6 y (LS)/4.3 y (Fourier). See SI Appendix, Fig. S2 for Fourier and LS spectra of the incidence time series, as well as a wavelet analysis which shows changes in the dominant frequency over time.

In the next two sections, we consider the role that stochasticity may play in sustaining the observed multiannual cycles. Apart from stochasticity being a fact of nature, the introduction of a stochastic model is further motivated by the observation that the fit in Fig. 1A displays low uncertainty compared with the variability in the observed data, suggesting that noise is needed to explain the data.

Demographic Stochasticity

Spurred by the observation that the deterministic model will always eventually decay into annual dynamics, it is natural to ask how a departure from determinism affects the dynamics. Demographic stochasticity is a ubiquitous source of noise in population dynamics (including epidemiological dynamics) (49, 50). This type of stochasticity occurs in finite populations and is fundamentally due to two factors. 1) Individuals being discrete (as opposed to the real numbers $[eqn]$ , $[eqn]$ , and $[eqn]$ of a typical differential equation-based model), and 2) The presence of probability rates such as $[eqn]$ , $[eqn]$ , $[eqn]$ and $[eqn]$ , which implies stochasticity in the timing of birth/death, infection, recovery, and immune waning in finite populations.

We implement demographic stochasticity by considering an SIRS model (with vital dynamics) in which individuals are discrete and may transition between compartments at rates given by Eq. 6 (Materials and Methods). Following the approach of refs. 51?–53, based on the van Kampen system size expansion, we derive an expression for the PSD $[eqn]$ of noise-induced fluctuations in the (nonseasonal) SIRS model with demographic stochasticity:

$$[eqn]$$

Here, $[eqn]$ , $[eqn]$ , $[eqn]$ , and $[eqn]$ are all functions of the parameters $[eqn]$ , $[eqn]$ , $[eqn]$ and $[eqn]$ . Their expressions are given in SI Appendix along with a derivation of Eq. 1. In Fig. 2A, we plot the PSD with parameters from the 2010 to 2025 fit, and note that while the peak is at 5 y, a range of frequencies are excited by demographic noise.

Demographic stochasticity excites multiyear cycles in smaller populations. (A) Theoretical power spectral density (PSD) of multiannual fluctuations excited by demographic noise (without seasonality), and avg. power spectrum of 50 stochastic simulations (with seasonality). See SI Appendix, Fig. S3 for plots of the dependence of the theoretical PSD on individual parameters. (B) Stochastic simulation of a seasonally forced SIRS model parameterized with M. pneumoniae parameters as determined by Bayesian fit to 2010 to 2025 data. In small populations (initial population size N0=20,000), demographic stochasticity leads to noisy multiyear cycles. (C) In intermediate populations (N0=250,000), multiyear cycles are prominent and fairly regular. (D) In larger populations (N0=2×106), annual dynamics dominates as the impact of demographic stochasticity is smaller.

By simulating the seasonally forced stochastic SIRS model using the Doob–Gillespie algorithm (details in Materials and Methods) (54), we find that multiyear cycles may indeed be sustained by demographic noise in sufficiently small populations, suggesting that the system exhibits a stochastic resonance (46). However, the effects of demographic stochasticity are weaker in larger populations, and clear multiannual cycles can only be sustained by this mechanism in populations of less than approximately 2 million individuals, see Fig. 2. Simple demographic stochasticity alone is thus unlikely to be the source of sustained multiannual cycles in M. pneumoniae, at least as far as well-mixed populations are considered. In SI Appendix, we consider a metapopulation version of the SIRS model with demographic stochasticity and show that multiannual cycles may persist in larger, incompletely mixing populations (SI Appendix, Fig. S11).

Environmental Stochasticity

We go on to consider the effects of extrademographic or environmental stochasticity. Random fluctuations in the transmission rate $[eqn]$ will be the primary source of environmental stochasticity considered, but we will also explore the effects of more generic exogenous shocks to the system. Such perturbations can have many different origins, including changes in contact and mobility patterns (such as during holiday periods or due to large events), superspreading of biological and behavioral origin (55??–58), as well as phenotypic variation in the pathogen itself due to mutations. We implement a version of the seasonal SIRS model with stochastic transmission rates based on the method given in ref. 59 (Materials and Methods). This SIRS model with environmental stochasticity is then used for parameter estimation as well as direct simulation. The parameters obtained by fitting this model to 2010 to 2025 data can be found in Table 1, labeled “stoch.”

As shown in Fig. 3, random fluctuations in the transmission rate tend to stimulate approximately 5-y cycles for a wide range of moderate noise levels. When stochasticity is very slight, the annual cycle dominates (e.g., $[eqn]$ in Fig. 3), and when environmental fluctuations are very strong, the epidemic dynamics lacks predictable patterns and becomes dominated by large but rare outbreaks, appearing as a Fourier spectrum with all frequencies represented and no discernible peaks (see, e.g., $[eqn]$ in Fig. 3). It should be noted that even a subdominant multiannual peak in the Fourier spectrum (with, e.g., half of the peak power of the annual signal) leads to a prevalence time series with a clear multiannual component. The averaging process of Fig. 3 leads to broader and flatter multiannual peaks than those found in individual spectra, but clearly shows the range of excited periodicities.

Moderate environmental noise stimulates multiannual cycles, especially approximately 5-y cycles. Average Fourier spectrum across 40,000 realizations (for each noise level indicated in the legend) of the prevalence I(t) of the seasonal SIRS model with parameters from the fit to 2010 to 2025 data from Denmark. The noise source is random fluctuations in the transmission rate β(t), implemented as Gamma-distributed multiplicative noise with mean 1 and variance σ2 (Materials and Methods and ref. 59). 80 y of each simulation time series were analyzed, and a 100 y “transient” was discarded.

While the 5-y cycle is the most typical multiannual cycle at moderate noise levels (given M. pneumoniae parameters), other cycles are excited as well, as shown in Fig. 3, with approximately 3-y as well as 6 to 7 y cycles being somewhat common at moderate-to-high noise levels ( $[eqn]$ ). In SI Appendix, Fig. S4, we show how the addition of continuous environmental noise to an initially purely annual system can rapidly excite the multiannual cycles.

The immediate impact of a singular perturbation (a “kick”) to the system (60) is explored in Fig. 4. The perturbation is implemented as either a discrete change in the infected population (Fig. 4A–C) or the susceptible population (Fig. 4D). Since just a single “kick” is given, the system will only transiently settle into a multiannual cycle and we thus measure the strength of the 5-y cycle during the 20 y immediately following the perturbation. In Fig. 4A, we show the trajectory of the system in the $[eqn]$ configuration space immediately before and after a single kick to the system. In this case, the perturbation consists of a sudden reduction in prevalence $[eqn]$ , but a sudden change in the susceptible population can have a similar effect (Fig. 4D). The annual attractor, where the system resides prior to the perturbation, appears as the central ellipse, while the resulting (transient) 5-y cycle encloses it. The rates at which $[eqn]$ and $[eqn]$ change is reflected in the spacing between the colored markers which are placed 1 wk apart. The colors of the markers themselves indicate the time of year. As is evident, the 5-y cycle can be quite complex. The perturbation that initiated the 5-y cycle is visible as the approximately vertical line connecting the annual and multiannual cycles, and the state of the system immediately prior to the perturbation is indicated by the red cross. Fig. 4B shows that the timing of the perturbation matters, with a kick during the low-transmission part of the season most likely to excite the multiannual cycle. In other words, a kick of a certain size (whether in the positive or negative direction) has a larger effect if it occurs when circulation is already low. However, Fig. 4C reveals that this seasonal dependence is completely accounted for by the relative magnitude of the perturbation, $[eqn]$ , i.e., the change in incidence divided by the current incidence at that time. If plotted as a function of this relative magnitude, the strength of the excited cycle shows only a very faint seasonal dependence. Fig. 4D is analogous to Fig. 4B, except that the perturbation is applied to the susceptible fraction $[eqn]$ . Here, seasonal dependence is fairly minimal, since the susceptible fraction fluctuates much less during a season (relative to its mean value) than the prevalence $[eqn]$ does.

One-off exogeneous perturbations (“kicks”) transiently excite multiannual cycles. (A) System trajectory depicted in (S,I) configuration space before and after a perturbation. The marker color denotes the time of year, and markers are 1 wk apart. The approximately vertical line starting from the red cross shows an exogenous “kick” (a discrete change in prevalence, in this case) which causes the system to transition from the annual (thick yellow line) to 5-y cycle (thin purple line). (B) After a discrete change ΔI in the prevalence, the system may temporarily settle into a 5-y cycle, as in A. The plot shows how the size of the perturbation and the time of year affect the strength of the resulting cycle. The color bar indicating the power of the 5-y cycle in the Fourier spectrum (relative to the highest peak) is shared between B, C, and D. Spectra were based on the first 20 y post-perturbation. The red cross indicates the perturbation implemented in A. (C) If the perturbation ΔI is proportional to the preperturbation prevalence I0, the prominence of the 5-y cycle ceases to depend on time of year. The red cross indicates the perturbation implemented in A. (D) Effects of perturbations in S(t), rather than the prevalence.

The mechanism of sustained cycling demonstrated here is fundamentally that of a quasicycle, a mechanism of population cycling proposed in ref. 61. It arises from an otherwise stable, underdamped system being subjected to noise which excites an intrinsic cycle that would otherwise decay. In other words, noise continually disrupts the system’s deterministic approach toward equilibrium. In SI Appendix, Fig. S8, we demonstrate that the same phenomenon occurs in the absence of seasonal forcing. In the ecological literature, there is a fundamental distinction between phase-remembering and phase-forgetting quasicycles (62?–64). In the former case, the cycle retains phase information over time, such that even temporally well-separated fluctuations exhibit a substantial correlation, typically due to an external periodic forcing. In the phase-forgetting case, this long-term synchronization is eventually lost. In SI Appendix, Fig. S9, we measure the autocorrelation of the prevalence signal $[eqn]$ in the fitted model and find that the multiannual cycles are of the phase-forgetting kind. If one considers only the autocorrelation of the raw prevalence time series (SI Appendix, Fig. S9A), one could be led to believe that the system is phase-remembering, since a substantial autocorrelation remains at large lags $[eqn]$ . However, as we show in SI Appendix, Fig. S9B, this residual autocorrelation is solely due to the annual modulation, and any phase information in the multiannual cycle is lost over time. We find that the autocorrelation of the multiannual cycle decays approximately exponentially with a half-life of 6 to 7 y.

The finding that environmental stochasticity can sustain the multiannual cycles of the system is significant because, unlike demographic stochasticity, this mechanism is independent of population sizes.

Aside from explaining the existence of multiannual cycles, the stochastic model also explains the temporary disappearance of the multiannual cycles between the late 1970s and mid-1980s, observed by Lind et al. (2). As shown in the wavelet analysis of model results in SI Appendix, Fig. S10, the stochastic model exhibits mode-hopping, i.e., spontaneous and temporary reversions to the annual attractor, similar to what can be seen in the wavelet transformed empirical data of SI Appendix, Fig. S2. The model fit of Fig. 1B also partially captures this phenomenon, with small infrequent changes in the transmission rate in lieu of continuous noise.

Beyond the Post-COVID Rebound of M. pneumoniae

Although the deterministic seasonal SIRS model captures the 2010 to 2025 dynamics well, predicting or forecasting M. pneumoniae incidence may be challenging. We know that long-term predictions using the deterministic model are limited by its tendency to revert to annual dynamics. We have also established that phase information in the multiannual cycles is lost over time when the system is subject to stochasticity, which the real-world system invariably is. Both of these findings suggest that a deterministic model is fundamentally limited in terms of long-term predictive power. In the short term, however, a deterministic model may still be appropriate for forecasting. Fig. 5A shows the trajectory obtained by continuing the deterministic seasonal SIRS model for 10 y beyond the 2010 to 2025 data (with parameters fitted to said data). However, the low uncertainty of the prediction offered by the deterministic model may be an underestimation. Fig. 5B shows the fit of the seasonal SIRS model with environmental stochasticity. While the 2010 to 2025 data are still captured well, prediction beyond a year or two is associated with high uncertainty. The noise level was estimated at $[eqn]$ using this model, a degree of stochasticity capable of sustaining pronounced multiannual cycles, as we saw in Fig. 3. We note that this level of stochasticity is substantial, in that $[eqn]$ is comparable in magnitude to the estimated amplitude $[eqn]$ of seasonality at approx. $[eqn]$ .

Prediction is hard. (A) Prediction using the deterministic seasonal SIRS model, parameters fitted to 2010 to 2025 data. (B) Prediction using the seasonal SIRS model with environmental stochasticity (noise in transmission, σ=0.12, fitted value). Noise can significantly affect long- and medium-term predictions, but moderate to high M. pneumoniae activity is predicted in the 2025/2026 season. The 95% CI is shown in light purple, and the 85, 75, and 50% intervals are shown in successively darker shades.

The sensitivity to noise can be measured via the Lyapunov exponents of the system. A Lyapunov exponent $[eqn]$ measures sensitivity to initial conditions—this also includes the tendency of the system to exponentially amplify ( $[eqn]$ ) or buffer ( $[eqn]$ ) perturbations/noise (63). SI Appendix Fig. S1 explores the local/instantaneous Lyapunov exponents (LLEs) (65) of the deterministic as well as stochastic model and shows that, while perturbations tend to be transiently amplified, all LLEs eventually become negative, indicating the eventual decay of the perturbation.

While sensitivity to stochastic effects may complicate medium-term predictions (note the spread in predictions in Fig. 5B from 2027 onward), some short-term features are relatively robust even at this noise level. Predictions agree that the 2025 to 2026 season is likely to be a moderate-to-high incidence year, exemplifying how deterministic models yield valuable short-term forecasts even in noise-driven epidemiological systems.

While long-term predictions using the deterministic model are limited by its decay into annual dynamics, the presence of stochasticity limits extended prediction in general, not only by causing variability in the incidence during any given year but also by enabling mode-hopping, i.e., transitions between annual and multiannual cycles (50, 66), which appear in the empirical data (SI Appendix, Fig. S2) as well as in model simulations (SI Appendix, Fig. S10). That the predictions for the years immediately following the large NPI-induced perturbation agree well between the stochastic and deterministic models fits well with the general observation that restorative forces are stronger far from equilibrium, leading the stochastic model to (temporarily) behave more like the deterministic model (46).

Limitations

In fitting the 2010 to 2025 time series, we used the PCR positivity rate as an incidence proxy, while we used the number of positive tests for the 1958 to 1995 data (total number of tests were not on record). Neither of these measures are perfect incidence proxies, and they suffer in distinct ways from biases introduced by varying testing intensities. If testing intensity is increased, this leads to an increase in the number of positive tests, while it leads to a decrease in the test positivity rate (since the denominator—the total number of tests—is increased). It may be possible to enhance fitting further by accounting for fluctuations in testing intensity when such data is available, something we have not attempted to do, beyond simple variance stabilization. Furthermore, the fidelity of the 1958 to 1995 data series is affected by both the (somewhat imperfect) test characteristics of the cold agglutination and complement fixation tests [since the samples considered positive are those which were positive in both of these tests (2)]. We have not taken age structure into account (67), nor have we considered detailed spatial effects which might enhance the impact of demographic stochasticity (68).

Discussion and Conclusion

The origins of nontrivial cycles in ecological systems and the extent to which stochasticity plays a role in shaping them are long-standing questions in ecology, including epidemiology (60, 69, 70). While the concept of quasicycles—periodic population fluctuations resulting from the interaction of noise with decaying intrinsic cycles—was proposed by Nisbet and Gurney (71), the theory has evolved since. Significant analytic advances for SIR-type systems (52) and predator–prey systems (63, 72) were made two decades ago, and quasicycles were proposed as a model of seasonal fluctuations in influenza incidence around the same time (73), albeit without empirical verification. In ref. 53, demographic stochasticity was considered as a potential driver of recurrent outbreaks of avian influenza in North American wild waterfowl and the observed cycles were quantified through wavelet analysis. The authors used a stochastic SIR-type model, extended to include environmental transmission (but not environmental stochasticity), and found that it could explain the occurrence of multiannual outbreak cycles, with the period being largely determined by the intensity of environmental transmission.

The effects of demographic stochasticity have been most widely considered, while environmental noise features less prominently in the epidemiological literature, and studies that consider both are rare. Pertussis (whooping cough) is another example of a disease where multiannual cycles have been studied in detail. In ref. 74, the authors used wavelet analysis to study the multiannual cycles of pertussis across 12 countries and noted that, while a 3 to 4 y periodicity was typically found [a result echoed by Metcalf et al. (75)], it is transient in nature. The multiannual cycles of pertussis were previously shown to likely originate from the system’s response to dynamical noise (76), and a study found that stochasticity is essential to even qualitatively capture the cycles of pertussis (60).

For M. pneumoniae, we find that noise-perturbed intrinsic dynamics is sufficient to explain the complex cycles, offering a parsimonious model for the dynamics. Importantly, the model also explains the temporary disappearance of the multiannual cycles observed in the Danish data from the late 1970s until the mid-1980s as a consequence of mode-hopping, i.e., the system reverting to the annual attractor due to a perturbation. This temporary disappearance of the multiannual pattern was previously noted by Lind et al. (2), who proposed changes in daycare attendance as a possible driver.

Earn et al. (77) demonstrated the existence of a rich landscape of stable cycles in the related, seasonally driven SEIR model at high transmission rates (basic reproductive numbers in excess of approx. 4). In the current case of M. pneumoniae however, multiannual cycles are unstable and require stochasticity or perturbations to be sustained. A similar systematic mapping of possible cycles in the seasonal SIRS model—taking stochasticity into account and including subharmonic resonances (71, 78)—would be highly useful.

With the COVID-19 NPIs phased out, another major—if less abrupt—change in M. pneumoniae epidemiology is on the horizon. Macrolide resistance has been growing rapidly in recent years, albeit with large geographical differences (79). It is currently unclear how this will affect epidemiological dynamics, but some studies have found macrolide-resistant strains to be associated with increased virulence and longer disease duration (80, 81), leading us to speculate that average transmissibility could be influenced as well. Extending the current model to incorporate observational data on macrolide resistance would be highly worthwhile, and holds the potential to estimate any effects of macrolide-resistant strains on M. pneumoniae dynamics.

Since the main objective of this study was to determine the mechanism of population cycling in M. pneumoniae, we opted for a fully mechanistic model. This has the further advantage of enabling straight-forward prediction of the overall dynamics in parameter ranges that differ from those obtained from the Danish data. Some parameters are likely intrinsic to the pathogen–host interaction (rates of waning and recovery), while others (transmission rates and birth/death rates) are likely to vary geographically, leading to different dominant periods according to Eq. 5. For adaptive forecasting, however, the formulation of a hybrid semimechanistic machine learning model would likely be advantageous. In a recent study (82), a machine learning model was shown to perform well when tasked with short- and medium-term forecasting of historical measles data. The authors note that “[t]he primary focus of mechanistic models has been understanding and characterizing the natural history of transmission [...], statistical and machine learning techniques in infectious disease modeling have primarily focused on improving forecasting accuracy without the explicit aim of inferring transmission dynamics.” While beyond the scope of this study, we encourage a synthesis of the two approaches for M. pneumoniae forecasting.

In this study, we have considered both demographic and environmental stochasticity as potential drivers of multiannual outbreak cycles. However, we have not attempted any detailed attribution of stochastic effects to behavioral or demographic sources. A future study could conceivably pinpoint the sources of stochasticity that most strongly shape M. pneumoniae dynamics by extending the current work to incorporate (time-varying) population structure (in terms of age and contact patterns) and mobility data.

In summary, we have offered a parsimonious explanation for the multiannual cycles of M. pneumoniae as a type of quasicycle, and provided a computational and analytic framework for similar analyses of other pathogens and locations. There is a need for comparative studies of the factors governing the diversity of complex cycles observed globally for pathogens such as M. pneumoniae and parvovirus B19, across countries and pathogens. Such comparisons may also enable coestimation of environmental sources of stochasticity, and the framework presented here is well suited to facilitate this.

Materials and Methods

The Seasonal SIRS Model.

The seasonally forced SIRS model with vital dynamics and constant population size is given by

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

where the instantaneous transmission rate $[eqn]$ is a nonnegative and periodic function with a period of 1 y, i.e., $[eqn]$ . The remaining parameters are the birth/death rate $[eqn]$ , recovery rate $[eqn]$ and immunity waning rate $[eqn]$ . Noise in the transmission rate $[eqn]$ is implemented as described in ref. 59; when numerically integrating the above set of differential equations, we draw multiplicative noise from a Gamma distribution with shape $[eqn]$ and scale $[eqn]$ . The parameter $[eqn]$ is the noise level referred to in Fig. 3.

Linearization of the autonomous system.

The analogous autonomous system to Eqs. 2–4 is obtained by substituting a constant transmission rate $[eqn]$ in place of $[eqn]$ . If perturbed away from equilibrium, $[eqn]$ performs an underdamped oscillation while approaching the equilibrium configuration $[eqn]$ , where:

$$[eqn]$$

Linearizing the system around this equilibrium and applying a small perturbation, the decay rate, and characteristic frequency are given by the real and imaginary parts of the eigenvalues of the Jacobian matrix:

$$[eqn]$$

The angular frequency of the intrinsic oscillation is then given by $[eqn]$ , where $[eqn]$ and $[eqn]$ . Since $[eqn]$ is much greater than $[eqn]$ and $[eqn]$ , it holds that $[eqn]$ , although we choose to keep $[eqn]$ in our calculations. The corresponding period is then given by

$$[eqn]$$

With the parameters fitted to the Danish 2010 to 2025 time series, the eigenvalues are $[eqn]$ and $[eqn]$ . Since all real parts are strictly negative, the equilibrium is asymptotically stable.

Bayesian Fitting of a Deterministic Model to 2010 to 2025 Data.

We fit the seasonal SIRS model as described above using the probabilistic programming language Stan (83, 84), with a sinusoidal seasonally forced transmission rate $[eqn]$ , to positivity rates given in the 2010 to 2025 dataset. $[eqn]$ is the NPI multiplier which is fixed to 1 until March 2020 and otherwise estimated. The priors used were

• $[eqn]$ where $[eqn]$ indicates a Gaussian (normal) distribution with mean $[eqn]$ and SD $[eqn]$ , while $[eqn]$ is the recovery time (2.5 wk).
• $[eqn]$
• $[eqn]$
• $[eqn]$
• $[eqn]$ such that $[eqn]$ (the stick-breaking method, ensuring $[eqn]$ )
• $[eqn]$ where $[eqn]$ indicates the uniform distribution with support on the interval $[eqn]$
• $[eqn]$
• $[eqn]$
• $[eqn]$
• $[eqn]$

Here $[eqn]$ and $[eqn]$ refer to the measured incidence proxy (positivity rate, arbitrarily rescaled) and the modeled incidence, respectively. Note that for the parameters $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ and the NPI multipliers, only the part of the normal distribution with positive support is used, since the parameters in question are necessarily nonnegative. The fit to the 1958 to 1995 MpCF-based data follows similar principles and is detailed in the appendix. For the purposes of fitting, the model (Eqs. 2–4) is implemented using the efficient discretization scheme described in ref. 25, with a time-step of $[eqn]$ . This discretization was used for Figs. 1 and 3–5 and SI Appendix, Figs. S4, S5, S7, S8, and S10, whereas the calculations of SI Appendix, Figs. S1 and S9 were performed using standard Euler integration of Eqs. 2–4. Convergence was assessed based on the absence of diagnostic warnings from CmdStanPy (84), including the absence of divergent chains, no exceedances of tree-depth, sufficient E-BFMI ( $[eqn]$ 0.3), and satisfactory $[eqn]$ ( $[eqn]$ 1.05) and effective sample size (bulk $[eqn]$ ) for all parameters.

Bayesian fitting of a stochastic model.

We fit a seasonal SIRS model with environmental stochasticity in the form of Gamma-distributed multiplicative noise in the transmission rate, as used in Fig. 5. The model set-up, using Stan (83), closely follows that of the deterministic model. A new fitted parameter, the noise level $[eqn]$ , is introduced which controls the magnitude of noise in $[eqn]$ , with prior $[eqn]$ (positive-support part of normal distribution only). The individual perturbations are then implemented by letting $[eqn]$ , where $[eqn]$ with the notation $[eqn]$ referring to a Gamma distribution with shape $[eqn]$ and scale $[eqn]$ . For computational efficiency, we only draw a new value of $[eqn]$ every other week.

Doob–Gillespie Simulations.

The stochastic simulations of Fig. 2 are based on the following system of transition rates:

$$[eqn]$$

with the time-dependent transmission rate given by $[eqn]$ . Note that here, $[eqn]$ , $[eqn]$ , and $[eqn]$ denote the populations of susceptible, infected, and recovered individuals (as integers), not the population fractions. Denoting the above rates by $[eqn]$ , $[eqn]$ , the Doob-Gillespie algorithm (54) thus prescribes iteratively propagating the state of the system by drawing a uniformly distributed random number $[eqn]$ , determining the time to next event $[eqn]$ and the nature of the event by drawing another random uniform number $[eqn]$ and determining the smallest $[eqn]$ for which $[eqn]$ . Event (or “reaction”) $[eqn]$ is then performed and the time is incremented by $[eqn]$ . The procedure is then reiterated until the desired final time is reached. The metapopulation simulations of SI Appendix, Fig. S11 use the same algorithm, but with modified rates. Details are given in SI Appendix.

Supplementary Material

Appendix 01 (PDF)