Improved Estimation of Human Lipoprotein Kinetics with Mixed Effects Models

Context Mathematical models may help the analysis of biological systems by providing estimates of otherwise un-measurable quantities such as concentrations and fluxes. The variability in such systems makes it difficult to translate individual characteristics to group behavior. Mixed effects models offer a tool to simultaneously assess individual and population behavior from experimental data. Lipoproteins and plasma lipids are key mediators for cardiovascular disease in metabolic disorders such as diabetes mellitus type 2. By the use of mathematical models and tracer experiments fluxes and production rates of lipoproteins may be estimated. Results We developed a mixed effects model to study lipoprotein kinetics in a data set of 15 healthy individuals and 15 patients with type 2 diabetes. We compare the traditional and the mixed effects approach in terms of group estimates at various sample and data set sizes. Conclusion We conclude that the mixed effects approach provided better estimates using the full data set as well as with both sparse and truncated data sets. Sample size estimates showed that to compare lipoprotein secretion the mixed effects approach needed almost half the sample size as the traditional method.


Introduction
Complex traits (e.g. metabolic pathways and fluxes) may often only be quantified indirectly using mathematical or statistical models. A common approach is to make a parameterized mathematical model describing the system on an individual level, and then use statistical methods to make inference from the estimated parameters on the population level. Thus, any uncertainty in the individual estimate is affecting the population estimate. The analysis is further complicated by the great variability seen in biological systems.
There are at least two conceptually different approaches to treat variability in mathematical models of biological systems. In the standard two stage (STS) approach [1,2,3], each individual is treated separately: measured data for an individual are used to estimate parameters that are assumed to be truly specific for that individual (the first stage). The obtained individual parameters are then treated with traditional statistical methods (the second stage), such as correlations with other covariates or comparison of parameters between groups.
Contrary to this, in a mixed effects model one assumes that each parameter in the model is either essentially the same for all individuals in the group, a function of some directly quantified individual-specific traits (covariates), or random with a probability distribution that is determined by covariates and parameters that are common for the whole group. One can then use all measured data to make parameter estimations for the groups.
An advantage of the STS approach is that it does not assume any a-priori information on the distribution of the parameters and on treatment effects (i.e. which model parameters are affected by treatment) or group effects (i.e. which model parameters are different between groups). A disadvantage is that, to adequately estimate parameters for each individual, large data sets are needed and to estimate treatment or group differences, large sets of individuals are needed. In case of complex and costly experiments, a modeling approach that utilizes smaller data sets is desirable to reach sufficient study power while minimizing the cost of the experiment.
The mixed effects, including nonlinear mixed effects (NLME), approach is widely described and used in the literature [3,4,5,6,7,8]. It has been shown that better estimates of population characteristics are obtained for mixed effects models compared to STS models in many systems. Sheiner and Beal showed already in 1980 that the STS method can produce biased and imprecise estimates of the variability between individuals [4]. The STS and mixed effects approaches have been compared for many systems [2,9,10,11].
Abnormal concentrations of blood lipids and lipoproteins are key risk factors for cardiovascular disease and are often observed in metabolic diseases such as diabetes mellitus type 2 (DM2) [12]. In DM2, the lipid abnormalities primarily include elevated plasma triglycerides and decreased high density lipoprotein (HDL) cholesterol [12] and it is believed that these changes are the cause of increased risk in cardiovascular disease in DM2 and related diseases [12,13]. Many therapies are targeted to normalize blood lipid levels, and potentially decrease the risk for developing cardiovascular disease.
The non-soluble lipids, such as cholesterol esters and triglycerides, are transported in the core of specific vesicles called lipoproteins. On the surface of these particles different proteins are attached. During fasting, these particles are synthesized in the liver and transport the lipids to peripheral tissues. As the lipids are delivered the particles become smaller and denser as the relative amount of protein (which is denser than the lipids) increases. The catabolic processes are combinations of lipases (lipoprotein lipase and hepatic lipase), which remove triglycerides from the core, and removal of whole particles due to specific uptake mechanisms. Thus, elevation of plasma triglyceride in disease states may be attributed to increased secretion, defined by the secretion rate (SR), or increased clearance.
To understand disease development and the role between the underlying metabolic defects and blood lipids levels the dynamics of lipid and lipoproteins are studied using stable isotope infusions and mathematical modeling [14]. Using such models, quantification of both the catabolic processes and the productions is possible [14].
The models allow for the quantification of the secretion rate (SR) of lipoproteins (either protein or triglyceride content) (mg/kg bw/day), the fractional clearance rate (FCR) representing the total turnover of the lipoproteins (pools/day), the fractional transfer rate (FTR) representing the fractional conversion of large particles to smaller particles due to the action of lipases (pools/day) and the fractional direct catabolic rate (FDCR) representing the direct catabolism/ removal of large particles (pools/day).
Data acquisition is labor intensive, in particular when the system is described by several sub-fractions. Improved methods for data analysis may therefore allow for a reduction of either the number of study subjects or the amount of data needed for each subject while maintaining study power.
In this paper we compare the traditional STS approach and a new mixed effects model approach. We use experimental data from fifteen healthy control subjects and fifteen DM2 subjects. DM2 is typically associated with elevated plasma triglyceride, which is known to be caused by increased production of lipoprotein particles from the liver [24,25,26,27]. The aim of the study was to investigate if the mixed effects approach improved parameter estimation compared to the STS approach for small sample sizes and for sparse data. We base the analysis on the established model of lipoprotein kinetics [15,19,23].
We show that, when shrinking the sample size, group properties are preserved using the mixed effects approach while the STS model fails to replicate significant differences between the groups. Moreover, individual estimates are better when using the mixed effects approach while using less data for each individual.

The Model
In this paper we develop a mixed effects population model of the very low density lipoprotein subclass (VLDL 1 and VLDL 2 ) kinetics. The model is based on our previous work on individual compartmental models [15], see Fig 1. The details of the model are described in the Materials and Methods section. We distinguish between model parameters (e.g. transfer rates in the model) and key parameters (e.g. biologically relevant parameters, often composed of several model parameters).
The data set consists of 15 healthy control subjects and 15 patients with DM2, as characterized in S1 Table. Estimated Parameters Using All Subjects Three data sets are used. The full data set (all data) and two reduced data sets: the first 4 hours (truncated data) and every second data point of the VLDL tracer data (reduced data).
The estimated key parameters using both the STS and the mixed effects approaches and the three data sets are summarized in Table 1.
All data. Using the full dataset the mixed effects approach typically produced smaller variations within groups compared with the STS approach ( Table 1). Both approaches confirmed previous findings of greater secretion of VLDL 1 particles (VLDL 1 SR), resulting in a larger VLDL 1 pool, a greater flux from VLDL 1 to VLDL 2 (expressed as indirect secretion of VLDL 2 ), and a larger VLDL 2 pool in the DM2 patients compared with the control subjects ( Table 1).
With both approaches the total VLDL 1 FCR was lower in DM2 patients compared with control subjects. However, the difference was not significant. In contrast to the STS approach, the mixed effects model showed significantly higher VLDL 1 FTR and lower VLDL 1 FDCR in DM2 patients compared with the control subjects ( Table 1).
Both approaches produced excellent fit to the data as seen in the residual plots (Fig 2A-2F). The STS approach generally produced slightly better fit to the data when quantified using the root mean square error (RMSE) of the residuals for all individuals.
Truncated data. The mixed effects model replicated all significant differences between the DM2 patients and the control subjects compared with using the full data set ( Table 1). Using the truncated data, the STS model still showed higher VLDL 1 SR in DM2 patients although with lower significance, but did not detect a higher VLDL 1 to VLDL 2 transfer for DM2 patients compared with control subjects (Table 1). In general the mixed effects model generated data with smaller CV ( Table 1), indicating a better coherence between the results using the mixed effects approach.
Comparing the mixed effects model using all data with the truncated data showed significantly lower VLDL 1 FTR in both control subjects (p<0.05) and in DM2 patients (p<0.05) using truncated data. In addition, the VLDL 1 to VLDL 2 transfer estimates were lower with truncated data in DM2 patients (p<0.01) ( Table 1).
Using truncated data, the STS estimates of VLDL 1 SR were significantly higher (p<0.05) and the VLDL2 FCR were significantly lower (p<0.05) in control subjects, compared with using all data ( Table 1).
The mixed effects model produced better fit to the data even when extrapolating the model curves to the excluded data (S1A-S1F Fig). This indicates that, in at least some individuals, the information about the trajectories is captured by the information in the first 4 hour data and that the mixed effects approach shares this information to provide better fit for all individuals.
Reduced data. After reduction of data by removing every second data point both methods still reproduced the major findings of study, however with higher significance for the mixed effects model compared with the STS model ( Table 1). Using the reduced data sets no estimate Model constraints, k 2,1 = k 1,2 , k 3,4 = 0.1 k 4,3 , k 9,8 = k 7,5 and k 0,9 = k 0,7 were used. Compartment 1 represents the free leucine in the plasma, compartments 3 and 4 represents leucine recycling in non-hepatic tissue and compartment 2 is the intrahepatic leucine that feeds the apoB synthesis compartment represented as a delay (D 1 -D 7 ). Newly synthesized apoB enters the plasma as VLDL 1 (large particles, compartment 5) or VLDL 2 (small particles, compartment 8). VLDL 2 may also be produced through conversion of VLDL 1 via compartment 6. Particles may leave the system through compartments 6, 7, 9 or 10. Healthy control subjects (Control) and type 2 diabetic patients (DM2) are compared using the two approaches and using three sets of data: The complete data set (All), the truncated data set (up to 4 hours, Truncated), and the reduced data set (every second data point removed, Reduced). The estimated parameters are presented as geometric mean (coefficient of variation were significantly different using neither the mixed effects model nor the STS model ( Table 1).
The fit to the curves were similar for the two methods (S2A- S2F Fig).

Reduced Study Sample Size
As the mixed effects model produced smaller variations we next tested how this affected the sample size needed to detect a difference between the groups. The (log-transformed) means and variances estimated using the full data set were used to calculate the sample size needed to detect a difference between the means with a power of 80%. To further test this we also used a resampling approach to resample the study population into smaller group sizes of 4, 6, 8, 10, and 12 subjects from each group. Both the mixed effects model and the STS model were applied to the re-sampled subgroups and the p-values for the comparison between the control subjects and the DM2 patients were recorded for 30 different re-sampled groups for each group size. The estimated sample sizes were generally smaller for the mixed effects approach compared with the STS approach, both using the sample size calculation and the resampling approach ( Table 2).
Using the mixed effects approach the increased VLDL 1 secretion in DM2 patients were correctly identified in all 30 resampled iterations using the mixed effects model for a study group size of 12 and 10. The difference was correctly identified in more than 80% of the iterations down to a group size of 6 subjects. In contrast the STS approach reached an 80% correct identification rate at a group size of 10. This indicates that to have an estimated study power of 80% in detection of the difference in VLDL 1 SR it is sufficient to have only 6 subjects from each group using the mixed effects approach but 10 using the STS approach ( Table 2). Similar Residual plots of enrichment data using all data. The residuals (model fit minus measurement data) for the three enrichment data sets (A and B, plasma leucine; C and D, VLDL 1 ; E and F VLDL 2 ) were plotted for the two methods (STS and NLME) and the two groups (A, C and E, Control; B, D and F, type 2 diabetes mellitus (DM2)). Both methods produced good fits to the data. The NLME approach used a sequential procedure; the plasma leucine were fitted first and the VLDL 1 and VLDL 2 curves were fitted using the leucine results. This may explain the worse fit for the plasma leucine in the STS approach. Lines, mean of mixed effects approach (red) and STS approach (black); Areas, mean ± SD for mixed effects approach (red) and STS approach (black). results were achieved for the VLDL 1 to VLDL 2 transfer and VLDL 1 pool ( Table 2). The sample size was also estimated directly using power calculations ( Table 2) with comparable results.

Computational Power
A major concern with mixed effects modeling is the extended computational time needed to estimate the likelihood function. For each calculation of the likelihood function, the individual parameters need to be estimated. Some steps can therefore be parallelized. Using our implementation in Matlab, a mixed effects model with 8 DM2 and 8 control subjects could still not be run on a traditional PC. On an 8 core cluster the computational time was between 24h and 42h in the 30 repeated sample size estimates, compared to less than 10 minutes for the STS model.

Discussion
The main result in this study is the novel application of mixed effects models to experimental data from human lipoprotein kinetic studies, comparing the lipid metabolism in healthy control subjects and DM2 patients. Compared to the traditional STS approach our mixed effects approach provided stronger results using the full data sets. The results also provide a basis for the use of smaller sample sizes when data are analyzed with the mixed effects approach. The mixed effects approach produced estimates with less variation which translated into the possibility to use fewer subjects with sustained study power as shown by direct calculations as well as simulations using resampled populations. Furthermore we show that significant results may be obtained using smaller experimental datasets, as given either by either shorter sampling times or less frequent sampling.
Studies of lipoprotein kinetics are labor intensive and require generation of large sets of data for each individual in order to obtain good parameter estimates. Often subjects are studied for several different tracers simultaneously which adds to the complexity. It is therefore of great interest to derive novel techniques that enables better estimates of kinetic parameters.
Many metabolic diseases are linked to disturbances in lipid metabolism. The traditional analysis of plasma concentrations cannot answer which underlying mechanism that is responsible for the observed disturbance [14]. In contrast lipoprotein kinetic studies allows for decomposition of the concentration into several components including secretion and clearance. In the case of type 2 diabetes, there is a large body of evidence that the major defect resulting in elevated plasma triglycerides is an increased secretion of triglyceride rich lipoproteins ( [24,25,26,27,28]). The rate of clearance, however, is usually lower in type 2 diabetic patients but the results are inconsistent with studies showing either non-significant (although with a trend) [24,26,28] or significant [25] differences. Our results confirm these results and may further contribute the analysis by providing better estimates of variability of the clearance parameters.
We tested the effectiveness of the mixed effects model using two different approaches to reduce data. First we reduced the data by either truncating the time series or by removing every second data point. Truncating the data results in loss of information about the tail shape of the curves, which usually lead to poorer estimates. The mixed effects model provided better predictions of the excluded data on an individual basis. Both the mixed effects approach and the STS approach showed significant differences in the estimated means compared to the full data sets (Table 1) for some parameters indicating that the absolute quantifications were poorer using the truncated data. However both approaches consistently recreated the significant differences between the DM2 patients and the control subjects as observed with all data, with p-values being smaller using the mixed effects approach. When data were instead reduced by 50% by removing every second data point, the absolute quantifications were not different compared with those obtained using all data. The mixed effects approach provided smaller p-values for the significant differences between the two groups.
The second approach to compare the two methods was to perform a power calculation and to estimate the sample size needed to detect differences with a power of 80%. The mixed effects approach needed smaller sample size compared with the STS approach, both when the sample size were calculated directly or when estimated using resampling.
An interesting observation was that when truncating the data at 4 hours, the mixed effects approach predicted the measurements between 4 and 8 hour much better than the STS approach. Typically subjects with normal plasma triglycerides does not present a "tail" of enrichment in VLDL 1 but subjects high plasma triglycerides do [19]. To account for these slowly turning over lipoproteins the compartments 7 and 10 are added to the model. It appears as the mixed effects model may capture this behavior in some individuals already before 4 hours and share this among all individuals. This is consistent with the fact that the STS approach estimates both the VLDL 1 SR and VLDL 1 FCR higher when truncated data is used, suggesting an inability to capture a slowly turning over pool.
We tested the mixed effects approach on an observational study. For interventional studies the statistical models for the parameters in the population model needs to be defined according to the study design.
Our mixed effects model was implemented in MATLAB. The computational time for the mixed effects approach is several magnitudes slower compared to the STS approach. This may be attributed to the fact that in each calculation of the maximum likelihood function individual parameters needs to be estimated for each individual. It is therefore not today feasible to run this on a traditional workstation PC. However as the estimation of the individual parameters are independent, the calculations may be done in parallel threads and may therefore be optimized for cluster or multiple processor workstations. However, the additional costs in computation are likely minor compared to the possible savings in study population size, study power and/or laboratory work.
In conclusion, we shown that mixed effects modeling is an effective tool to estimate kinetic parameters in human lipoprotein kinetic studies. Either sample size or the amount of data may be greatly reduced while retaining study power.

Study Population
Data used in this study are previously published [15,24,29]. In this study, we randomly selected 15 healthy control subjects and 15 patients with type 2 diabetes with similar age and BMI (S1 Table) from our study cohorts. The diagnosis of DM2 was based on glucose tolerance test results in the diabetic range according to World Health Organization (1999) criteria or on the use of oral diabetes medication. The original study designs were approved by the Helsinki University Central Hospital ethics committee, and each subject gave written, informed consent as described in the original publications [15,24,29]. All samples were collected in accordance with the Helsinki Declaration.

Data
The details of the experimental protocol, sample analysis etcetera are explained elsewhere [15,24]. In short, N C = 15 healthy control subjects and N D = 15 patients with type 2 diabetes with similar age and BMI were studied (S1 Table). Following an overnight fast a bolus infusion of [ 2 H 3 ]-leucine was infused. Blood was collected and the amounts of labeled and unlabeled leucine were measured as free in plasma and in the major protein component, apolipoprotein B-100 (apoB), in VLDL 1 and VLDL 2 by GC/MS. The tracer-to-tracee ratio (i.e. labeled-to-unlabeled leucine) was calculated for free leucine in plasma (sampled at 2,4,6,8,10,12,15,20,30, and 45 minutes and 1, 2, 3, 4, 6, and 8 hours after the injection) and as leucine in VLDL 1 and VLDL 2 (sampled at 0.5, 0.75, 1, 1.25, 1.5, 2, 2.5, 3, 4, 5, 6, 7, and 8 hours). Pool sizes of VLDL 1 and VLDL 2 are also measured at 0, 4, and 8 hours after the injection, and recalculated to leucine levels in apoB in VLDL 1 and VLDL 2 as explained in [15].
Three data sets are used. The full data set and two reduced data sets: the first 4 hours (truncated data) and every second data point (reduced data). Furthermore reduced group sizes (i.e. 12, 10, 8, 6 and 4 subjects) are tested.
The full data set is available as an online supplement in file S1 File Data.csv with data format given in S2 File Data Format.docx.

The Individual Model
The individual model is based on the compartmental model formulation developed by Demant and Packard [19,23]. The model has been used to determine the kinetics of apoB-100 in the VLDL 1 and VLDL 2 subclasses in individuals with a wide range of conditions (e.g. [15,18,19,23,24,29]). The model has a 4 compartment catenary chain that models the free leucine in plasma and the intra-hepatic leucine pool before entering the synthesis module consisting of a seven compartment delay. Each sub-system of VLDL 1 and VLDL 2 have three compartments (5, 6 and 7 for VLDL 1 and 8, 9 and 10 for VLDL 2 , Fig 1), two compartments are representing a delipidation chain (5 and 6, and 8 and 10 respectively) resulting in transfer to denser particles and one compartment (7 and 9 respectively) represents particles being cleared from the circulation by other mechanisms as summarized in Fig 1. We denote the mass of the tracer (labeled leucine) for an individual i, in a compartment j, at time t by q i j ðtÞ and the constant mass of tracee (unlabeled leucine) by Q i j . Thus, free leucine is denoted by q i 1 ðtÞ and Q i 1 , VLDL 1 leucine in apoB is denoted by q i V 1 ðtÞ ¼ q i 5 ðtÞ þ q i 6 ðtÞ þ q i 7 ðtÞ and Q i V 1 ðtÞ ¼ Q i 5 ðtÞ þ Q i 6 ðtÞ þ Q i 7 ðtÞ, and VLDL 2 leucine in apoB is denoted by q i V 2 ðtÞ ¼ q i 8 ðtÞ þ q i 9 ðtÞ þ q i 10 ðtÞ and Q i V 2 ðtÞ ¼ Q i 8 ðtÞ þ Q i 9 ðtÞ þ Q i 10 ðtÞ. To have comparable results we use the same model of error and method of parameter estimation for both the mixed effects model (below) and the individual model. The output variables are assumed to be log-normally distributed as this result in relative errors. Therefore the data are log transformed before analysis and the output variables for each sample time point t k are defined by the logarithm of the corresponding variables with the added normally distributed error term as: We assume that the errors are of the same magnitude for VLDL 1 and VLDL 2 , i.e. s i Since only a small amount of tracer material have been incorporated into the lipoproteins at 0.5 hours (the time point for the first measurement of tracer kinetics of VLDL 1 and VLDL 2 ), the uncertainty is larger for measurements taken after only half an hour. Therefore, we assume that for a parameter ϕ i ! 1. In principle, separate values of s i (t k ) could be used for each k. The state variables are given by solving the following system where q i and Q i are state vectors for the tracer-and tracee system respectively, U i is the inflow vector of tracee material into the system (which is only through compartment 1 for this model), and K i is the compartmental matrix of the model (see, e.g., [30,31]). The steady-state assumption implies that the left hand side in Eq (1)

The Population Model
For the STS approach each parameter in the model is estimated individually (except ϕ i , which is set to 10). Key variables are then derived from the model parameters as is explained below. For the mixed effects approach the model parameters for all individuals are estimated simultaneously by assuming that the individual model parameters are members of probability distributions determined by the whole population. For this mixed effects model we separate the leucine model to the rest and first estimate the parameters in the leucine model as carefully explained in [32]. For the VLDL kinetics part of the model the individual model parameters are where each Z i x $ Nð0; o 2 x Þ, except Z i 0;6 which is allowed to vary uniformly between 0 and 5, since we could not obtain good fits to the data for all individuals with a normally distributed Z i 0;6 . The indicator functions 1 i C and 1 i D are used to indicate if individual i belongs to the control-or DM2 group. Observe that we assume that all ω x are group independent. The reason for using fixed population parameters for k 10,8 is that this parameter is hard to estimate when only VLDL 1 and VLDL 2 data are available. Better fits to the data were not obtained for individual-specific k 10,8 (for the STS models the parameters for each subject were re-estimated with k 10,8 fixed to the group values obtained from the mixed effects model). In fact, ω 10,8 tended to zero when it was included in the mixed effects model. The variability parameters in the residuals are assumed to be identical for all individuals, Eq (2). A final remark is that we assume that the covariance matrix O of Z i ¼ ½Z i p 1 ; . . . ; Z i 0;10 is diagonal to keep a reasonable dimension of the parameter space.

Geometric Mean and Coefficient of Variation of a Log-Normal Distribution
Here we explain how the population values of the parameters in the tables and figures are computed. Let θ i , i = 1,. . .,N be the estimated values for a parameter for the subjects in one of the groups. Almost all key parameters are skewed over the population (for both the STS-and mixed effects models). Therefore, we assume that they are log-normally distributed and compute the sample geometric mean which is the maximum likelihood estimator for the median of a log-normal distribution [33]. The coefficient of variation (in percent) is given by 100 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi expðô 2 Þ À 1 p ; whereô 2 is the unbiased estimate of the sample variance of the logarithm of the θ i 's. When statistical tests are performed we compare logarithms of the parameters and use the conventional t-test theory. This can be done since the logarithm of the geometric mean in Eq (3) is

Parameter Estimation
For parameter estimation of the mixed effects model, we apply the maximum likelihood framework. Here we quickly present how the likelihood function of the population parameters is approximated. For a more thorough analysis and further references we refer to [32,34]. The population likelihood function is a marginal likelihood function, where the vector Z i 2 < r of all individual-specific parameters (r = 9 for the VLDL model) are integrated out. Let θ be the vector of population parameters and let Y = [Y 1 ,. . .,Y N ] the complete set of measurement data for the N individuals. Assuming that the subjects are picked independently in the population the population likelihood function takes the form where l i = l i (η i ) = l i (η i ;Y i ,θ) is the log-likelihood function for subject i. Let the residuals at time t k be denoted ε i k and let R i k be its covariance matrix, hence R i k ¼ diagðsðt k Þ 2 ; sðt k Þ 2 ; S 2 ; S 2 Þ for the VLDL model. The individual log-likelihood function is given by where m is the number of outputs (thus, m = 1 for the free leucine model and m = 4 for the VLDL model). The integrand, exp(l i ), in Eq (4) is hard or impossible to integrate exactly and must therefore be approximated. The Laplace approximation method is based on a second order truncation of the Taylor expansion of l i around a stationary pointẐ i (see, e.g., [32,34] for the derivations). Eq (4)  where H Z ðl i ðẐ i ÞÞ is the Hessian matrix of l i with respect to η i . In the first order methods second order derivatives are disregarded. In the model considered in this paper the covariance matrix R i k does not depend on η i and the Hessian becomes In the first order conditional estimation (FOCE) method the mode of the individual loglikelihood function is used forẐ i , hencê in contrast to the first order (FO) method that usesẐ i ¼ 0.
For the STS model, the individual log-likelihood function in Eq (5) with the last three terms removed is maximized for each subject. In that case the uncertainty parameters s i and S i are estimated separately for each individual and ϕ i is set to 10. It should be noted that k 10,8 is individual-specific for both models. To find the minima of the negative log-likelihood functions the Broyden-Fletcher-Goldfarb-Shanno (BFGS) optimization method is used (see [35]).

Parameters
We distinguish between model parameters (e.g. transfer rates in the model) and key parameters (e.g. biological relevant parameters, often composed of several model parameters). The calculation of the key parameters VLDL 1 secretion rate (SR), VLDL 2 direct secretion rate (DSR), VLDL 2 indirect secretion rate (IndSR) (i.e. the transfer from VLDL 1 to VLDL 2 ), VLDL 1 and VLDL 2 fractional clearance rate (FCR), VLDL 1 fractional direct catabolic rate (FDCR), VLDL 1 (4b) fractional transfer rate (FTR) is defined as The VLDL 1 and VLDL 2 pool sizes are defined as Q i V 1 and Q i V 2 . Although the mixed effects approach produces direct estimates of the distribution of model parameters, the key parameters are evaluated for each individual separately.

Implementation
Both the mixed effects model and the STS model were implemented in MATLAB (The Math-Works Inc., Natick, MA).

Statistical Analysis
Log-normally distributed variables are presented as geometric mean and coefficient of variation and were log-transformed before being compared using t-tests assuming equal variance. Mann Whitney test was used if the transformed data were not normally distributed. A p-value below 0.05 was considered significant. For power analyses a desired study power of 80% was used to calculate sample size.
Supporting Information S1 Fig. Residual plot of enrichment data using reduced data. The residuals (model fit minus measurement data) for the three enrichment data sets (plasma leucine, VLDL 1 and VLDL 2 ) were plotted for the two methods (STS and NLME) and the two groups (Control and type 2 diabetes mellitus (DM2)). Both methods produced good fits to the data. Lines, mean of mixed effects approach (red) and STS approach (black); Areas, mean ± SD for mixed effects approach (red) and STS approach (black). (DOCX) S2 Fig. Residual plot of enrichment data using truncated data. The residuals (model fit minus measurement data) for the three enrichment data sets (plasma leucine, VLDL 1 and VLDL 2 ) were plotted for the two methods (STS and NLME) and the two groups (Control and type 2 diabetes mellitus (DM2)). The NLME produced good fit to the data, even when extrapolating the curves between 4 and 8 hours. The STS approach in the other hand fails to produce a good fit for the extrapolated data. Lines, mean of mixed effects approach (red) and STS approach (black); Areas, mean ± SD for mixed effects approach (red) and STS approach (black). (DOCX) S1 File. Data for 15 healthy individuals and 15 DM2 subjects. (CSV) S2 File. File format description for S1 File Data.csv file. (DOCX) S1 Table. Basic subject characteristics. Data from 15 healthy control subjects and 15 type 2 diabetes (DM2) patients were used in the study. As expected DM2 patients had higher plasma glucose and insulin concentrations and showed a typical dyslipidemia with elevated plasma triglycerides (TG), low high-density lipoprotein (HDL) cholesterol and smaller low-dense lipoprotein (LDL). Data is mean ± SD unless otherwise stated. a median (IQR); Ã , p<0.05 vs Control; ÃÃ , p<0.01 vs Control; ÃÃÃ , p<0.001 vs Control. (DOC)