IDENTIFIABILITY OF INFECTION MODEL PARAMETERS EARLY IN AN EPIDEMIC

2.1. The deterministic and stochastic SEIR models.

The deterministic version of the SEIR model [10, 14] that we will consider is

where the variables S, E, I, and R represent the populations of susceptible, exposed, infected, and removed patients, respectively, and $N=S+E+I+R$ denotes the total population. Time is measured in days. We use the simplest, or SEIR-without-vital-statistics, model, which assumes that N is constant with no births and deaths. There are more complex versions with additional parameters, but the identifiability issues we want to describe occur even for this simplest model. The sole nonlinearity is the βIS/N term, which moves patients from the susceptible compartment to the exposed compartment according to transmission rate coefficient β.

We will interpret the model in the following way. The parameter α is the time constant of movement from exposed to infected; thus we assume that on average, the patient spends 1/α days as exposed before transitioning to infected, where we assume viral shedding begins. We will also make the assumption that symptoms are present in patients in the I compartment, so that the case can be observable for the first time. After 1/γ days in the I compartment, on average, the patient is removed from the population and does not return to the susceptible class.

Our principal interest is in determining what information can be inferred from measured reports of infected cases I(t). We address two obvious limitations of these assumptions. First, perhaps not all infected cases are reported. Thus, the true infected number may be c₁I instead of I. Second, a portion of the infected cases may be asymptomatic and are not reported for this reason. Thus, the true infected number may be c₂c₁I. In either case, the true number of infected may not be knowable. If the true number of infected is proportional to the reported I, the meaning of the contact transmission parameter β will be changed. However, many of the purposes of using the model, such as to forecast future I(t), may still proceed unaffected.

In addition to the deterministic version, we will also consider the SEIR model as a set of stochastic differential equations with Poisson noise. In this version, we will calculate trajectories as follows. For each time step, the right-hand side of the equations will be evaluated by selecting from a Poisson distribution and then integrated using an Euler method step. In other words, the values

are chosen to represent the contribution of the right-hand side at each step, i.e.,

This version treats the SEIR model as a stochastic system for greater fidelity. However, our main conclusions about identifiability will be relevant for both the deterministic and stochastic versions.

2.2. Parameter estimation.

Parameter estimation is customarily achieved by locating, implicitly or explicitly, the optimum of some auxiliary function that measures the fitness of the parameters. In some methods, the likelihood or marginal probability is maximized, while in others, an error or loss function is minimized.

In one method for estimating parametersfrom daily reports of the single observable I(t), we will choose a particular loss function based on data assimilation and explicitly minimize it. This approach will be useful for illustrating the geometry of the minima of the loss function in two different parts of the SEIR trajectory. Our choice for the loss function will be the data assimilation error in I(t) incurred while using the proposed set of parameters to optimally reconstruct the trajectory $(S(t),E(t),I(t),R(t))$ from the observed I(t). The use of data assimilation to reconstruct unobserved variables is the basis of modern numerical weather prediction, and it has started to appear in epidemic modeling [6, 8, 12]. For the deterministic SEIR, we employ a standard EnKF (ensemble Kalman filter) [23, 22] to reconstruct the dynamics. For the stochastic SEIR, we use an EnKF tailored to Poisson noise instead of the standard Gaussian assumption. The EnKF used for this purpose is based on the PKF (Poisson Kalman filter) from [7].

Data assimilation gives a way of reconstructing all variables of a differential equations model from partial observations, for example, by measurements of one key variable. For SEIR model (2.1), if the parameters β, α, and γ are known, the observable I(t) or, alternatively, the daily changes $\text{Δ}I(t)=I(t)-I(t-1)$, are in general sufficient for reconstructing the other three variables S, E, and R. Figure 1(a) shows a trajectory of a stochastic SEIR model (2.1) with parameters β = 1.1, α = 0.2, and γ = 0.5, and with initial conditions $S={10}^6$, $E={10}^2$, $I=0$, $R=0$. The inputs to the data assimilation algorithm are the model, the exact parameters, and the daily reports of new infections $\text{Δ}I(t)=I(t)-I(t-1)$. The assimilation algorithm uses the initial condition S = 10⁶, E = 0, I = 0, R = 0. That is, it is allowed to know the (constant) total population but no information about the initial caseload. The EnKF is used to estimate the most likely values of S(t), E(t), I(t), and R(t) given the reports $\text{Δ}I(t)$. Figure 1(b) shows the resulting reconstructed trajectory, a reasonably accurate version of the original.

If the parameters are not known and incorrect parameters are used in the model, the reconstruction in general will be farther from the original. This leads to a convenient loss function to consider for the purposes of parameter estimation. Let $L(\beta ,\alpha ,\gamma )$ denote the mean squared difference between the observed $\text{Δ}I(t)$ and the reconstructed $\text{Δ}I(t)$ from the EnKF, over a time interval [T₁, T₂]. Then minimization of L as a function of the parameters should lead to the correct, or generating, parameters.

To begin, we carried out this idea on the deterministic SEIR model (2.2) with a standard simplex minimization algorithm [17]. We started the simplex algorithm with 1000 starting guesses for the parameters β, α, and γ that varied from the exact values by about 50%. Figure 2(a) shows the cumulative results of the minimization procedure for a trajectory of length 100 days, using two different intervals of observations, [T₁, T₂] = [0, 50] or [50, 100], with 1000 realizations of starting parameter guesses. There is a dramatic difference, depending on whether the time interval [0, 50] or [50, 100] is used for the input I(t). The red dotted curve is a histogram of approximate parameters using $\text{Δ}I(t)$ from the interval [0, 50]. The black histogram uses the interval [50, 100]. While the histogram shows no identifiability on [0, 50], on the interval [50, 100] the method finds the correct parameters with less than 0.1% error on over 95% of the 1000 attempts.

The success of this simple approach to parameter estimation on [50, 100] (or the complete interval [0, 100], not shown) is due to the fact that the SEIR model (2.1) is structurally identifiable from I(t), as long as the peak of the epidemic can be observed. However, one can see that this approach fails on the outbreak part of the epidemic, as shown by the histogram in red. On the time interval [0, 50], the input I(t) is not sufficient for constraining the three parameters.

Figure 2(a) also shows a test of a completely different approach to parameter estimation. We applied MCMC to sample the posterior density of the parameters given the observations, namely $P\left(\beta ,\alpha ,\gamma \mid \text{Δ}{I}_{\text{obs}}(t)\right)$ for t in the same intervals as above. In the deterministic SEIR (used for the MCMC computation of the posterior) the likelihood $P\left(\text{Δ}{I}_{\text{obs}}(t)\mid \beta ,\alpha ,\gamma \right)$ is a product of Poisson densities, which allows easy sampling of the true posterior $P\left(\beta ,\alpha ,\gamma \mid \text{Δ}{I}_{\text{obs}}(t)\right)$. In Figure 2(a) we show the three marginals of the posterior. We notice similar qualitative behavior for this estimator, namely that the parameters are identifiable from the second half [50, 100] of the epidemic(blue curve) but are almost completely unidentifiable from I(t) during the first half [0, 50] (green curve).

Figure 2(b) returns to minimization of the data assimilation error $L(\beta ,\alpha ,\gamma )$ as above but applied to the stochastic SEIR model and using a Poisson-based EnKF. The histogram shows the variation over 1000 different realizations of Poisson noise. For the interval [50, 100], the variation is increased for stochastic SEIR in comparison to the deterministic SEIR, but the estimates are unbiased around the correct parameter settings. For [0, 50], no meaningful estimation occurs.

In summary, for both deterministic and stochastic versions of the SEIR model, both data assimilation-based and MCMC-based algorithms are able to identify the three parameters given I(t) from the time interval [50, 100], and fail on the interval [0, 50]. The intervals [0, 50] and [50, 100] are chosen to be representative of intervals for which identifiability fails and succeeds, respectively. Similarly chosen intervals show the same results; that is, early in an epidemic before the peak is reached, there is a structural reason why the parameters will not be identifiable. We address that reason in the next two sections.

3.1. Asymptotic behavior of linear models.

Consider a linear initial value problem consisting of a vector differential equation $\dot{x}=Ax$, satisfying initial conditions $x(0)=x_0$, where $x=\left[x_1,\dots ,x_n\right]$. Assume A has distinct real eigenvalues.Then solutions are of the form

where ${\lambda }_1>{\lambda }_2>\cdots >{\lambda }_n$ are the eigenvalues of A. Because of the exponential form of the solutions, as t moves away from zero, the solutions begin to closely approximate

Assuming $c_{11}\ne 0$, this means that for each $i,x_i(t)\approx c_i{x}_1(t)$ for some constant c_i.

Example. Consider the linear initial value problem

which we write as $\dot{x}=Ax,x(0)={\left[\begin{array}{c}{E}_0{I}_0\end{array}\right]}^T$, where

Let $A=PD{P}^{-1}$ be the diagonalization, where the columns of P are eigenvectors of A. The diagonalization exists because $\alpha ,\beta ,\gamma >0$ implies that A has distinct real eigenvalues ${\lambda }_1>{\lambda }_2$. The solution is

We can consider separate cases, depending on the constants and B₂. Although the B_i have no particular physical significance, they are formally significant because they represent linear combinations of E₀ and I₀ that grow exponentially with exponent λ_i, respectively. Thus if one of the B_i is zero, the solutions ) and I(t) will evolve exactly proportionally. If both are nonzero, they will still behave asymptotically proportional to each other, with exponent λ₁, the larger eigenvalue.

To be more precise, in what we will call the exactly proportional case, one or both of the B_i are zero. If $B_1=B_2=0$, the solution is identically zero. If one of the B_i = 0 or, equivalently, the $e^{{\lambda }_{i}t}$ term of the solution, is absent, then $I(t)=cE(t)$ for some constant c and for all t.

In what we call the approximately proportional case, both $B_i\ne 0$ and the solution will be

meaning that $I(t)\approx cE(t)$ asymptotically, where $c=I_1/E_1$. Note that in all cases, $I(t)\approx cE(t)$, with the approximation improving exponentially in time.

3.2. Identifiability in linear systems.

A general approach to assessing identifiability in linear systems is suggested in [27]. To search for alternative solutions to (3.1) with the same output I(t) but different E(t) and different parameters $\left({\alpha }^{\prime },{\beta }^{\prime },{\gamma }^{\prime }\right)$, define the coordinate change z = Sx for a nonsingular matrix

Specifically, we seek an S that satisfies

for some F. The new variable z will reproduce I(t) as its second entry, using a “dynamically compensating” F(t) as its first entry, with a different set of parameters, determined below.

This equation is expressible as $\left[\begin{array}{c}\begin{array}{cc}0& 1\end{array}\end{array}\right]Sx=\left[\begin{array}{c}\begin{array}{cc}0& 1\end{array}\end{array}\right]x$. From (3.3), this constraint is

Transposing yields

for all t. Now we split the argument into two cases, depending on the initial conditions (see (3.3)).

Case 1 (approximately proportional). In this case, $B_1\ne 0$ and $B_2\ne 0$. Then for two different times t₁, t₂, the rows of the leftmost matrix in

are linearly independent. Since all matrices on the left side are nonsingular, s₂₁ = 0 and s₂₂ = 1, and therefore,

With this change of coordinates, we can consider the alternative system to (3.1) as $\dot{z}=S\dot{x}=SAx=SA{S}^{-1}z$, where

and where we have set $s_{12}=s_{11}-1$ to match the desired form (3.1). This gives a family of alternative solutions of (3.1) sharing I(t), but with different parameters and different E(t), that are indexed by the single parameter s₁₁. The revised E(t) is $F(t)=s_{11}E(t)+\left(s_{11}-1\right)I(t)$. These solutions exactly match I(t) for all t ≥ 0 and satisfy

The approximately proportional case provides a one-dimensional family of alternative solutions. As promised in [27], these alternative solutions show that in the approximately proportional case, the parameters of (3.1) are unidentifiable from I(t). That is, on the basis of I(t) alone, one cannot distinguish among the infinite set of solutions of (3.5). If our information about the system (3.1) or its parameters is to be inferred from I(t), the existence of multiple solutions consistent with the observations of I(t) will make recovering the parameters effectively impossible.

Case 2 (exactly proportional). Now assume that either B₁ or B₂ is zero. Then $I(t)=cE(t)$ for all t.

The proportionality constant c can be calculated from the equations and depends only on the parameters α, β, γ. Keeping the approximation $S\approx N$ and substituting $I=cE$, we get

which implies

The largest solution c of this quadratic equation is real and positive, assuming that $\alpha ,\beta ,\gamma >0$.

4.1. Unidentifiability in prepeak SEIR.

The SEIR model (2.1) is a coupled set of nonlinear differential equations, but at the beginning of the epidemic, S ≈ N.As the first cases of exposed individuals begin to transition into the infected class, note that the second and third equations approximate a linear system

This approximation was exploited in [13] to derive a formula $R_0=1+(L+D){\lambda }_1+LD{\lambda }_1^2$ for the reproductive number $R_0=\beta /\gamma$ in the case when β is unknown but the latent and infectious periods $L=1/\alpha$ and $D=1/\gamma$ and the exponential growth rate λ₁ from (3.3) can be independently estimated.

According to the previous section, we will observe the asymptotics of the approximately linear dynamics,

for some c as t moves away from 0. In fact, this behavior is apparent in Figure 3(b), which is a magnification of panel (a). The trace of I(t) appears to be a constant proportion of E(t), and this is confirmed in Figure 3(c), where the ratio is plotted versus time.

Figure 4 shows the results of a parameter estimation computation using the data from Figure 3, which sets β = 1.1, α = 0.2, and γ = 0.5. We run data assimilation on the time interval [0, T] using only the daily case numbers $\text{Δ}I(t)$ as input, for various choices of T. To simplify the situation, we will fix the parameter γ = 0.5 to be the exact value and attempt to estimate β and α. We accomplish this by minimizing L(β, α, 0.5) as described in section 2.2.

The function L(β, α, 0.5), sampled at 10,000 random values, is displayed in Figure 4, projected onto the β and α axes, respectively, for ease of analysis. For “prepeak” values of T, the parameters β and α are not well estimated. As T increases and approaches the epidemic peak 60 < T < 80, the parameter estimates gradually become quite accurate. This corroborates our finding in Figure 2 that parameter estimation fails to isolate correct parameters early in the epidemic.

The lesson from Figure 4 is that as the proportion of susceptibles $S(t)/N(t)$ decreases from 1, the error bounds on the parameter estimates will grow. The parameters are identifiable for [0, T] for T well above 50 due to the fact that $S(t)/N(t)<1$, and the parameter estimation will degrade continuously as T is decreased. This degradation is shown explicitly in Figure 4.

4.2. The unidentifiability manifold.

The dynamical compensation results of the previous section explain the phenomenon seen in Figure 4. The unidentifiability manifold, in this case a surface, is plotted in Figure 5. The red dots identify the parameter points $\left({\beta }^{\prime },{\alpha }^{\prime },{\gamma }^{\prime }\right)$ whose evaluated loss function computed on the time interval [0, 50] is in the lowest 1% of points (among 10,000 random points sampled ). The points lie extremely close to the unidentifiability surface (3.8). The wide distribution of the points shows the impossibility of estimating the generating parameter set (β = 1.1, α = 0.2, γ = 0.5) with any accuracy. The color shading on the surface corresponds to reproductive number $R_0={\beta }^{\prime }/{\gamma }^{\prime }$. We note that R₀ is not significantly constrained by the parameters with minimal loss function.

The MCMC approach introduced in section 2.2 shows a similar story. In this case, we use observations of the deterministic model (2.2) and apply MCMC using a single realization of I(t) in the time interval [0, 50] as observable. The true parameters lie inside the envelope of the posterior, as shown in Figure 5(b). The Metropolis–Hastings algorithm within MCMC rejects thousands of proposals that do not lie on the surface and accepts only those that do.

Since the unidentifiability surface is a two-dimensional set, we conclude that even if one of the parameters is known, the other two are not identifiable—the set of possible parameters will only be reduced to a one-dimensional curve. For example, with fixed γ, the data assimilation error on the interval [0, 50] has a poorly defined minimum as a function of (β,α). To illustrate this, we fix $\gamma ={\gamma }^{\prime }=0.5$ in the unidentifiability manifold equation (3.8) to yield the curve ${\alpha }^{\prime }=\alpha (\beta -\alpha /c)/\left({\beta }^{\prime }-\alpha /c\right)$. This curve is plotted in blue in Figure 6(a). The plotted red points are the 1% of (β, α) pairs with smallest values of the loss function. Instead of a localized ball near the true value (β, α) = (1.1, 0.2), there is a curve of pairs equally fitting the observed data, which are therefore indistinguishable from the loss function. These pairs form the flat minima of the loss function seen in Figure 4 for times T preceding the epidemic peak.

Similarly, if we fix a different parameter, we see the same phenomena when trying to estimate the other two parameters. For example, fixing ${\alpha }^{\prime }=\alpha =0.2$, the slice through the unidentifiability manifold (3.8) is ${\gamma }^{\prime }=\gamma +\alpha \left({\beta }^{\prime }-\beta \right)/(c\beta )$, a line. Figure 6(b) shows the line in blue, with the near-minimal pairs of the loss function shown as red dots. Finally, fixing ${\beta }^{\prime }=\beta =1.1$ yields the curve ${\gamma }^{\prime }=\gamma +(\beta -\alpha /c)/(1+\beta c/\left({\alpha }^{\prime }-\alpha \right))$ from the manifold (3.8), shown in Figure 6(c).

On the other hand, fixing the first two parameters on the unidentifiability surface implies that the third can be determined. That is, if we have knowledge of the true α and γ, setting ${\alpha }^{\prime }=\alpha$ and ${\gamma }^{\prime }=\gamma$ in (3.8) implies that ${\beta }^{\prime }=\beta$, so there is a unique solution with those parameter settings. Thus even on the prepeak interval [0, 50] in the example, if α and γ are known, then β is structurally identifiable from the observations of I(t).

Of course, there are many other figures of merit that could be minimized to determine the parameters from the observed I(t); these could be based on data assimilation errors, maximization of likelihood, or some other probabilistic measure. However, during the prepeak part of the epidemic, they will all be susceptible to the alternative solutions that are equally compatible with I(t), implicit in dynamical compensation.

A perhaps more intuitive, if less geometric, view of the unidentifiability surface is that it is the set of parameters for which the leading eigenvalue λ₁ of the resulting system is equal to the λ₁ (see (3.3)) of the underlying system that generated I(t). (In fact, this leads to an alternate derivation of (3.8).) Thus, if we trust the parameter estimation algorithm to return a parameter set that is at least on the unidentifiability surface, then it will have the correct λ₁. Even if the parameters are wrong, this fact can be exploited for uncertainty quantification purposes, as we discuss in the next section.

4.3. Uncertainty quantification.

The unidentifiability surface (3.8) is useful for theoretical reasons to show the impossibility of isolating the original parameter set p from the infinity of other systems that approximately share I(t) during the beginning portion of an epidemic. Next, we suggest that it may be useful in practice for uncertainty quantification.

It turns out to be a helpful fact that the unidentifiability surface generated by an arbitrary parameter set p indexes the set of parameter sets that share the observed I(t). Assume that we use a parameter estimation algorithm with input I(t), and estimate the parameter set as p′ and that lies on the surface. The roles of p and p′ are symmetric, so we can also consider that p lies on the unidentifiability surface generated by p′. That means we can reverse the roles: switch the primed and unprimed variables in (3.8), noting that c must be replaced by c′ computed from (3.7) with unprimed variables replaced with primed variables.

For instance, assume the correct parameters are $p=(\beta ,\alpha ,\gamma )=(1.1,0.2,0.5)$ but that a parameter estimation algorithm instead returns, for example, an estimate $p^{\prime }=\left({\beta }^{\prime },{\alpha }^{\prime },{\gamma }^{\prime }\right)=(0.852,0.25,0.4)$ that lies on the unidentifiability surface. The set p′ given here is just for illustration; in this case it was chosen by making an arbitrary choice of α′ and γ′ and then computing the corresponding β′ lying on the surface (3.8). Next, we ignore the origin of p′ and consider what we can infer from it. In Figure 7(a), we produce 30 trajectories of the stochastic SEIR by perturbing p′ 10% to new values $p^{″}=\left({\beta }^{″},{\alpha }^{″},{\gamma }^{″}\right)$. We have overlaid as a yellow curve the original trajectory that produced I(t), generated by the parameters p. There is a large amount of variability in the 30 trajectories.

Figure 7(b) shows trajectories of 30 stochastic SEIR systems, where we have randomly changed α′ and γ′ by 10% to α″ and γ″, but this time we computed the″ that lies on the surface. We reiterate that the surface, being the corresponding β unidentifiability surface of p′, can be computed from p′ and is therefore known to us, even if the original p is unknown. The ensuing trajectories are much more faithful to the original system, given that they share the leading dynamical eigenvalue λ₁. Thus, even starting with a mildly incorrect parameter set p′, by querying nearby points p″ on its unidentifiability surface, we see reasonable facsimiles of the underlying dynamics generated by the original parameters p.

Note that there are limitations on how far the incorrect parameters p′ can be from the original parameters p in order for the trajectories produced in this way to be representative of the original systems. In particular, the constant c in the proportionality $I(t)\approx cE(t)$ is in general different for the new system, and so its trajectories will be different. Our informal observation is that if the alternative parameters are within about 20% of the originals, the approximating trajectories may still be useful for uncertainty quantification.

This observation opens up the possibility of using the unidentifiability surface for uncertainty quantification purposes by studying the spread of nearby solutions as a function of uncertainty in the parameters. If an uncertainty in the estimate can be determined from the algorithm generating the estimate, bootstrapping techniques can be used to move along the surface (3.8) and quantify the variance of key aspects of the family of nearby trajectories. We leave a more complete analysis of this phenomenon, and its possible utility to forecasting, to future investigation.

5.1. The SEUIR model.

The epidemic model in [18] was used to represent populations in a specific city and included extra external inputs from other cities.The underlying SEIR-style model is

with constant total population $N=S+E+U+I+R$, where $0<\alpha <1$. The new variable U represents unreported infected cases, while I is reserved for reported infected cases. For SEIR, we will consider I(t) as the observable variable.

For simplicity, we rewrite the parameters as $z=1/Z$, $d=1/D$, $w=\alpha /Z$ to arrive at the equivalent but more user-friendly system

where $N=S+E+U+I+R$, with parameters $\beta ,z,w$, and $d,0<w<z$, which we call the SEUIR model. Figure 8 displays a sample trajectory of the SEUIR model.

5.2. Unidentifiability in SEUIR.

Again consider the prepeak portion of the epidemic, where $S\approx N$. Then there is an approximating linear system

which will exhibit dynamical compensation. Given our experience with SEIR, we consider solutions of (5.3) where E, U, and I are proportional, say $U(t)=bE(t)$ and $I(t)=cE(t)$. One checks that if $E(t),U(t),I(t)$ are such solutions, then $E(t)=E_0{e}^{[\beta (b+c)-z]t}$, where

It will be convenient in proving the lemma below to note the identities