Bariatric surgery improves postprandial VLDL kinetics and restores insulin-mediated regulation of hepatic VLDL production

Dyslipidemia in obesity results from excessive production and impaired clearance of triglyceride-rich (TG-rich) lipoproteins, which are particularly pronounced in the postprandial state. Here, we investigated the impact of Roux-en-Y gastric bypass (RYGB) surgery on postprandial VLDL1 and VLDL2 apoB and TG kinetics and their relationship with insulin-responsiveness indices. Morbidly obese patients without diabetes who were scheduled for RYGB surgery (n = 24) underwent a lipoprotein kinetics study during a mixed-meal test and a hyperinsulinemic-euglycemic clamp study before the surgery and 1 year later. A physiologically based computational model was developed to investigate the impact of RYGB surgery and plasma insulin on postprandial VLDL kinetics. After the surgery, VLDL1 apoB and TG production rates were significantly decreased, whereas VLDL2 apoB and TG production rates remained unchanged. The TG catabolic rate was increased in both VLDL1 and VLDL2 fractions, but only the VLDL2 apoB catabolic rate tended to increase. Furthermore, postsurgery VLDL1 apoB and TG production rates, but not those of VLDL2, were positively correlated with insulin resistance. Insulin-mediated stimulation of peripheral lipoprotein lipolysis was also improved after the surgery. In summary, RYGB resulted in reduced hepatic VLDL1 production that correlated with reduced insulin resistance, elevated VLDL2 clearance, and improved insulin sensitivity in lipoprotein lipolysis pathways.

In the multicompartmental model, the transfer rate from compartment i to j is given by non-negative coefficient k(j,i).If there is no transfer from compartment i to compartment j, then k(j,i) is set to 0. Consequently, the rate of change of material in compartment Qi is given by; where, the first sum represents the total flux into the compartment Qi from other compartments and second sum represents outflux from the compartment Qi.This system can be represented by the following system of differential equations; where, U is the input, Q is the vector of model compartments or the state variables, and K is the matrix that holds transfer rate coefficients.The elements on the main diagonal of K are the negative sums of the elements in the corresponding column.Being defined this way, the diagonal element gives the negative of the total outflux rate from the corresponding compartment.

Modeling tracer injection, liver and plasma modules
Tracer injection, liver and plasma modules are represented by a 28 compartments system (Fig. S2), which is adapted from a previous model (Adiels et al., 2005) and adjusted to account for the postprandial regulations.Compartments 1 and 8 represents free plasma leucine and glycerol, respectively, where the labeled metabolites are injected.Plasma glycerol is converted to TG using the molar weight ratio of 885/92 g/g (Adiels et al., 2005) and the fractional leucine content of apoB is taken as 0.01212g/g (Demant et al., 1998).Leucine and glycerol are taken up by the liver and secreted as VLDL apoB and TG, respectively.The hepatic processing in the leucine and glycerol pathways are represented by 7 and 5 consecutive transient delay compartments, respectively (compartments 12,13,18,19,20 for glycerol and compartments 21-27 for leucine).The transfer coefficients between leucine and glycerol delay compartments are set to 7/ ) and 5/ * , respectively.After the delay, produced apoB is secreted into the plasma either as VLDL1 or VLDL2 (compartments 4 and 6, respectively).The fraction of apoB secreted into the compartment 4 is d4_27, and thus the fraction that is secreted into the compartment 6 is (1-d4_27).
In the model, each plasma apoB compartment is coupled to a TG compartment.If the i th compartment represents an apoB pool, then the i+10 th compartment is chosen to represent the associated plasma TG compartment for ease of illustration.The fractional transfer coefficients between plasma VLDL apoB compartments are defined in the classical way.This is, coefficient for apoB transfer from compartment i to j is defined as k(j,i).The coefficient for the corresponding TG transfer from i+10 th compartment to the compartment j+10 is k(j+10, i+10), which is proportional to k(j,i).During this transfer, due to the TG hydrolysis, only a certain fraction of the TG that leaves compartment i reaches to the compartment j.If this fraction is given by fj,i, then k(j+10,i+10)=k (j,i).fj,i.The fraction of the TG that is removed from the circulation during this transfer is (1-fj,i).
TG flows into the compartment 16 is either from compartment 20 or from compartment 14.The particles transferred from compartments 20 and 14 into the compartment 16 should have the same apoB/TG ratio.If the flux from compartment i to j is defined by F(j,i)=k(j,i)Q(i), this constraint can be imposed by setting F(16,14)/F(6,4)=F (16,20)/F(6,27) at the steady state.This leads to the following distribution term; Hence, the fraction of the TG flux from compartment 20 to the compartments 14 and 16 are given by d14_20 and (1-d14_20), respectively.The compartment 28 is the 'exit' compartment, which has no physiological meaning.
In order to calculate the enrichments, the same liver, tracer injection and plasma modules are used for tracer and tracee subsystems.This way, each tracee compartment coupled to a tracer counterpart.This raises the total number of compartments for liver, tracer injection and plasma modules to 56.The first 28 compartments represent the tracee subsystem as shown above (Q) and remaining 28 compartments are used to represent the tracer subsystem (q) in the following form, where q is the vector of tracer compartments, K is the matrix that holds transfer rates that is defined above and u is the external input, which represents the labeled leucine and glycerol injections.This way, the enrichment of the metabolite in compartment i is calculated as the ratio of q(i) and Q(i).

Gastrointestinal module
Two hours after the tracer injection, patients were asked to consume a liquid mixed meal.Nutritional content of the mixed meal is given in Table S1.In the model, digested lipids are absorbed and incorporated into chylomicrons (CM) in the form of TG.The gastrointestinal module is represented by an 11 compartments system (Fig. S2), where the compartments are labelled with the latter D. Compartments D1-3 represent stomach, D4-9 represent small intestines.Lipid absorption and TG secretion into the plasma is only allowed from intestinal compartments.Colon is represented with a single closed compartment (D10) and used to calculate the amount of unabsorbed lipids.Compartment D11 represents plasma chylomicron TG concentration, which links gastrointestinal module to the hepatic module through competition for clearance.Two hours after the tracer injection, patients were asked to consume a liquid mixed meal, which consisted of two servings of Fresubin® Protein Energy Drink enriched with a palmitate shake (Table S1).Due to the reduced stomach size, it took longer for patients to consume the meal after the surgery.Therefore, the meal input into the first stomach compartment D1 is represented with a piecewise continuous square wave function in the following form; where mTG is the approximate amount of TG that can be produced from consumed lipids, tc is the average time it took for patient to consume the meal in hours and it is set to 0.25 and 0.5 hours for pre-and post-surgery conditions, respectively.For converting consumed lipids to TG and calculating mTG, it was assumed that 100gr TG contained 95gr fatty acids, as suggested before (Ratnayake & Galli, 2009).
In the post-surgical configuration, the stomach is reduced to a single compartment (D1) to account for the reduced stomach size, and this compartment is directly attached to the second intestinal compartment to represent the bypassed duodenum.In Figure S2, bypassed stomach and intestinal compartments after the surgery are shown with faded colors, where biliopancreatic limb is only shown for illustration purposes.The rate of change of the amount off lipids in gastric compartment Di is given by the following differential equation, Half gastric emptying time (T1/2) following a liquid test meal in healthy humans was measured as 80 minutes on average, with no correlation with age, sex or BMI (Hellmig et al., 2006;Mishima et al., 2009).Taking this as a reference, and solving system of equations given by S8 for all three gastric compartments simultaneously, parameter ks is estimated as 2.005 h -1 .Furthermore, using identical transfer and absorption rate coefficients for intestinal compartments (kf and ke, respectively) and assuming a known malabsorption coefficient, the relation between intestinal lipid absorption and transfer rates, ke and kf, can be established in the following form; where, ma is the malabsorption coefficient and dc is the number of intestinal compartments for pre and post-surgery configurations.The malabsorption coefficient is calculated by the ratio of the lipids in the column to the total amount of consumed lipids.We used the values estimated by Odstrcil et al.
for pre-and post-surgery lipid malabsorption coefficient as, 8% and 27% respectively (Odstrcil et al., 2010).Same absorption coefficient is used for all intestinal segments, except for the first intestinal compartment in the post-surgical configuration (D5, the roux-limb).A separate parameter was employed for the lipid absorption from this compartment.However, the estimated value of this parameter approached to 0 during the model fit.Therefore, it was fixed to 0 in the final model.
The diagram above shows a simple two compartment system to illustrate the transfer dynamics between model compartments.In the model, the rate of change of material in the compartment Vj is given by the following differential equation; Here, kj,i is the transfer rate of material from compartment Vi to compartment Vj and defined as follows; where,  G $,& is the constant transfer rate and pj is the affinity of lipoprotein species Vj to the LPL.The binding probability of lipoprotein species Vj to LPL depends the ratio of its size to the average lipoprotein particle size in circulation.Since VLDL1, VLDL2 and chylomicrons are triglyceride rich lipoproteins, we assume that the size of each particle is proportional to its TG content.In the preprandial state, plasma concentrations are constant.Therefore, the affinity for each lipoprotein species (pj) should be fixed to 1 and the clearance rate should be solely determined by  G $,& , which is estimated from the data.In the postprandial state, the affinity becomes a dynamically changing function of the ratio of the size of the corresponding particle to the average lipoprotein particle size in circulation.Consequently, increased particle size increases the affinity for the corresponding lipoprotein species, where increased total plasma TG concentration reduces its affinity through competition.This way, the clearance rate for each lipoprotein species depends on the concentration of other lipoproteins as well.Under these assumptions, the affinity of the particle Vj to LPL is given by, where, VTGj and VapoBj are the amount of TG and apoB in the Vj fraction, respectively.The subscript b is used for baseline values for the corresponding metabolite.By normalizing variables over their baseline values, pi becomes a dimensionless parameter that is 1 during the steady state.The total plasma apoB concentration virtually remains the same in the post-prandial state.Therefore, the baseline and post-prandial total plasma apoB terms cancel out in the denominator.Finally, since the relation between volume and surface area is nonlinear, the term is raised to na th power, which is left as a free parameter to be estimated from the data.

Modeling insulin action
Insulin module consists of 12 transient delay compartments, which are labelled with the letter I in Figure S2.Compartment I1 represents plasma insulin concentration.Since the relation between plasma insulin and the components of the lipoprotein metabolism is captured as a unidirectional link, the plasma insulin is used as a forcing function in the model.In order to obtain a smooth and continuous plasma insulin curve, insulin input is approximated with the following bell-shaped Gaussian function; where  is time,  determines the height of the peak,  is the peak location and parameter  determines the width of the curve.When used as input, this curve is flexible enough to generate the measured plasma insulin concentration.Using Eq.S13 as input, the rate of change of plasma insulin concentration is given as; Here Ib is the basal insulin concentration in pmol/L and  '&(6  > enforces insulin concertation to remain at the basal level when there is no input from f. Since lipoprotein metabolism does not feedback onto the insulin module, the parameters of the insulin input function f(t) (, , ) are separately estimated for each patient and used in the model.Estimated values for ,  and  for the population average are given in parameter table below (Table S2).Parameter k1ins is fixed to the population average 2 h -1 to ensure better control over continuous insulin input curve.This way, input to the first delay compartment is fixed and the delay associated with insulin action is controlled by the free delay parameter that was estimated in the model.Resulted plasma insulin concentration for each patient is given in Figure S3.As it seen from the Figure S3, method produces a good fit for each patient.
Although, the method leaves out some of the dynamic properties, such as the biphasic responses produced by some patients, it must be noted that these small deviations would not have a potent effect in the response compartment due to the long delay.

Insulin mediated LPL stimulation
Insulin regulates lipoprotein lipase at transcriptional and posttranslational levels (Fried, Russell, Grauso, & Brolin, 1993;Sadur & Eckel, 1982), which is compromised in insulin resistant patients (Panarotto, Rémillard, Bouffard, & Maheux, 2002).In the model, plasma insulin concentration is transmitted to the insulin response compartment (Ir) through a delay, which is given by 10 consecutive delay compartments (I2-11).This delayed insulin response is then incorporated into the corresponding transfer coefficient given by Eq.S11, which reads as follows: where Rmax is the maximal response that can be produced by insulin and it is fixed to 0.25 based on the approximate value obtained during a hyperinsulinemic clamp study (Sadur & Eckel, 1982).Kins is the insulin responsiveness parameter and Irbase is the response produced by the basal insulin concentration.Therefore, the term in the brackets in Eq.S15 is equal to 1 in the fasting state for all patients and it dynamically responds to the insulin change over basal level in the postprandial state.
The potency of the insulin-mediated postprandial response depends on the insulin response parameter Kins, which represents plasma insulin concentration over the basal level that is necessary to evoke the half maximal response,  7?@ /2.Kins is left a free parameter and estimated for each patient from the data, where 1/ &(6 was added to the cost function as a regularization term to prevent Kins to converge arbitrarily small values.A smaller Kins value corresponds to a higher insulin responsiveness. On the other hand, for very large Kins values the term in the square brackets in Equation 1 approaches to 1, which indicates that there is no detectable dynamic insulin-mediated lipoprotein lipolysis stimulation in the model.Hence, the reciprocal of Kins is a measure of insulin-mediated stimulation of lipoprotein lipolysis.We introduce lipoprotein lipolysis insulin sensitivity index (lipoprotein lipolysis ISI) in arbitrary units as follows; where, Ks is a scaling factor and set to 20 to contain the calculated indices between 0 and 1 for convenience.
Insulin mediated suppression of the hepatic apoB production has also been incorporated into the model as a delayed forcing signal generated by the portal vein insulin concentration, where portal vein insulin signal was derived from the plasma insulin data.However, insulin-mediated effects on the

Parameter Estimation and calculating production, fractional catabolic and transfer rates
Kinetic rate parameters were estimated by employing MATLAB's optimization toolbox.For each patient, model parameters were estimated by minimizing a cost function, defined as the weighted sum of squared differences between model simulation and the data.For the model simulations reported in Figures 1 and 2 reported in the results section in the main text, the computational model was calibrated using population average data and obtained parameter values are given in Table S2.
During the parameter estimation process, multiple searches were initiated from randomly selected points in the parameter space.The parameter search resulting in the minimum error was selected as the optimal parameter set.
For model calibration and analysis following data is used; plasma VLDL1,2 apoB and TG concentrations, plasma VLDL1,2 apoB and TG enrichments, plasma chylomicron TG concentration, plasma leucine and glycerol enrichments, plasma insulin concentration.For each patient, parameter estimation is carried out by minimizing a cost function, which is defined as the sum of squared differences between model simulation and the data.The model consists of a system of ordinary differential equations that governs the rate of changes of internal state variables in the following form: where  ⃗ is the vector of the internal state variables,  _⃗ is the vector of external inputs and  ⃗ is the vector of model parameters (e.g.rate constants, transfer coefficients, initial conditions and scaling factors), and  ⃗() is the measurements.In the model, the first 28 entries of the vector  ⃗() represents the tracee subsystem (Q) given with the hepatic, plasma and tracer injection modules in Figure S1.Q represents unlabeled biochemical concentrations.State variables x29-56 represent the tracer subsystem (q), where variable q(i)=x(i+28) corresponds to the tracer counterpart of the tracee compartment Q(i).State variables x57-68 represent gastrointestinal tract (D).Finally, variables x69-80 represent insulin subsystem (I).
Since the system in hand is partially observed, measurements  ⃗() can be performed on a limited number of state variables and their certain combinations.Namely, plasma VLDL1,2 apoB and TG concentrations, plasma VLDL1,2 apoB and TG enrichments, plasma chylomicron TG concentration, plasma leucine and glycerol enrichments, plasma insulin concentration.In eq.S18, g is the mapping that associates internal state variables to the observables  ⃗() with scaling parameters in  ⃗ .The plasma biochemical concentrations are the sum of the corresponding compartments in the tracee subsystem.For instance, plasma VLDL2 TG concentration is given by Q( 16)+Q( 17).In the other hand, the enrichment of the corresponding metabolite is given by the tracer to tracee ratio (TTR), which is the ratio of the sum of the tracer and tracee compartments.For instance, the TTR for the plasma VLDL2 TG is calculated by the following; The model parameters ( ⃗ ) are estimated by minimizing the difference between model simulations and experimental data, which is measured by the cost function ^ ⃗ a that calculates the weighted sum of squared differences between model simulation and measured data; where  5,$ ( P )  the data for the j th observation measured at the k th time point and  5,$,P is the measurement error associated with this data point.In the parameter estimation process, the parameter space Q is searched for the parameter values that minimize the cost function, ^ ⃗ a.For normally distributed error with zero mean this corresponds to the maximum likelihood estimate (MLE).During the parameter estimation process, multiple searches were initiated from randomly selected points in the parameter space generated by Latin Hypercube Sampling, and MATLAB's nonlinear optimization function lsqnonlin was used.The parameter search that was resulted in minimum error was selected as the optimal parameter set,  ⃗ * .In order to evaluate the uncertainty associated with the estimated parameters, asymptotic confidence intervals were calculated by using the covariance matrix C=F -1 , where F is the Fisher Information matrix, and given by the negative hessian of the log-likelihood function calculated at the optimum point.This provides information about the curvature of ^ ⃗ a in the neighborhood of the optimum point  ⃗ * , where a flat profile in the direction of a certain parameter indicates unidentifiability and results in an infinite confidence bound.We approximated the hessian of the log-likelihood by the first order model sensitivities given by the Jacobian matrix that is calculated during the parameter estimation, H=J T J (Vanlier, Tiemann, Hilbers, & van Riel, 2013).The parameter confidence interval for the i th parameter is given by; where,  T,!3 is the  quantile of the t-distribution with  degrees of freedom.Estimated model parameters and calculated confidence intervals for the population data are provided in Table S2.2,84±1,24 2,5±0,06 d 0,62±0,49 0,23±0,09 (Pre: pre surgery estimated value ± confidence bound, Post: post-surgery estimated value ± confidence bound) After model parameters were estimated for each patient, we calculated the production/secretion (PS), fractional catabolic (FCR) and fractional transfer rate (FTR) parameters by utilizing calculated pool sizes and kinetic transfer rates.For calculating each kinetic parameter, the sum of fluxes from associated compartments was divided by the sum of corresponding pool sizes.Flux from compartment i to j is calculated as F(j,i)=k(j,i)Q(i).Afterwards, PS, FCR and FTR parameters were calculated by employing fluxes and pool sizes by the following equations; Pre-surgery (Table S3) and post-surgery (Table S4) kinetic rates are calculated for each patient.Since the postprandial pool sizes and transfer rates dynamically change over time, postprandial kinetics were calculated for every 30 minutes and their averages over their respective time-courses were represented as the final value.Physiological parameters for the population averages were calculated and reported as mean ± standard deviation (sd).

Figure S3 :
Figure S3: Pre-(red) and post-surgery (blue) plasma insulin concentrations simulated by the time dependent input function given in Eq.S13 (N=24).The last panel is the population average.
pathway remained undetectable due to the uncertainty associated with portal vein insulin concentrations and resulting parameter estimates.Therefore, insulin mediated suppression of the hepatic apoB production was removed from the final version of the model.

Figure S5 :
Figure S5: (A) Pre-and post-surgery lipoprotein lipolysis insulin sensitivity index (lipoprotein lipolysis ISI).(B) Pre-and post-surgery insulin mediated peripheral lipolysis suppression per insulin increased over basal level (D(%)insulin) from the clamp studies vs lipoprotein lipolysis ISI.(C) Pre-and postsurgery insulin mediated endogenous glucose production (EGP) suppression per D(%)insulin vs lipoprotein lipolysis ISI.(D) Pre-and post-surgery insulin mediated glucose rate of disappearance (Rd) stimulation per D(%)insulin vs lipolysis ISI.In each panel, the regression lines for the pre-and postsurgery data are shown in red and blue, respectively.Correlation coefficients (r) and associated pvalues are reported in each panel.

Table S1 :
Content of the mixed meal.

Table S2 :
Pre and post-surgery estimated parameter values for population average data.