## Reconstruction of disease transmission rates: applications to measles, dengue, and influenza

Alexander Lange

Institute of Thermodynamics and Thermal Process Engineering, University of Stuttgart, Germany

E-mail: lange@itt.uni-stuttgart.de

## Abstract

Transmission rates are key in understanding the spread of infectious diseases. Using the framework of compartmental models, we introduce a simple method that enables us to reconstruct time series of transmission rates directly from incidence or disease-related mortality data. The reconstruction exploits differential equations, which model the time evolution of infective stages and strains. Being sensitive to initial values, the method produces asymptotically correct solutions. The computations are fast, with time complexity being quadratic. We apply the reconstruction to data of measles (England and Wales, 1948-67), dengue (Thailand, 1982-99), and influenza (U.S., 1910-27). The Measles example offers comparison with earlier work. Here we re-investigate reporting corrections, include and exclude demographic information. The dengue example deals with the failure of vector-control measures in reducing dengue hemorrhagic fever (DHF) in Thailand. Two competing mechanisms have been held responsible: strain interaction and demographic transitions. Our reconstruction reveals that both explanations are possible, showing that the increase in DHF cases is consistent with decreasing transmission rates resulting from reduced vector counts. The flu example focuses on the 1918/19 pandemic, examining the transmission rate evolution for an invading strain. Our analysis indicates that the pandemic strain could have circulated in the population for many months before the pandemic was initiated by an event of highly increased transmission.

## Keywords

infectious disease modeling; disease transmission; differential equation models

## 1 Introduction

Essential for modeling an infectious disease epidemics is the knowledge of the transmission rates - the rates at which susceptibles become infected by contagious individuals [1, 2]. These rates are determined by the contact behavior of the involved hosts as well as the risk of transmission during contact, they are specific to

the pathogen and its transmission route [3]. Transmission rates fluctuate and often systematically change over time as we will examine in this paper for data of three different infections.

Combined with basic medical and demographic information, transmission rates are the most natural reference in predicting the time evolution of disease prevalence. Public health policies rely on such predictions, referring to the control of endemic infections as well as to the design of measures against emerging diseases, pandemics, or bio-terroristic threats [4-6]. In practice, however, there is no straightforward way these rates are obtained from epidemiological data. The methods known are rather complicated and/or require lots of computing power [7-13], others are less developed at this stage [14-17]. Addressing these issues, we will present a very simple method for reconstructing transmission rates from incidence or mortality data and illustrate its wide-ranging applicability.

There are lots of computing and modeling strategies in infectious disease epidemiology, including probabilistic simulations [9, 10], networks [18, 19], and compartmental models [1, 20], to mention a few. The conceptually simplest and mathematically most tractable ones fall into the last category. In compartmental models, transmissions between hosts are assumed to happen at random, incorporated through versions of the mass-action principle [21]. Formulated in terms of differential equations, these models keep track of expectation values that represent susceptible, infective, and recovered individuals (SIR), as well as other relevant compartments of the host population [20]. The transmission mechanism is implemented by a (timedependent) coefficient, the transmission rate β . Normalized with respect to an infection-specific removal rate λ , the resulting unit-free number R 0 = β/λ defines the basic reproduction of infections [22] - a parameter well-known because of its importance for disease control and pathogen evolution [1,3].

Whereas disease incidence and mortality define quantities of major interest to public health [23], transmission rates offer the more natural parameter in characterizing how a disease system is changing over time. Due to the internal dynamics, disease prevalence often behaves strangely, for example, shows biennial time patterns even if infection is forced by an annual period (known for measles [24,25]) or increases even if risk factors are decreasing (reported for dengue [26]). Prevalence is sensitive to initial conditions and time resolution, it even shows chaotic behavior [27]. The kind of information that can be captured by transmission rates will be illustrated here for three examples: measles during the pre-vaccination era, the currently re-emerging dengue epidemic, and the Spanish flu pandemic.

The measles example is used to develop the method. It offers a test for the reconstruction results and its systematic errors. Data from England and Wales have been studied extensively, and much is known about the temporal pattern of transmissions [8,10,24].

Epidemiological data are tainted with reporting errors, including missing or false diagnoses. To compensate these errors the concept of the reporting proportion has been introduced (e.g., [8, 24]). In addition to the reporting proportion known from the literature, we introduce a second one, enabling the definition of an effectively constant population, useful when demographic information is insufficient.

Besides focusing on technical questions - regarding parameters, errors and extensibility to new compartments - we study epidemiological questions and illustrate how the reconstruction can be applied to obtain conceptual results. Here dengue represents a generic example. Increasing cases of dengue hemorrhagic fever (DHF) in Thailand need to be explained based on decreasing transmission rates [26]. Our method allows us to re-evaluate former explanations based on strain interaction [28] and demographic transitions [29]. The third application is the influenza pandemic 1918/19, which we study based on mortality data from the U.S. Here we use the reconstruction method to investigate irregularities prior to the pandemic peak, referred to as herald waves [30].

## 2 Methods

In this section we develop a procedure for reconstructing transmission rates from time series data. Measles infections are used as the main example. For other diseases the methodology requires adjustments, as we illustrate for two other infections in Results.

## 2.1 Basic model

One of the simplest settings for an infectious disease epidemic is given by a system of three first order differential equations,

$$S ^ { \prime } = \Gamma [ S, I, R ] - \beta \, \Sigma [ S, I ]$$

$$I ^ { \prime } & = \beta \, \Sigma [ S, I ] - \lambda I$$

$$R ^ { \prime } & = \gamma I - \rho R \,.$$

These equations determine the time evolution of compartments, S, I , and R , which in respective order quantify the average numbers of susceptible, infective , and recovered individuals [20]. Their sum,

$$N = S + I + R \,,$$

represents the size of the host population . Usually, this number differs from population data published in demographic reports (e.g., [31]); it only includes parts of the population. Most obviously this applies to sexually transmitted infections where the host population represents the sexually active part.

Utilized as balance equation, N ′ = ν -µN , the derivative of (2) determines the rate,

$$\Gamma = \nu - \mu S + \delta R + \Theta [ I, R ] \,,$$

at which individuals enter, leave, and (possibly) re-enter the susceptible compartment - through birth ν , death µS , and (e.g., for flu-like infections) after recovery δR . The remaining term,

$$\Theta = ( \lambda - \gamma - \mu ) I + ( \rho - \delta - \mu ) R \approx 0 \,,$$

obtained by adding the three equations (1) and comparing the result to (3), incorporates disease-induced mortality rates specific to the compartments I and R . For many epidemics these rates can be neglected. The recovery rate and the decay of immunity are then given by γ = λ -µ and δ = ρ -µ , respectively.

↦

The number of transmissions is determined by the rate at which susceptibles become infected. This rate, β Σ, which defines removal in (1a) and incidence in (1b), involves a functional and a coefficient. The functional, ( S, I ) → Σ[ S, I ], represents the contacts between infected and susceptible individuals; usually it is modeled via mass-action [21],

$$\Sigma = S I / N \,.$$

Refined models contain fractional powers of S and I [33]. The coefficient, β , defines the transmission rate - the parameter we intend to reconstruct from time series data.

## 2.2 Beta reconstruction

The starting point is a version of (1b),

$$I ^ { \prime } = W - M / q \,,$$

which models the time evolution of disease prevalence in terms of quantities that are observable and thus stored in epidemiological data bases (e.g., [23]). W , given by W = β Σ, represents incidences of infection (number of new infections per time unit). M , given by M = qλI , represents incidences related to a focal condition that infected individuals are affected by; the parameter q models the associated likelihood. If, for example, the infection causes deaths then M can be utilized to model the induced mortality. Here the case-fatality ratio c determines the likelihood, q = 1 c λ -µ λ ≈ 1 c . Another example for such a condition is being recovered . Due to (1b) and (1c), recovery from infection happens with likelihood q = γ/λ at rate M = γI .

By comparing the incidence terms in (6) and (1b) we obtain a formula for the transmission rate,

$$\beta ( t ) = \frac { W ( t ) } { \Sigma [ S ( t ), I ( t ) ] } \,.$$

Incidences W are observed but contacts Σ need to be calculated. To determine S t ( ) and I ( ) we distinguish t two cases, proposing that either incidence W t ( ) or mortality M t ( ) is given. For now, we assume that these data are smooth functions over time.

First case: W t ( ) is given; I ( ) t and S t ( ) need to be calculated. The time-evolution of the infectives is reconstructed by solving a first order linear differential equation,

$$I ^ { \prime } + \lambda I = W \,,$$

which combines (1b) and (7). The solution (see, e.g., [34]),

$$I ( t ) = I ( 0 ) e ^ { - \lambda t } + \int _ { 0 } ^ { t } e ^ { - \lambda ( t - r ) } W ( r ) \, d r \,,$$

involves a convolution integral and an initial value, which without additional information we approximate by the steady state solution, I (0) = W (0) /λ . The time-evolution of the susceptibles is reconstructed based on the recovered individuals (1c). This requires to solve the same type of equation once again. Relying on the balance equation (2), the result reads

$$S ( t ) = N ( t ) - I ( t ) - R ( t ) \,,$$

where similar to (9a),

$$R ( t ) = \left ( N ( 0 ) - I ( 0 ) - S ( 0 ) \right ) e ^ { - \rho t } + \gamma \int _ { 0 } ^ { t } e ^ { - \rho ( t - r ) } I ( r ) \, d r \,.$$

The initial values are approximated by steady state solutions of (1). That is, I (0) = W (0) /λ and, provided mass-action (5) applies, S (0) = N (0) / R 0 , where R 0 = β (0) /λ defines the basic reproduction at t = 0. This number can be estimated by the ratio of 'births' (expressed by the mean age at attack [35], a) and 'deaths' (mean life expectancy, 1 /µ ), R 0 = 1 /a / µ .

Second case: M t ( ) is given; I ( ) t and W t ( ) need to be calculated. Prevalence follows by I = M/ qλ ( ), incidence W I [ ] is obtained from (8) by differentiation. The reconstruction of β t ( ) is now reduced to the first case.

Figure 1: Data and transmission rates. The red curve in Panel A shows the weekly incidences of measles infections in England and Wales [23]. Estimates of the corresponding disease prevalence (yellow) are computed according to (9a), assuming an infectious period of 9 days. In Panel B, the reconstruction of β t ( ) is shown for constant (black) and time-dependent mortality (green). The constant mortality rate is estimated via the mean lifetime (1 /µ = 75years, cf. Fig. 2D). Numerical differences between the reconstructed betas (blue) are very small, plotted in the lower part of the diagram. Time resolution has been increased from weeks to days.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0640v1-figure-c1906284-56b2-462b-8b51-089b40fefc99.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0640v1-figure-c35febe0-7885-49d1-a9be-616d66a9cdd2.jpg)

The reconstruction is robust with respect to the initial values. As (9a) and (9c) reveal, the influence of the initial values vanishes exponentially with rates λ and ρ , respectively. That is, for large time ( t glyph[greatermuch] 1 /λ and glyph[greatermuch] 1 /ρ ), I and R are determined exclusively by the convolution integrals.

In (7), contacts Σ are not restricted to mass-action (5). After specifying the functional Σ[ S, I ], one may derive an expression for β t ( ) that only depends on W t ( ) or M t ( ) and on an initial value (cf. [14-16]).

## 2.3 Including demographics

So far we have assumed that the mortality rate µ and the removal rates, λ and ρ , do not change over time. The reconstruction formulas with time-dependent rates are very similar to the ones derived above, only the exponents of (9a) and (9c) are replaced by integrals [34]. Applying our initial conditions, we obtain

$$I ( t ) = \frac { W ( 0 ) } { \lambda ( 0 ) } e ^ { - \int _ { s } ^ { t } \lambda ( s ) } d s + \int _ { 0 } ^ { t } e ^ { - \int _ { r } ^ { t } \lambda ( s ) d s } W ( r ) \, d r \,, \quad \,.$$

$$R ( t ) = \left ( \stackrel { \cdot \cdot } { N } ( 0 ) \left ( 1 - \frac { 1 } { \mathcal { R } _ { 0 } } \right ) - \stackrel { \cdot \cdot } { \frac { W ( 0 ) } { \lambda ( 0 ) } } \right ) e ^ { - \int _ { 0 } ^ { t } \rho ( s ) } d s + \gamma \int _ { 0 } ^ { t } e ^ { - \int _ { \tau } ^ { t } \rho ( s ) d s } I ( r ) \, d r \,.$$

For measles, however, the improvement of the reconstruction resulting from these refined formulas is virtually invisible (Fig. 1B).

Unlike mortality, the birth rate ν (= N ′ + µN ) and the population size N have not been regarded as constant over time. In contrast to the previous result, the impact of a changing population size (discussed below) is clearly visible (Fig. 2C).

## 2.4 Scaling invariance

Provided contacts are modeled by mass-action (5), our setting is invariant with respect to systematic errors glyph[epsilon1] involved in the data. That is, the basic equations (1) can be interpreted as a model for only a proportion

of incidences W eff = ηW , susceptibles S eff = ηS , infectives I eff = ηI , and recovered individuals R eff = ηR in an effective population N eff = ηN , where η = 1 -glyph[epsilon1] defines the reporting proportion with respect to the actual but unknown values. This is easily confirmed multiplying each additive terms in (1) by η .

The transmission rate obtained for such an effective description is the same as for the actual values, β eff = β . This can be seen by expanding the basic equation (7),

$$\beta = \frac { N \ W } { S \ I } = \frac { \eta N \ \eta W } { \eta S \ \eta I } = \beta _ { \text{eff} } \, ;$$

the additional ' η ' cancel.

The reporting proportion is obtained from (1a) - by substitution of (3) and (7), S ′ = ν -µS + δR + Θ -W , by subsequent integration over time (up to some large t , so that the left-hand side becomes small in comparison with the terms on the right-hand side), 0 ≈ ∫ t 0 S ′ ( r ) dr = ∫ t 0 ( ν r ( ) -··· -W r ( )) dr , and by some infection-specific simplifications. For measles (i.e., δ ≈ 0, Θ ≈ 0, R glyph[greatermuch] 0 1), the integral simplifies to 0 ≈ ∫ t 0 ( ν r ( ) -W eff ( r /η ) ) dr so that the ratio between cumulated incidences and births (cf. [7]),

$$\eta _ { \nu } ( t ) & = \frac { \int _ { 0 } ^ { t } W _ { \text{eff} } ( r ) \, d r } { \int _ { 0 } ^ { t } \nu ( r ) \, d r } \,,$$

as well as moving averages η ν ( ) t of this ratio, yield estimates of the reporting proportion η . In numerical simulations we utilize the harmonic mean to calculate η . If, for example, the average is taken over the whole time period then η ν = n /∑ i ≤ n ν t ( i ) /W eff ( t i ), where n is the number of data points t i .

For measles, the scaling behavior is of particular importance. For the pre-vaccination data from England and Wales the reporting proportion is known to be significantly smaller than one [7, 24]. If we therefore apply the reconstruction (9) to raw data, the method does not work.

## 2.5 An effective population of constant size

Whenever our knowledge about the demography is insufficient, we propose that the size of the population is constant over time, N t ( ) = N . Even for simplicity we could make such an assumption. However, it is not trivial to estimate an appropriate constant effective size. Scaling invariance together with (12a) does not work. For our measles example, the resulting size of η ν N ≈ 0 5 . × 46 × 10 6 = 23 × 10 6 is far too low.

If the population size is chosen too low, the reconstruction produces unrealistically high beta values and eventually collapses (as happening with N eff = 23 × 10 , cf. 6 orange curve in Fig. 2C). This is not surprising when looking at the reconstructed susceptibles S eff = N eff -I eff -R eff in the denominator of (11). By lowering N eff the number of susceptibles S eff approaches zero and eventually turns negative. If, on the other hand, the population size is chosen too large, the reconstructed beta decreases below reasonable values (cf. green curve in Fig. 2C) and, as (11) reveals, tends towards the trivial solution, β = W /I eff eff ≈ λ . The sketched behavior is quite sensitive and can therefore be used to narrow down an appropriate size for the effective population. Assuming that β t ( ) ≈ const, we obtain N eff ≈ 33 × 10 . 6 This number is much smaller than the number of inhabitants (46 × 10 ) in England and Wales at the time [31]. 6

When proposing a constant population size a straightforward application of (12a) is problematic. First , Equation (12a) contains the birth rate ν but the reconstruction equations (9) rely on the death rate µ . Second , A constant population requires that N ′ = 0 (i.e., ν = µN ), even if in reality the size changes. These inconsistencies are solved, however, when replacing ν in (12a) by µN . The effective population is then given

Figure 2: Parameters. For each of the three infections studied in the article, Panel A shows the cumulative distribution function (CDF) of the lower bound λ ∗ of the removal rate (13) with respect to the data points. The mean infectious periods used in our numerical simulations are indicated by dashed lines. The Panels B-F refer to measles. In B, the two kinds of reporting proportions - based on natality (green) and on mortality (blue) - are plotted over time. The red line indicates the reporting proportion that corresponds to a constant effective population (of size 32 8 . × 10 ). 6 In Panel C, the reconstruction results (2-year moving averages) are shown for the population sizes plotted in Panel D: the actual size (black), a slowly varying size (blue), and three constant sizes (green, red, orange). In respective order, the reconstructions are performed based on reporting proportions η ν ( ), t η µ ( ), and t η µ ; the mortality rate (dashed) is plotted in Panel D too. In the two remaining panels, E and F, the mean infectious period is varied: 6 days (orange, light gray), 9 days (red, gray), 12 days (green, dark gray). Panel E shows reconstruction results (2-year moving averages) for a constant (colored) and a variable population size (black; reporting errors are indicated gray). Panel F shows reconstruction results (1-month moving averages) for a variable population size; amplitudes are higher and more regular for longer infectious periods.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0640v1-figure-38c2a27f-fa4b-447c-9ba9-0d2db5901e8a.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0640v1-figure-8be94a46-8166-48fd-b433-a1c4c75bb357.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0640v1-figure-450712ed-fe5a-4fa6-8ac1-582de0e49f21.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0640v1-figure-941f7f59-1180-4270-9915-d37147ade070.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0640v1-figure-55890bd7-7882-4676-b1e7-f63041db46cc.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0640v1-figure-ab71dddc-1615-4290-b6f0-e0469f4dfff0.jpg)

by

$$N _ { \text{eff} } ( t ) = \frac { \int _ { 0 } ^ { t } W _ { \text{eff} } ( r ) \, d r } { \int _ { 0 } ^ { t } \mu ( r ) \, d r } \,.$$

glyph[negationslash]

Here, η µ = N eff /N (= η ν ) defines the reporting proportion, but again one must employ suitable averages. For our measles data we obtain N eff = 32 8 . × 10 6 when proposing a constant effective size during the period 1948-67; see the reconstruction results (red curves) in Figures 2C and 3.

Alternatively, one might start with a moving average N eff ( ) t of (12b) and apply the reconstruction to raw data W t ( ). Such a procedure works, although it does not seem to be robust (blue curve in Fig. 2C).

## 2.6 Numerical simulations

When modeling real data, the functions W t ( ) and M t ( ) are given at a certain resolution in time (e.g., t n +1 -t n = 1week). Therefore, it must be explained how continuous-time formulas are applied to functions defined an a discrete time line { t n } n ∈ N . The derivative of a function F t ( ) at t = t n is approximated by F ′ ( t n ) = F t ( n +1 ) -F t ( n -1 ) 2∆ t . This definition works well as long as infectious periods are short and removal rates λ are large. Admissible values are given by

$$\lambda \geq \lambda ^ { * } = \sup _ { n } \frac { | I ^ { \prime } ( t _ { n } ) | } { I ( t _ { n } ) } \,,$$

imposing that incidence rates are non-negative; see (8). The values (1 /λ = 9 9 7 days) utilized for the , , infections in this article (measles, dengue, flu) fulfill this condition for nearly all data points ( &gt; 99%), and λ glyph[greatermuch] λ ∗ for most of the data points (Fig. 2A).

Numerical simulations indicate that the specific choice of λ is not important on a large time scale (cf. the 2-year moving average of beta in Fig. 2E). This result, which exempts us from applying other methods to determine precise infectious periods, is confirmed by the following approximate equation and its λ -independent right hand side,

$$\frac { \beta ( t ) } { \lambda } \approx \frac { N ( t ) } { N ( t ) - N ( 0 ) \left ( 1 - \frac { 1 } { \mathcal { R } _ { 0 } } \right ) e ^ { - \rho t } - \int _ { 0 } ^ { t } e ^ { - \rho ( r - t ) } W ( r ) d r } \,.$$

The equation is derived from (11) and (9c) by incorporating the estimates: W ≈ λI , λ ≈ γ , I glyph[lessmuch] N .

The reconstruction formulas, (9a) and (9c), contain integrands of the form e κr F r ( ). Their discrete versions require particular care, especially when time resolution is small (i.e., κ ∆ t glyph[greatermuch] 0) and the exponential term dominates. We utilize the Gaussian quadrature rule,

$$\int _ { 0 } ^ { t } e ^ { \kappa r } F ( r ) \, d r = \sum _ { n } w _ { n } e ^ { \kappa t _ { n } } F ( t _ { n } ) \,,$$

with an index independent weight w n = 1 -e -κ ∆ t κ &lt; ∆ . The weight accounts for the exponential contribution t in between observations, imposing that ∫ t n t n -∆ t e κr F t ( n ) dr = w e n κt n F t ( n ) , ∀ n .

Usually integrals are approximated well by summation (15) with the trivial weight, w n = 1. Here this would lead to huge errors: the integral in (9a), for example, must be corrected by w n ≈ 2 / 3, if employing weekly data and an infectious period of one week (where κ ∆ t ≈ 1). This might be a reason why the method has not been proposed in the literature before. Another similar reason might be the non-trivial scaling

behavior discussed in the two previous subsections. In [14,17], the scaling problem has not been addressed. By computing ' β/N ' the authors avoid the problem of adjusting N ; they only need to estimate the initial value β (0) /N (0). While being apparent in numerical simulations, scaling is absent in theoretical calculations and therefore unnoticed in [15,16].

## 2.7 Computations

The eight figures in this article can be re-computed in less than 20 min. On a usual laptop, our Mathematica [32] implementation

```
                                                                                                                                                                                                        (*)                                                                                                                                                                                                       
                                                                                                                                                                                                        "*)                                                                                                                                                                                                      

                                                                                                                                                                                                       "                                                                                                                                                                                                        #                                                                                                                                                                                                        )                                                                                                                                                                                                        }                                                                                                                                                                                                        :                                                                                                                                                                                                        |                                                                                                                                                                                                        ]                                                                                                                                                                                                        \
                                                                                                                                                                                                       |
                                                                                                                                                                                                      
	                                                                                                                                                                                                       '                                                                                                                                                                                                       -                                                                                                                                                                                                        +                                                                                                                                                                                                        <
```

takes about 2 sec to reconstruct beta for 1,000 data points (Figs. 3A, 4A) and about 100 sec for 7,000 points (Figs. 3C, 4B). Because of the convolution integrals in (9), the computation time increases with the square of the data size.

## 2.8 Error estimates

For our measles data, the transmission rate β t ( ) is known to have an annual period [24]. The reconstruction, however, is showing frequency modes with periods longer than one year, especially when the time resolution is low (cf. Figs. 3B, 4AC). Interpolation of data successfully removes the extra modes, suggesting that these modes are caused by an insufficient data resolution.

Another major source of error is the reporting proportion, which we calculate by a version of (12a). For both sources of error we determine the error of beta after the reconstruction (Fig. 3). For simplicity we only consider the case where I ( ≈ W/λ ) is given and R is computed according to (9c). We assume that η ν ( ) t ≈ η ν

Figure 3: Errors. Transmission rates are reconstructed for a variable (black) and a constant population size (red, N eff = 32 8 . × 10 ). 6 The resulting curves are very similar. For the original (1-week) resolution, shown in Panels A and B, the deviation stays within the error bounds. Errors due to reporting (gray) and time resolution (blue) - both estimated for a variable population size - are shown in the panels on the right hand side. Panels C and D show reconstructions for an increased resolution (days); interpolation is done via cubic polynomials. Here the deviation exceeds the tiny error bounds (in Panel D, barely visible). The reconstruction at high resolution produces more uniform looking curves, β t ( ), with higher amplitudes than the reconstruction at the original resolution. The reconstruction at high resolution has an annual instead of a biennial period.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0640v1-figure-6b881cfc-e283-4e58-8af3-9c093d1806ea.jpg)

and | I ′ ( ) t | /I t ( ) glyph[lessmuch] λ . Taylor expansion of (7) then yields

$$| \Delta \beta ( t ) | & = \, \frac { | 1 - R ( 0 ) e ^ { - \rho t } + R ( t ) | } { \lambda \, \overline { \eta } _ { \nu } N ( t ) } \, \beta ( t ) ^ { 2 } \, | \Delta \eta _ { \nu } ( t ) | + \mathcal { O } ( | \Delta \eta _ { \nu } ( t ) | ^ { 2 } ) \\ & \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cd.$$

$$+ \left ( 1 + \frac { \gamma } { \rho } \left ( 1 - e ^ { - \rho t } \right ) \right ) \left | 2 - \frac { \lambda } { \beta ( t ) } - \frac { 2 \, R ( 0 ) } { e ^ { \rho t } N ( t ) } \right | \frac { \bar { \eta } _ { \nu } \beta ( t ) ^ { 3 } } { \lambda ^ { 2 } N ( t ) } \left | \Delta \eta _ { \nu } ( t ) \right | \left | \Delta I ( t ) \right | \quad \ \ ( 1 6 b )$$

$$+ \left ( 1 + \frac { \gamma } { \rho } \left ( 1 - e ^ { - \rho t } \right ) \right ) \frac { \bar { \eta } _ { \nu } \beta ( t ) ^ { 2 } } { \lambda N ( t ) } \left | \Delta I ( t ) \right | + \mathcal { O } ( | \Delta I ( t ) | ^ { 2 } )$$

as | ∆ η ν ( ) t | → 0 and | ∆ ( ) I t | → 0. The error of the reporting proportion is estimated by the deviation from the mean, | ∆ η ν | = | η ν -η ν | . The error of the data resolution is estimated by the error related to the time resolution, | ∆ I | = | I ′ | ∆ + t O ((∆ ) ), where ∆ t 2 t n = ( t n +1 -t n ) / 2.

## 3 Results

For three sets of infectious disease data we reconstruct beta as a function of time. Each set represents an epidemiological scenario for which the method can provide new insights into the transmission dynamics. Epidemiological results are discussed immediately, methodological aspects are reviewed in the Discussion.

## 3.1 Measles

Transmission rate reconstructions have taught us about the epidemiology of measles. The first result [24] was obtained with data from the pre-vaccination era in England and Wales (Fig. 1A). Biennial patterns of incidence were related to transmission events with an annual period. In particular, the three school terms and the Christmas holidays were identified as periods of high transmission (cf. our reconstructions in Fig. 4CD). Subsequent work [8-10] introduced probabilistic methods - to achieve better parameter estimates and, being able to control errors, to put the beta reconstruction on a more trustworthy foundation. The more recent work [11-13] combines nonlinear optimization techniques and compartmental models, which can significantly shorten the computation time [12].

Our approach is similar to the classical one [24]. In that it is based on SIR dynamics (1) and requires the reconstruction of the involved compartments. In contrast to the discrete-time approach of [24], we use continuous-time expressions and work with differential equations, similar to [14,15]. For our numerical simulations we define a discrete calculus, including a refined representation (15) of the convolution integrals, which model the time evolution. Initial conditions of (9a) and (9c) are chosen in accordance with steady state solutions, where (9c) employs the initial value of the basic reproduction. In practice, we choose the initial value to match the basic reproduction at the end of considered time period (end of 1966), where the mean age at attack and the mean life expectancy are available most accurately [35].

The reconstruction requires knowledge of the mortality distribution, although knowing the life expectancy turns out to be sufficient (Fig. 1B). The mean infectious period (13) is another input parameter. The reconstruction, however, does not rely much on its precision (Fig. 2EF). If demographic information is available, as being the case for England and Wales [31], the distribution of births is used to estimate the reporting proportion (12a). This parameter is needed for calibration; it represents the epidemiological data - which usually suffer from systematic errors - as effective numbers (11). Otherwise, if demographic information is incomplete, one assumes that the population is of constant size, determined by the mean life expectancy (12b) or, even simpler, by the requirement of convergence. The size of the effective population

Figure 4: Beta reconstruction and data reproduction. Panels A,C show β t ( ) reconstructed from weekly incidences of measles in England and Wales (cf. Fig. 1A). The infectious period is assumed to be 9 days. Panels B,D illustrate the reconstruction for a 1-day resolution, obtained by interpolation of the original data via cubic polynomials. Panels C,D zoom into Panels A,B so that the three school terms of high transmission and the drop in summer can be seen. The peak shortly after the Christmas holidays has been smoothed out by the plotting resolution of 1 month. At the original (1-week) resolution, the biennial period - induced by the pattern of incidence - is visible. For a 1-day resolution, this is different: the beta curves of odd and even years are similar. Errors due to insufficient time resolution are indicated blue; these errors vanish at high resolution. In Panels E-H, prevalence data I ( ) are reproduced from the reconstructed t β t ( ). Reproduction results profit from higher resolution (Panels E,F) and from shorter infectious periods (Panels G,H). For a 1-day resolution and an 8-day infectious period (Panel H), the data are reproduced almost identically.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0640v1-figure-6128b71a-d607-4d14-9e1d-2e9ea4ad74b5.jpg)

happens to be very sensitive (Fig. 2C), allowing for such a criterion. Both our reconstructed betas -based on a variable (black) and on a constant population size (red) - resemble the curves known from the literature. The deviation between our two curves is small (Fig. 3) and stays within the error bounds of the original resolution (Fig. 3B).

The periodic behavior of β t ( ) tells us about another parameter, the infectious period. We find that for short infectious periods (1 /λ &lt; 10 days) the reconstructed transmission rates show biennial frequency modes (Fig. 4A): in 1962 and 1964 the beta curves are similar (beta decreases on average), in 1963 and 1965 the curves are similar too (but beta increases), whereas in 1962 and 1963 the curves look different (Figs. 3B, 4C). Even if prevalence of measles infections is known to show various frequency modes [25], modes of the transmission rate that are larger than one year are unrealistic (at least for the U.K. [36]). Therefore, the result must represent an artifact, which we assume is caused by the poor resolution of the original data. In fact, errors are large every second year (Figs. 3B, 4C). If, by interpolation, we increase the resolution from weeks to days, we observe that biennial modes change into annual modes (cf. Figs. 3D, 4D). The same effect, however, can be achieved by employing the original data with longer yet more unrealistic infectious periods.

How about the ability to reproduce the original data? We perform the reproduction by solving the basic equations (1a, 1b) employing the reconstructed betas together with the parameters used for the reconstruction (Fig. 4E-H). Here, shorter infectious periods turn out to be in favor (Fig. 4G vs. 4H); the quality of the reproduction in Figure 4H is striking. Similar to the reconstruction, reproduction results are satisfactory only at high resolution (Fig. 4E vs. 4F). The poor time resolution in [24] might thus explain why the reproduction of data (from betas computed in [24]), investigated in [37], did not work so well. The same problem might occur when reproducing data based on betas obtained in [14]; their beta curves seem to have a low resolution as well.

We conclude that there is a trade-off between the quality of the reconstruction and the quality of the reproduction balanced by the length of the infectious period. Working with higher resolution - obtained by interpolation of the original data - is beneficial for both, reconstruction and reproduction. This is not surprising, as the assumption underlying the applicability of interpolation is smoothness, being consistent with our methodology relying on differential equations. Good estimates of the infectious period, however, are relevant only for the short-term behavior of beta (Fig. 2F). Long-term, the length of the infectious period is not important (Fig. 2E).

To widen the applicability of the method we start to extend the reconstruction to models with multiple infective stages. For measles, it is common to also introduce a compartment of exposed (infected and not yet infective) individuals E with removal rate ε . Accordingly, we replace Equation (1b) by

$$E ^ { \prime } = \beta S I / N - ( \varepsilon + \mu ) E$$

$$I ^ { \prime } & = \varepsilon E - \lambda I \, ; & ( 1 7 b )$$

incidence is then given by W = εE . In doing so we reconstruct beta according to an SEIR model, which relies on more realistic assumptions than an SIR model. Consequently, the result (Fig. 5) is more precise, it resembles the beta curve of the SIR-model (Fig. 3CD) reconstructed at higher time resolution.

Figure 5: SIR versus SEIR reconstruction. The diagrams show β t ( ) reconstructed for SIR (black) and SEIR systems (green) with mean incubation period of 1 /ε = 5 days. The deviation between the beta values is plotted (blue) in the lower part of the diagrams. For both models, the mean infectious period is 9 days. Computed at original (1-week) resolution, the SIR reconstruction shows an artifactual two year period, whereas the SEIR reconstruction does not.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0640v1-figure-5b2dbe35-4da4-4a28-bd97-ff6662127b72.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0640v1-figure-134db5e0-f7a3-4a57-9520-aa287e8dc175.jpg)

## 3.2 Dengue

Another infection that involves multiple infective stages is dengue. This is a vector-born infection caused by an RNA virus that occurs in four serotypes with oscillating prevalence in the human host-population. Primary infections are acute and only last for a couple of weeks. Individuals are protected against reinfection with the same serotype strain but live at risk to become severely ill and develop symptoms referred to as dengue hemorrhagic fever (DHF) when reinfected with another serotype. This phenomenon is known as an antibody-dependent enhancement (ADE) effect [38].

Recently, dengue has been reported to re-emerge in some regions of the world, although transmission rates have successfully been reduced by vector-control measures [26]. The apparent contradiction has been explained by two competing mechanisms: clinical cross-protection [28] and demographic transitions [29]. The first approach focuses on host-pathogen interactions, the second one on changes in the host population. But which of two offers the correct explanation? This is the question we try to answer in this section utilizing the beta reconstruction. Different from both these approaches, our reconstruction relates DHF data and transmission rates via a formal procedure. This may falsify one of the suggested explanations.

Three compartments in the human host-population are relevant for the dengue epidemic: susceptible, infected, and exposed individuals. A model based on the contact functional Σ[ S, I ] is given by

$$S ^ { \prime } & = N ^ { \prime } + \mu ( N - S ) - \beta \Sigma [ S, I ] & ( 1 8 \text{a} )$$

$$I ^ { \prime } & = \beta \Sigma [ S, I ] + \beta k \Sigma [ E, I ] - \lambda I$$

$$E ^ { \prime } & = F [ I ] - \mu E - \beta k \Sigma [ E, I ] \,,$$

where strain-specific information is encoded by the number k . The functional F represents the replenishment of exposed individuals. The exposed compartment refers to individuals who were infected at least once and have therefore been exposed to one or more serotypes. These individuals are also susceptibles, even if becoming reinfected with a new serotype strain is less likely than getting infected for the first time; the parameter k ( &lt; 1) models the likelihood of reinfection. We do not need to estimate k via probabilistic arguments or extra data, the implementation is achieved through consistency requirements explained below.

In [28], removal of infected individuals into the exposed compartment has been proposed to be timedelayed and reduced by cross-protection. This mechanism can be modeled by two additional compartments,

$$R ^ { \prime } = ( \lambda - \mu ) I - ( \rho + \mu ) R - \beta l \Sigma [ R, I ]$$

$$P ^ { \prime } & = \beta l \Sigma [ R, I ] - \mu P \,,$$

where R represents temporarily and P permanently immunized individuals. The contact term βl Σ, which corrects for multiple serotypes, accounts for events of permanent immunization, the term

$$F [ R ] & = \rho R$$

quantifies the rate F R I [ [ ]] at which the exposed compartment is replenished.

Similar to the parameter k , we do not need to estimate l or ρ to reconstruct beta. We only estimate the likelihood that relates DHF rates M and incidences of secondary infection,

$$q = \frac { M } { \lambda H } \,,$$

where H = hI (with h &lt; 1) represents the infected individuals with prior exposure. It is sufficient to consider secondary infections here, as DHF data rarely account for cases resulting from primary infection: p λ I 1 ( -H ) + p λH 2 ≈ p λH 2 , as h glyph[greatermuch] p /p 1 2 . Due to ADE, the risk of developing DHF after re-infection is p /p 2 1 = 16 9 times higher than without prior exposure; see legend to Figure 3 in [28] and Table 1. .

The unknown DHF-likelihood q is determined by the global behavior of β t ( ). For dengue in Thailand, beta is known to decline (blue curve in Fig. 6D), corresponding to an increasing mean age of the individuals suffering from DHF [28].

It may surprise that q is all that is needed to reconstruct the correct transmission rates; our beta curve (black) successfully fits the (blue) reference curve in Figure 6CD. We computed β t ( ) following the second case in Methods while using mass-action (5). The two sets of dengue data (in Fig. 6) are taken from [28], representing DHF cases in Thailand (right) and a simulated model scenario (left).

To understand how the method works one must realize that β t ( ) is not reconstructed literally - via S and E , as defined in (18). The method keeps track of the exposed individuals E and reconstructs a truncated version of (18b),

$$H ^ { \prime } = \beta E K / N - \lambda H \,,$$

where K = kI models the compartment of individuals that infect exposed individuals. When expressing replenishment of exposed individuals (18f) by a multiple of natality,

$$F [ N ] & = \eta \, \mu N$$

(given by mortality, as the reconstruction is done for a constant effective population), one obtains the relevant subset of (18),

$$\frac { E ^ { \prime } } { \eta } = \mu N - \mu \frac { E } { \eta } - \beta \, \frac { E } { \eta } \frac { K } { N }$$

$$K ^ { \prime } = \beta \, \frac { E } { \eta } \, \frac { K } { N } - \lambda K \,.$$

L

0

R

unitsof

H

b

L

0

R

unitsof

H

b

12

10

8

6

4

12

10

8

6

4

40

50

60

70

80

90

100

1985

1985

1990

1990

1995

1995

B

2000

D

2000

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0640v1-figure-4822a33b-07fe-4663-b7ac-b01c5a5b4ea6.jpg)

1985

1990

1995

2000

10

9

8

7

6

10

9

8

7

6

35

30

25

20

15

10

5

0

Year

Figure 6: Dengue. The diagrams on the right hand side represent data from Thailand, the ones on the left hand side a model scenario. Red curves show DHF data, blue curves the transmission rates. These curves are taken from Figure 4 in [28] (at a resolution of 220 × 160 and 700 × 160 sq. pixels, for real and model data, respectively). Panels A and B show the reconstruction results for scaling factors q deviating from the optimal values ( q/p 2 = 2 55 and 2 19, resp.) . . by ± 5% (orange/green). In accordance with [28], the risk of developing DHF is assumed to be p 2 = 0338, and the mean infectious period 1 . /λ = 9days. The black curves in Panels C and D are reconstructed from DHF data (red curves in Panels E and F). The result is nearly perfect, except for the beta value in 1996 (Panel D). The blue curve in D is computed based on the mean age of the infected individuals suffering from DHF [28] (i.e., the blue curve is not related to the red curve in Panel F by computation). The reconstructed betas are plotted gray, their half-year moving average black. DHF cases (black curves in Panels E and F) are reproduced by employing the reconstructed betas following the exposed SIR-dynamics (explained in the text). The data curve (Panel F) is almost identically reproduced (red spots are barely recognizable). For the model scenario (Panel E), reproduction is good for the first 10-20 years.

40

50

60

70

80

A

90

C

100

F

Table 1: Parameter estimates. Values (in adequate powers of weeks) are based on the optimal choice of q = p η 2 (cf. Fig. 6AB) and the most likely duration of cross-protection (Fig. 7). The computation uses steady state solutions of (18) and (22), where β is approximated by the mean value.

|      | model data     | real data      |
|------|----------------|----------------|
| η    | 2.19           | 2.55           |
| 1 /ρ | 52             | 25             |
| β/λ  | 7.54           | 8.02           |
| h    | .738           | .803           |
| k    | .337           | .315           |
| l    | .580           | .542           |
| α    | 6.54           | 7.02           |
| ν 0  | - 11 . 0       | - 9 . 56       |
| ν 1  | .0174          | .0201          |
| ν 2  | 2 . 4 × 10 - 7 | 2 . 3 × 10 - 7 |

Comparison with the basic SIR equations (1a, 1b) reveals the following correspondences, E/η → S and K → I , where

$$k = \frac { h } { \eta } \quad \text{and} \quad \eta = \frac { q } { p _ { 2 } } \,.$$

These definitions demonstrate that the DHF-likelihood q is the only parameter required as input. Numerical simulations (Fig. 6AB) determine its value at 2-3 times the risk of developing DHF after re-infection (cf. Table 1).

The other parameters ( h, k, l ; cf. Table 1) are calculated with the help of steady state solutions ( S ′ = E ′ = I ′ = R ′ = 0) and mass-action (5): from (22) we extract expressions for E and I , plug them into the steady state versions of (18), and solve for N ′ . The result maps the mean duration of cross-protection 1 /ρ to the time-derivative of the population size N ′ ; first order expansion offers a good approximation,

$$N ^ { \prime } ( 1 / \rho ) = \nu _ { 0 } + \nu _ { 1 } \left ( 1 + \alpha / \kappa \right ) / \rho + \underbrace { \nu _ { 2 } \left ( 1 + \alpha / \kappa \right ) ^ { 2 } / \rho ^ { 2 } + \dots } _ { \approx 0 \,, \, \text{cf. Table 1} } \,.$$

glyph[negationslash]

The formula contains the transmission rate ratio κ = k/l . For our model data, where N ′ = 0, this ratio can be calculated: κ 0 = 581. . The situation is different for real data from Thailand, as these data involve demographic changes [29] (i.e., N ′ = 0). Nevertheless, we use κ = κ 0 as an estimate; this value is close to the one ( κ = 1 1 85) mentioned in [29] (legend to Fig. 6A). / .

As the parameter ρ does not enter the truncated dynamics (22) explicitly, we model clinical crossprotection by shifting the time line and estimate its duration by multiple reconstructions. We calculate likelihoods for all possible durations (Fig. 7AB). Our model data show a peak at one year (Fig. 7A), which is in perfect agreement with the simulation setting in [28]. Real data from Thailand show several peaks the latest and highest one at half a year (Fig. 7B). A half-year period corresponds to a decreasing effective population (cf. the dashed line in Fig. 7D) and thus confirms the conclusions in [29].

Therefore, our results confirm that clinical cross-protection and demographic transitions - suggested to explain the discrepancy between DHF counts and decreasing transmission rates - represent possible and non-exclusive mechanisms. Without many changes we were able to extend our method to an infection that

Figure 7: Clinical cross-protection versus demographic transitions. Similar to Figure 6, the diagrams on the right hand side represent dengue data from Thailand, the ones on the left hand side the model scenario in [28]. Panels A and B show the likelihoods for the various mean durations of cross-protection - obtained by shifting the time line underlying the reconstruction and minimizing the squared distance between the timeshifted reconstruction and the original (blue) beta curve (in Fig. 6). The size of the time-shift represents the duration of cross-protection, the frequency of the distance minimums measures the corresponding likelihood. Panels C and D illustrate the (nearly linear) relationship between the mean duration of cross-protection 1 /ρ and the change of the effective population N ′ with respect to the transmission rate ratio κ = k/l (bluegreenish color). The model data incorporate cross-protection for a period of one year (cf. legend to Fig. 4 in [28]), corresponding to the likelihood maximum (red) in Panel A. For real data (Panel B), the most likely period of cross-protection is slightly below half a year. When utilizing the transmission rate ratio from the model data (red dot), the Thailand data show a declining effective population ( N ′ = -12 9 . / month) as proposed in [29].

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0640v1-figure-7fd9a630-91dd-4f2e-b059-c7f8bb19e71b.jpg)

involves several stages. We learned that the method is relying on an SIR-structure similar to the one of the measles example but with exposed individuals taking over the role of susceptibles. Multiple application of the reconstruction procedure, as required in our third example, will represent another extension.

## 3.3 Influenza

We apply the method to multiple pathogenic strains, examining mortality data [23] related to the 1918/19 influenza pandemic in the U.S. Here we try to understand how the pandemic strain has replaced the residing strain after entering the host population. To this end, we determine the strain-specific infections for periods covering the pandemic peak and reconstruct the corresponding transmission rates. Case-fatality ratios enable us to express the death rate as a linear combination over the strains,

$$M = \lambda _ { 1 } I _ { 1 } / c _ { 1 } + \lambda _ { 2 } I _ { 2 } / c _ { 2 } \,.$$

Transmission rates do not superpose in a linear fashion, glyph[negationslash]

$$\beta \neq \beta _ { 1 } + \beta _ { 2 } \,,$$

which makes the reconstruction problem a non-trivial task.

Before we start, we like to emphasize that our modeling approach is again minimalistic. Unlike [39], for example, we do not incorporate extra information (about school terms, temperature changes, or changes in human behavior) in addition to the mortality data around 1918/19. We only propose two standard SIR models - one for the residing strain (1) and one for the invading (supposedly pandemic) strain (2).

The time evolution of the SI-compartments of the residing strain is given by

$$S _ { 1 } ^ { \prime } & = \rho _ { 1 } ( N - S _ { 1 } - I _ { 1 } ) - \Delta _ { 2 } S _ { 1 } - \beta _ { 1 } S _ { 1 } I _ { 1 } / N & ( 2 7 \text{a} )$$

$$I _ { 1 } ^ { \prime } & = \beta _ { 1 } S _ { 1 } I _ { 1 } / N - \lambda _ { 1 } I _ { 1 } \,.$$

Replenishment of susceptibles is dominated by the decay of immunity, ρ 1 ≈ δ (with 1 /δ = 5-10years, as common for flu). The equations include two further assumptions. First. Previously infected individuals are partly cross-protected against the other strain. That is, the invading strain, which infects S 2 ( t 0 ) -S 2 ( ) = t -∫ t t 0 S ′ 2 ( τ ) dτ individuals up to time t &gt; t 0 , reduces the number of individuals susceptible to the residing strain approximately at the rate

$$\Delta _ { 2 } S _ { 1 } & = \chi \frac { S _ { 1 } } { N } S _ { 2 } ^ { \prime } \,.$$

The parameter χ ≤ 1 defines the amount of cross-protection. Second. We reconstruct the transmission rate β 1 ( ) from data prior to the invasion, when only strain 1 is present, and - assuming that strain 1 does not t change phenotypically - extrapolate its averaged annual shape (red curves on the left-hand side panels in Figure 8) beyond the pandemic peak.

Motivated by the neutral network mechanism [40], the invading strain is modeled by an almost identical copy of the residing strain (27),

$$S _ { 2 } ^ { \prime } & = \rho _ { 2 } ( N - S _ { 2 } - I _ { 2 } ) - \Delta _ { 1 } S _ { 2 } - \beta _ { 2 } S _ { 2 } I _ { 2 } / N & ( 2 8 \text{a} )$$

$$I _ { 2 } ^ { \prime } & = \beta _ { 2 } S _ { 2 } I _ { 2 } / N - \lambda _ { 2 } I _ { 2 } \,.$$

Figure 8: Exchange of strains. The four sets of diagrams (Panels AB-GH) illustrate invasion scenarios, which differ by the time t 0 the pandemic strain is introduced into the human host population (cf. figure labels). The diagrams on the right hand side show the time evolution of infectives and susceptibles, I 1 (red), I 2 (blue), I 1 + I 2 (light gray), S 1 (orange), S 2 (cyan) - the diagrams on the left hand side the corresponding transmission rates, β 1 (red), β 2 (blue), β (gray). Infectious periods of 1 week = 1 /γ 1 = 1 /γ 2 ensure that the scale 'weeks/N' corresponds to units of R 0 . The other parameter values ( N = 1 5 . × 10 , 8 c 1 = 1%, R 0 = 2) are chosen so that the moving average of β 1 ( ) t is constant before the pandemic peak t p . These choices define a realistic model scenario [41]. The transmission rate of the resident strain β 1 ( ) (= t β t ( ) for t &lt; t 0 ) is continued beyond t 0 by a combination of the annual shape β t ( ) prior to the peak and a copy of β t ( ) that is stopped right after the peak, β 1 = (1 -r β ) + rβ . Ratios r &gt; 0 support the assumption that I 1 contributes to the peak; here we chose r = 1 4. / Independent of this choice, I 1 ( ) quickly declines after the t peak. At values below 1 (Panels B,D) and below 10 -2 (Panels F,H) the resident strain is removed from the system.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0640v1-figure-7cb45314-06af-46fb-9ef3-7625c79a7e6e.jpg)

Initially, susceptibles to the invading strain S 2 are only replenished by birth, ρ 2 ≈ µ glyph[lessmuch] δ . That is, the replacement of S 2 is small. Therefore, (28a) is approximated by S ′ 2 = -∆ 1 S 2 -λ I 2 2 -I ′ 2 , where due to (28b) the contact term β S I /N 2 2 2 is expressed by λI 2 + I ′ 2 . According to (25), we further replace I 2 and I ′ 2 by M and I 1 .

The number of susceptibles to the invading strain S 2 changes around the pandemic peak t p . Before and after the peak, S 2 ( ) can be approximated by a constant. That is, the cross-protection term reads t

$$\Delta _ { 1 } S _ { 2 } ( t ) = \chi \left ( \frac { \theta ( t - t _ { p } ) } { \mathcal { R } _ { 0 } ^ { 1 } } + \frac { \theta ( t _ { p } - t ) } { \mathcal { R } _ { 0 } ^ { 2 } } \right ) S _ { 1 } ^ { \prime } ( t ) \,,$$

where θ t ( ) denotes the unit step function and R k 0 an estimate for the basic reproduction of strain k .

The number of susceptibles to the residing strain S 1 ( ) t is nearly constant, even in the vicinity of the pandemic peak. That is, S ′ 1 ( ) t ≈ 0. One may therefore neglect cross-protection for the invading strain (29) and only consider contacts in (28a). The outcome barely changes.

Either way, S ′ 2 is defined in terms of M M , ′ , I 1 , I ′ 1 , and S ′ 1 . As cross-protection ∆ 2 S 1 [ S , S 1 ′ 2 ] is the only remaining term containing the index 2, the system of strain 1 (27) can now be solved. Finally, the transmission rate of the invading strain is calculated by

$$\beta _ { 2 } = \frac { I _ { 2 } ^ { \prime } + \lambda _ { 2 } I _ { 2 } } { S _ { 2 } I _ { 2 } } \, N \,,$$

where again S 2 , I 2 are expressed by M I , 1 and their derivatives.

Following the instructions outlined above we perform numerical simulations where we vary the time t 0 &lt; 1918 at which the invading strain is introduced into the population. The blue curves in the left-hand side panels of Figure 8 show the evolving transmission rates β 2 ( ) in four example scenarios. t Although these scenarios vary largely with respect to the choice of t 0 , the beta reconstruction produces very similar results. For up to three years after introduction, the transmission rate of the invading strain ( I 2 glyph[greatermuch] 0) fluctuates around values of 1 (in units of R 0 ) until suddenly it becomes singular, resembling a Dirac delta at the pandemic peak t p . A closer look reveals that β t ( p ) never exceeds values 5-7 (in units of R 0 ). Then β 2 ( ) t increases and approaches the values of the residing strain β 1 ( ) (red) as t t →∞ .

The gray curves in the left-hand side diagrams of Figure 8 illustrate the reconstruction performed without differentiation between the two strains. One observes a huge increase of β t ( ) at the pandemic peak, which only stops after a couple of years. We confirm that β is not the result of a linear superposition of β 1 and β 2 . Except for the singularity at the pandemic peak, the transmission rates of the invading strain β 2 ( ) behave t very regularly and do not seem to resemble a pandemic at all. There is no need to model a higher-thannormal case-fatality ratio [30] and/or a shorter infectious period for strain 2. In reality, parameters might have changed back to normal values (characteristic for annual flu) very quickly, explaining why small (or no) changes of these parameters accomplish a valid effective theory (cf. Fig. 8CD).

The scenario in Figure 8G, where the pandemic strain is introduced at the earliest time t 0 (week 5 / 1915), shows strong fluctuations with singularities. Even more irregular behavior of β 2 ( ) occurs in scenarios of yet t earlier t 0 , although certain parameter combinations ( c 2 , χ ) produce singularities also if strain 2 is introduced later. Scenarios with singularities of β 2 ( ) are unrealistic. t The singularity at t p is an exception, if interpreted as the initial cause of the pandemic. We cannot decide if changes of the virus (mutation or reassortment) are responsible for the singularity of β 2 ( ) in 1918. t Our computations only suggest that the pandemic could have been initiated by any kind of singularity in the transmission rates. Sudden changes in the host-contact

behavior, for example, induced by extended world-wide migration at the end of the war, represent another plausible cause.

Herald waves could have started as early as 1915/16. The reconstruction scenarios in Figure 8 indicate the possibility of three epidemic waves prior to the pandemic peak - periods during which both strains could have been present. The proportion of the two strains depends on our model parameters (e.g., χ ) but in no obvious way. Even if the prevalence of the invading strain is far above zero, the corresponding transmission rate prior to the pandemic peak never exceeds 1 much. However, the beta reconstruction always predicts an exchange of strains - at or immediately after the pandemic peak, where the prevalence of the residing strain reduces to nearly vanishing amounts.

The precise dynamics underlying the exchange of strains around 1918/19 was certainly more complicated. Nevertheless, the similarity of our results (Fig. 8, right) to those in [40] (Fig. 3A) demonstrates that the modeling approach via the reconstruction of transmission rates can unveil features inherent to the antigenic shift of influenza.

## 4 Discussion

We have developed a methodology that enables us to formulate models of disease transmission by extracting information from time series of incidence or mortality. Based on the assumption that the epidemic is governed by SIR dynamics (i.e., infection related compartments and some version of mass action) we reconstruct transmission rates as a function of time.

The reconstruction is given by explicit formulas (9, 10), with data represented by functions defined on a continuous time scale. Mathematically these formulas are easy to handle and straightforwardly provide us with error estimates (16). In numerical computations we replace the continuous-time expressions by carefully chosen discrete versions - critically important for the integrals involved; see (14). Doing so, we are able to reconstruct transmission rates for data of measles, dengue, and influenza. These examples, which cover very different types of infections [3], illustrate the method's wide range of application.

Besides mathematical and conceptual flexibility, the involved computations are fast. The core expression is a convolution integral, which determines the time complexity of the reconstruction to be quadratic. It takes less than two minutes to compute Figure 2B on a usual laptop. Using a probabilistic sampling method (e.g., [9, 10]), it may take a day to accomplish the reconstruction for a data set of similar size; 20 hrs is the reference given in [10].

The large difference in computation time is easy to understand. Sampling methods offer probabilistic analysis with confidence intervals whereas our method only determines the error resulting from error estimates of the input variables. Moreover, sampling provides intrinsic parameter estimates, which a priori our approach is lacking. Therefore, the comparison is not entirely fair. The more recent methods [11-13], which utilize nonlinear optimization techniques, are more similar to our compartmental approach. Yet, determining beta values and parameters at once, these methods need much more computing power and, concluding from the parameter behavior, might not gain so much.

Namely, when performing the reconstruction we learned that most of the parameters - removal rates and initial values - are not required to be precise. Different infectious periods, for example, produce similar outcomes, at least on a large time scale (cf. Figs. 2EF, 4E-H). Mortality is sufficiently defined by mean values (Fig. 1B). Initial values can be estimated using steady state solutions and external data. Their influence decreases exponentially - with decay of immunity or inverse lifetime - so that that asymptotically, for

large time, the reconstruction returns the correct transmission rates.

Only one parameter must be taken from a small range of values as otherwise the reconstruction yields unrealistic results or even collapses (Fig. 2C). For measles and flu, this parameter is the effective population size (if assumed to be of constant size), for dengue, it is the DHF-likelihood of secondary infections. Suitable parameter values are determined by the long-term behavior of beta. That is, a posteriori our method provides necessary parameter estimates.

Previous work on pre-vaccination data of measles infections in the U.K. [7-10, 24] offered a successful test of the method. We reconstructed the expected annual patterns of high and low transmission. More importantly, the reconstructed betas reproduced the original data. For dengue, we immediately re-obtained the data (Fig. 6D) from the reconstructed betas (Fig. 6F). For measles, reproduction did only work for a few years (Fig. 4E). Here we had to increase the time-resolution, which turned out to have a huge impact on the error of the computations (cf. Fig. 4AC). After interpolation, the error of the reconstructed betas disappeared (Figs. 4BD). With those betas we were then able to correctly reproduce the original data (Fig. 4F-H).

We also tested a simplified version of the method, which proposes an effective population of constant size -useful if demographic data are not available. The reconstruction results turned out to be quite similar to those based on a variable population size (Fig. 3).

Elaborate probabilistic simulations are often difficult to verify. Here our reconstruction can serve as a tool, allowing to cross-check assumptions in a simple SIR setting. The dengue application is an example where the method was used to re-evaluate two seemingly contradictory hypotheses, which claim to explain controversial observations related to the re-emerging dengue epidemic in Thailand [26]. Our reconstruction has shown that the two hypotheses, the one based on a medical [28] and the other on a demographic mechanism [29], are both plausible and not contradicting.

The beta reconstruction can be combined with other methods and approaches, leading to novel hypotheses relying on the time evolution of transmission rates. This has been demonstrated here with the influenza example. Having to differentiate between multiple strains [40, 42], flu was the most advanced application out of the three. We performed the reconstruction based on mortality data and simultaneously modeled the evolution of an invading strain. The resulting transmission rates suggest that the invasion could have happened many month ahead of the pandemic (Fig. 8), which is telling us about the origin of herald waves.

## Acknowledgements

This work was supported in part by a grant to David Earn, Troy Day, Jonathan Dushoff, and Junling Ma from the Canadian Institutes of Health Research (CIHR). The author thanks these researchers for valuable discussions during his stay at McMaster University.

## References

- [1] Anderson, RM, May RM (1991) Infectious Diseases of Humans: Dynamics and Control. Oxford Univ. Press.
- [2] Grassly NC, Fraser C (2008) Mathematical models of infectious disease transmission. Nat Rev Microbiol 6: 477-87.

- [3] Lange A, Ferguson NM (2009) Antigenic diversity, transmission mechanisms, and the evolution of pathogens. PLoS Comput Biol 5(10): e1000536.
- [4] Ferguson NM, Cummings DA, Cauchemez S, Fraser C, Riley S, Meeyai A, Iamsirithaworn S, Burke DS (2005) Strategies for containing an emerging influenza pandemic in Southeast Asia. Nature 437: 209-14.
- [5] Fraser C, Donnelly CA, Cauchemez S, Hanage WP, Van Kerkhove MD, Hollingsworth TD, Griffin J, Baggaley RF, Jenkins HE, Lyons EJ, Jombart T, Hinsley WR, Grassly NC, Balloux F, Ghani AC, Ferguson NM, Rambaut A, Pybus OG, Lopez-Gatell H, Alpuche-Aranda CM, Chapela IB, Zavala EP, Guevara DM, Checchi F, Garcia E, Hugonnet S, Roth C; WHO Rapid Pandemic Assessment Collaboration (2009) Pandemic potential of a strain of influenza A (H1N1): early findings. Science 324: 1557-61.
- [6] Legrand J, Egan JR, Hall IM, Cauchemez S, Leach S, Ferguson NM (2009) Estimating the location and spatial extent of a covert anthrax release. PLoS Comput Biol 5(1): e1000356.
- [7] Finkenstadt BF, Grenfell BT (2000) Time series modelling of childhood diseases: a dynamical systems approach. J R Stat Soc C 49: 187-205.
- [8] Bjornstad, ON, Finkenst¨dt, BF and Grenfell, BT (2002) Dynamics of measles epidemics: a Estimating scaling of transmission rates using a time series SIR model, Ecological Monographs, 72: 185-202.
- [9] Morton A, Finkenstadt BF (2005) Discrete time modelling of disease incidence time series by using Markov chain Monte Carlo methods. J R Stat Soc C 54: 575-594.
- [10] Cauchemez S, Ferguson NM (2008) Likelihood-based estimation of continuous-time epidemic models from time-series data: application to measles transmission in London. J R Soc Interface 5: 885-97
- [11] Hooker G, Ellner SP, De Vargas Roditi L, Earn DJD (2011) Parameterizing state-space models for infectious disease dynamics by generalized profiling: measles in Ontario. J R Soc Interface 8: 961-974.
- [12] Word DP, Cummings DA, Burke DS, Iamsirithaworn S, Laird CD (2012) A nonlinear programming approach for estimation of transmission parameters in childhood infectious disease using a continuous time model. J R Soc Interface 9: 1983-97
- [13] Word DP, Young JK, Cummings DA, Iamsirithaworn S, Laird CD (2013) Interior-point methods for estimating seasonal parameters in discrete-time infectious disease models. PLoS One 8(10): e74208
- [14] Pollicott M, Wang H, Weiss H (2011) Recovering the time-dependent transmission rate from infection data via solution of an inverse ODE problem. J Biol Dyn 6: 509-523
- [15] Hadeler K (2011) Parameter identification in epidemic models. Math Biosci 229: 185-189
- [16] Hadeler K (2012) Parameter estimation in epidemic models: simplified formulas. Can Appl Math Q
- [17] Mummert A (2013) Studying the recovery procedure for the time-dependent transmission rate(s) in epidemic models. J Math Biol 67: 483-507
- [18] Read JM, Keeling MJ (2003) Disease evolution on networks: the role of contact structure. Proc Roy Soc Lond B 270: 699-708.
- [19] Keeling MJ, Eames KT (2005) Review Networks and epidemic models. J R Soc Interface 2: 295-307.

| [20]   | Kermack WO, McKendrick AG (1927) A contribution to the mathematical theory of epidemics, Pro- ceedings of the Royal Society of London Series A 115: 700-21.                                                                                                                     |
|--------|---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| [21]   | Heesterbeek H. (2005) The law of mass-action in epidemiology: a historical perspective. In: Cuddington K, Beisner B, editors. Ecological paradigms lost: routes of theory change. Elsevier; pp. 81-105.                                                                         |
| [22]   | Diekmann O, Heesterbeek JAP, Metz JAJ (1990) On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations. J Math Biol 28: 365-82.                                                                   |
| [23]   | International Infectious Disease Data Archive. McMaster University. http://iidda.mcmaster.ca/ .                                                                                                                                                                                 |
| [24]   | Fine, PEM and Clarkson, JA (1982) Int J Epidemiology, 11: 5-14.                                                                                                                                                                                                                 |
| [25]   | Earn DJ, Rohani P, Bolker BM, Grenfell BT (2000) A simple model for complex dynamical transitions in epidemics. Science 287: 667-70.                                                                                                                                            |
| [26]   | Thammapalo S, Nagao Y, Sakamoto W, Saengtharatip S, Tsujitani M, Nakamura Y, Coleman PG, Davies C (2008) Relationship between Transmission Intensity and Incidence of Dengue Hemorrhagic Fever in Thailand. PLoS Negl Trop Dis. 2: e263.                                        |
| [27]   | Earn DJD, Rohani P, Grenfell BT (1998) Persistence, chaos and synchrony in ecology and epidemiology. Proc R Soc B 265: 7-10.                                                                                                                                                    |
| [28]   | Nagao Y, Koelle K (2008) Decreases in dengue transmission may act to increase the incidence of dengue hemorrhagic fever. Proc Natl Acad Sci USA 105: 2238-43.                                                                                                                   |
| [29]   | Cummings DAT, Iamsirithaworn S, Lessler JT, McDermott A, Prasanthong R, Nisalak A, Jarman RG, Burke DS, Gibbons RV (2008) The Impact of the Demographic Transition on Dengue in Thailand: Insights from a Statistical Analysis and Mathematical Modeling. PLoS Med 6: e1000139. |
| [30]   | Taubenberger JK, Morens DM (2006) 1918 influenza: the mother of all pandemics. Emerg Infect Dis 12: 15-22.                                                                                                                                                                      |
| [31]   | Hicks, J, Allen, G (1999) A Century of Change: Trends in UK statistics since 1900. House of Commons Library, Research paper 99/111.                                                                                                                                             |
| [32]   | Wolfram Research, Inc. (2007) Mathematica, Version 6, Champaign, IL.                                                                                                                                                                                                            |
| [33]   | Liu W, Hethcote HW, Levin SA (1987) Dynamical behavior of epidemiological models with nonlinear incidence rates. J. Math. Biol. 25: 359-380                                                                                                                                     |
| [34]   | Boyce WE, DiPrima RC (2005). Elementary Differential Equations and Boundary Value Problems. Wiley.                                                                                                                                                                              |
| [35]   | Griffiths DA (1974) A catalytic model of infection for measles. Appl. Statist. 23: 330-9.                                                                                                                                                                                       |
| [36]   | Gomes MC, Gomes JJ, Paulo AC (1999) Diphtheria, pertussis, and measles in Portugal before and after mass vaccination: a time series analysis. Eur J Epidemiol 15: 791-8.                                                                                                        |

- [37] Mollison D, Din SU (1993) Deterministic and stochastic models for the seasonal variability of measles transmission. Math Biosci 117: 155-77.
- [38] Ferguson N, Anderson R, Gupta S (1999) The effect of antibody-dependent enhancement on the transmission dynamics and persistence of multiple-strain pathogens. Proc Natl Acad Sci USA 96: 790-4.
- [39] He D, Dushoff J, Day T, Ma J, Earn DJD (2013) Inferring the causes of the three waves of the 1918 influenza pandemic in England and Wales. Proc. R. Soc. B 280: 20131345
- [40] Koelle K, Cobey S, Grenfell B, Pascual M (2006) Epochal evolution shapes the phylodynamics of interpandemic influenza A (H3N2) in humans. Science 314: 1898-903.
- [41] Mills CE, Robins JM, Lipsitch M (2004). Transmissibility of 1918 pandemic influenza. Nature 432: 904-6.
- [42] Smith DJ, Lapedes AS, de Jong JC, Bestebroer TM, Rimmelzwaan GF, Osterhaus ADME, Fourchier RAM (2004) Mapping the antigenic and genetic evolution of influenza virus. Science 305: 371-6.
Alexander Lange
Reconstruction of disease transmission rates: applications to measles, dengue, and influenza
Transmission rates are key in understanding the spread of infectious
diseases. Using the framework of compartmental models, we introduce a simple
method that enables us to reconstruct time series of transmission rates
directly from incidence or disease-related mortality data. The reconstruction
exploits differential equations, which model the time evolution of infective
stages and strains. Being sensitive to initial values, the method produces
asymptotically correct solutions. The computations are fast, with time
complexity being quadratic. We apply the reconstruction to data of measles
(England and Wales, 1948-67), dengue (Thailand, 1982-99), and influenza (U.S.,
1910-27). The Measles example offers comparison with earlier work. Here we
re-investigate reporting corrections, include and exclude demographic
information. The dengue example deals with the failure of vector-control
measures in reducing dengue hemorrhagic fever (DHF) in Thailand. Two competing
mechanisms have been held responsible: strain interaction and demographic
transitions. Our reconstruction reveals that both explanations are possible,
showing that the increase in DHF cases is consistent with decreasing
transmission rates resulting from reduced vector counts. The flu example
focuses on the 1918/19 pandemic, examining the transmission rate evolution for
an invading strain. Our analysis indicates that the pandemic strain could have
circulated in the population for many months before the pandemic was initiated
by an event of highly increased transmission.