Bayesian Multi-level Mixed-effects Model for Influenza Dynamics

Abstract

Influenza A viruses (IAV) are the only influenza viruses known to cause flu pandemics. Understanding the evolution of different sub-types of IAV on their natural hosts is important for preventing and controlling the virus. We propose a mechanism-based Bayesian multi-level mixed-effects model for characterising influenza viral dynamics, described by a set of ordinary differential equations (ODE). Both strain-specific and subject-specific random effects are included for the ODE parameters. Our models can characterise the common features in the population while taking into account the variations among individuals. The random effects selection is conducted at strain level through re-parameterising the covariance parameters of the corresponding random effect distribution. Our method does not need to solve ODE directly. We demonstrate that the posterior computation can proceed via a simple and efficient Markov chain Monte Carlo algorithm. The methods are illustrated using simulated data and a real data from a study relating virus load estimates from influenza infections in ducks.

1 INTRODUCTION

Influenza A viruses (IAV) can spread beyond species barriers and adapt rapidly to new hosts and environmental conditions and thus represent a major threat to human and veterinary health. IAV are divided into many sub-types and a pandemic can occur when a new and very different virus emerges that both infects people and has the ability to spread efficiently among people. Wild aquatic birds, particularly certain wild ducks, are the natural hosts for most IAV. Understanding the evolution of different sub-types of IAV in their avian hosts is an important prerequisite to better surveillance and preparedness for pandemic threats.

Our study is motivated by the virus load data from influenza infections in ducks for seven different IAV strains (Figure 1). The seven virus isolates were obtained from an ongoing long-term surveillance of IAV circulation in wild birds in Minnesota, USA (Keeler et al., 2013; Wilcox et al., 2011). Five viruses were isolated from dabbling ducks. Two additional viruses were isolated from surface lake water. For six of the influenza strains, measurements for five ducks were collected, for one strain only (H6N2) data for four animals is available. There are total 34 ducks. For each duck the virus loads at different time points after inoculation are measured. Virus load was quantified based on the Matrix gene copy number in cloacal and oropharyngeal samples estimated by real-time polymerase chain reaction. Further details related to the experimental design and methodology are available from the previously published studies (Lebarbenchon et al., 2011, 2012; Handel et al., 2014).

The genetic diversity of IAV naturally arose from genetic re-assortments and the seven virus strains displayed in Figure 1 were identified on the basis of full-genome sequencing. We aim to investigate the effects of IAV genotype on the within host transmission dynamics in wild duck populations. The mechanism-based ordinary differential equations (ODE) models have been widely used to describe influenza viral dynamics due to their transparent structures. Parameters in ODE often have biologically interpretations and can be used to investigate the effects of possible interventions, such as treatment strategies. When the dynamical process are measured for multiple subjects as shown in Figure 1, the hierarchical mixed-effects ODE models appear to be reasonable for modelling the infection dynamics at population level. The mixed-effects models include random-effects to account for within-subjects correlation and allow for the viral dynamic processes to share certain similar patterns between subjects while still having distinct individual characteristics. For the study described in Figure 1, the random effects can come from both strain and individual levels.

Estimating mixed-effects ODE models is an important but challenging problem because most of the ODEs do not have closed form solutions due to their nonlinear features. Several methods and software have been developed to estimate mixed-effects ODE models by solving ODE numerically. Both the commercial software package, MONMEM (Beal & Sheiner, 1980) and the free popular software, MONOLIX (Lixoft, 2012) use EM-type algorithm such as stochastic approximation EM algorithm (Kuhn & Lavielle, 2005) to estimate non-linear mixed effect models. Huang et al. (2006) and Huang and Wu (2006) use a Bayesian MCMC approach to estimate mixed-effects for a HIV infection ODE model.

The problem with all these methods is that they need to numerically solve the ODE equation for every updating of the parameters and thus are very intensive in computation. For estimating ODE parameters from a single subject, the two-step-based approaches have been proposed in the literature to reduce this computational complexity. Examples include Campbell and Steele (2012), Dondelinger et al. (2013) and Macdonald et al. (2016) among many others. For multiple subjects, based on the profiling approach proposed in Ramsay et al. (2007), Wang et al. (2014) developed a semi-parametric method to estimate mixed-effects ODE models. They use spline basis functions to approximate the dynamic process and estimate the spline coefficients, the ODE random effects and fixed effects in three nested levels of optimisation. As mentioned in Wu et al. (2012), this spline smoothing-based parameter cascade approaches may not be flexible enough to deal with the complicated local features of the data. In addition, they need more efficient optimisation techniques and complicated iterative computation algorithms to obtain an estimator.

Huang et al. (2020) propose a Bayesian model that directly combine the smoothness, system observations and ODE together. In this model, the ODE constraints are replaced with a probability expression and combined with non-parametric regression into a joint likelihood function. The benefit of this Bayesian approach to parameter estimation in ODEs is that, for some special ODE formulations, the closed form conditional distribution for ODE variables can be obtained which substantially facilitates the convergence of MCMC process. Moreover, this approach can be applied to situations where only partial components of the state variables are observed. In this study, we incorporate the hierarchical mixed-effect model into this framework to investigate the variability among different levels. We combine all subjects together by including random effects to ODE parameters to account for the variations among different strains and individuals.

An important practical problem in multi-level mixed-effects model is to select the random effect components and decide which coefficients vary at different levels with non-zero variance. Assume there are $d$ potential ODE parameters which can include random effects, one can use standard model selection measures, such as Akaike's information criterion (AIC) or Bayesian information criteria (BIC), to compare all $2d$ models for random effect variable selection. However, when $2d$ is large, it is very computational intensive because estimating parameters in each ODE model is more complicated than dealing with regular statistical models. Miao et al. (2009) use traditional subset selection methods AIC and BIC to select an ODE model. Zhang et al. (2015) proposed a new ODE model selection method which employs statistical estimation of the full model followed by a combination of a least squares approximation and the adaptive lasso. On Bayesian model selection, a great amount of work exists for fixed effects; however, very little exists on selection of random effects. The main difficult is that it is not clear what penalty for model complexity is appropriate for comparing models with different numbers of random effects.

Chen and Dunson (2003) proposed an approach to parameterise the mixed model so that functions of the covariance parameters of the random effect distribution are incorporated as regression coefficients on standard normal latent variables. They allow random effects to effectively drop out of the model by choosing mixture priors with point mass at zero for the random effects variances. We incorporate this approach into our multilevel Bayesian ODE framework and perform selection on the random effects component in order to identify ODE parameters that vary across different virus sub-types.

Most closely related to the current paper are results by Liu et al. (2019) that proposed a Bayesian mixed-effects ODE model which also estimated the ODE solution with a linear combination of basis functions. However, the major difference is that they approximated the ODE solution by minimising a penalty term, while we replace the ODE constraint with a probability expression and implement a fully Bayesian estimation. Moreover, their method was only evaluated by a simple one variable ODE example without considering random effects variable selection.

The rest of the paper is organised as follows. Section 2 describes ODE models for population-level influenza dynamics. Section 3 outlines the Bayesian algorithm for ODE parameter estimation. Section 4 presents results from a simulation study. Section 5 illustrates the methodology via an application to the viral load data for influenza infections in ducks illustrated in Figure 1, and Section 6 summarises the results.

2 ODE MODEL FOR INFLUENZA DYNAMICS AT POPULATION LEVEL

2.1 ODE model for within-host influenza dynamics

The infection dynamics in avian hosts has been recognised and a mechanistic target-cell limitation model can be used to describe virus kinetics. The model tracks uninfected cells $U$⁠, infected cells $I$⁠, and free infectious virus $V$⁠. They are all functions of the time $t$⁠. Assume that cells become infected at rate $a$⁠, infected cells produce virus at rate $p$ and die at rate $b$⁠, and free virus is cleared at rate $c$⁠. The dynamics of the system can be described by the following differential equations

In this model, the state variables are $U(t)$⁠, $I(t)$ and $V(t)$⁠, and the unknown parameters are $a$⁠, $b$⁠, $p$ and $c$⁠. These types of models have been used in several analyses of influenza A virus within-host infection dynamics (see e.g. Nowak & May, 2000; Smith & Perelson, 2011 for reviews). This simple model can describe most data for influenza virus infections rather well, see for example, a study in Baccam et al. (2006). After an initial rise in virus load, uninfected target cells become depleted, leading to a subsequent virus decline and resolution of the infection. This so-called target cell limited model is basically equivalent to a simple epidemic model, which produces a single infectious disease outbreak in a susceptible population. This is also a typical inverse problem for ODE and our goal is to estimate the parameters using the measurements of the state variables. The observed state variables are not exactly the ODE solution because they are measured with noise.

In the experiment of influenza infections in ducks, only the virus load was observed among the three components involved in the dynamics. It is necessary to study parameter identifiability of ODE models with only $V$ measured before developing or applying statistical methods to estimate model parameters; otherwise, including unidentifiable unknown parameters in a model is likely to result in biassed estimates and thus misleading conclusions. Starting with Model (1) and recalling that only $V$ is experimentally measured, we can eliminate $U$ and $I$ from the model to obtain

which is independent of $p$⁠. Thus, we can let $p=1$ and the unknown parameters for ODE model only include $a$⁠, $b$ and $c$⁠.

2.2 Multi-level mixed-effects ODE model for influenza infections

Model (1) describes the infection dynamics at individual level. In order to combine all subjects together and determine how infection kinetics differ among strains and subjects, we consider mixed-effects model which includes both fixed effect components and random effect components. Motivated by the examples of influenza infections in ducks in Figure 1 which show strong variability of infections among both individuals and strains, we consider a two level random effect model in which the first level comes from the strain and the second level comes from the individual.

Denote the number of strains by $S$ and the number of the subjects in the sth strain by $ns$⁠. The total number of subjects is $n=∑s=1Sns$⁠. Define a vector of ODE parameters $θ=(a,b,c)$⁠. Let $θsj=(asj$⁠, $bsj$⁠, $csj)$ represent the parameters for the jth subject in the sth strain. For the jth subject in the sth strain, $msj$ virus load observations $ysj1,…,ysjmsj$ are obtained according to independent additive noise model

where the noise is Gaussian $ϵsji∼N(0,σsj2)$ and $Vθsj(ti)$ is the solution of virus load to ODE equation (1) at time $ti$ based on coefficients $θsj$⁠. In multi-level mixed-effect model, we include random effects at both strain and individual levels, and have the following formulation

Here $μ$ represents the fixed effect of the dynamical coefficients, $αs$ represents the random deviation for the sth strain, and $βsj$ represents the random deviation for the jth subject in the sth strain. $∑S$ and $∑I$ are $d$ dimensional covariance matrices, where $d=3$⁠. In hierarchical mixed-effect ODE model (3), it is of primary interest to estimate the fixed effects $μ$ and random effects $∑I$ and $∑S$⁠.

3 PARAMETER ESTIMATION

We estimate the multilevel mixed-effects ODE model by incorporating the hierarchical random effects into the Bayesian ODE framework developed in Huang et al. (2020). Section 3.1 gives a brief description of this Bayesian framework. Then in Section 3.2 we build a multi-level ODE model by incorporating (3) into this framework. A hierarchical Bayesian method for selecting the random effect components in ODE coefficients is also developed in this section. Section 3.3 summarises the full conditional posterior distributions for all the variables.

3.1 Bayesian ODE parameter estimation approach at individual level

In a general continuous time dynamical system, the evolution of $K$ states $x(t)≡{x1(t)$⁠, $…,xK(t)}$ is described by a set of $K$ ODEs

where $θ∈ℝd$ is the parameter vector of ODE and $f(·)={f1(·),…,fK(·)}$ is a vector of known appropriately smoothing functions. Without loss of generality, assume that the first $K0$ states are measured at time points $t1,…,tm$⁠. A common situation in practice is that only part of the system is measured, that is, $K0<K$⁠. Model (1) can be fit into this general framework with $d=3$⁠, $K=3$⁠, $K0=1$⁠, $x1=V$⁠, $x2=U$⁠, and $x3=I$⁠.

In Huang et al. (2020), $xk$ and its derivative were expressed in terms of a basis function expansion

where $ϕ(t)≡(ϕ1(t),…,ϕq(t))T$ is a vector of basis functions such as polynomial spline bases or B-spline bases. Denote $q×K$ matrix $C≡[c1,…,cK]$ and the observation matrix $Y≡[y(t1),…,y(tm)]$⁠, where $y(t)≡(y1(t),…,yK0(t))T$⁠. The likelihood function of the model is

where $σ2≡{σ12,…,σK02}$⁠. Here we assume that the noise variability $σk2$ does not depend on the observed state variable and is a constant across different time points. This homoscedastic assumption has been widely used in the literature of ODE parameter estimation. The consistency of the two-stage least square estimator has been established in Liang and Wu (2008) in cases where the true variability in the noise is not necessarily constant. Moreover, for ODE problem (1), if the virus load at a time point is below the threshold, it cannot be precisely measured and usually is assigned a tiny number close to 0. It will cause numerical problem in the analysis if we assume that the variability of noise is proportional to the observed virus load. Therefore, the homoscedasticity assumption we made can also facilitate numerical stability in the analysis. For each $σk2$⁠, we employ an inverse Gamma prior

for some shape and rate hyper-parameters $(ασ,βσ)$⁠.

The exact ODE constraint (4) was replaced with the following probability expression which models a normal distribution for the discrepancy between $x˙(t)$ and $f(x(t),θ,t)$⁠.

where $δ2≡{δ12,…,δK2}$⁠. The advantage of the probability (8) for ODE is that in situations where $f(x(t),θ,t)$ only depends on $xk$ up to the first order of power, the full conditional posterior distribution for $ck$ becomes multivariate normal. For each $δk2$⁠, we employ an inverse Gamma prior

for some shape and rate hyper-parameters $(αδ,βδ)$⁠.

3.2 Random effects selection in hierarchical mixed-effect ODE model

Combining all subjects and strains together, the observation model (6) and ODE requirement (8) lead to the following joint distribution of the whole observations over state parameters $Csj$⁠, ODE parameters $θsj$⁠, observations $Ysj$⁠, and the remaining parameters

where $P(θ)$ is determined by (3), and $P(σ2)$ and $P(δ2)$ are properly chosen priors for parameters $σ2$ and $δ2$ which are given in (7) and (9), respectively. For those components of $C$ which are not observed, they are not included in $POBS(Y|C,σ2)$ but can still be estimated because they appear in $PODE(C|θ,δ2)$⁠.

To select strain-level random effects in (10), we follow Chen and Dunson (2003) and use a modified Cholesky decomposition of $∑S$⁠: $∑S=ΛΓΓTΛ$⁠. $Λ$ is a positive diagonal matrix with diagonal elements $λ=(λ1,…,λd)$ proportional to the random effects SDs, so that setting $λl=0$ is equivalent to dropping the lth random effect from the model. $Γ$ is a lower triangular matrix with diagonal elements equal to 1 and off-diagonal elements describing the random effects correlations. In the case of independent random effects, $Γ$ is simply the identity matrix $Id$ and the diagonal elements of $Λ$ equal to the random effects SDs. Applying the covariance decomposition to (3) we have

which leads to

where

and $P(μ)$⁠, $P(Λ,Γ)$⁠, and $P(∑I)$ are priors for $μ$⁠, $Λ$⁠, $Γ$ and $∑I$⁠, respectively.

Motivated by practical considerations, we want to choose priors that facilitate posterior computation. For this reason, prior distributions that are conditionally conjugate are desirable. For $μ$ and $∑I$⁠, we choose the following conjugate priors

where $ν0>0$⁠, $η0$ is a $d$-dimensional vector, $Λ0$⁠, $Ω0$ are $d×d$ positive definite matrices. Here $η0,Λ0,Ω0,ν0$ are hyper-parameters and are known.

In choosing priors for $Λ$ and $Γ$⁠, we assume

where $γ={γ1,…,γd(d−1)/2}$ are elements of $Γ$⁠. Here, the indicator function $I(·)$ sets zero for the elements of $γ$ that correspond to zero $λ$⁠. Further specify

where $ZI-N+(pl0,ml0,sl02)$ denotes the density of a zero inflated half normal distribution consisting of a point mass at zero (with probability $pl0$⁠) and a $N(ml0,sl02)$ density truncated below by zero.

3.3 The full conditional posterior distributions

From (10) and (12), the full conditional distributions for all parameters can be obtained as

Using priors (13), the full conditional distributions of $bs$ and $μ$ are normal, $bs∼N(μ‾s,∑‾s)$⁠, $μ∼N(μ‾μ,∑‾μ)$⁠, where

The full conditional distribution of $∑I−1$ is Wisharts, $∑I−1∼Wishart(ΩI−1,n+ν0)$⁠, where

The full conditional distributions of $σsj2$ and $γsj2$ are inverse Gammas. The full conditional distribution of $γ$ given $Λ$ and $bs$ has a multivariate normal distribution. The full conditional distribution of $λl$ is zero inflated half normal. The full conditional distributions of $X$ and $θ$ depend on the explicit form of function $f(·)$⁠. For situations where $f(x(t),θ,t)$ only depends on each argument up to the first order of power, we can derive a closed form full conditional posterior distribution for $θ$ and each component of $X$⁠. For model (1), $f(·)$ is a linear function of both $θ$ and all the components of $X$⁠, the full conditional distributions for $θsj$ and the components of $Csj$ are all multivariate normal. Thus we can use closed form Gibbs sampling scheme in the MCMC process.

4 SIMULATION STUDY

In this section, we analyse data from simulated studies to illustrate our approach. The design of simulations is similar to the motivation experiment of influenza infections in ducks described in Section 1. In particular, we consider six strains and generate data for 36 subjects, six for each strain. For each subject, we assume that measurements of viral load are taken at every 0.5 time unit on the interval $[1,14]$⁠, that is, we have 27 virus load measurements for each subject. We let the population-level dynamic parameters $μ=(9,2/3,2/3)$⁠. The individual dynamic parameters $θsj$ for the jth subject in the sth strain are generated by (11) with the individual level covariance matrix $∑I=diag(0.9,0.1,0.1)$ and the strain-level covariance matrix $∑s$ being the product of $Λ$ and $Γ$⁠, where $Λ$ takes

and $Γ$ takes

Therefore, we have six different settings for $∑s$⁠. The random variation at strain level occurs on parameter $a$ for $Λ1$⁠, on parameters $b$ and $c$ for $Λ2$⁠, and on all parameters $a$⁠, $b$ and $c$ for $Λ3$⁠. The parameters $a$⁠, $b$ and $c$ are generated independently for $Γ1$⁠, but for $Γ2$⁠, the correlation is considered.

We use the initial condition $U0=11$⁠, $I0=0$⁠, and $V0=0.1$⁠. The trajectories of $U(t)$⁠, $I(t)$⁠, and $V(t)$ based on $μ=(9,2/3,2/3)$ and this initial condition is shown in Figure 2. Then, based on generated true parameters $θsj$⁠, the observations $ysji$ are generated by perturbing the solution of the differential equation (1) with a within-subject measurement error, that is, $ysji=Vsj(ti)+ϵsji$⁠, where $Vsj(ti)$ is the numerical solution of viral load to the differential equation (1) for the jth subject in the sth strain at time $ti$⁠. It is assumed that the within-subject measurement error $ϵsji$ is normally distributed with $N(0,0.12)$⁠. For each setting, we replicate the simulation 100 times. Thus our conclusion is based on the summary over analyses of 100 simulated data sets.

We applied the Gibbs sampling algorithm described in Section 3.3 to all the simulated data and test the performance by comparing the estimated dynamic parameters with the true dynamic parameters. We choose hyper-parameters $η0=0$⁠, $Λ0=Ω0=105I3$⁠, $ασ=αδ=βσ=βδ=ν0=10−5$ to reflect the non-informative priors for $μ$⁠, $∑I$⁠, $σ2$ and $δ2$⁠. The hyper-parameters for $λl$ are chosen to be $ml0=0$ and $sl02=100$⁠. The hyper-parameters for the elements of $Γ$ are chosen to be $γ0=0$ and $R0=I3$⁠. To study the effect of the prior probability that $λl=0$⁠, we use $pl0=0.2$⁠, 0.5, and 0.8, respectively.

To fit the non-parametric model, Cubic B-splines with 25 equally spaced knots in $[1,14]$ are used to approximate the ODE solution. To choose appropriate initial values for the elements of $θsj$ and $Csj$⁠, we solve the least square problem for RSS$sj$ = $∑i=1msj(ysji−Vsj(ti))2$ which are function of the ODE parameters $θsj$ and initial condition $U0,I0,V0$⁠. To ensure the identifiability of this optimisation problem, we fix $I0=0$⁠, $V0=ysj1$ and obtain $U0,θsj$ by minimising RSS$sj$⁠. Then the estimated values for $θsj$ are taken to be the initial values for $θsj$ and the estimated values for $U(ti),I(ti),I(ti)$ are used to derive the initial values for $Csj$ in our MCMC procedure.

For the analysis of each data set, the MCMC sampler was run for a total of 25,000 cycles. The first 5000 cycles were discarded as burn-in, and the remainder of the chain was thinned by keeping one out of every 10 samples, resulting in a total of 2000 samples for post-MCMC analysis. To compute the posterior probabilities of each of the $23$ models for strain level random effects, we simply add up the number of occurrences of each model and divide by the total number of iterations.

Table 1 shows the summary of posterior probabilities over 100 simulated data for eight different strain-level random effect models. There is overwhelming evidence that the true model is selected with the highest posterior probabilities. The only exception is for setting $∑s=Λ2Γ2$⁠, where the model $(a,b,c)$ has the largest posterior probability while the true specification, that is, the ODE mixed model with correlated strain-level random effects on $b$ and $c$⁠, has the second largest posterior probability. This conclusion holds regardless of the choice for the prior probability of $λl=0$⁠. One possible explanation is that the correlation between $a$ and $(b,c)$ causes the inclusion of parameter $a$ in the model and thus the higher probability for model $(a,b,c)$⁠. We further studied the marginal probabilities from Table 1 for the true parameters to be included in the model and obtained the consistent results. For example, under $pl0=0.5$⁠, the probabilities for including random effects of $(a)$⁠, $(b,c)$⁠, and $(a,b,c)$ in the model are 95.7%, 57.8% and 58.4%, respectively for settings $Λ1Γ1$⁠, $Λ2Γ1$ and $Λ3Γ1$⁠. This shows that the proposed method obtains good estimates for strain-level random effects.

To evaluate the performance of our method in different scenarios, we define the average relative estimation error (ARE) of a parameter $θ$ as

where $θ^i$ is the estimate of $θ$ in the ith simulation and $N$ is the number of simulation runs (here $N=100$⁠). Table 2 presents the posterior summaries for the fixed effect ODE parameters over 100 simulated samples including mean, SE, and ARE, which are properties associated with the population level. Here we choose non-informative prior $pl0=0.5$⁠. It can be seen that the estimated values are quite close to the true values, suggesting good performance of the algorithm. Parameter $a$ is slightly underestimated. The reason is that $a$ is mainly determined by the decay rate of uninfected cell $U$⁠. As shown in the upper panel of Figure 2, $U$ decays very fast and quickly becomes flat. Therefore, only the measurements at very beginning time are useful in the estimation of $a$⁠. We can increase the accuracy of parameter estimation for $a$ by including more measurements and choosing more knots at the beginning. But in practice, only virus load $V$ is observed, we can only choose knots according to the patterns of virus load changes as described in Section 5.

We have also conducted simulations to study how the performance of our algorithm is influenced by the number of subjects per strain, that is, $ns$⁠. We found that increasing $ns$ can reduce the bias in the estimation of population level ODE parameters. For example, under setting $Λ1Γ1$ with $pl0=0.5$⁠, if the number of subjects per strain $ns=20$⁠, we obtain that the AREs are 0.13, 0.08 and 0.05 for the fixed effect coefficients a, b and c, respectively, which are smaller than the corresponding values of 0.18, 0.09 and 0.07 obtained based on $ns=6$⁠. We also studied the effect of the number of virus load measurements $m$ for each subject and found that the conclusion is pretty similar, that is, the bias for the estimation of the fixed effect coefficients decreases with $m$⁠.

5 APPLICATION TO THE EXPERIMENT OF INFLUENZA INFECTIONS IN DUCKS

In this section, we apply the proposed methodology to the data set of viral load dynamics following influenza infections in ducks as shown in Figure 1. Cubic B-splines with 10 equally spaced knots between day 1 and day 4 and 10 equally spaced knots between day 4 and day 14 are used to approximate the ODE solution. The knots are chosen in this way in order to capture the patterns of virus load changes which increase sharply from day 1 to day 3 followed by a slow and continuous decrease until day 14. We have tried different number of knots in the analysis and obtained very similar results. For example, if we choose seven knots between day 1 and day 4 and five knots between day 4 and day 14, the conclusion is pretty similar for both visualisation and parameter estimation. Since no previous information is available about the strain level variation in duck flu, we choose non-informative prior $pl0=0.5$ on the set $λl=0$ for $l=1,2,3$ to reflect a balance among outcomes. For all the other hyper-parameters, we use the same values as in the simulation study to induce relatively diffuse priors on those parameters. We also follow the similar way as in the simulation study for choosing the initial values of $θsj$ and $Csj$⁠. We discard the first 5000 iterations as burn-in. Posterior probabilities are calculated based on 20,000 iterations collected thereafter. To compute posterior summaries of the parameters, we thin the chain by a factor of 10.

Figure 3 displays the individual fitted curves based on the solutions of the viral load $V$ to the ODE equations. The virus load of the jth subject in the sth strain is estimated from ODE using parameters $θ^sj$ which is taken as the average of the 2000 posterior samples for $θsj$⁠. We found that the model provided a good fit to the observed data for all of the subjects.

To assess sensitivity to prior model probability, the analysis was repeated for $pl0=0.2$ and $pl0=0.8$⁠, respectively. Table 3 lists the posterior probabilities of the eight sub-models. There is evidence that the ODE mixed model which includes between-strain random variations for the parameter $a$ alone is optimal. This conclusion holds for three choices of the prior probability for $λl=0$⁠. The marginal probabilities of including the strain level random effects for the parameter $a$ are 83%, 62% and 82%. In contrast, the corresponding marginal probabilities are 27.9%, 9.4% 29.85% for including the parameter $b$ and 24.75%, 8.6% 26.6% for including the parameter $c$⁠. This indicates that the cells are infected at different rates for different strains, but both the death rates of the infected cell and the clearance rates of the free virus show strong homogeneity across different strains.

In addition to investigating differences among strains, investigators are interested in assessing the overall population-level effect of the dynamic parameters. The population posterior means and the corresponding SE as well as the 95% highest density interval (HDI) for the three parameters are summarised in Table 4. This result is also robust to the various prior choices. Table 5 displays the posterior estimation of the dynamic parameters for each strain, that is, $μ+αs$ in (3). It can be seen that, for parameters $b$ and $c$⁠, there are no distinct between-strain differences. But for parameter $a$⁠, the between-strain differences are quite clear especially for sub-types H4N8 and H4N6* whose values are much larger and smaller, respectively, than the values of the remaining sub-types. This result is consistent with the posterior probability pattern shown in Table 3.

Since the seven sub-types are identified according to the genomic characterisation, our random effect model selection analysis shows that the genetic differences do have some effects on the transmission dynamics of IAV through changing the infection rate parameter $a$⁠. As our analysis was based on a limited number of tested viruses, such a finding should be reinforced with additional studies focusing on other virus genotypes, at other locations and in different avian hosts. This result nevertheless provides some insight to help biologist for further investigation.

6 DISCUSSION

In this paper, we propose a Bayesian multi-level mix-effects ODE model to describe the dynamics of influenza virus infection. Both the population and individual dynamic parameters can be estimated. In comparison with fixed-effect modelling, mixed-effects modelling has the advantages of obtaining more precise parameter estimation by pooling all data together. By relaxing the ODE constraint with a probability expression, our method does not need to solve ODE directly. The variations at different levels are simultaneously considered by incorporating both strain-specific and subject-specific random effects for each ODE coefficient. One advantage of the proposed method is that the closed form conditional posterior distributions for all variables and parameters can be obtained which substantially facilitate the convergence process. To identify the variation at strain level, we perform random-effects model selection based on a decomposition of the covariance matrix of the strain-specific random-effects distribution. The decomposition enables us to specify a prior distribution which can exclude one or more random effects by setting their variances to zero. This re-parameterisation of the covariance matrix results in straightforward and efficient posterior computation via a Gibbs sampler.

We have presented large-scale simulation examples and an actual long-term surveillance study of IAV circulation in wild birds to illustrate how the Bayesian procedures can be applied to influenza virus dynamics investigation. For the simulation studies, it was seen that the models provided fairly good fits to both ODE random effects and ODE fixed effects. The estimates for both population and individual dynamic parameters are also very close to the true values. For the real data, the proposed model fitted the observed viral load reasonably well for all subjects in our study. Using random effects variable selection, we identify that the cell infection rate $a$ varies among different strains while the other two ODE parameters are quite homogeneous across the strains.

We use a simplified model to explain the observed patterns and characterise the biological mechanisms of IAV infection with the main goals of retaining crucial features of influenza dynamics. We plan to generalise it to more complicated models with consideration of more infected cell and virus compartments in order to provide more accurate descriptions for long-term IAV dynamics. One challenge is the identifiability problem of model parameters due to the complexity of the models and the factor that only viral loads are observed among ODE variables. We need to consider the trade-off between the complexity and applicability of influenza dynamic models. It has been shown in Lebarbenchon et al. (2011) that the environmental factors, particularly temperature and humidity, appear to be important determinant of IAV persistence and therefore within-host fitness. We will consider modelling ODE parameters as function of environmental factors in order to integrate data and models at both the infection dynamics and the persistence of the virus in the environment within a combined ecological and within-host framework.

DATA AVAILABILITY STATEMENT

Data and code can be obtained from the github repository https://github.com/hwhuanghep/random_effect_ODE.

ACKNOWLEDGEMENTS

The author thanks the editor, associate editor and two referees for many helpful comments and suggestions which led to a much improved presentation. The author also thanks Andreas Handel for sharing the data set and Xiao Song for helpful discussion. This research is supported by Division of Mathematical Sciences (National Science Foundation) grant DMS-1916411.

REFERENCES

Author notes