The Parameter Houlihan: a solution to high-throughput identifiability indeterminacy for brutally ill-posed problems

TL;DR

The paper introduces the Parameter Houlihan, a novel method combining machine learning and data assimilation to select physiologically relevant model parameters, improving forecasting accuracy and addressing identifiability issues in clinical data modeling.

Contribution

It presents a new approach that integrates machine learning with data assimilation to optimize parameter selection, reducing identifiability problems in nonlinear physiological models.

Findings

The method improves forecasting accuracy over traditional parameter selection.

Houlihan-selected parameters generally minimize forecasting errors.

The approach offers a pathway to enhance physiological fidelity in models.

Abstract

One way to interject knowledge into clinically impactful forecasting is to use data assimilation, a nonlinear regression that projects data onto a mechanistic physiologic model, instead of a set of functions, such as neural networks. Such regressions have an advantage of being useful with particularly sparse, non-stationary clinical data. However, physiological models are often nonlinear and can have many parameters, leading to potential problems with parameter identifiability, or the ability to find a unique set of parameters that minimize forecasting error. The identifiability problems can be minimized or eliminated by reducing the number of parameters estimated, but reducing the number of estimated parameters also reduces the flexibility of the model and hence increases forecasting error. We propose a method, the parameter Houlihan, that combines traditional machine learning techniques with data assimilation, to select the right set of model parameters to minimize forecasting error while reducing identifiability problems. The method worked well: the data assimilation-based glucose forecasts and estimates for our cohort using the Houlihan-selected parameter sets generally also minimize forecasting errors compared to other parameter selection methods such as by-hand parameter selection. Nevertheless, the forecast with the lowest forecast error does not always accurately represent physiology, but further advancements of the algorithm provide a path for improving physiologic fidelity as well. Our hope is that this methodology represents a first step toward combining machine learning with data assimilation and provides a lower-threshold entry point for using data assimilation with clinical data by helping select the right parameters to estimate.

Tables

Table 1. Table 1: Demographic information and summary statistics are reported for the four participants whose retrospectively collected data are included in the study.

Data Summary
Participant ID	P1	P2	P4	P5
Age	$40-50$	$40-50$	$40-50$	$40-50$
Disease Status	T2D	T2D	No Diabetes	No Diabetes
Medications	metformin	metformin	—	—
Total $\#$ glucose measurements	$304$	$211$	$520$	$322$
Total $\#$ meals recorded	$124$	$76$	$370$	$184$
Total $\#$ days measured	$16$	$16$	$53$	$52$
Mean measured glucose	$113\pm 25$	$127\pm 32$	$92\pm 17$	$101\pm 16$

Table 2. Table 2: The rank ordering choice of the four parameter selection methods for the feature-metrics mean, μ 𝜇 \mu , and standard deviation, σ 𝜎 \sigma .

Rank-ordered parameters per selection method out of $21$ possible parameters
method	$1$	$2$	$3$	$4$	$5$	$6$	$7$	$8$	$9$	$10$	$11$
LASSO $\mu$	$a_1$	$C_1$	$V_p$	$t_p$	$R_m$	$C_3$	—	—	—	—	—
LASSO $\sigma$	$R_g$	$C_3$	$U_m$	$a_1$	$C_1$	$t_p$	$R_m$	$V_p$	–	—	—
Linear regression $\mu$	$a_1$	$C_1$	$C_3$	$R_m$	$t_p$	$V_p$	$U_m$	$R_g$	$C_4$	$U_b$	$U_0$
Linear regression $\sigma$	$R_g$	$C_3$	$U_m$	$a_1$	$C_1$	$R_m$	$V_p$	$t_p$	$k_{decay}$	—	—
Ridge regression $\mu$	$a_1$	—	—	—	—	—	—	—	—	—	—
Ridget regression $\sigma$	$R_g$	—	—	—	—	—	—	—	—	—	—
PCA $\mu$	$a_1$	$C_1$	$C_3$	$R_m$	$t_p$	$V_p$	$U_m$	—	—	—	—
PCA $\sigma$	$R_g$	$C_3$	$U_m$	$a_1$	$C_1$	—	—	—	—	—	—

Table 3. Table 3: The mean squared error (MSE) between forecast/smoothed and measured glucose. The machine-based methods, almost always selected the parameter set that achieved the MSE minimum, but for some individuals, certain hand-chosen parameters matter.

Rank-ordered parameters per selection method
	MSE for MCMC				MSE for UKF
parameter	P1	P2	P4	P5	P1	P2	P4	P5	method-feature-metric pairs
$a_1$	$822$	$1140$	$338$	$296$	$809$	$1270$	$304$	$356$	LASSO($\mu )$, LR($\mu$), PCA($\mu$)
$R_g$	$655$	$1180$	$475$	$288$	$672$	$1490$	$470$	$401$	LASSO($\sigma$), LR($\sigma$) , ridge($\sigma$), PCA($\sigma$)
$t_p$	$807$	$1020$	$448$	$349$	$788$	$1050$	$407$	$420$	by-hand, high-influence
$V_p$	$820$	$1120$	$332$	$320$	$805$	$1300$	$313$	$362$	by-hand, high-influence
$E$	$681$	$1250$	$655$	$500$	$721$	$1380$	$704$	$724$	by-hand, low-influence
$\alpha$	$501$	$1250$	$526$	$346$	$526$	$1580$	$528$	$394$	low-influence
$t_d$	$530$	$1080$	$730$	$674$	NaN	1260	NaN	$480$	low-influence
Rank-ordered parameter pairs per selection method
$(a_1,C_1)$	$570$	$1080$	$285$	$276$	$698$	$1290$	$258$	$330$	LASSO($\mu )$, LR($\mu$), PCA($\mu$)
$(R_g,C_3)$	$593$	$923$	$210$	$297$	$613$	$1260$	$215$	$385$	LASSO($\sigma$), LR($\sigma$) , ridge($\sigma$), PCA($\sigma$)
$(a_1,R_g)$	$578$	$1130$	$292$	$296$	$614$	$1400$	$269$	$343$	Union of rank $1$ over methods
$(\alpha ,E)$	$454$	$1174$	$518$	$345$	$483$	$1310$	$535$	$520$	low-influence
$(\alpha ,t_d)$	$432$	$993$	$525$	$347$	NaN	$1120$	NaN	NaN	low-influence
$(E,t_d)$	$462$	$1030$	$592$	$487$	$643$	$1190$	NaN	$490$	low-influence
Rank-ordered parameter $3$-tuple per selection method
$(a_1,C_1,V_p)$	$569$	$1060$	$284$	$276$	$663$	$1310$	$260$	$329$	LASSO($\mu$) ($1^{st}$)
$(a_1,C_1,t_p)$	$518$	$864$	$247$	$275$	NaN	NaN	$234$	$294$	LASSO($\mu$) ($2^{nd}$)
$(R_g,C_3,U_m)$	$590$	$922$	$190$	$294$	$618$	$1140$	$228$	$391$	LASSO($\sigma$), LR($\sigma$), PCA($\sigma$)
$(a_1,C_1,C_3)$	$431$	$1020$	$261$	$274$	$1330$	$1110$	$251$	$340$	LR($\mu$) PCA($\mu$)
$(C_2,E,\alpha )$	$442$	$1020$	$518$	$346$	$479$	$1250$	$535$	$515$	low-influence
$(t_d,C_2,\alpha )$	$432$	$894$	$525$	$347$	NaN	$1190$	NaN	NaN	low-influence
$(t_d,E,\alpha )$	$398$	$956$	$479$	$343$	NaN	$1120$	NaN	$520$	low-influence
$(t_d,E,C_2)$	$464$	$941$	$592$	$489$	$630$	$1190$	NaN	NaN	low-influence
Rank-ordered parameter $4$-tuple per selection method
$(a_1,R_g,C_1,C_3)$	$398$	$864$	$182$	$288$	$649$	$985$	$265$	$324$	Union of rank $2$ over methods
Full Houlihan for $\mu$ and $\sigma$
$(a_1,C_1,V_p,t_p,R_m,C_3)$	$414$	$862$	$217$	$229$	$661$	NaN	$236$	$291$	Lasso($\mu$)
$(R_g,C_3,U_m,a_1,C_1,t_p,R_m,V_p)$	$375$	$863$	$182$	$231$	$632$	$942$	$224$	$289$	Lasso($\sigma$)
Method with the lowest MSE
	Lasso	Lasso	Lasso/Union	Lasso	low-influence	Lasso	Lasso	Lasso

Table 4. Table 4: Full list of parameters for the ultradian glucose-insulin model Keener and Sneyd ( 2008 ) . Note that IIGU and IDGU denote insulin-independent glucose utilization and insulin-dependent glucose utilization, respectively.

Ultradian model parameters
Name	Nominal Value	Meaning
$V_p$	$3$ l	plasma volume
$V_i$	$11$ l	interstitial volume
$V_g$	$10$ l	glucose space
$E$	$0.2$ l min-1	exchange rate for insulin between remote and plasma compartments
$t_p$	$6$ min	time constant for plasma insulin degradation (via kidney and liver filtering)
$t_i$	$100$ min	time constant for remote insulin degradation (via muscle and adipose tissue)
$t_d$	$12$ min	delay between plasma insulin and glucose production
$k$	$0.5$ min-1	rate of decayed appearance of ingested glucose
$R_m$	$209$ mU min-1	linear constant affecting insulin secretion
$a_1$	$6.6$	exponential constant affecting insulin secretion
$C_1$	$300$ mg l-1	exponential constant affecting insulin secretion
$C_2$	$144$ mg l-1	exponential constant affecting IIGU
$C_3$	$100$ mg l-1	linear constant affecting IDGU
$C_4$	$80$ mU l-1	factor affecting IDG
$C_5$	$26$ mU l-1	exponential constant affecting IDGU
$U_b$	$72$ mg min-1	linear constant affecting IIGU
$U_0$	$4$ mg min-1	linear constant affecting IDGU
$U_m$	$94$ mg min-1	linear constant affecting IDGU
$R_g$	$180$ mg min-1	linear constant affecting IDGU
$\alpha$	$7.5$	exponential constant affecting IDGU
$\beta$	$1.772$	exponent affecting IDGU

Equations

$$||Y-\mathcal{F}(x,\theta )|{|}_Y^p$$

$$\underset{\Theta }{\min }(w_1(||Y-\mathcal{F}(x,\theta )|{|}_Y^p)+w_2(\#\{{\Theta }_i\}))$$

$$PCR(i)=\sum _j{\sigma }_{PC}(j)\ast P{C}_j(i).$$

$$\frac{d{I}_p}{dt}$$

$$\frac{d{I}_i}{dt}$$

$$\frac{dG}{dt}$$

$$\frac{d{h}_1}{dt}$$

$$\frac{d{h}_2}{dt}$$

$$\frac{d{h}_3}{dt}$$

$$f_1(G)$$

$$f_2(G)$$

$$f_3(I_i)$$

$$f_4(h_3)$$

$$\kappa$$

$$I_G(t)=\sum _{j=1}^N\frac{m_{j}k}{60}\exp (k(t_j-t));N=\#\{t_j<t\}$$

Full text

data assimilation; identifiability; machine learning; inverse problems; physiology; Markov Chain Monte Carlo; glucose-insulin

I Introduction

We want to use data and our understanding of the world to better manage health — we want evidence and understanding to guide clinical and personal health-related decisions. Of course at a high level this is generally what medicine is about: interventions are undertaken only when they are understood or predicted to improve an individual’s health. However, traditionally this prediction is done in a non-personalized manner, meaning that interventions treat the ”mean” person or patient. Personalized and precision medicine were conceptualized to relax this constraint by tailoring an intervention to a person. While genetics offers a path to personalizing treatment, we can also use data science machinery together with personal (Albers et al. (2017)) and population-scale data to better personalize treatment (Albers et al. (2012); Xu et al. (2016); Albers et al. (2014)). Specifically here, we want to leverage our knowledge encapsulated in mechanistic physiologic models and combine it with free living or clinical data to allow this knowledge and data to be used to make decisions related to health. In this context, computational problems related to personalized medicine can be broken into two broad categories: forecasting, where we make quantitative predictions about a patient’s future state that can be used by clinicians and patients to take corrective action, and phenotyping (Hripcsak and Albers (2017); Pathak et al. (2013); Pivovarov et al. (2015); Halpern et al. (2016); Rajkomar et al. (2018)), where we identify properties of macroscopic observables that can be used to classify patients into subgroups that can give clinicians and researchers actionable insight into commonly occurring treatment outcomes and biological phenomena.

The idea of using mechanistic models and data assimilation in biomedicine or healthcare is old, but what is new is attempting to integrate models with variable complexity with sparse, noisy free-living and clinically collected data. Many mathematical biology models were designed to have variable degrees of biological fidelity, fidelity that we do not necessarily want to eliminate or reduce, but fidelity that we generally need to constrain in the usual case where we cannot estimate all the parameters because of data limitations that always exist in practice. This problem poses a significant barrier to using data assimilation—enough of a barrier that often data assimilation is not even attempted because the models, given data are hopelessly poorly resolved. This paper poses a machine learning solution to this problem—by using machine learning to identify and rank-order which model parameters are the most necessary to estimate.

Returning to the more practical contexts of phenotyping and forecasting, both applications impose particular demands on certain aspects of computational machinery used to model data. The properties we focus on here are the selection of the model parameters to estimate and the ensuing identifiability of a model, or ability to uniquely solve for parameters that yield optimal solutions (Westwick and Kearney (2003); Ljung (1987); Schoukens et al. (2016)). Our goal is to strike a balance between identifiability and model fidelity in situations where a model is not fully identifiable if all or even sometimes when any model parameters are estimated, given the available data. The method we develop here can facilitate both forecasting and phenotyping studies, and we evaluate this method in the context of modeling glucose dynamics using mechanistic models, machine learning and data assimilation.

The Houlihan, or the Houlihan throw, is a lasso throw used for roping livestock, e.g., a horse. It is used often under difficult circumstances such as picking out, from a substantial distance, a single horse from among a crowd of horses standing close together. It is a particularly flexible technique that can be used in a variety of circumstances. In this spirit we intuitively define the Houlihan method(s) as a collection of methods that use for selecting the most productive model parameters to estimate; specifically, the collection of methods uses machine learning techniques applied to simulated model output under parameter variation subject to a set of features, e.g., the mean of a state.

II Background

The larger biomedical context of this work is the application of data science machinery used to personalize forecasts and phenotypes via a broadly defined regression. While there are many linear versions of regression that have been successfully applied to healthcare data (Levine et al. (2016, 2018); Hripcsak et al. (2011); Hripcsak and Albers (2013)), here we focus on a specific type of nonlinear regression—data assimilation—in an effort to take advantage of potentially important nonlinearities present in most biological systems. Nonlinear regression approaches such as deep learning and related methods (Perotte et al. (2015); Ranganath et al. (2016); Rajkomar et al. (2018); Lasko et al. (2013)) have seen some success in a number of biomedical applications thanks to their ability to approximate arbitrary, non-linear functions. While the flexibility of universal approximator approaches (Hornik et al. (1989, 1990)) is particularly useful when little is known about the system and data are plentiful, this approach does not always work well when data are sparse and non-stationary, leading to problems such as poor generalization to new or unobserved individuals, problems with quickly changing health conditions, and difficulties with fast, accurate prediction with very few, e.g., $20$, data points. Unfortunately, many health data and healthcare situations fit one or more of these data pathologies (Hripcsak and Albers (2012); Pivovarov et al. (2014)).

In order to exploit the complex yet rich quantity of available health data, it is natural to consider ways of constraining the search space for machine learning methods. One way to do this is to constrain the model search space in accordance with as much expert knowledge as possible. To achieve this here we turn to mechanistic models developed by mathematical biologists and physiologists Keener and Sneyd (2008), which are typically formulated as dynamical systems (Brin and Stuck (2004); Guckenheimer and Holmes (1983); Arrowsmith and Place (1990)), e.g., $x_{t+1}=f(x_{t},\theta)$, or differential equations (Arnold (1992, 1983)), e.g., $\frac{dx}{dt}=f(x,\theta,t)$, where $x$ are the time-varying states of the system and $\theta$ are the physiologic parameters that govern the process. For example, in the case of phenotyping type 2 diabetes one way of constraining the search space of a regression is to regress the data onto a nonlinear physiologic model Albers et al. (2017, 2018) instead of regressing the data onto a universal approximator Hornik et al. (1989, 1990) function space such as neural networks. The way this is done is using data assimilation.

Data assimilation (DA) is a collection of methods (Jazwinski (1998); Lorenc (1988); Law et al. (2015); Ash et al. (2016); Reich and Cotter (2015); Haug (2012); Ristic et al. (2004); Candy (2009); Evensen (2009, 2003)) concerned with performing the types of non-linear regressions we describe for dynamical systems, and centers itself around forecasting and inferring mechanistic states under available observations; it solves both forward and inverse problems (Stuart (2010); Zenker et al. (2007); Banks et al. (2014)). There have been many successful applications of mechanistic modeling and data assimilation in biomedicine (Albers et al. (2018); Hirata et al. (2010a, b); Schiff (2011); Dukic et al. (2012); Miao et al. (2011); Chee and Fernando (2007); P et al. (2003); Fabietti et al. (2007); Kovatchev et al. (2009); Coutinho et al. (2003); Mirowski et al. (1980); Thabit et al. (2015a); Leelarathna et al. (2013); Cobelli et al. (2009); Thabit et al. (2015b); Glass and Courtemanche (2000); Christini and Glass (2002); Hall et al. (1997); Mackey and Glass (1977); Zenker et al. (2007); Parker et al. (20xx, 2001); Donnet and Samson (2013); Bonate (2005); Sadean and Glass (2009); Kristensen et al. (2005); Sedigh-Sarvestan et al. (2012); Selgrade et al. (2009); Sedigh-Sarvestani et al. (2012); Lin et al. (2004, 2011); Albers et al. (2017); Levine et al. (2017)). However, mechanistic models that are typically developed in biological laboratory settings are often not designed to interface with health data collected in the process of delivering care or in free-living situations—in particular, the physiologic models often model macroscopic states that are observable from routinely collected data but are governed by a composition of unobservable mechanisms. While these models capture the dynamics we are interested in and constrain the regression to a smaller class of functions, their high-fidelity creates issues of identifiability and ill-posedness, problems for which this paper develops a practical, machine-learning-based work-around.

To understand how identifiability works for these machines, consider a trivial case of identifiability for the model $\frac{dx}{dt}=abx$. If we assume that $a$ and $b$ are unknown parameters, they cannot both be identifiable without another equation that could uniquely determine one of them. This topic and the the associated methods for handling this situation are too old and wide ranging to give complete background (Ljung (1987); Westwick and Kearney (2003); Hines et al. (2014); Schoukens et al. (2016)). We can, however, give a broad sketch of how identifiability has been traditionally approached. Identifiability analysis generally follows one of three pathways: analytical methods, e.g., showing algebraically that all parameters can be uniquely solved for (Simens et al. (2017); Eisenberg and Hayashi (2014); Walch and Eisenberg (2016)); numerical methods (Schoukens et al. (2016); Westwick and Kearney (2003); Brouwer et al. (2018)); and heuristic, knowledge-driven sensitivity analyses where certain parameters are chosen based on computational experiments or knowledge of the system. In many complex, non-linear mechanistic physiologic models algebraic methods and linear computational methods are not tractable or applicable. In these situations nonlinear methods can be applied, but nonlinear methods usually have to be constructed for a particular situation (Chis et al. (2011)), and, much like nonlinear optimization, generally do not have clearcut or simple resolutions (Chis et al. (2011); Hines et al. (2014)). These problems pose a significant roadblock to parameter inference in the context of DA. Nevertheless, there exist methods for working to remain within a traditional identifiability framework, e.g., Hines et al. (2014) uses Bayesian inference to determine when parameters can be made identifiable.

The usual way of addressing these issues focuses on making sure the model is identifiable or finding ways of making it more identifiable (Hines et al. (2014); Schoukens et al. (2016); Westwick and Kearney (2003)). This work is often performed using substantial intuition about the important features encoded in the model, and parameterizing and grouping sub-processes. However, this creates silos of expertise and prevents wide-spread dissemination and evaluation of mechanistic models in potential application domains. Therefore, to progress toward understanding complex physiology via model refinement and selection, and to provide solutions in clinical situations that come with constraints of time-sensitive solutions, we must find a robust way of coping with brutally ill-posed problems and accept certain impurities and inaccuracies.

Here, we develop and evaluate a method for rank-ordering mechanistic parameters based on their ”influence” on important dynamical features, in order to improve forecast accuracy and help determine which models most faithfully represent a given system. This provides a starting point from which to estimate parameters, prune the model, etc., that can be automated.

III Conceptual construction of the Houlihan approach

III.1 Conventional operational use of data assimilation with ill-posed problems

The standard method of applying data assimilation (DA) or control in generic situations follows roughly the following steps: (i) select a model, (ii) work out identifiability, (iii) select a filter or inference method, (iv) find an optimal solution for states and parameters. This requires very careful experimental constructions, generally dense data streams, can be expensive, and requires relatively simple models, all situations that lie outside of what is possible in applications and even many basic science settings. The approach for applying DA in operational, complex, high-dimensional settings where accurate real-time forecasts are imperative is to: (i) select or develop a model, (ii) tune and fix parameters offline, often by hand or using a combination of by-hand and numeric tuning that allows the model to reconstruct or forecast states within some tolerance, (iii) select an inference scheme, and (iv) estimate states only and make a forecast. This is a tried and true method and is used in situations such as weather and climate forecasting (Hamill (2006); Evensen (2009, 2003)). Neither of these approaches apply to biomedical situations that, by comparison, have a different set of constraints and problems, including: (i) the models are smaller, so they can be simulated faster and estimated faster, allowing for potentially many models to be used simultaneously; (ii) there are less data relative to the number of unknown parameters, so while parameter estimation is necessary Albers et al. (2017) not all parameters can be estimated; (iii) models are not generated from first principles and their application to given individuals is potentially highly variable necessitating the use and potentially the averaging of many models; and (iv) tuning would have to be done for millions of people frequently, e.g., every patient in every ICU potentially every day, a process that is not likely to be practically possible. Because of these reasons, choosing which parameters to estimate is a significant barrier to the adoption and use of DA in biomedical situations.

III.2 Houlihan approach to ill-posedness

Here we are operating under a different situation from the more canonical DA application setting, one more heavily constrained by imperfections of data that will never disappear because the data are collected in the process of managing health instead of data collected in a controlled manner explicitly for the DA. In particular we assume: (a) we do not know the right model but we have some models we can try, (b) we do not know whether a given model is identifiable and that we do not have enough data to estimate all model parameters well anyway but that we have enough data to estimate at least one parameter, (c) for a given model, we don’t know what parameters are the most useful to estimate, given that we cannot estimate all of them. Given this situation we develop a method for rank ordering which parameters to estimate, subject to features we want to capture, when we have no idea how to choose which parameters to estimate or when we must choose parameters in a more high-throughput setting where we are using many models at once.

This solution involves stacking machine learning on top of DA: machine learning is applied to simulated model output to select the important parameters to estimate to best synchronize the model with the data, and then we use DA restricted to estimation of the parameters chosen by machine learning. In this way, the method will scale to a high-throughput setting and can be applied to many different models with high dimensional parameter spaces more easily. And while we know that this method may not lead to a unique solution in function, parameter, or initial condition space, the set of solutions will be reduced to a workable set of solutions that allow forward progress to be made.

Conceptually, we are proposing to: (i) assume a model, (ii) simulate the model under discrete parameter variation creating a grid in parameter space for which at every point we have simulated data from the model (i.e., the instance of one attractor of the model for a given set of initial conditions at the parameter grid point), (iii) select features, e.g., the mean, of the attractors that are important for estimating the physiologic system, and then (iv) use a machine learning algorithm to identify the parameters that have the greatest impact on the features. While for the authors the geometric intuition of this method originates from bifurcation theory—we will discuss this in a later section Kuznetzov (1998)—one useful way to think about the problem is in the inverse problems context. As was the case for the bifurcation theory context, this discussion is allegorical; we are not proposing a formal inverse problems regularization framework. From a high level, given data, $Y$, and a model, $\mathcal{F}$ with a state space $x$, the task is to find a set of parameter values, $\Theta_{i}$, of which there may be many if the system is not identifiable, that minimize:

[TABLE]

for some $p$, $p=2$ being the commonly applied least squares minimization. The core of the identifiability issue is that, for complex models, and especially given sparse data, there may be many sets of parameters $\Theta_{i}$ that minimize the distance between the model and the data. In this case a goal might then be to balance the number of potential minimizing parameter sets, the number of $\Theta_{i}$’s, against the distance between the model and the data via an optimization algorithm, e.g.,

[TABLE]

where the $w_{i}$’s are continuous functions. This framework, a formal regularization methodology, has many advantages, but can induce many complexities that increase rather than decrease the barrier to using data assimilation in more data-poor environments. Moreover, this relatively complex methodology may not be applicable in more high-throughput situations where, e.g., many models are used in a model averaging context. Therefore, motivated by the goal of an imperfect but practical solution, we postulate that if we carefully select the right parameters that maximize the parameter subspaces that can be explored relative to a set of desired features, we can often, effectively but imperfectly, solve the optimization problem. Effectively but not rigorously, we are regularizing a priori, by selecting and reducing the parameter set to be estimated before we go about estimating the parameters given data. Given the framework above, such a solution may be well handled by a tool from sparse machine learning such as lasso Hastie et al. (2015) because it uniquely rank-orders parameters by their predictive power, but it is easy to imagine using other methods. But, it is important to be clear that we are hypothesizing that the parameter subspaces that allow maximal exploration of dynamics relative to a given feature, e.g., the mean, will contain sets of parameters, $\Theta_{i}$ that also find relatively good minima of Eq. 1. In our evaluation we will see cases where this hypothesis fails, but we will also see that this hypothesis generally holds true in our data set, and this conclusion is the point of the paper.

In short, here we are assuming a problem is ill posed and a system that is likely not identifiable, and given this situation, we are trying to cope. Therefore, we are not really solving an identifiability task because we are not trying to find the best or most representative model that admits unique parameter estimates; rather we are solving a problem more akin to, but not literally, a regularization task. We are starting from a point where the problem is both brutally ill-posed and likely non-identifiable, and where investigating identifiability using analytic methods, or even many numeric-by-hand methods are intractable. In this case we are assuming there will be a few different parameter combinations that represent reasonable parameter estimates. In this situation each combination of parameter represents a hypothesis for how the system works. More importantly, the method we present here is a flexible entry-point for using data assimilation with a complex nonlinear model and data collected in an uncontrolled environment rather than directly solving an identifiability problem.

IV Data cohort

We test and evaluate the Houlihan methodology in the context of modeling and forecasting blood glucose collected in a free-living setting — via a type 2 diabetes self-management moblie application. The dlood glucose and nutrition data used here were collected retrospectively from four participants, two with type 2 diabetes and two without diabetes, using custom-designed mobile applications for capturing self-monitoring data (Mamykina et al. (2015)). These data are summarized in Table 1. We acquired two types of data: 1) fingerstick blood glucose measurements taken at the discretion of each of the 4 participants (roughly 3-10 times per day) and 2) estimates of carbohydrate consumption over time (roughly 1-5 meals per day) determined by a certified dietitian’s analysis of the daily meal logs (with photos and descriptions) reported by each participant.The data are documented more completely in (Albers et al. (2017); Levine et al. (2017)) and are available on PhysioNet upon request.

V Methods

V.1 Glucose-insulin physiologic model

The Houlihan method was conceived in the context of DA with a mechanistic model, and while it could be used in any nonlinear regression context, this paper will be restricted to the setting where we begin by projecting data onto a mechanistic dynamical system and then work to decide which parameters of that dynamical system we should estimate to represent the data. The mechanistic model is more formally either a dynamical system when time is discrete or a system of ODEs when time is continuous. Explicit versions of such systems form parameterized families of functions that are physically meaningful but generally do not satisfy nice function space properties such as completeness and are not universal approximators. The more general theory of dynamical systems can be found in many books (Brin and Stuck (2004); Guckenheimer and Holmes (1983); Arrowsmith and Place (1990)), but here we will restrict our use of these details to an absolute minimum. We will assume that the systems we use have at least one invariant density; the invariant density is likely defined relative to a SRB-measure (Young (2002); Ruelle (1976); Sinai (1972); Bowen (1975)) rather than Lebesque measure, but the point is that for a given set of parameter values and initial conditions, the states have a probability density function associated with them denoted $\Lambda$. This invariant density can potentially depend on both the parameters and the initial conditions for a set of parameters.

As previously noted, we want to use DA to model the glucose-insulin system of a human being. We begin with a particular mechanistic glucose-insulin model, here the ultradian model that has been detailed in Sturis et al. (1991); Keener and Sneyd (2008); Albers et al. (2012, 2014, 2017), and has $6$ states and $21$ parameters; its details can be found in the appendix A. The model has unknown identifiability properties, especially when only glucose is measured, but we have strong evidence that at least some of the model parameters and states are not identifiable (Levine et al. (2017)). The Houlihan method rests on quantifying how the invariant densities of the synthetic data sets and their properties vary as parameters of the mechanistic model(s) vary. Specifically, the Houlihan method decides which parameters to estimate by varying the parameters of the ultradian model, observing how the invariant densities and their properties vary, and then using this information to select parameters to estimate by ranking ordering their importance using statistical inference or machine learning. The synthetic data used to select parameters to estimate will be generated by solving the ultradian model using an adaptive version of Runga-Kutta, ode23 in Matlab and will consist of $10^{5}$ simulated data points.

V.2 Stochastic filtering and inverse problems methods

We use two previously documented data assimilation formulations, an unscented Kalman filter (Julier and Uhlmann (2004); Julier et al. (1995); Want and Merwe (2000); Gove and Hollinger (2006); Want et al. (2000); Wan and Merwe (2001)) (UKF) whose details can be found in Albers et al. (2017) and a Metropolis-within-Gibbs Markov Chain Monte Carlo (MCMC) method whose details can be found in Levine et al. (2017); Cotter et al. (2013). As previously mentioned, these DA methods are used with the ultradian model (Sturis et al. (1991)) for performing the DA tasks. We only use these methods over the course of evaluating the Houlihan methods; the exact implementation of the DA methods can be found in Albers et al. (2017); Levine et al. (2017).

V.3 Analytical construction and intuition for throwing the Houlihan around the right parameters

While the approach we are proposing is new, the allegorical geometric intuition motivating this approach comes from bifurcation theory and in particular the bifurcation sets defined in the 1970’s (Sotomayor (1973)) and the analytic geometry vision of bifurcation theory and singularities in parameter space Arnold (1983). Bifurcation sets are the low-dimensional sets or manifolds that denote transition/bifurcation surfaces between topologically equivalent invariant sets, partitioning the parameter space into a set of equivalence classes. It is this idea of partitioning the parameter space into equivalence classes that differently impact dynamical featuers we care about is they key motivational insight. In our context we want to partition the parameter space by influence on some feature or set of features, denoted the feature-metric, of the dynamics. Feature-metrics are calculated from the time-series of the simulated model (dynamical system), e.g., a mean. We do not want to be as rigid as requiring topological equivalence as was defined in the bifurcation sets framework, or necessarily strict classes, but we do want to partition the parameter space according to how parameters influence a dynamical feature we care about. The over-arching idea is that the subsets of parameter space that have the highest influence on the feature-metric are the parameters that will be the most useful to estimate to minimize Eq. 1. And, knowing the most useful parameters to estimate provides a systematic way of choosing the parameters to estimate until the system is either identifiable or identifiable enough to be serviceable; in practice serviceable might mean that the errors are within desired tolerances, that parameter estimates are unique, or that the parameter estimates have few enough equilibria or minima that they can be made useful. To make this more precise, begin with the following terms, which are functions of a parameter vector, $p$.

Feature metric: the feature of the dynamical system we wish to influence, denoted $g(p)$; feature metrics are estimated from the time-series of the simulated model output and vary with parameter variation.

Influence: the amount that a parameter influences the feature-metric, denoted $F_{i}(g(p))$ for the i-th parameter.

Influence equivalence: a rule that defines equivalence of influence, e.g, all parameters $i$ such that $a_{j}\leq F_{i}(g(p))<a_{j+1}$. This allows for us to introduce a partition over influence, $\{a\}_{j\in J}$ called an impact set, which represents the transitions or boundaries between influence equivalence classes.

Parameter influence sets: the sets of parameters with equivalent influence according to the influence equivalence rule.

Demonstrative example: Begin by defining the dynamical system $f$ with state variables $x_{i}$ and parameters $p_{i}$ assuming at least four parameters. Next define the feature-metric as the mean of a single state variable $x_{*}$, $\mu_{x_{*}}$ (i.e. we are interested in how each parameter ”influences” the state’s mean). Set the influence function to be the absolute linear correlation, $|\beta_{i}|$, between the feature-metric, $\mu_{x}$, and values of the parameter $p_{i}$. In this example, the influence function is a vector-valued function, with a scalar metric (linear correlation between parameter and the state’s mean it induces) corresponding to each parameter. The influence per parameter defines a probability mass function (PMF) with support $[0,1]$ with values $\frac{|\beta_{i}|}{\sum_{j}{|\beta_{j}|}}$. Finally, we define influence equivalence as membership in a given quartile of the PMF defined by the influence function. Note that the impact set is defined by the PMF quartiles, and the influence sets are the parameters in respective quartiles of the PMF. Depending on the separation observed in the impact sets, we could ultimately choose to estimate parameters only from the upper equivalence class(es); i.e. the set of parameters with $|\beta_{i}|$ in the upper quartile. $\Box$

This example takes a narrow interpretation of the flexible construct we develop for identifying equivalence classes of parameter influence. However, even the above example allows for wild topological variation within a given equivalence class. For example, within a given equivalence class one would easily imagine there being many topologically distinct invariant sets due to both parameter variation and initial condition variation. Presumably there are other similar equivalence class violations such as ergodicity properties (Pugh and Shub (1989); Burns and Wilkinson (2010)), $k-LCE$ stability (Albers and Sprott (2006); Albers et al. (2006)), etc. These issues can all be addressed by defining the various properties, e.g., the influence function, differently, or more restrictively such that we end up with increasingly more restrictive constructions such as the original notion of bifurcation sets. This flexibility in equivalencies is the point of this construction: we can, depending on our goals, data, etc., have substantial flexibility in how we set up how to choose what parameters to estimate all while explicitly acknowledging what we know we do not know we are preserving. For example, if we define the feature-metric to be the mean, we know we are allowing the system to explore or have many different coexisting invariant densities as long as they have a mean that lies within a given equivalence class.

Visual example: Figure 1 shows three cases of the outcome of the Houlihan analysis. The left-most plot in fig. 1 shows the case where the rank-ordering of influence is on a by-coordinate basis; meaning, the equivalence classes were collections entire coordinates, here where each equivalence class has a single member. The middle plot in fig. 1 shows a case where the influence can be portions of different coordinates, but still there are is not joint dependence between variables. The right-most plot in fig. 1 demonstrates an example where the influence equivalence includes joint coordinate relationships. In this paper we will only address the first of these cases, leaving the more complex situations for later work.

V.4 Computational moving parts for throwing the Houlihan around the right parameters

The computational task of selecting parameters to estimate involves defining the equivalence-like classes, finding their boundaries, rank ordering the parameters by importance and has, broadly, five moving parts. First, select the feature-metric(s), $g(p)$, e.g., mean. Second, formulate the representation of the space of parameters and their variation, including (i) parameter grid resolution, (ii) parameter perturbation range, (iii) parameter variation type, e.g., joint versus individual by-parameter parameter variation. Third, choose an influence function that defines how to model the parameterized variation of the feature-metric variation with parameter variation. Fourth, choose a method for rank ordering these parameterizations by influence. Sometimes steps three and four can be done using a single method, e.g., linear regression with a $L_{1}$ regularization or by using lasso with cross validation, and sometimes it is done in two steps, e.g., linear regression with a threshold on the $\beta$’s, partitioning the $\beta$’s into equivalence classes. And fifth, decide which parameters to keep or which equivalence classes, or which impact sets are important.

Feature metrics

We use two feature metrics, mean and standard deviation of the invariant density generated by mechanistic model with set parameter values and initial conditions.

Parameter grid

We begin with the nominal parameters (Sturis et al. (1991); Keener and Sneyd (2008); Albers et al. (2017)), and then vary them in intervals of $\log_{2}$ over $10$ decades in both directions. For example, for parameter $i$ the parameter grid point for the $k^{th}$ decade was set as $p_{i}(nominal)2^{k}$. We did not consider joint-variation of parameters, but varied parameters independently while holding all other parameters fixed at their nominal values.

V.4.1 Parameter selection methods: Influence functions, impact sets, and ranking

Given a feature metric as a covariate or input vector, e.g., the means of attractor densities for a set of parameter values, we use several methods for selecting the best set of parameters to estimate in a DA context. Some of these methods are stock—linear regression with lasso—some are standard practice—parameter selection using knowledge of the model—and some are modifications of existing methods—see PCA-lariat below. We will see that the method for selecting the parameters matters, although not as much as the feature metric, and it is clear that sophisticated machine learning methods could be useful in this context.

Covariates or input vectors

All of the methods below take a covariate matrix as input. The covariates correspond to vectors: one dimension of the covariate matrix corresponds to a feature metric calculated at every point along the parameter variation, e.g., the mean of a simulated attractor at every point along a one-dimensional parameter curve.

By hand selection parameter selection — parameter selection using knowledge

In our previous work we selected parameters to estimate by hand as they were tied to certain dynamical features, physiologic knowledge want to fit something in particular to solve a problem, e.g., phenotyping. We selected $E$ and $V_{p}$ because they seemed to have an impact on the mean (Albers et al. (2014)) and $t_{p}$ because it was related to liver function; the results can be found in Albers et al. (2017).

Automatic parameter selection using linear regression

A basic method for determining influence is the linear dependence between the feature metric and parameter variation. In this setting we perform a linear regression between the feature metric and the parameters and we keep all $\beta$’s for which $\beta_{i}>(\beta_{1})(\kappa_{LR})$. Here we set $\kappa_{LR}=20\%$ or $0.2$, meaning that we keep all the parameters that have a regression coefficient that explains at least $20\%$ of the regression coefficient with the highest influence.

Automatic parameter selection using Lasso and cross validation

A natural way of reducing the number of parameters in a model is to select parameters that have a lot of power explaining the feature metric while simultaneously being non-redundant. One way of achieving this is to use lasso, or $L_{1}$ regularization to enforce a sparse representation of the parameter system (Hastie et al. (2015); Tibshirani (1996); Breiman (1995); Zou and Hastie (2005)). We use the standard lasso formulation (Hastie et al. (2015)) with cross validation to determine the rank-ordering of parameters; the optimal value of $\lambda$, or the optimal number of parameters, is set using a cutoff of one standard error. Lasso automatically and uniquely rank orders parameters. We keep the parameters within one standard error of the minimum mean squared error (MSE) ensuring a sparse representation of the model.

Automatic parameter selection using elastic net approximation of ridge regression

In addition to lasso regularization, we also use ridge regression, or $L_{2}$ regularization (Hastie et al. (2015); Zou and Hastie (2005); Jaggi (2014)). We compute the ridge regression selected parameters using an elastic nets formulation with $\alpha$ set to $0.0001$ where elastic nets formulation approachs $L_{2}$ regularization, and select the number of parameters using cross validation in the same way as is done in the lasso setting. We keep the parameters within one standard error of the minimum mean squared error (MSE) ensuring a sparse representation of the model.

Automatic parameter selection using PCA-lariat with a single metric

To add diversity to the set of methods for selecting parameters beyond linear regression-based methods, we devised a principle component analysis (PCA) (Pearson (1901); Hotelling (1935); Jolliffe (2010)) based algorithm for computing an influence function, then implement a rank-ordering scheme for defining influence equivalence. The method we develop, PCA-lariat, follows seven steps. First, estimate the PCs for the feature-metric, $g(p)$, taking care to de-trend the summary. Second, estimate the percentage of the variance captured by the $i-th$ PC, $\sigma_{PC}(i)$. Third, identify the important PCs, or the PCs that explain variance above a threshold, $\kappa_{PC}$; we use $5\%$. Fourth, for each important PC, rank-order the contribution of each parameter or coordinate to the PC. Fifth, collect all the coordinates for all the important PCs that contribute proportionally to a given PC above a set threshold, $\kappa_{C}$, $PC_{j}(i)>\kappa_{C}$; we use $10\%$. Sixth, for the important parameters for the important PCs, estimate the contribution per parameter:

[TABLE]

And seventh, rank order the important parameters by $PCR$ and select the parameters above a given threshold, $\kappa_{I}$; we use $0.1$, or $10\%$.

Multi-directional parameter wrangling

Combining models, or model averaging can be very useful for improving results (Hoeting et al. (1999); Raftery et al. (2005); Claeskens and Hjort (2008); Albers et al. (2017)), especially when you either know you want to adjust to multiple feature-metrics, or you do not know what feature metrics are important. Here, we only consider using set operations over methods, and consider three cases. First, we take the union of: (number of rank-ordered parameters, feature-metric, influence function) using one parameter per influence function, two feature metrics, mean and standard deviation. Second and third, we take the union of: (number of rank ordered parameters, feature-metric, influence function) using one parameter per influence function and one feature metric, either mean or standard deviation.

V.5 Evaluation scheme

The evaluation of the Houlihan methods is done in four steps. First, we apply the Houlihan methods to the ultradian model to select parameters to estimate and compare the parameter selections as the method is perturbed. Second, we use both the UKF and the MCMC DA methods to estimate these Houlihan-selected parameters for the four people in our cohort and calculate the mean squared error (MSE) between the data and the model state estimates (MCMC methods) and forecasts (UKF methods). Third, we use both the UKF and the MCMC DA methods to estimate parameters for both parameters that were previously chosen by hand in previously published work (Albers et al. (2017)) and parameters that the Houlihan methods determined were low-influence parameters and again calculate the MSE between the data and model state estimates and forecasts. Fourth, we compare the MSE for the variously selected parameter sets.

VI Results

The results come in two stages. First, we present the rank-ordered parameters selected by different methods in order to demonstrate: (i) which parameters the methods selected, (ii) that the methods selected some but not all parameters, (iii) how the parameter selection varied across methods, and (iv) the rank-ordering of parameters by method. Second, we evaluate the methods by using the parameters selected in each method to forecast glucose with the UKF and smooth glucose with MCMC; methods are compared via the MSE between measurements and predictions.

VI.1 Parameter selections by method

Table 2 shows the rank-ordered parameters selected by each parameter selection method. The methods were sensitive to the feature metric; the mean and standard deviation-based methods did not select the same parameters as important.

For a given feature metric, all selection methods identified the same top two parameters — all methods ranked $a_{1}$ and $C_{1}$ as the top influencers of the mean, and ranked $R_{g}$ and $C_{3}$ as the top influencers of the standard deviation. However, the entire influence sets differed substantially. This indicates that influence set structure, as defined (upper quartile of influence), is sensitive to choices of influence functions and influence equivalence definitions.

Interestingly, the equivalence classes of high and low parameter influence are preserved under perturbations to the influence function. Fig. 2 shows how the $l_{1}$, $l_{2}$ and PCA-based methods rank-order parameters according to how they influence the mean. While lasso is expected to preserve the ordering with different $\lambda$ (it fits one-at-a-time), ridge regression also remains robust to variations in the regularization term, $\lambda$, adding parameters one at a time.

Most methods find only $5-6$ influential parameters out of $21$, greatly reducing the dimension of the parameter space. In all cases, the methods gave an entry point for which parameters to begin estimating; the next question, then, is whether using the Houlihan approach helps to reduce forecasting errors and improve convergence of parameter estimates.

Redundancy and influence

Our goal is to select parameters to estimate during forecasting and smoothing tasks. We aim to facilitate this goal by identifying small parameter sets that have significant, minimally redundant influence over important dynamical features. Accomplishing this can minimize problems in identifiability, multiple coexisting invariant sets, etc. Fig. 3 visualizes variation of the feature-metric, mean and standard deviation of the invariant density with parameter variation, as well as how the methods partitioned parameters into a high and low-influence equivalence class. It is clear that some variations in some parameters create large shifts in the mean and variance (e.g. $a_{1}$), whereas the mean and variance features are far less sensitive to other parameters, like $E$ and $t_{d}$.

While the mean and standard deviation are not always influenced by the same parameters, the methods select parameters that have both high influence and relatively orthogonal influence; e.g., in the case of the mean the methods generally select $a_{1}$ and $R_{g}$ first. The low influence parameters, by comparison, are not able to move the mean or standard deviation appreciably and are therefore not able to fully explore the space. Similarly, the low influence parameters are relatively redundant. Following this logic one might predict that estimating alpha and $C_{2}$ would lead to the most accurate model estimates while estimating $E$ and $t_{d}$ would lead to the least accurate model estimates.

Comparison with by-hand selection

In our previous work (Albers et al. (2017)) we selected parameters to estimate by hand based on our desire to estimate certain parameters related to physiologic function, e.g., $t_{p}$, and because of their obvious influence on parameters, e.g., $V_{p}$ as could be deduced from other previous work (Albers et al. (2014)) to influence the mean state. The automated methods selected $V_{p}$ and $t_{p}$ as high influence parameters, but not $E$, a parameter the methods determined was a low-influence parameter.

VI.2 Parameter selection method evaluation

To evaluate the effectiveness of the machine-selected parameters compared to low-influence parameters as characterized by their $\beta_{i}$’s, and the by-hand-selected parameters we used in our previous work, we compare the mean squared error (MSE) between the data and the forecasts for the various parameter combinations as shown in table 3. Fig 4 provides a visual summary of the results in table 3—the plots are calculated directly from table 3—for the MCMC smoothing setting, and demonstrates that all Houlihan-based parameter sets (of any size) noticeably out-performed the by-hand and low-influence parameter sets. Moreover, we see that most Houlihan-based methods achieve similar overall accuracy for parameter sets of cardinality $\leq 3$. In addition, Houlihan-based methods that selected parameter sets with 4 or more parameters achieved the best performance, and there is a general trend of improved fit with more parameters—this contrasts sharply with the by-hand parameter selections, whose performance tapered with more than 3 parameters (probably due to unforeseen issues of identifiability).

In particular, lasso chose parameters with the lowest MSE between forecasts and measurements in $7$ of $8$ cases. In one case, taking the union over methods shared the same MSE with lasso. And, in one case, the lowest MSE was observed with a pair of low-influence parameters. In this case it was the parameter-pair combination, $\alpha$ with $E$, that mattered. This result implies that generally low influence parameters may, for some people, be physiologically important and explore particular pathophysiology necessary to synchronize to the individual. We also know that as the number of parameters increased to $3\geq$, some of the MCMC parameter estimates with the lowest MSE found multiple, competing equilibria, were not unique, and sometimes did not fully converge. For example, Fig. 5 shows parameter estimates of two different parameters—one that converges and one that does not—for two parameter sets for P1 with standard deviation as the feature-metric. When lasso-selected parameters are restricted to two parameters for P1, then both parameters, $R_{g}$ and $C_{3}$ converge producing a MSE of $600$; $C_{3}$ is shown in Fig. 5. In contrast, lasso restricted to the one standard error minimum selects eight parameters, has a lower MSE of $375$ but at least one of the parameters, $t_{p}$, does not converge well as shown in Fig. 5. This means that as we increased the flexibility, we lowered the MSEs but possibly came at the expense of physiology or convergent parameter estimates.

VII Discussion

Summary Our most broad conclusion is that the machine-selected parameters work better than hand-selected parameters and that the Houlihan methods are a scalable method for selecting which parameters of a mechanistic model to estimate using DA methods. This means that stacking machine learning techniques on top of, or together with, DA is a helpful strategy, especially when models are complex and data are sparse, as in our glucose modeling example.

Houlihan methods We intuitively define the Houlihan method(s) as a collection of methods for selecting the most productive model parameters to estimate with machine learning techniques applied to simulated model output under parameter variation subject to a set of features, e.g., the mean of a state.

using machine learning to

Feature metric selection matters: For all methods, the feature metric (mean or standard deviation) was the first-order driver of differences in parameter rank orderings. This choice is highly problem-dependent. In some biomedical applications, sensitivity of the mean to parameter perturbation is not especially important for a good fit; e.g., there are physiologic systems where variation in the mean across people is small, but excursions, peaks, number of peaks, location of peaks, etc., may be a more important types of features to capture.

The cutoff matters: The cutoff for influence has a substantial impact on the ability to estimate parameters. For example, lasso-selected parameters usually minimized MSE, but the induced MSE and MCMC convergence were both sensitive to the influence cutoff. All the methods had this sensitivity, and estimating optimal cutoffs automatically would be beneficial.

The selection method sometimes matters: For the high-ranked parameter choices, the feature-metric was the primary difference between selected parameter sets. However, as the number of parameters included was increased, the methods diverged. We suspect that as the complexity of feature metrics and ranking methods increases, e.g., using nonlinear regressions, there will be more sensitivity of the parameter selections to the methods.

Physiology matters: We know from carefully considering the convergence properties of the MCMC chains that some of the lowest MSEs for the runs with three or more parameters didn’t converge well. Meaning, as we increased the flexibility, we lowered the MSE but possibly at the expense of physiologic fidelity or convergent parameter estimates. For pure forecasting applications this may or may not matter, but when we want the parameters to be meaningful, we need the parameter estimates to converge, not necessarily to a unique set of parameters, but to distinctly different parameter estimates that can be treated as hypotheses. Another problem that can arise because of physiology is that different people with different physiology can be sensitive to different parameters. For example, the physiological feature that is important to personalize the model for a particular person may not be related to the properties captured by the feature metric, e.g., the mean, and in this circumstance parameters identified as low influence relative to the feature-metric will not be estimated. A potential example of this is P1, for whom estimating $\alpha$ and $E$ achieved the lowest MSE despite $E$ and $\alpha$ being low-influence parameters relative to both the mean and standard deviation.

Effective parameter space exploration: Abstractly, a mechanistic model is a parameterized family of functions whose parameters, depending on the model, have varying degrees of independence. From this perspective, the goal of the Houlihan methodology is to find a way to explore the maximal amount of the parameter space while minimizing the redundancy between parameters. The feature-metrics and the influence functions define which subsets of the parameter space are most useful to open for exploration, which in turn defines which dynamics can be explored. For example, focusing on variations of the mean may close off other dynamical features such as amplitude variations or any feature that is not uniquely defined by the variation of the mean. We do not yet have a good method for understanding how a feature-metric may influence other, potentially valuable explorations. We acknowledge that understanding and quantifying how limited feature-metrics influence the effective parameter space of a model is an important, unexplored problem.

Computational complexity: We consider only the case here where we vary any one single parameter while leaving all other parameters fixed at their nominal values; this means that the dimension of the input for regression used to select the most useful parameters scales linearly in the number of parameters. If we were to co-vary parameters, meaning if varied all parameters at once, depending on how one choose to partition the parameter space, the computational complexity would explode. In this way, the framework we present here does not solve the computational complexity problem of exploring parameter space. Instead, the results in this paper show that even by only considering feature-metric variation along one-dimensional subspaces of the full parameter space we can gain substantial insight into which parameters have the most impact on the features we are interested in approximating. Moreover, we can also see the limitations of this approach — we do observe synergy between parameters where combinations of some low-influence parameters for some people can end up having a high influence on the model fit.

Obvious extensions: In this paper, we stack machine learning on top of DA, which has many potential extensions. Feature-metrics could be generalized to be multi-dimensional both over states and over types of feature-metrics. Feature selection methods could be developed or employed to select feature metrics. The estimates of influence could be calculated to include jointly varying parameters—this would be computationally expensive and would require computational innovation in high-dimensional settings, cf the computational complexity discussion above. Moreover, this problem is not necessarily a simple extension because the parameter spaces of mechanistic models are not likely to form a basis for the model space, in contrast to the parameters of the space of polynomials which do form a basis. Of course this lack of a basis structure is part of the problem—parameters of mechanistic models and likely the physiology they represent are redundant, likely for biological reasons such as robustness. We use linear regression and PCA-based machine learning methods; it is likely that more sophisticated machine learning methods e.g., full elastic nets, support vector machines, deep learning, sparse machine learning (compressed sensing), Bayesian methods, model averaging and ensemble learning could all be used and would likely improve the parameter selections. Similarly, further stack of machine learning techniques on top of the Houlihan methods would likely be productive. For example, greedy, Gibbs-sampling-like rotation between sets of parameters that are identifiable and explore different subsets of the parameter space could minimize both model errors and identifiability issues. And finally, feature-metrics could be made substantially more sophisticated, insightful and tailored to circumstance or physiologic knowledge, such as preserving power in certain frequency bands. More sophisticated feature metrics could also be used to gain insight into potentially meaningful constraints on parameters for use in operational DA.

VIII Conclusion

We devised a methodology for rank-ordering parameters of a mechanistic model and using this rank-ordering to select an effective subset of parameters to estimate when projecting biomedical data onto the model via data assimilation. This methodology specifically targets parameter sets that avoid issues of model identifiability and parameter-estimation convergence problems, improving forecasting and phenotyping performance of data assimilation methods that use mechanistic biological models. Using machine learning to select parameters to estimate worked: the machine-chosen parameters reduced the mean-squared error between estimates and forecasts and data in nearly all cases by factors as large as three. These results imply that combining mechanistic and non-mechanistic machine learning could be a particularly productive direction of future research and could greatly aid in our ability to use computational machinery to both help deepen our physiologic understanding and help clinicians achieve more positive outcomes in clinical settings.

Appendix A Ultradian model

The model is comprised of a set of six ordinary differential equations; the model is non-autonomous because it has an external, time-dependent driver, consumed nutrition. The six dimensional state space made up of three physiologic variables and a three stage filter. The physiologic state variables are the glucose concentration $G$, the plasma insulin concentration $I_{p}$, and the interstitial insulin concentration $I_{i}$. The three stage filter $(h_{1},h_{2},h_{3})$ which reflects the response of the plasma insulin to glucose levels Sturis et al. (1991). The model was designed to capture ultradian oscillations missing in previous models. The ordinary differential equations that define the model are Keener and Sneyd (2008):

[TABLE]

The state variables include physiologic processes that have been parameterized, including: $f_{1}(G)$ represents the rate of insulin production; $f_{2}(G)$ represents insulin-independent glucose utilization; $f_{3}(I_{i})G$ represents insulin-dependent glucose utilization; $f_{4}(h_{3})$ represents delayed insulin-dependent glucose utilization. These functions are defined by:

[TABLE]

The nutritional driver of the model $I_{G}(t)$ is defined over $N$ discrete nutrition events Albers et al. (2014), where $k$ is the decay constant and event $j$ occurs at time $t_{j}$ with carbohydrate quantity $m_{j}$

[TABLE]