Structural identifiability analysis of age-structured PDE epidemic models.

Computational and mathematical models rely heavily on estimated parameter values for model development. Identifiability analysis determines how well the parameters of a model can be estimated from experimental data. Identifiability analysis is crucial for interpreting and determining confidence in model parameter values and to provide biologically relevant predictions. Structural identifiability analysis, in which one assumes data to be noiseless and arbitrarily fine-grained, has been extensively studied in the context of ordinary differential equation (ODE) models, but has not yet been widely explored for age-structured partial differential equation (PDE) models. These models present additional difficulties due to increased number of variables and partial derivatives as well as the presence of boundary conditions. In this work, we establish a pipeline for structural identifiability analysis of age-structured PDE models using a differential algebra framework and derive identifiability results for specific age-structured models. We use epidemic models to demonstrate this framework because of their wide-spread use in many different diseases and for the corresponding parallel work previously done for ODEs. In our application of the identifiability analysis pipeline, we focus on a Susceptible-Exposed-Infected model for which we compare identifiability results for a PDE and corresponding ODE system and explore effects of age-dependent parameters on identifiability. We also show how practical identifiability analysis can be applied in this example.

Introduction

Identifiability analysis addresses the question of whether the parameters of a mathematical model can be uniquely identified, or otherwise, given the system observations. There are two general types of identifiability analysis: structural (assuming data to be noiseless and arbitrarily fine-grained) and practical (using real data). These analyses are crucial to interpreting biologically relevant predictions from computational and mathematical models that rely on parameter values that are estimated from experimental datasets. Previous studies have shown that unidentifiable models can lead to vastly different predictions for different parameter values that are able to produce the same observed data (Kao and Eisenberg 2018). In the context of biological models, parameter identifiability gives increased confidence in model predictions that are based on estimated parameters. This is particularly relevant for predictions applied to topics such as diseases, treatment outcomes, vaccination efficacy, and future infection dynamics. Specifically, the use of epidemiological models has become increasingly necessary and useful for predicting disease outcomes, as we have seen in the case of COVID-19, and these efforts should be as accurate and confident as possible to produce the most useful predictions. Thus, we focus our work in the context of epidemic models due to their use in many diseases and to compare with corresponding parallel work previously done for ODEs regarding structural identifiability.

Many theoretical results and algorithms for structural identifiability of linear and non-linear ordinary differential equations (ODEs) models have been derived (Miao et al. 2011; Audoly et al. 2001; Ljung and Glad 1994; Pohjanpalo 1978; Chappell et al. 1990; Cobelli and DiStefano 1980; Chis et al. 2011) and have been applied to study an array of diseases. There are models for the spread of cholera (Eisenberg et al. 2013), vector-borne diseases (Kao and Eisenberg 2018; Zhu et al. 2018; Tuncer et al. 2016), and general infectious diseases (Tuncer and Le 2018; Evans et al. 2005, 2004). Methods for structural identifiability are often analytical in nature (e.g. the commonly used differential algebra approach (Audoly et al. 2001; Ljung and Glad 1994)), but may also be numerical, or some blend of the two (Jacquez and Greif 1985; Jacquez and Perry 1990; Raue et al. 2009). Numerical approaches to structural identifiability typically provide insight at the local, rather than global, level.

ODE models are not always sufficient for representing disease dynamics when, for example, age is an important factor. In many diseases, different age groups may respond differently to diseases or disease interventions based on differences in susceptibility (Gardner 1980), mortality (Castillo-Chavez et al. 1991), and contact rates (Mossong et al. 2008; Glasser et al. 2012). Age preferences have been shown to have a significant effect on mixing patterns that drive disease spread (Valle et al. 2007). This type of information can be incorporated into a mathematical model for epidemics, for example, by explicitly representing age-structured distributions of the susceptible and infected subpopulations. Age-structured models have been used to evaluate and optimize age-targeted vaccination strategies for several diseases (Harris et al. 2019; Keeling and White 2011; Hao et al. 2019; Castillo-Chavez and Feng 1998; Shim et al. 2006), and it has been shown that population heterogeneity and non-random mixing can have a substantial effect on the prediction of disease outbreaks and herd immunity (Glasser et al. 2016; Feng et al. 2015).

In a continuous setting, age structure is introduced by creating state variables that are functions of both age and time, and also by introducing transport terms to represent constant aging of a population. These model formulations result in PDEs, rather than ODEs, making them more biologically accurate; however, such models also typically more complex to analyze. There are other PDE formulations of epidemic models such as age of infection models (Thieme and Castillo-Chavez 1993; Inaba and Sekine 2004) or spatially explicit models (Bertuzzo et al. 2009; Wu 2008). The focus of this paper will be specifically on age-structured models. While there are alternative approaches to incorporating age structure, such as to use compartmentalization to stratify populations (Chow et al. 2011; Ferguson et al. 1996) or to use a stochastic discrete representation such as an agent-based model (Ajelli et al. 2010; Bian 2004), we focus here on the PDE formulation. Importantly, the issues of both structural and practical parameter identifiability for these types of models still needs to be elucidated.

While structural identifiability analysis has been extensively applied in the context of ODE epidemic models, it has not been widely explored for age-structured PDE models (Evans et al. 2004, 2005). PDE models present additional difficulties due to the increased number of variables and derivatives as well as the presence of boundary conditions. Identifiability of age-structured PDE models was first considered by Perasso et al. (2011) for a simple SI-type model and was again considered by Perasso and Razafison (2016) for a population dynamics model including only birth and death. These analytical results have not yet been extended to more comprehensive PDE epidemic models, though some identifiability results have been shown numerically (Tuncer et al. 2016). For these types of problems, modelers have relied on practical identifiability analyses which are dependent on either available data or data that has been synthesized for the purpose of analysis.

In this paper, we derive analytical identifiability results for more complex age-structured epidemic models using a differential algebra framework that is simpler to apply than the previous (aforementioned) previous approaches. We apply this methodology to a compartmental epidemic Susceptible-Exposed-Infected (SEI) model both with and without age structure explicitly modeled, and we compare the identifiability results under different conditions such as the presence of immigration in the system or for age-dependent model parameters. To corroborate our structural identifiability results, we also perform a practical identifiability analysis for one version of the model. While we focus on age-structured models, the methods we present herein can be applied to other PDE formulations, with certain assumptions about smoothness and uniqueness of solutions.

Methods

Identifiability framework

To consider parameter identifiability of a model, we must first define the model structure. From an identifiability or parameter estimation perspective, a model is composed of two pieces: (1) a state model or set of state equations, which describe the underlying system of interest, and (2) a measurement model or set of output equations, which describe the observation or measurement of interest (i.e., the quantities of interest for which one either has or plans to collect data). In this study, we assume the following general first-order partial differential equation (PDE) model structure, which is commonly used in age-structured PDE models in biology and health:where is a vector of state variables describing the system of interest, is the vector of parameters, t and a are two independent variables. The model may include initial and boundary conditions that are generic (unknown), written in terms of input parameters, or known. In the following definitions, we assume generic initial and boundary conditions. Non-generic initial conditions can be addressed as detailed in Pia Saccomani et al. (2003), and we assume that non-generic boundary conditions can be addressed analogously. Here we consider time and age as the independent variables, although the approach given here will also work for more general variables, as well as for additional derivatives of further independent variables on the left hand side. Additionally, represents the vector of model outputs (measured variables) and represents any known inputs (e.g. forcing functions) to the system, if they exist, which are assumed to be arbitrarily differentiable (or at least sufficiently differentiable for determining the input-output equations, described below). Finally, f and g are rational functions of the model variables and parameters such that the model can be written equivalently as differential polynomials by clearing denominators.

In structural identifiability analysis, it is assumed that the inputs and outputs are known (observed) perfectly – in this case, we assume that u, and y are known for all times and all ages and that there is no noise in the data. Similarly, the independent variables t and a are assumed to be known. Thus, the resulting identifiability results can be interpreted as the upper limit of what can possibly be identified with perfect data collection methods. However, structural identifiability is a necessary condition for identifiability in the realistic data case, making it an important first step in the parameter estimation process.

A model such as (1) can be thought of as a map from parameter space (e.g. for ) to output space (the space of trajectories for ), often termed the model map . With this framing, the structural identifiability question becomes that of whether or not the model map is injective. In particular, we define structural identifiability as follows:

Definition 1

For a model of the form (1), an individual parameter p in is uniquely (or globally) structurally identifiable if for almost every value and almost all initial conditions, the equation implies . A parameter p is said to be non-uniquely (or locally) structurally identifiable if for almost any and almost all initial conditions, the equation implies that p has a finite number of solutions.

Definition 2

Similarly, a model of the form (1) is said to be uniquely (respectively, non-uniquely) structurally identifiable for a given choice of output if every parameter is uniquely (respectively, non-uniquely) structurally identifiable, i.e. for almost every value and almost all initial conditions, the equation has only one solution, (respectively, finitely many solutions). Equivalently, a model is uniquely (respectively, non-uniquely) structurally identifiable for a given output if and only if the map is injective (respectively, finite-to-one) almost everywhere, i.e. if there exists a unique set of parameter values which yields a given trajectory almost everywhere.

Definition 3

A model of the form (1) is said to be structurally unidentifiable for a given choice of output if it is not uniquely or non-uniquely identifiable, namely if at least one parameter is structurally unidentifiable, i.e. for a set of and initial conditions of nonzero measure, the equation has infinitely many solutions. Equivalently, a model is structurally unidentifiable for a given output if and only if the map is infinite-to-one on a set of positive measure.

We note that because some specific parameter values or initial conditions may render an otherwise identifiable model unidentifiable (e.g. if all parameters or initial conditions are zero), we define structural identifiability generically, i.e., for almost all parameter values and initial conditions (Audoly et al. 2001; Pia Saccomani et al. 2003). Because the definitions apply for almost all parameter values and initial conditions, there may be some initial conditions of , or trajectories of (corresponding to a measure zero set of parameters and initial conditions) for which an otherwise identifiable model may be unidentifiable. We additionally note that the system’s initial and boundary conditions may also depend on unknown parameters.

Overview of the differential algebra approach

To determine structural identifiability of a model with defined outputs, we extend the differential algebra approach (Ljung and Glad 1994; Audoly et al. 2001) to the age-structured PDE case. This approach relies on converting the model system into a system of input-output equations, that is, a set of monic polynomial equations in terms of the known variables , , and their derivatives, with rational coefficients in the parameters . The input-output equations represent an implicit form of the model map, and as a set of differential equations, they are input-output equivalent to the original system. This means that for the same inputs and parameter values, the equations generate the same output trajectories (Eisenberg 2019). We then consider the map from the input parameter space to the coefficients of the resulting monic polynomial.

A key insight of the differential algebra approach in the rational ODE case is that injectivity of the model map can be evaluated by considering the coefficient map, i.e. the map from parameters to coefficients of the input-output equations (Audoly et al. 2001; Ljung and Glad 1994). In the rational ODE case, the coefficient map typically provides complete identifiability information for the model. If this map is injective, then the model is uniquely structurally identifiable. If the model is unidentifiable, this map determines the identifiable combinations of parameters. For example, while some model parameters may not be independently identifiable, the product or sum of two or more model parameters may be identifiable. This may enable the model to be reparameterized in a way that allows the model to become structurally identifiable (Chappell and Gunn 1998; Cole et al. 2010; Meshkat et al. 2009, 2011, 2014; Meshkat and Sullivant 2014). It should be noted that a set of identifiable combinations is not necessarily unique – there may be many sets of parameter combinations that are equivalent. For example, the sets of identifiable combinations and are equivalent since .

The differential algebra approach for ODE models typically uses characteristic sets to calculate the input-output equations (Audoly et al. 2001; Ljung and Glad 1994). Characteristic sets are a differential analog of the Gröbner basis (Ritt 1950) based on differential polynomial psuedodivision. For more details on the characteristic set approach, the reader is referred to (Audoly et al. 2001; Ljung and Glad 1994; Pia Saccomani et al. 2003; Hong et al. 2020; Eisenberg 2019).

This methodology has previously been applied to ODE systems, where derivatives are taken with respect to only one variable. In this work, we generalize this methodology to PDE systems of the form in (1); in particular, we consider age-structured systems where derivatives are taken with respect to both time and age. We show that the same framework can be applied to these systems, with the additional consideration of boundary conditions. In age-structured systems, a boundary condition must be specified at , usually in the form of a Dirichlet boundary condition. Thus, in addition to performing identifiability analysis on the model equations, one must also consider the equations evaluated at to explicitly include information at the boundary. This is analogous to methods for determining identifiability of initial conditions in ODE systems (Pia Saccomani et al. 2003).

Extending the differential algebra approach to age structured PDE models

In this section, we discuss our generalization of the differential algebra approach to the PDE case. For PDE models, the primacy of the differential ideal generated by the model and measurement equation differential polynomials is not necessarily guaranteed (at least, to our knowledge it has not been shown to be the case), meaning that the characteristic set may not be unique. Thus, instead of a characteristic set based approach, we will use a simpler, substitution-based approach based on previous work by Eisenberg (2019) (similar to that used in Eisenberg et al. (2013)). In this approach, we use basic substitution (which, as long as one is careful not to divide by quantities that could be zero, etc., does not change the underlying solutions to the system) to eliminate the unobserved variables, resulting in a monic, reduced system of differential equations which is input-output equivalent to the original system (Eisenberg 2019). This input-output equivalent system has equivalent identifiability properties to the original, allowing us to use it to evaluate identifiability of a certain class of PDE systems.

As is typical in ODE systems, once we have a set of input-output equations (i.e. monic differential polynomials only in terms of the known variables and their derivatives with respect to either independent variable and the parameters ), we assess identifiability based on the coefficients of the monomial terms. As in the ODE case, this requires that the input-output equations be in some sense minimal or fully reduced, with the monomials of the input-output equations linearly independent (Ovchinnikov et al. 2021; Bearup et al. 2013). Given these reduced input output equations, we assume that our observed trajectory contains sufficiently many algebraically independent points that we can successfully identify the coefficients of the input-output equations. In other words, from a given observed solution of over age and time, we assume the observed solution contains an arbitrary number of distinct points at which the input-output equations can be evaluated to solve for the input-output coefficients (in particular, more distinct points than the number of coefficients). In practice most trajectories meet this criterion, although some care must be taken to avoid systems or measure zero sets of initial conditions for which the solutions are at steady state, or have key parameters equal to zero, etc. (Hong et al. 2020). With the coefficient-identifiability assumption (Eisenberg 2019), the rest of the differential algebra approach follows as with the ODE case, and is described below.

Procedurally, the process of performing identifiability analysis for PDEs is almost identical to the ODE case, though substitutions may be made more complicated by the presence of an additional derivative. To obtain input-output equations from model equations, a typical first step is to solve the output equations for state variables and, using substitution, eliminate these state variables and their derivatives from the model equations. For simplicity, all nonzero terms in the model equations are moved to the right-hand side. Remaining state variables can then be eliminated by solving one or more of the model equations. The particular substitutions required will of course vary depending on the model and output equations. The goal is to obtain a system of polynomial equations in terms of the outputs and their derivatives, with all state variables eliminated, whose coefficients are in terms of the model parameters. In the case that square roots arise, for example, they must be eliminated by multiplying by the conjugate.

Once input-output polynomial equations are obtained, we divide each equation by the coefficient of its leading term to obtain a monic polynomial (see below for details about the ranking of terms). From these monic polynomials, we obtain the set of polynomial coefficients that defines the coefficient map. This map, as in the ODE case, can be used to determine complete identifiability information for the model equations.

Note that thus far, we have not considered any information about the initial or boundary conditions—an important consideration of identifiability, In particular, analysis of the boundary conditions is a key component of identifiability analysis for PDEs that does not occur in the ODE case. In the case of age-structured systems, there is a single boundary condition at . Thus, we treat the boundary condition as similar to an initial condition for ODEs (see e.g. the methods in Pia Saccomani et al. (2003) for the ODE case). If information is known about the boundary conditions, such as values or functional forms, this information can be substituted into the model equations evaluated at to obtain boundary condition equations. Since the values of model outputs are known at , the same identifiability analysis as above can be applied to boundary condition equations. This may reveal additional identifiability information about the model or any parameters that are present in the boundary conditions.

We note that, as in the initial condition case, if identifiability is evaluated without considering boundary conditions, this implicitly assumes that the boundary conditions are generic and unknown. However, there are some important aspects of the boundary conditions which are more subtle than the ODE case of initial conditions, and may warrant further exploration, such as that boundary conditions still may vary over time. Exploring how different boundary condition types may affect identifiability would be an important direction for future work.

Our procedure for performing structural identifiability analysis, which can be applied to both ODE and age-structured PDE systems, is outlined in Fig. 1. In our structural identifiability analyses herein, we utilize Mathematica to perform the required algebraic manipulations and evaluate map injectivity. All Mathematica notebooks are provided in our code repository at https://github.com/epimath/age-pde-identifiability.

Fig. 1

Flowchart describing the methodology for identifiability analysis of ODE and age-structured PDE models based on using a differential algebra approach

Practical identifiability analysis

Practical identifiability was evaluated using the Fisher Information Matrix (FIM) and profile likelihoods (Jacquez and Perry 1990; Raue et al. 2009; Venzon and Moolgavkar 1988; Eisenberg and Hayashi 2014). Here, we follow the methodology laid out in (Eisenberg and Hayashi 2014).

For ODE models when evaluating identifiability, the FIM is computed by , where X is the sensitivity matrix defined by for parameters and time points . This form is a simplified form of the FIM for the case of normally distributed measurement error, but it is frequently used for evaluating structural identifiability due to its easily interpretable form and clear connection to the sensitivity matrix X (the two matrices share the same rank, and identifiable parameter combinations can be naturally understood through the sensitivity matrix as compensation between parameters will be reflected in dependences in the columns). For an age structured system, we define where the sensitivity matrix is defined by where for are pairs of age and time points that form a grid. Here, y denotes the observable model output. For multiple model outputs (e.g., and ), the matrix X can be defined by concatenation (i.e., ). The rank of F indicates the number of identifiable parameters and parameter combinations. If F is singular, the model is structurally unidentifiable; if F is close to singular (meaning a large condition number or a determinant near zero, e.g. below machine precision), the model may be practically unidentifiable. The FIM can also be used to determine which parameter subsets are identifiable or unidentifiable (Eisenberg and Hayashi 2014).

We compute profile likelihoods by fixing a particular parameter and fitting all remaining parameters to maximize a likelihood function. We define the likelihood by assuming constant, normally distributed measurement error with standard deviation equal to 10% of the mean of the data, and mean equal to the model trajectory. Maximizing the likelihood function is equivalent to minimizing a cost function based on the negative log likelihood (in this case equivalent to least squares). This is performed at many values across the range of values for . The profile likelihood for is given by the maximum values of the likelihood function across the range of values for . If the profile likelihood is flat, is structurally unidentifiable; if it is sufficiently shallow, is practically unidentifiable. For more details on profile likelihood approaches to identifiability, the reader is referred to (Raue et al. 2009; Venzon and Moolgavkar 1988).

In our analyses herein, we solved the age-structured PDEs numerically using Method of Lines, using a first order finite difference approximation of the spatial derivative, and using Matlab’s ODE solver ode45.

SEI model for TB epidemiology

In several case studies below, we consider a simple mathematical model of the epidemiology of tuberculosis (TB). TB is transmitted via the aerosolized bacterium Mycobacterium tuberculosis from person to person via lung inhalation. TB is unique from many infectious diseases in that most people develop a latent infection, meaning that they present no symptoms and cannot transmit the disease, and the infection can remain latent for an extended period of time. Thus, our model consists of three mutually exclusive compartments: Susceptible (S), Exposed/latent (E), and Infectious (I). Through contact with infectious individuals, susceptible individuals (called susceptibles) may progress to the exposed class (i.e. latent infection) or the infectious class (i.e. primary infection). This is slightly different from the standard SEI model, in that some individuals bypass the exposed class and immediately become infectious, but is similar to what has been used to model tuberculosis epidemics previously (Ozcaglar et al. 2012; Guzzetta et al. 2011; Renardy and Kirschner 2019; Vynnycky and Fine 1997). Individuals who are exposed/latent may progress to the infectious class where they can transmit bacteria. We assume all individuals are born into the susceptible class, and each compartment has a compartment-specific death rate.

A schematic for this model is shown in Fig. 2. This schematic can be translated into a mathematical or computational model by representing it as a system of ordinary or partial differential equations (as has been done by, e.g., Guzzetta et al. (2011); Renardy and Kirschner (2019); Vynnycky and Fine (1997)), a stochastic model, or an individual-based model (as has been done by, e.g., Guzzetta et al. (2011); Renardy and Kirschner (2020)). Here, we will first consider the model as a system of ODEs and then extend the results to age-structured PDEs.

Fig. 2

SEI diagram and parameter definitions for a mathematical model of TB epidemiology. S represents the susceptible population, E represents the exposed/latent population, and I represents the actively infectious population. Arrows describe flow between mutually exclusive compartments

Case study 1: density-dependent transmission

For simplicity, we first demonstrate the identifiability analysis framework on a model with density-dependent transmission, i.e. transmission is proportional to SI rather than . In this case, the per-individual rate of contact increases as the population density increases. This is not necessarily a realistic model for TB, since TB is often spread only among close contacts, but simplifies the differential algebra for the purposes of demonstration.

In order to perform parameter identifiability analysis, we must first define what the outputs (i.e., the quantities for which we could obtain reliable, biologically relevant time series data) of the specific model are. For all of the following examples, we will assume that we can observe a fraction of the exposed individuals, i.e. where is the reporting/detection rate. This type of data could be obtained by, for example, contact tracing to identify individuals who have been exposed to TB disease, or testing a proportion of the population for exposure. This is a useful practice even in areas with a low TB burden such as the US, since detection and treatment of latent infection is crucial for TB elimination (LoBue and Mermin 2017). Since latent disease is asymptomatic, it is highly unlikely that all exposed/latent individuals will be identified, and thus we assume instead that only a fraction of the exposed individuals are identified.

In addition, we may be able to observe all or some of either all active infections (prevalence) or new active infections (incidence). The relevant output function depends on the reporting and surveillance data for the region of interest. For example, if all infectious cases are monitored and reported then these data would represent total prevalence. If not all cases are detected due to misdiagnosis, for example, then the data would represent only a fraction of the total prevalence. If new infectious cases are reported upon diagnosis but then lost to follow-up, then incidence is an observable output while prevalence is not. Again, since it is possible that not all infected individuals will be diagnosed, we may only be able to observe a fraction of incidence. We present here results for observing prevalence over time, but the identifiability analysis of model parameters for incidence would follow similarly.

We assume that , i.e. that we are able to obtain data for a fraction of the total number of infections at any given time. We assume no other information about the system parameters or initial conditions. The case of incidence is more complex due to nonlinearities, but the methodology presented here can still be applied, and we provide such an example in our code repository.

ODE SEI model with density-dependent transmission

Though identifiability analysis for ODE systems is well established, we include the technical details here for completeness and for comparison when we move to the PDE system. The model can be written as a system of ODEs:We consider the case of generic (unknown) initial conditions. We do not intiially assume a fixed population, and thus the birth rate and death rates are independent. We begin by rewriting the system of Eq. (2) in terms of the outputs and . For ease of notation, we adopt the notation . Using the substitutions and , we obtain the following equations:It now remains to eliminate S from these equations so that we obtain equations in terms of solely the model parameters and outputs. This can be done by solving for S in the equation for :We note that we divide by in this example for simplicity; however the same input-output equations can be obtained without dividing, using Ritt’s pseudodivision (Ritt 1950) (just with more involved algebraic computations required). This expression can then be substituted into the remaining two equations to obtain two input-output equations that are rational equations in terms of the outputs and their derivatives. By multiplying by a common denominator and moving terms to one side, these input-output equations can be rewritten as the following polynomial equations.Dividing the first equation by and the second equation by to form monic polynomials, we get the following set of polynomial coefficients.We consider this set of coefficients as a function of the model parameters, and we seek to determine whether this map is injective. If it is not, then we consider whether it is injective as a function of some parameter combinations. We find that and are structurally identifiable. We then simplified the above set by removing duplicates up to sign (e.g. and are redundant), subtracting off known quantities where we can (e.g. if is identifiable and and are identifiable, then so is , a simpler identifiable combination), and identifying coefficients that are combinations of previously established known quantities (e.g. can be found by dividing coefficients and ). Through simplification and substitution of the above set, we find that other identifiable parameter combinations areIf we further assume that , then and also become identifiable while , c, , and remain unidentifiable. However, the following combinations are identifiable: . Consequently, if any one of , c, , or is known, then the rest become identifiable.

Age-structured PDE model with density-dependent transmission and constant parameters

Now that identifiability results are established for the ODE representation of the SEI model 2, we extend these results to a more complex PDE model that captures the same epidemiological setting (SEI), but includes age structure. We incorporate age structure into the model by adding a transport term to capture population aging. To simplify the model, we assume that transmission occurs between individuals of the same age, so that the transmission term is proportional to S(t, a)I(t, a) for any given age a and time t. While this is a significant simplification, data has shown that people tend to interact primarily with others in their age group (Mossong et al. 2008). The model thus becomesfor and . Note that all state variables are now functions of both time and age. The boundary conditions are , , and we begin by assuming initial conditions are generic. Note that the birth rate is now incorporated into the boundary condition rather than into the differential equation for S. We begin by assuming that model parameters are constants, independent of age.

Identifiability of model equations with generic boundary conditions

We first consider identifiability of the model parameters without explicitly accounting for the presence of boundary conditions, meaning that we assume there is no a priori knowledge about the boundary conditions and that they cannot be observed. We call this case generic boundary conditions. To perform the identifiability analysis, we utilize the same framework as we did for the above ODE system analysis. That is, we rewrite the model system as a set of polynomial equations in terms of the outputs and their partial derivatives, and then consider the map from parameter space to the polynomial coefficients.

We adopt the notation . Using the substitutions and , as in the ODE case, we obtain the following equations:We now eliminate S to obtain equations in terms of only the model parameters and outputs. This can be done by solving for S in the third equation:We assume here that generically to avoid division by zero, although one can also obtain the same input-output equations without dividing by by solving instead for and substituting that into the equations (but this requires more extensive algebraic calculations and multiplication by for the first equation). As before, this expression can then be substituted into the remaining two equations to obtain two input-output equations that are rational equations of the outputs and their derivatives. By multiplying by a common denominator and moving terms to one side, these input-output equations can be rewritten as the following polynomial equations.We choose, without loss of generality, derivatives with respect to time to be higher ranked than those with respect to age for our ranking of terms within the polynomial, and let be higher ranked than , yielding the ranking . Dividing the first equation by and the second equation by to form monic polynomials, we get the following set of polynomial coefficients.We find that, as was the case for the ODE model, , and are structurally identifiable. Through simplification and substitution, we find that other identifiable parameter combinations areIf then and become identifiable as well as the combinations .

Incorporating non-generic boundary conditions

In our model, the boundary conditions are in fact not generic, but rather are known to be , where c represents the birth rate. Here, the boundary condition contains an unknown model parameter, c. To determine the identifiability of the boundary condition parameter c, we perform the same analysis as above but with the model equations evaluated at (boundary condition equations). However, to avoid reducing equations to zero, we let the boundary conditions be set to an extremely small but assumed known quantity that should not affect the model simulations significantly but enables us to evaluate structural identifiability in this case (e.g. reflecting a small amount of congenital tuberculosis in the population (Cantwell et al. 1994)). Then, noting that S(t, 0) is constant and thus , the relevant equations areBy transforming the system (4) to a set of polynomial equations we obtain the following polynomial coefficients:Since is identifiable from the model equations, we thus find that is an identifiable combination, which is consistent with the ODE model. Thus, the full set of identifiable parameters and combinations for this model iswhich is equivalent to the set of identifiable combinations for the ODE model.

Identifiability of initial conditions

With generic initial conditions, we do not gain any additional identifiability results. However, if we assume that certain information is known, such as the total population size N at the initial time, then further model parameters can become identifiable. Here, we assume that the initial population size as well as its derivatives with respect to time and age are known for all ages (noting that the age-specific population could potentially be determined from census data, and derivatives perhaps roughly estimated from multiple years of census data). Under these conditions, we may solve for the initial susceptible population S(0, a) as , where N represents the total population size. Similarly, since we assume the derivatives of N are known at , we may solve for the derivatives of S at in terms of the derivatives of E and I. It is thus straightforward to solve for S and its derivatives at in terms of the outputs and .

To determine parameter identifiability from the initial condition with this additional knowledge, we evaluate the model equations at and perform identifiability analysis analogous to the previous analysis for the boundary conditions. We obtain the following set of polynomial coefficients:From these coefficients, we find that and are all identifiable. Combining this with the previously established identifiability of from the boundary equations, we find that all parameters are identifiable if the initial population size and its derivatives are known.

Case study 2: frequency-dependent transmission

While density-dependent transmission is convenient for simplicity, it is not generally a realistic model for TB. Rather, most models for TB epidemiology use frequency-dependent transmission models wherein the transmission term is proportional to . In this section we present identifiability results for the SEI model with frequency-dependent transmission, assuming that our output functions are and . In addition to the ODE and constant-parameter PDE models, we will also consider variations including immigration and age-dependent death rates. The full details of the analysis are cumbersome and thus are not provided here, but can be found in the code repository.

ODE SEI model with frequency-dependent transmission

For completeness, we include here identifiability results for the ODE version of the model. This model can be written as the following system of equations:where . We find that the set of identifiable parameters and parameter combinations for this model is

Age-structured PDE SEI model with frequency-dependent transmission and constant parameters

with boundary conditions , . From the model equations, we find that the parameters and combinations are identifiable. We further find from the boundary equations that is identifiable. Thus, the full set of identifiable parameters and combinations for this model iswhich is consistent with the ODE model. If the initial age distribution of the total population and its derivatives were known (as in Sect. 4.2.3 for the previous case study), then all parameters in the model would be identifiable.

SEI model with frequency-dependent transmission and immigration

The models considered thus far have included influx of new individuals only via birth. However, it may be necessary to also consider immigration when studying epidemics in fluid populations. If immigration is included in the model, with the assumptions that the immigration rate is constant over time and all immigrants are susceptible, the ODE system becomeswhere m is the immigration rate. If we assume that the age distribution of immigrants is known, then the corresponding age-structured PDE system becomeswhere is some known function such that and , with boundary conditions and .

The identifiability results for the ODE system (7) remain the same as before but with c replaced by , i.e., the identifiable parameters and combinations areFor the PDE system (8), we get the following identifiable parameters and combinations from the model and boundary equations:Thus, in this case the identifiability results for the PDE model differ from those for the ODE model in that the combinations and are both identifiable, whereas in the ODE model only the sum is identifiable. In other words, if all immigrants are susceptible then in an age-structured model it is possible to differentiate between birth and immigration, whereas in an ODE model it is not.

Age-structured PDE SEI model with frequency dependent transmission and age-dependent parameters

A major benefit of age-structured models is the ability to incorporate age-dependent parameters such as death rates, contact/transmission rates, and reactivation rates to make the model more realistic. Here, we consider the case of age-dependent death rates. There are several ways to introduce age dependence to the death rate: one could assume a specific functional form (e.g., piecewise-constant, linear, exponential), that the function belongs to a particular family (e.g., polynomials, splines), or that the function is completely arbitrary with some conditions on smoothness. In this section, we present identifiability results for four cases: (1) piecewise-constant death rate functions for all three compartments, (2) a single exponential death rate function for all compartments, (3) a single polynomial death rate function for all compartments, and (4) a single arbitrary function for all compartments. For simplicity, we will consider the SEI model without immigration.

PDE SEI model with piecewise-constant death rates case

We first consider the case that the death rates are piecewise-constant functions of age with discontinuities at ages , i.e. for and . This represents the case where death rates are known independently for different age groups. Then the system (6) becomes a sequence of PDE systems obeying the equationsfor and , , with boundary conditions , andfor . We will assume that the locations of the discontinuities are known, and that the sub-intervals are wide enough that the usual structural identifiability results hold.

Since the model equations remain essentially unchanged, structural identifiability of the model equations has already been established in Sect. 5.2. Thus, from the model equations we obtain thatare identifiable for all i. From the boundary equations at we again get that is identifiable, as in the constant parameter case. We gain no additional information about the parameter values from the boundary conditions at for .

PDE SEI model with exponential death rate case

One may also wish to use a continuous function to represent age-dependent death rates. One possibility is to use an exponential function so that death rate increases exponentially with age. For simplicity, we assume that death rates are the same for all three compartments. Then the model system can be written asfor and , with boundary conditions , , and . We write the system in this way rather than directly inserting so that we maintain a differential polynomial form.

To obtain an input-output equation from the model equations, we use the I equation to solve for S and the E equation to solve for (note that these choices are somewhat arbitrary). Then by noting that and and substituting into the S and equations, we obtain rational equations in terms of the outputs and their derivatives. Thus, by multiplying by a common denominator, we obtain differential polynomials. We find that are the identifiable parameters and parameter combinations.

By performing similar analysis on the boundary equations, we find that and are identifiable. Thus, the full set of identifiable parameters and combinations for this version of the model isThis shows that both parameters from the exponential death rate are identifiable. The identifiability of the remaining parameters is consistent with the results for the ODE model in Sect. 5.1.

PDE SEI model with polynomial death rate

Another option for a continuous death rate is to use a polynomial, . For example, a polynomial may be fit to the available data for death rate by age group. Then the model equations becomefor and , with boundary conditions and . The age a is considered to be observable. By using the I equation to solve for S and substituting into the remaining equations, we obtain a set of input-output equations that are differential polynomials. We find that the corresponding monic polynomial contains the monomial terms for all . Thus, all coefficient parameters in are identifiable and can be uniquely recovered. We also find that are identifiable from the model equations and is identifiable from the boundary equations.

PDE SEI model with arbitrary death rate function

The final and most general option we consider is to allow to be an arbitrary function of a, i.e. the model equations becomefor and . In this case, we use the I equation to solve for S and then use the E equation to solve for . We obtain the following expression for in terms of the outputs and parameters:By substituting (9) into the remaining equation for S and performing the identifiability analysis, we again find that are identifiable. We note that the second derivatives of and appear in the input-output equations, and hence must be at least twice differentiable. We then find that the right hand side of (9) consists entirely of observable quantities and identifiable parameter combinations; hence is identifiable. Note that while we use the term “identifiable” here, what we actually mean is that can be determined from observable quantities and identifiable parameter combinations. Further, we obtain from the boundary equations that is identifiable.

Identifiability results for all cases of the ODE and age-structured PDE model systems are summarized in Table 1. We note that in the cases where the same death rate is used for all compartments, the death rate function is always identifiable. In cases where the compartments have different death rates, the rates for the susceptible and infected compartments are identifiable. These results hold under the assumption that a fraction of the total prevalence and a fraction of all exposures can be observed at all times and ages. It has been shown previously that other output functions may lead to non-identifiability of death rates even in simple population models (Perasso and Razafison 2016).

Table 1

Summary of identifiable parameters and combinations for the ODE and age-structured PDE systems where the output functions are and

Case	Identifiable parameters and combinations
Density-dependent ODE (Sect. 4.1)	μ_S,μ_I,βk_I,βcϵ,δk_Ik_E,δ(ϵ-1)ϵ,μ_E+δ
Density-dependent PDE with constant parameters (Sect. 4.2)	μ_S,μ_I,βk_I,βcϵ,δk_Ik_E,δ(ϵ-1)ϵ,μ_E+δ
Frequency-dependent ODE (Sect. 5.1)	β,δ,ϵ,μ_S,μ_E,μ_I,ck_I,k_Ek_I
Frequency-dependent PDE with constant parameters (Sect. 5.2)	β,δ,ϵ,μ_S,μ_E,μ_I,ck_I,k_Ek_I
Frequency-dependent ODE with immigration (Sect. 5.3)	β,δ,ϵ,μ_S,μ_E,μ_I,(c+m)k_I,k_Ek_I
Frequency-dependent PDE with immigration (Sect. 5.3)	β,δ,ϵ,μ_S,μ_E,μ_I,ck_I,mk_I,k_Ek_I
Frequency-dependent PDE with piecewise-constant death rates (Sect. 5.4.1)	β,δ,ϵ,ck_I,k_Ek_I∪⋃i=0n-1{μSi,μIi,μEi}
Frequency-dependent PDE with exponential death rate (Sect. 5.4.2)	β,δ,ϵ,κ,μ_0,ck_I,k_Ek_I
Frequency-dependent PDE with polynomial death rate (Sect. 5.4.3)	β,δ,ϵ,ck_I,k_Ek_I∪⋃i=0n{κ_i}
Frequency-dependent PDE with arbitrary death rate function (Sect. 5.4.4)	β,δ,ϵ,ck_I,k_Ek_I,μ(a)

Practical identifiability

The above structural analysis determines identifiable parameters and combinations under the assumption that data are noiseless and are available for arbitrarily many time and age points. In reality, epidemiological data are often noisy and are collected only at discrete time points and for discrete age groups. Thus, it is also necessary to perform practical identifiability analysis to determine which parameters and parameter combinations are still identifiable given more realistic data. The practical identifiability of model parameters given realistic data may differ from their structural identifiability.

As a proof of concept, we explore the practical identifiability of the constant-parameter model from Sect. 5.2 using profile likelihoods. We assume that data are collected monthly for 5-year age groups from age 0 to 100 over a span of 20 years. We synthetically generate data by choosing a baseline parameter set, simulating the model, and adding noise that is normally distributed with a standard deviation of 5%. Using these synthetic data, we perform parameter estimation using the Matlab function fminsearch, and evaluate practical identifiability of the model parameters as detailed in Sect. 2.4.

We compute the Fisher Information Matrix (FIM) using the outputs and at the above-specified ages and time points. We find that for a wide variety of parameter sets within reasonable ranges, the rank of the FIM is eight, implying that there should be eight identifiable parameters and parameter combinations. This is consistent with our structural identifiability analysis performed in Sect. 5.2, which revealed a set of eight identifiable parameters and combinations. From the structural analysis, we see that all parameters except , and are structurally identifiable. However, from profile likelihoods we observe that the practical identifiability of these parameters does vary across the parameter space. Here, we demonstrate two baseline parameter sets that lead to qualitatively different identifiability profiles (Fig. 3). Practical identifiability varies across the parameter space, from clear minimum consistent with the structural identifiability results, to practical unidentifiability as the likelihood surface flattens out around the best-fit value. For both parameter sets, and are not practically identifiable.

Fig. 3

Practical identifiability of the constant-parameter model from Sect. 5.2. Cost functions (black lines) based on likelihoods are shown for two sets of parameter values, (panel A) and (panel B). Practical identifiability varies across the parameter space, from clear minimum consistent with the structural identifiability results, to practical unidentifiability. Note that all model parameters are naturally bounded below by zero, and thus a lower bound on the confidence interval is only meaningful if it is greater than zero. The x axes in the plots above extend to zero for all parameters that do not have such a lower bound. Red dashed lines denote 95% confidence thresholds, and red dots denote local minima (corresponding to maximum likelihoods). See Sect. 2.4 for more details of the analysis

Using the parameter set , we find that certain structurally identifiable parameters are only practically identifiable on one side, that is, we can establish a reasonable lower or upper bound but not both. While the profile likelihoods do have a unique minimum, the profiles are flat and confidence intervals for these parameters are thus wide. In particular, the transmission rate cannot be reasonably bounded from above, and and cannot be distinguished from zero.

Using a slightly different parameter set , in which reporting rates and death among the exposed population are increased, we find that all structurally identifiable parameters are also practically identifiable. While this parameter set is not necessarily biologically realistic (since , for example), this demonstrates that there are regions in the parameter space where these parameters can be well estimated.

Summary and discussion

Creating mathematical and computational models can greatly impact the biological community through predictions that can be used and tested. Parameterization of models is a key factor in connecting them to the biology, and thus identifying parameters accurately is an important step of model development. Parameters are typically estimated from data; however, data may be sparse and noisy, and even perfect data does not guarantee unique estimability of parameters. Understanding if model parameters are identifiable lends credibility to the model system. If parameters are not identifiable, then the applicability of those models becomes more theoretical. While this may add value to understanding, it limits the potential for predictions to be translated in a meaningful way toward impacting the biological system. To this end, we focus on parameter identifiability in both ODE and PDE models. While there is a solid literature on structural identifiability analysis for ODE models, there is not a unified framework for studying identifiability in PDEs. Here, we have presented a methodology for performing structural identifiability analysis on age-structured epidemic models using a differential algebra framework. We have applied this methodology to explore the structural identifiability of a simple age-structured Susceptible-Exposed-Infected (SEI) epidemic model under various assumptions.

We found that when the age-structured model uses constant parameters and does not include immigration, then the identifiability results are identical to those for the corresponding ODE model. In the presence of immigration, however, differences in identifiability arise: the rates for birth and immigration appear only as a sum in the identifiable combinations of the ODE model, while they appear separately in those for the age-structured PDE model formulation. Further, we found that age-dependent death parameters in the PDE model are able to be uniquely recovered in at least the following three cases: piecewise-constant death rates, a single exponential death rate for all compartments, and a single polynomial death rate for all compartments.

These findings imply that identifiability results for corresponding ODE and PDE models are likely equivalent when parameters are constant, but that additional information can be gained from PDE models when parameters vary with the non-temporal variable(s). Rigorous derivation of necessary and sufficient conditions for different models to have identical identifiability results remains to be done, however based on our results herein, we conjecture that identifiability for the ODE and PDE forms of an age-structured model should be equivalent when there are not age-dependent parameters present in the model. Assuming data is available for all time and age points, we also expect age-dependent parameters to be identifiable in the PDE model if the corresponding constant parameters are identifiable in the ODE model. Additionally, future work to investigate how discretization of the PDE (e.g. converting age dynamics into linear ‘chains’ of compartments for age groups (Hurtado and Kirosingh 2019; Hurtado and Richards 2020)) impacts identifiability of the model would be a useful direction given the common use of compartmental representations of aging in biological models. Since the PDE system can be viewed as the limit of age-structured compartmental ODE models as the number of age compartments approaches infinity, we expect the identifiability results for the PDE system to also be a limiting case of those for compartmental ODE systems.

More generally, this work underscores some of the directions that will need to be explored to extend the differential algebra approach to the PDE case. The framework presented here could be extended to other PDE models, such as reaction-diffusion and advection-diffusion equations, in a variety of application areas. In particular, the assumption used here that age-interactions are only local (i.e. individuals interact with only people their own age) may be reasonable for reaction-diffusion models. But for more general models, this method will need to be extended to handle integrals across some variables—either to account for measurement processes that aren’t age-structured or spatially resolved, or to account for non-local interactions. However, since this methodology relies on the ability to express the input-output equations as differential polynomials, there is a limit to the amount of complexity that it can handle. Extension of concepts such as Gröbner bases and characteristic sets (which enable automated generation of the input-output equations for ODEs) to the PDE case would be an important step in this direction. However, highly nonlinear models and/or output functions may require the use of other tools such as numerical methods. One such method is the FIM, which we have demonstrated in the context of an age-structured model with preferential mixing informed by real data. Additional tools for examining identifiability of ODE models, e.g. based on Lie derivatives and observability methods (Villaverde et al. 2018; Yates et al. 2009; Raue et al. 2014) may also be able to be translated to the PDE case. In applying this methodology to other classes of PDE models, care will need to be taken regarding required existence and uniqueness properties. Rigorous derivation of the necessary and sufficient conditions for this methodology to be applied remains a direction for future research. Similarly, expanding the types of PDE models considered will also require further consideration of how boundary conditions (particularly age or spatially varying boundary conditions) may affect model identifiability—either in the case where the boundary conditions are known and may improve model identifiability or in the case where they involve unknown parameters to be estimated. Finally, the approaches demonstrated here using PDEs may also open up the possibility of examining identifiability of stochastic systems that can be characterized using PDEs (e.g. via Kolmogorov equations). We have previously begun to explore this possibility for a class of continuous time Markov models used in cancer modeling (in which case the Kolmogorov equations could be written as ODEs) (Brouwer et al. 2017).

A significant limitation of structural identifiability analysis is that it assumes the availability of perfect data. In practice, data may contain noise and bias, and may be sampled only at discrete time points. Thus, practical identifiability analysis must also be performed to gain a complete understanding of the parameter identifiability of a model given the available data. As demonstrated in our examples, practical identifiability of a model may vary not only with the quality of the data but also with the corresponding values of model parameters. Practical identifiability thus needs to be handled on a case-by-case basis for models and their corresponding parameter sets. Taken together, practical and structural identifiability provide rigor to the modeling process, adding to the applicability of model predictions.

A further limitation of this methodology for structural identifiability analysis is that it relies on the ability to eliminate system variables to obtain input-output equations. The required algebra to do so may be intractable for sufficiently complex models and/or outputs, and may necessitate extending concepts such as characteristic sets and differential Gröbner bases (Boulier 2006) to the PDE identifiability setting. Alternatively one might need to rely on numerical approaches to identifiability. For example, an age-structured model with non-local mixing would require integral terms in the state equations, making this approach substantially more difficult. Similarly, other concepts from the ODE case, such as the idea of “persistently exciting” inputs (Ljung and Glad 1994), will need to be translated to the PDE case as well. Structural identifiability analysis of such problems is a direction for future work.