Fitting dynamic models with forcing functions: Application to continuous glucose monitoring in insulin therapy

The artificial pancreas is an emerging technology to treat type 1 diabetes (T1D). It has the potential to revolutionize diabetes care and improve quality of life. The system requires extensive testing, however, to ensure that it is both effective and safe. Clinical studies are resource demanding and so a principle aim is to develop an in silico population of subjects with T1D on which to conduct pre-clinical testing. This paper aims to reliably characterize the relationship between blood glucose and glucose measured by subcutaneous sensor as a major step towards this goal. Blood-and sensor-glucose are related through a dynamic model, specified in terms of differential equations. Such models can present special challenges for statistical inference, however. In this paper we make use of the BUGS software, which can accommodate a limited class of dynamic models, and it is in this context that we discuss such challenges. For example, we show how dynamic models involving forcing functions can be accommodated. To account for fluctuations away from the dynamic model that are apparent in the observed data, we assume an autoregressive structure for the residual error model. This leads to some identifiability issues but gives very good predictions of virtual data. Our approach is pragmatic and we propose a method to mitigate the consequences of such identifiability issues. Copyright © 2011 John Wiley & Sons, Ltd.


Introduction
Type 1 diabetes (T1D) is a chronic autoimmune disorder characterized by dysregulated blood-glucose (BG) levels due to an inability of the pancreas to produce insulin, the hormone that promotes uptake of glucose by cells [1]. Persistent exposure to high glucose levels (hyperglycaemia) causes long-term diabetes complications and organ dysfunction [2]. The standard therapy is based on multiple insulin injections, using a combination of short and long acting insulin analogues, informed by frequent BG selfmonitoring [3]. Treatment by continuous subcutaneous insulin infusion (CSII [4]) is on the rise and uses a portable electromechanical pump to mimic nondiabetic insulin delivery, infusing insulin at preselected rates-basically a slow basal rate with patient-activated boosts at mealtimes. However, intensive insulin therapy aiming to achieve near-normal glucose control is associated with an increased risk of low BG levels (hypoglycaemia), potentially leading to seizures, unconsciousness, brain damage and even death [5]. Optimization of insulin therapy is confounded by large day-to-day and diurnal variability in insulin requirements influenced by factors such as exercise, stress, and recurrent illness [6--8].
Self-monitoring of BG offers only a snapshot, each time, of the underlying glucose excursion, thus, making for considerable uncertainty in determining the right treatment decision to achieve and Facility, Addenbrooke's Hospital, University of Cambridge, UK [24]. The sample size was not based on any power calculation as the study concerned (APCam01 in [24]) was exploratory. Participants and, as appropriate, their carers gave informed consent/assent. The study was approved by the Cambridge Research Ethics Committee (REC Ref 06/Q0108/350).
A glucose sensor was fitted to each participant at least 24 h prior to the study, and, following a run-in period and calibration as suggested by the manufacturer, the Guardian RT CGM system took SG readings every 5 min. BG was measured every 15 min by collecting samples via a venous cannula. The study ran from 17:00 until 12:00 the following day, giving m = 77 BG and n = 228 SG measurements for each individual. There are a small number of missing SG measurements, which we shall treat as unknown parameters in order to retain a balanced data set. Two self-selected meals were eaten at 18:00 and 08:00 the following morning to maintain a normal carbohydrate intake, each meal containing a mean (SD) of 87 (23) g carbohydrates. Prandial insulin boluses were given with the meals. During the study, the Guardian RT was calibrated, using recent BG measurements, shortly after 17:00 and every 6 h thereafter, thus splitting the study period into five distinct 'calibration periods' for each individual.

Background
The model that we will describe in subsequent sections is essentially an extension of the non-linear regression where y j and x j denote response and independent variables, respectively, denotes a set of regression parameters, and j ∼ N(0, 2 ), say. As a simple example, let us consider a situation in which the regression function g() represents exponential decay and x is elapsed time: Note that (1) is the unique solution to dg dx = − 2 g, g( , and so (1) and (2) are equivalent specifications of the same regression function. Now suppose that we can only express our regression function in terms of differential equations (with the corresponding initial conditions) as we do not know the analytic solution. If we know that a unique solution exists, however, then we know that solution is simply a deterministic function of the inputs, and x, albeit of unknown form. If we can find a way to evaluate the solution, we may thus treat it as we would any other deterministic function. We may then exploit standard graphical modelling theory [29,30] to evaluate the full conditional distributions of any unknown inputs, e.g. 2 . In this paper, we make use of the BUGS software [31,32] with WBDiff interface [33] installed (to allow specification of differential equations). The differential equations, described in the following subsection, are solved numerically by the software using a Runge-Kutta algorithm [34], and Metropolis-Hastings samplers [35,36] are typically used for sampling the unknown inputs. Now suppose that the differential equations depend on some additional quantity, such as the ambient temperature or pressure, say, whose evolution through time is driven by external factors. This happens in many settings. For example, wind stress may be a factor in modelling ocean circulation [37], whereas light intensity, temperature, availability of food/nutrients, and wind speed may all be important in ecological modelling [38]. In our case, the equations depend on BG concentrations but in other areas of diabetes research insulin concentrations may be used, e.g. [39]. We cannot usually model such quantities but may be able to observe their values over a series of times. If we interpolate between the observations, then we can approximate the relevant quantity at any time within the observation period. If this observation period envelopes the time-frame over which we wish to evaluate our regression function, then solving the differential equations is still possible, and the interpolated series is referred to as a forcing function. WBDiff has not been designed with the specification of forcing functions in mind, but we show how they may be accommodated in Section 3.5.
As we shall see later, our regression function (also referred to herein as the 'dynamic model') is somewhat imperfect-there are clear fluctuations away from the fitted model apparent in the observed data. For the purposes of testing the artificial pancreas' control algorithm, it is important that we are able to predict realistic sensor data, and so we consider an AR model for the residual errors ( j above)-see Section 3.3. We fit the AR and dynamic models simultaneously, to fully account for uncertainty and also to capture posterior correlation between them. However, we find that they are somewhat confounded and identifiability issues arise unless the degree of autocorrelation in the AR process is constrained. Our approach to this is pragmatic and involves exploring various ways of limiting the extent to which the two models may interact, as discussed at the end of Section 4 and in the discussion.

Glucose kinetics
The CGM sensor measures glucose concentrations in the interstitial fluid, which can be related, mathematically, to the BG concentration via a compartmental model [40,41]: where IG denotes interstitial glucose. Hence, IG increases at a rate proportional to BG but is 'used up' according to a first-order process-the more there is, the faster it disappears. The sensor does not measure IG directly but, instead, measures electric current in the interstitial fluid and maps this to a scaled measure of IG via an assumption of proportionality. When the system is calibrated, the appropriate scale is chosen by equating scaled current with recent measures of BG. To account for this calibration we transform (3) to the same scale, by defining normalized interstitial glucose NIG = IG and choosing = p 1 / p 2 so that NIG is equal to BG at steady state: Let SG ij denote the jth measured SG concentration for individual i(i = 1, . . ., N = 12, j = 1, . . . , n = 228). Similarly, let BG il ,l = 1, . . ., m = 77, denote the lth measured BG concentration for individual i. Further, denote the times at which SG ij and BG il were measured by t ij and s il , respectively. A simple model for fitting individual i's data is then where NIG ij is the solution to (4) at time t ij . This is a deterministic function of three unknown inputs: (i) the value of p 1 ; (ii) the initial condition NIG(t = 0); and (iii) the form of BG(t). We assume that each individual has a distinct, but unknown, value of p 1 , which we denote by p 1 i . Often the initial conditions will be known, but in general they are not, and so these may also be treated as unknown parameters; in this case denoted by NIG 0 i , i = 1, . . . , N . Regarding the form of BG(t), we assume that linearly interpolating ‡ between the observed values for each individual, BG il ,l = 1, . . ., m, provides a satisfactory approximation, although see later for further discussion. Denoting the forcing function for individual i by BG i (t), we then have

Calibration
As demonstrated in Figure 1(b), the simple model above can perform poorly. With some careful thought, however, as to the nature of the underlying calibration mechanism, we can do much better. We stress, though, that the details of calibration are actually unknown to us-implementation details of the calibration procedure are proprietary and comprise guarded know-how by the respective CGM-system manufacturers to retain competitive advantage. In what follows, we make basic assumptions about how the process might work in order to construct a reasonable model. We first assume that SG is given by sensor, which we assume is subject to some error, , such that I m = I + , where I is the true current. We also assume that the true current I is related to interstitial glucose IG through I = IG/S + I B , where S denotes current sensitivity and I B represents a baseline current that is present even in the absence of IG. Hence, Let ik denote the kth calibration time for individual i(k = 1, . . ., K = 4). For convenience, we also define i0 = t i1 and i(K +1) > t in . We can then write down the calibration period to which each SG ij belongs as (In this paper, the calibration times are assumed known, but, in general, we may wish to acknowledge some uncertainty regarding their values.) A more realistic model for SG observations is then where B ik and F ik (k = 1, . . . , K +1) are unknown, individual, and calibration-period-specific parameters referred to as the calibration shift and calibration scale-factor, respectively. We consider two different models for the residuals ij . Assuming that the only source of 'error' is the measured current, then a homoscedasticity assumption leads to ij ∼ N(0, 2 i P[i, j] ), i = 1, . . ., N , j = 1, . . . , n, whereas an AR model gives Here i is an unknown individual-specific parameter controlling the degree of autocorrelation among the residuals. Note that in order to obtain the term, which is equal to one except at the calibration times, when it 'adjusts' the AR process, we make the assumption that current sensitivity, S, and do not change over time for a given individual. We cannot specify the autoregressive model, as given by (6)-(7), in BUGS, however, since a logical relationship for the response variable is not allowed. One way to get around this is to assume that each individual's SG-series arises from a CAR distribution [42,43]. A more flexible and intuitive approach, though, is to reexpress (6)-(7) as To define i1 , we could simply specify i1 = SG i1 −CIG i1 . However, note that the AR model does not penalize large s, and so unless we control their size, through an informative prior on i1 , say, they can become large and force the underlying 'model fit' {CIG ij , j = 1, . . ., n} away from the observed data, leading to implausible parameter estimates (see later for discussion). Note, though, that as the initial condition for the differential equation is unknown, we are free to choose the time to which it relates. If we choose the time of the first sensor reading t i1 , then we may express the initial condition deterministically: Hence, specifying a prior for i1 means that there is no need to model the initial conditions. Figure 2 shows a graphical representation of the full model in the case of autoregressive errors.

Priors
Calibration shifts and scale-factors may be correlated; in addition, scale factors must be positive whereas shifts may be negative. We therefore define C ik = log F ik B ik , i = 1, . . . , N , k = 1, . . . , K +1, which we assume arise from a bivariate normal 'population' distribution. If we believe that calibration parameters reflect characteristics of the individual (and/or sensor), then we may wish to assume individual-specific means i and an intra-individual covariance : C ik ∼ MVN 2 ( i , ). We may then wish to assume that the i s also arise from a bivariate normal distribution, with unknown 'global' mean and interindividual covariance : If, on the other hand, we believe that there is no correlation among calibration vectors for the same individual, then we might assume that it is the C ik s that are drawn from the right-hand side of (8) instead.
Directed acyclic graph (DAG) corresponding to 'calibrated' model with AR(1) process for the residual errors. For simplicity, the case in which there is only one calibration period for each individual is depicted. Each variable in the statistical model corresponds to a node and links between nodes show direct dependence. The graph is directed because each link is an arrow; it is acyclic because by following the arrows it is not possible to return to a node after leaving it. Square nodes denote known constants whereas circular nodes represent either deterministic relationships (i.e. functions) or stochastic quantities, i.e. quantities that require a distributional assumption. Stochastic dependence and functional dependence are denoted by solid and dashed arrows, respectively. Repetitive structures, such as the 'loop' from i = 1 to i = N , are represented by 'plates', which are nested if the model is hierarchical. The 'plate' in light-type on the right-hand side is shown to indicate the nature of dependence between successive observations. Nodes and BG i denote the entire set of population parameters, and the set of observed blood glucose concentrations for individual i, {BG il ,l = 1, ..., m}, respectively.
We assume fairly standard, vague (but proper) priors for , , and : bivariate normal and inverse-Wishart, centered at our best a priori guess with large variance. Throughout, the prior standard deviation specified for vague normal priors is 100 whereas Wishart priors are made as vague as possible by setting the degrees of freedom equal to the dimension, two in this case. The remaining parameters, , p 1 i , and the initial conditions NIG 0 i or initial residuals i1 (i = 1, . . . , N ), are transformed appropriately and assumed to arise from normal population distributions. Except for the initial residuals, these population distributions have unknown means and log-standard deviations with vague normal priors. The population mean initial residual is assumed to be zero and the population standard deviation is assigned an informative uniform prior on (0, 0.5), where the upper bound ensures that initial residuals greater than one are unlikely. The transformations applied are logarithmic for the residual standard deviations, initial conditions and p 1 parameters, and logistic for the i parameters (no transformation is required for the i1 s).

Implementation issues
Dynamic models in BUGS are 'packaged' in one of two ways [33]. One option is to specify the differential equations using the BUGS language and pass these as arguments to a generic ordinary differential equation (ODE) solver. The alternative is to edit and compile a template module for 'hardwiring' the ODE system into the software. In so doing, we create a new logical function in the BUGS language, which provides access to the numerical solution. We pass any parameters required to define the ODE system as arguments to the new function.
BUGS relies heavily on graphical modelling theory [29,30]. However, it is important to note that in graphical modelling terms, forcing functions are non-standard nodes. At first glance we might think they are logical nodes, since they are deterministic in nature. However, logical nodes are deterministic functions of other nodes in the graph: when they are defined in terms of 'time' they are functions of specific times, whereas forcing functions are defined for all times, in particular, they are functions of the dummy variable of integration, which, technically speaking, does not belong in the graph (Figure 2). In cases where the ODE-system is defined using the BUGS language, however, the integration-dummy is accessible within the model description, since it is required to define the differential equations. Hence, we can also use this to define forcing functions, linearly interpolating between the 'forcing data' via the BUGS language. One very significant advantage of hard-wiring the ODE-system, however, is that it is potentially much faster to compute. In this case, the integration-dummy is only available within the hard-wired module, and so we need to pass the forcing data, along with the times to which they relate, as parameters to the new function, and perform the interpolation within-see the Appendix.

Results
We begin by fitting the basic model given by (5), with the individual-specific parameters, p 1 i , NIG 0 i , and i , assumed to arise from log-normal population distributions with unknown means and standard deviations assigned vague normal and log-normal priors, respectively. We run 10 000 iterations in WinBUGS 1.4.3 [31,32] with WBDiff interface [33] installed. WinBUGS code for the main models considered in this paper is given in the Appendix.
Primarily, we work with 'hard-wired' systems of equations by writing and compiling specialized BUGS-modules, but, in this case, we also considered specifying the model entirely via the BUGS language. The former took around 5 min on a 2.13-GHz machine (when coded efficiently-see the Appendix), whereas the latter approach was substantially slower, taking ∼ 50 min. Point and interval estimates for the population parameters are presented in Table I, and a typical model fit is shown in Figure 1(b)-individual '10' was chosen as their data best illustrate the incremental benefit of increasing the model complexity.
Visual inspection of the model fits confirms that they are generally poor. Of primary importance in this paper is our ability to predict new data, and, to this end, the basic model is clearly inadequate. We might still wonder, however, whether it provides meaningful parameter estimates. From a clinical perspective, we are interested in the time delay that exists between glucose appearing in the blood and it then showing up on the sensor. This can be seen by looking at the relative positions of BG and SG peaks and/or troughs in Figure 1  get a rough idea of its size purely from the data, we performed a crude correlation analysis in which the SG values were lagged by 0, 5, 10 min, etc., and the correlation between BG and lagged-SG was calculated. The largest correlation coefficient (0.937) was obtained at a lag of 15 min, suggesting that a model-based estimate of the population median delay (from Table I) of exp(3.58) ≈ 36 min could be somewhat inaccurate.
To fit the 'calibrated' model given by (6) with ij ∼ N(0, 2 i P[i, j] ), we need to choose between the available exchangeability assumptions for the calibration parameters C ik , i = 1, . . . , N , k = 1, . . . , K +1. In particular, we wish to explore whether or not the calibration parameters reflect characteristics of the individual (and/or sensor), that is, whether to include individual-level means for these parameters. To address this we fit the model with and without individual means, and assess performance by looking at the posterior mean deviance, as a measure of model fit, and the Deviance Information Criterion, which penalizes the former by adding a penalty equal to the 'effective number of parameters' [44]. We find that there is no support for individual-level means as the mean deviance is virtually the same for both models, regardless of the effective number of parameters. Hence, we proceed with a three-level, as opposed to a four-level, model.
We ran 300 000 iterations of WinBUGS for the three-level calibrated model in a little over 3.5 h. To reduce the amount of computer memory required to store the output, we retained only every fifth sample for each parameter. The final 50 000 of the resulting 60 000 samples were then used for inference. The increased run-length here is due to a high level of autocorrelation in the output for the newly introduced C ik vectors. Parameter estimates for the population parameters are shown in the fourth and fifth columns of Table I. Now the population median time delay is 16.3 min, which is consistent with our crude empirical estimate. The inter-individual variability is relatively low, corresponding to a coefficient of variation, on the time-delay scale, of around 18 per cent. This is most usefully expressed, however, in the form of a prediction interval for new individuals' time delays, which accounts for uncertainty in the population mean and variability estimates: a 95 per cent interval is given by ( Although the model fit is much improved, we would still like, for prediction purposes, to be able to track the considerable fluctuations away from such fits that are apparent in the data. As can be seen in Figure 1(c) the residuals are serially correlated, and so we attempt to model them via the autoregressive (AR) process (7). However, a problem arises when we run the MCMC simulation. Recall that there is no penalty for large ij s, only their stochastic components ij (see (7)) need to be small. A good fit to the data, then, can often be obtained by setting i = 1 and choosing other parameters such that the underlying 'model fit' {C I G ij , j = 1, . . . , n} lies a roughly constant (with respect to time) distance away from the data. Then the residuals are all similar, consistent with i = 1, and each requires only a small stochastic component. We have some control over this phenomenon in specifying an informative prior for the initial residual i1 . However, it still occurs for some individuals unless we constrain the value of , and while the resulting model fits well in terms of { ij , j = 1, . . ., n}, the underlying CIG ij values are often implausible.
To impose the required constraint we assume that logit( i / max ), as opposed to logit( i ), arises from some normal population distribution (with unknown parameters), where max is assigned a specific value. However, now all of the individual i s are estimated equal to max , whatever value of max we choose. We address this by choosing the maximum value possible (in increments of 0.05) that still leads to plausible CIG-series for all individuals. Note that in so doing we find that there is no support for individual-specific i s, and so we also set i = ∀i, where logit( / max ) is assigned a vague normal prior. The value chosen for max was 0.8, and WinBUGS was again run for 300 000 iterations, retaining only every 5th sample. This took around 11.5 h, and point and interval estimates for the population parameters from the final 50 000 samples are presented in Table I. These are in good agreement with results from the previous model (without the AR process). The calibration shift B seems a little higher with less variability but we would expect from these figures that the underlying CIG-series are similar to before, as illustrated in Figure 3(b) for individual 10. Note that the population median residual standard deviation is now around half its previous value, at 0.118. This corresponds to the stochastic component of the residuals ij = SG ij − ij , indicating that the -series offer a substantial improvement over the CIG-series, as we would hope, and as is demonstrated in Figure 3 Finally, we wish to acknowledge some uncertainty regarding our choice of upper bound max . Choosing a value for max , it seems, is tantamount to fixing at that value. We would like to acknowledge that could lie between 0.75, which gives inferior model fits { ij , j = 1, . . . , n}, and 0.85, which leads to implausible CIG-series. But there is no point trying to acknowledge this uncertainty via a prior distribution, since we know that the posterior will be concentrated on the upper boundary (we may as well set = 0.85). Instead we specify as a 'distributional constant'-a fixed distribution as opposed to a fixed value (see [45], for example). We specify ∼ Unif(0.75, 0.85) but we prevent learning about from the likelihood-in graphical modelling terms, acts as a parent of {SG ij , i = 1, . . . , N , j = 1, . . ., n} but the SG ij s are not considered to be children of . WinBUGS code for 'cutting feedback' from the likelihood in this way is presented in the Appendix, and further comments on the use of such techniques are given in the discussion below. Results from this model give point estimates for population parameters identical to those obtained with max = 0.8, modulo Monte Carlo error, except for a very slight difference in the population median , which is still given by 0.118 to three significant figures.

Discussion
We have characterized the relationship between BG and SG, in children and adolescents, for the Guardian RT CGM system. Various hierarchical models have been explored to determine the most appropriate model, which assumes that calibration parameters for the same individual/sensor are not correlated. Our model can be viewed as a hierarchical, nonlinear regression, where the regression function is given by numerically solving a differential equation with accompanying initial condition, forcing function, and unknown parameter. While the forcing function is not naturally accommodated within the BUGS software, due to it representing a new class of graphical node, we have illustrated how it can be incorporated.
Even without modelling autocorrelation among the residuals our model fits remarkably well, given its simplicity. However, there are clear, unexplained fluctuations away from such fits, which we have accounted for by simultaneously fitting the nonlinear regression and an AR(1) process for the residuals. § While this has presented several practical challenges, we emphasize that the resulting model fits are very satisfactory (the majority of residuals-nearly 80 per cent-correspond to percentage differences between data and model fit of less than 2 per cent), and predictions from the model faithfully reproduce features present in the observed data. Sudden deviations away from the fitted curve, which are more prominent with other CGM systems, might be better handled by assuming heavier tailed, t-distributed stochastic residuals ( ij ) in the AR process, but we have not explored this yet. Note that Breton and Kovatchev [28] find the (unbounded) Johnson family of distributions [46] useful for the FreeStyle Navigator TM system. Another possible modification to the model for the residuals is to allow for different behaviour in hypoglycaemic, euglycaemic, and hyperglycaemic glucose ranges. However, any differences that could reliably be ascertained from our analyses were small, and it was thought that adapting the model would be of little practical benefit, especially considering that increasing the complexity of the residual model may exacerbate the apparent identifiability problem.
The motivation behind this work is to accelerate artificial pancreas development by providing a means of simulating large quantities of realistic sensor data, which can be fed into the control algorithm to test its response [17]. To this end, it is important to account for all sources of variability in glucosesensor data. In particular, we have estimated the population distribution of regression parameters, including their means/medians, their inter-and intra-individual variabilities, and any correlations that may exist between parameters. This allows us to generate 'virtual patients' by simulating realistic parameter sets from the population distribution. In addition, the Bayesian nature of our model enables full acknowledgement of the uncertainty associated with each population parameter estimate when simulating the virtual patients. Each virtual patient's parameters can be input, along with observed or simulated BG data, into the nonlinear regression function to derive CIG-values. However, the regression function is parsimonious and simply adding Gaussian white noise to this is insufficient for the purposes of generating realistic sensor data. Hence, characterizing the residuals has been a vital step in our analyses. The resulting parameter estimates can be used, again fully acknowledging their uncertainty, to simulate AR processes to be added to the derived CIG-series.
It is somewhat disappointing that we have had to constrain the degree of autocorrelation in the AR process and that the posterior distribution for is then always concentrated on the artificial boundary. This is mainly due, we think, to an inability to encapsulate within our prior distribution our 'common sense' knowledge as to the relationship between the observed data and the model-predicted CIG-series. Basically, we believe that the CIG-series should fit the data reasonably well and that the AR process should then account for mild fluctuations around that fit. But the unconstrained posterior is located near = 1 with CIG-series that often lie, implausibly, some (roughly) constant distance from the data. The resulting individual-level parameters, B, F, and p 1 are clearly inappropriate for the given individual, but do not necessarily have outlying values in terms of the population distribution of those parameters. Hence, it is not possible to circumvent the problem by constraining B, F, and p 1 . We have to remember that our model represents a gross simplification of the underlying process. To put too much emphasis on fitting the data when the model is known to be 'wrong' would be a mistake, in our opinion-the model is designed to provide meaningful parameter estimates and reasonably realistic predictions, which we believe it achieves.
One approach to avoiding the confounding problem would be to perform the analysis in two stages. We could first fit the calibrated model without the AR process and use this for inference on the parameters. We could then apply an autoregressive model to the residuals in order to characterize any fluctuations away from the deterministic model. However, this would ignore uncertainty in the model fit and prevent the model fit being adjusted, even slightly, to accommodate the different error structure. Hence, our efforts to fit nonlinear regression and AR models simultaneously. However, fixing max is tantamount to fixing and we would prefer to acknowledge some uncertainty regarding the latter. We achieve this by specifying as a 'distributional constant', as opposed to a fixed value, by placing a valve in the graphical model that allows information to flow from prior to likelihood but not vice versa. This allows us to be uncertain about without the model fit being 'tweaked' inappropriately. Point estimates provided by this approach are identical to those obtained with max = 0.8, but we would prefer to use the former for prediction as the level of uncertainty would be more realistic. Cutting the feedback from one or more sources of likelihood in such a way is growing in popularity-see, for example, [45, 47--54]. The motives vary, but it is typically used for 'multiply imputing' missing data or combining different sub-models that might otherwise be somewhat inconsistent, due to misspecification, say.
Other ways to cut feedback in the model begin with duplicating the SG data. A homoscedastic model could be specified for the first set of data and the same B ik , F ik , and p 1 i parameters could then be used to define the CIG ij s needed for fitting an AR model to the second set of data. Without appropriate valves/cuts in the graph, this would lead to excessive precision due to using the data twice. If we cut the feedback from the second set of data to the B ik , F ik , and p 1 i parameters, however, then  Table I). This is a Bayesian analogue of the two-stage approach described above where uncertainty in the model fit is now acknowledged. Further research is required to address whether or not this approach is preferable to specifying as a distributional constant. Another option would be to also create duplicate B ik , F ik , and p 1 i parameters, a set for each model fit, and to cut feedback from the AR set to the common population parameters (see [45], for example). This does not help for our data, however, as the parameters for the AR fit become implausible with unconstrained .
One area that we have not considered in this paper is uncertainty in the forcing function. Here we have simply linearly interpolated between a series of observed values, but it is natural to think that those observations might be subject to some error, resulting in a somewhat jagged forcing function. This is likely at odds with our prior beliefs as we might expect a biological process to be largely smooth. Moreover, if there is noise in the time-series then, presumably, we would rather use the underlying 'true' values. Hence, we may wish to specify a separate sub-model for the BG data. This is the subject of ongoing research and will form the basis of a future report. Another methodological issue is the numerical stability of ODE solving algorithms, which can be sensitive to the values of the input parameters, in some more complex settings precluding the use of vague prior distributions, say. The works of Ramsay et al. [55] and Campbell [56] offer an alternative approach that circumvents this problem. However, it is, as yet, unclear to us how this might be implemented in a flexible modelling framework such as BUGS or JAGS [57]. More robust solving algorithms for BUGS are currently under investigation.

Appendix A: WinBUGS/WBDiff code
WinBUGS code for the calibrated model (without AR errors) is presented below. Most of the code is self-explanatory but some notes, pertaining to the line numbers given in the right-hand margin, are provided below for clarity. (ii) the vector of times at which the solution is to be evaluated; (iii) a set of parameters required to fully specify the differential equation-these are defined on lines 11-15 (see below); (iv) the time to which the initial condition applies; and (v) the numerical tolerance to be used by the solving algorithm in determining whether or not the solution is sufficiently accurate-a value of 10 −6 is typical. Pseudo-code for the hard-wired function is given in Appendix B. • Lines 11-15: In order that our new, hard-wired component can evaluate the forcing function, we must supply it with the forcing data, and the times to which they relate. More generally, we may also need to supply the number of forcing data but here this is the same for all individuals, and so this information can be hard-wired instead of being passed as a parameter. The only other parameter required to fully specify the differential equation is p 1 i . • Line 19 onwards: Throughout, mean.x and prec.x denote the unknown population mean and precision of appropriately transformed x. • Lines 47-48: zero[], prior.prec.C [,], and prior.mat.C [,] are specified in the data set: they are given by a two-dimensional vector of zeros, 0.0001× I 2 , and I 2 , respectively, where I 2 denotes the 2×2 identity matrix.
The above model is adapted to incorporate autoregressive errors (with a common autocorrelation parameter for all individuals) by first making the j-loop on line 6 run from 2 to n rather than 1 to n and by replacing CIG[i,j] on line 7 with phi[i,j]. The following code is then inserted into that j-loop. Here the cut(.) function makes a 'copy' of the argument rho.dc. This always has the same value as the argument, and so uncertainty is propagated into the graph, but the 'cut' acts as a valve in the graph, preventing information from flowing through rho in the other direction (back to rho.dc). Hence, feedback from the likelihood is prevented and no learning about rho.dc can occur.

Appendix B: Pseudo-code for hard-wired differential equation with forcing function
Here we present pseudo-code for evaluating the differential equation(s) in our 'hard-wired' module.
(The reader is referred to the WBDiff documentation [33] for more general information.) The input parameters, including the forcing data, are available via a vector named theta and we are required to evaluate the equation(s) at some arbitrary time t. Note that the elements of theta are arranged as defined by the BUGS code in Appendix A, i.e. m forcing times, followed by m forcing data, followed by p 1 . Note also that the current solution at time t is also made available by the solving algorithm-this is denoted by NIG in the pseudo-code below. The main difficulty is in finding which elements of the forcing data are relevant to time t, i.e. which successive pair of observations lie either side of t, so that we can interpolate between them. We should be aware that the code will be called many, many times in solving the equation(s), and so it may be prudent to think about efficiency. Here we realize that each time the code is called, there is a very good chance that the 'forcing interval' will be the same as that most recently used or the subsequent interval. Hence, each time we evaluate the equation(s), we store, in a global variable named prev, the index of theta corresponding to the end of the forcing interval, so that we can refer to it next time.