A Mathematical Model of Dietary Lipid Absorption and Postprandial Chylomicron Dynamics
1Department of Biomedical Engineering, Linköping University, Linköping, Sweden
2Department of Experimental Medical Science, Lund University, Lund, Sweden
3SUND sound medical decisions; Linköping, Sweden
4Department of Health, Medicine and Caring Sciences, Linköping University, Linköping, Sweden
*Corresponding author: Gunnar Cedersund, Email: gunnar.cedersund@liu.se. Address: Department of Biomedical Engineering, Linköping University, 58185 Linköping, Sweden. Phone: +46-702-512323Abstract
Obesity and related conditions such as dyslipidemia impose an increasing burden on healthcare systems worldwide. These conditions are associated with altered postprandial chylomicron (CM) metabolism, the elusive and critical first step in lipid metabolism. This step remains elusive because it is governed by large interindividual variations and a complex set of intestinal processes. In particular, the second meal effect (SME) implies that enterocytes release previously stored fat during subsequent meals. To deal with this complexity, CM and lipid metabolism have previously been explored using mathematical modeling. However, existing models primarily describe TAG dynamics following a single meal or are too complex for practical personalization across datasets. Herein, we address these limitations by presenting a small-scale mathematical model of CM dynamics that incorporates the SME. The presented model successfully describes data from six clinical studies of both single and repeated meal interventions. Model performance was further evaluated by predicting independent datasets using a BMI-dependent calibration. Finally, to demonstrate model applicability, we simulated full-day responses consisting of three sequential meals in individuals with varying BMI values, with qualitative agreement to clinical observations. This work supports our understanding of the SME, person-specific CM postprandial responses, and mechanisms underlying obesity.
Article notes
Competing Interest Statement
The authors have declared no competing interest.
Introduction
Obesity and related metabolic conditions, such as insulin resistance, atherosclerosis, and type 2 diabetes (T2D), are placing an increasing burden on our healthcare systems. In T2D research, mathematical models of insulin and glucose regulation have substantially advanced our understanding of disease progression and helped facilitate the development of treatment strategies (1,2). In contrast, fewer mathematical models have focused on postprandial lipid metabolism. Dyslipidemia, characterized by elevated circulating lipid levels, is an important feature of metabolic disease progression and may in part be driven by the postprandial appearance of chylomicrons (CM) (3). As CM dynamics are highly complex, incompletely understood, and subject to substantial inter-individual variation, the development of robust mathematical models of CM appearance and postprandial metabolism may be important for advancing our understanding and treatment of metabolic diseases.
A mathematical model of postprandial CM dynamics needs to describe key biological mechanisms. One such mechanism is the intestines’ ability to store triglycerides (TAG) for later release, which contributes to a lipid absorption phenomenon known as the ‘second meal effect’ (SME) (4,5). The SME describes how subsequent meals will result in higher CM-TAG release compared to the first meal, even when the meals are identical (6). Several studies have observed the SME in humans (4,7–11), demonstrating complex relationships between CM-TAG release, prior meal composition, and the ingestion of other macronutrients. For example, Evans et al. showed that a low-fat meal, after a high fat breakfast, still triggered a CM-TAG release (4). Moreover, CM-TAG appearance differs between individuals, with reported differences related to sex, insulin sensitivity, and lipid-form (emulsion or solid) (12,13).
Only a handful of mathematical models have been developed to describe CM-TAG release. One example of such a model is the model by Leohr et al. (14), a semi-physiological model that was used to quantify the CM-TAG response to a high-fat meal in both obese and healthy individuals. While this model describes CM-TAG data for a single meal, it does not aim to capture the SME. Other models describe TAG dynamics after a single meal (15), but these do not explain CM-TAG dynamics. Sosa et al. (11) focused primarily on CM release and investigated differences in CM appearance and the SME between insulin-sensitive and insulin-insensitive individuals, through both experimental data collection and systems biology analysis. Their analysis focused on a single study of sequential meals, and employed a mathematical model, supporting the hypothesis that stored enteral TAG bypasses re-synthesis and lipolysis, enabling rapid secretion. However, there is still a need for a simple but physiologically meaningful model, integrating data from different studies, capable of reproducing the complex dynamics of the SME, while remaining suitable for individualized parameterization.
Herein, we present such a mathematical model designed to describe CM plasma appearance from repeated mixed meals of several independent studies (Fig. 1A), supporting our understanding of the SME and mechanisms underlying obesity.
Results
In this study, we have iteratively developed and evaluated a mathematical model for plasma CM-TAG dynamics following both single and multiple meals (Fig. 1A). As a first step, the model was trained individually to six different studies (Fig. S1), where two model parameters exhibited substantial interpopulation variability (Figs. S2–S3). As a second step, the mathematical model was trained to all studies jointly while allowing group-specific variability (study and intervention) in the two parameters identified with substantial variability (Fig. 1B, model training). Both parameters were found to be negatively associated with the mean body mass index (BMI) reported in the corresponding studies, using linear regression. As a third step, this relation between BMI and model parameters was used to predict CM dynamics of an independent dataset not used in model training (Fig. 1B, model validation). Finally, the established model was used to investigate repeated sequential meals in individuals with varying BMI values (Fig. 1B, utility).
Individual model agreement with CM dynamics for single and multiple meals
The mathematical model of CM postprandial dynamics can describe six different studies individually (Fig. S1). Details regarding each intervention and cohort in the training dataset can be found in the supplementary material (Table S1). The agreement between model simulations and data was evaluated using χ2-tests with 95% confidence interval (see material and methods). The model passed the χ2-test for all studies individually (Table 1).
A parameter analysis was performed to evaluate the interstudy differences (Figs. S2-S3). We found that two parameters displayed substantial interstudy variability: kBreakdown (governing the rate of CM clearance), and kL1 (governing the rate of transport in intestine and enterocytes). This assessment was performed by comparing the parameter boxplot distributions (Fig. S2) and parameter-to-parameter correlations (Fig. S3). For the boxplot distributions, the confidence intervals (CIs) for both kL1 and kBreakdown appeared broader for intervention including single meals, in comparison to intervention including multiple meals. For the parameter correlations, the analysis indicates that kL1 and kBreakdown correlate, and thus both could potentially be described with a joint description. Together, these insights indicate that intervention and study variability could be described by allowing freedom in these two identified parameters.
By allowing freedom in the parameters kBreakdown and kL1, the model can simultaneously describe all estimation data (Fig. 2), confirmed by a χ2-test with 95% confidence interval (90.1<Tχ2=148.8). Additional tests were performed to evaluate if only one parameter was sufficient to represent interstudy differences (Fig. S4, freedom only in kBreakdown; Fig. S5, freedom only in kL1). We found that freedom in only one of these parameters could pass the χ2-test (kBreakdown, 144 < Tχ2=148.8; and kL1, 132.4 < Tχ2=148.8). However, both cases were not accepted, as they did not sufficiently capture the SME and qualitative assessment of the meal responses (Figs. S4-S5). Therefore, to more accurately represent interstudy variability, model development proceeded with allowed freedom in both parameters (Fig. 2).
Using the model to predict independent validation data of CM appearance
Values of the two parameters describing study- and population-differences (kBreakdown and kL1) have to be assumed to be able to make predictions of new individuals or meal interventions (Fig. 3). To address this, the relationship between the study-specific parameters and known study covariates, e.g., BMI, age, and Homeostatic Model Assessment for Insulin Resistance (HOMA-IR) (Table S1), was assessed via a linear-regression analysis (Fig. 3A-B). This was performed using the parameter-value (in loge) corresponding to the best model agreement with data. We found that the highest correlation was observed between BMI and kBreakdown with a Pearson correlation coefficient of R = 0.63 (Fig. 3A). A slightly smaller correlation coefficient was observed between BMI and kL1 (R = 0.57) (Fig. 3B). Using these regression equations, we could now get a plausible value of the calibration parameters based on the known BMI of an individual.
We used the BMI-dependent parameter calibration to make predictions of an independent study (18) (Fig. 3C-E). In this study, plasma CM was measured after a single meal in three different cohorts: women (Fig. 3C; mean BMI of 22), insulin sensitive men (Fig. 3D; mean BMI of 26), and insulin resistant men (Fig. 3E; mean BMI of 25). The model simulations agreed with the data from the women cohort (Fig. 3C) and the data from the cohort of insulin sensitive men (Fig. 3D), whereas the model underpredicted the CM for the insulin resistant men (Fig. 3E).
Model predictions of daily CM variations for different individuals and meal patterns
To highlight the usability of the constructed model, we investigated CM appearance during three sequential meals (breakfast, lunch, and dinner containing 15 g, 30 g, and 30 g of fat respectively) during a day (Fig. 4). The simulations were performed using different BMI inputs to visualize the relation between CM dynamics and metabolic dysregulation (Fig. 4A). Here, we can observe that increasing BMI (used herein as a proxy for metabolic dysregulation) is connected to higher peaks of CM during the day. To highlight this behavior further, we simulated separate CM dynamics for individuals with different BMI (Fig. 4B) and mean levels of circulating CM (Fig. 4C). Together, these three figures (Fig. 4A-C) indicate increasing CM peak value with higher BMI and higher mean CM throughout the day. More specifically, the model predicts postprandial chylomicron exposure approximately 2-fold higher for people with a BMI of 30 compared to a BMI of 20. This is in qualitative agreement with literature, where obese individuals (BMI 30 or above) have been observed to have a 2- to 3-fold higher postprandial chylomicron exposure than lean participants (19). Similar observations have been reported in other works (20,21).
We also simulated how meal patterns affect the mean level of CM during the day (Fig. 4D-E). The model predicted that increasing the time between meals reduces the SME, with a substantial reduction between having a 2-hours gap and a 5-hours gap between meals (Fig. 4D, blue compared to orange). Importantly, we found that mean levels of circulating CM also decreased with longer time periods between meals (Fig. 4E). In this meal-timing investigation, a SME is predicted when a meal is consumed 5-hours post the first meal (Fig. 4D, orange), which is in agreement with clinical observations (22). Furthermore, the model predicts a diminishing SME effect between 5-hours and 8-hours (Fig. 4D, orange compared to yellow), which is in qualitative agreement with reports of chylomicron remnants clear after approximately 6- to 8-hours after eating (23).
In conclusion of the meal intervention predictions, the model indicates qualitative agreement to many clinical observations in circulating CM levels and provides a small-scale mathematical model capable of simulating CM postprandial dynamics.
Discussion
We have constructed a mathematical model describing plasma CM appearance (Fig. 1A) and evaluated its ability to describe data and utility in several steps (Fig. 1B). First, we found that the model could sufficiently describe six different clinical studies individually (Fig. S1; confirmed using individual χ2-test). Second, we found that the most important parameters to describe group variability (study and intervention) were kL1 (governing enterocyte transport of lipid) and kBreakdown (CM clearance from plasma) (Figs. S2-S3). By allowing freedom in these two parameters and restricting variability in all other parameters, the model was able to simultaneously fit all estimation data (Fig. 2), confirmed by a χ2-test (90.1 < 148.8 = Tχ2). A relationship between the intervention-specific kL1 and kBreakdown to study-specific BMI was found using linear regression (Fig. 3A). Third, the established BMI-specific calibration was used to predict new independent validation data (Fig. 3B). Finally, we investigated how the model could be used to predict postprandial chylomicron dynamics of various meal interventions for people with various BMIs (Fig. 4). In summary, we have constructed a small-scale mathematical model of CM appearance and the SME, providing a foundation for creating and extending models describing postprandial lipid metabolism.
The model’s ability to predict SME under various conditions needs to be further evaluated. When the model was trained using only the Barrows et al. (9) dataset, the model could produce a sharp second peak (Fig. S1L), but when trained simultaneously to all studies, the peak is less prominent (Fig. 2B). This is mainly because the SME peak magnitude varies from study to study, e.g., peak in Vors et al. (13) compared to Barrows et al. (9) (Fig. 2L compared to Fig. 2B). These discrepancies could be due to several factors not accounted for in the mathematical model and in many cases are not reported in studies herein, such as differences in how CM was measured, the form of consumed fats (oil, emulsion, or spread), or interpopulation variability not captured by BMI (e.g., insulin sensitivity). While BMI and insulin sensitivity have been reported to correlate significantly (24), we found weak agreement with the insulin resistant group in the Hodson et al. study despite using the BMI-dependent calibration (Fig. 3E). We advise future works to implement insulin sensitivity in their model predictions to potentially overcome this shortcoming. Despite these listed model limitations and potential future works, the mechanisms incorporated in our small-scale model are sufficient to describe differences in the SME (Figs. 2-4; as confirmed by individual χ2-test and qualitative agreement to literature).
There are several biological mechanisms that were omitted from the presented model. We do not include the production of CMs from plasma-derived non-esterified fatty acids (NEFA) and glycerol, nor the conversion of dietary nor plasma derived carbohydrate conversion into CM. Furthermore, there are other aspects to expand upon: gastric emptying crosstalk interactions between different macronutrients and meal volumes, and the fact that different fatty acids may have different metabolic rates. We have recently published an ethanol rate of appearance model, with a gastric emptying module (25), and future work could involve integration of such gastric emptying mechanisms to better describe lipid absorption, or interaction with other foods. Furthermore, as mentioned earlier, in a recent work by Sosa et al. (11) a comprehensive model is presented, which includes several of these pathways. Thus, in future work, we will evaluate if one could make use of the description of more detailed mechanisms described in the Sosa et al. model, to achieve a better simultaneous agreement to all data.
The relationship between covariates such as BMI and HOMA-IR in relation to their effect on chylomicron dynamics needs to be evaluated. There exists a known relationship between metabolic dysregulation and elevated levels of circulating TAGs, where mechanisms are e.g. poor clearance of CM-TAGs, or increased Very Low Density Lipoprotein (VLDL)-TAG secretion (26). The model analysis indicated the clearance of CM from plasma (kBreakdown) as significant for describing interstudy variability, and we show that this parameter has a tentative correlation with BMI (Fig. 3A). The perceived correlation with BMI and model parameters could be reasoned to be explained by the decreased clearance due to a change in insulin sensitivity. Yet, when we tested this correlation, we did not appear to find a similar correlation between the breakdown parameter and HOMA-IR (Fig. S6). A larger dataset with more covariate differences would be beneficial to understand these relationships further. Nevertheless, doing model personalization as presented herein could possibly be a strong approach if further developed, connecting rate parameters with individual anthropometrics, after the model training (not during).
There are several additional key aspects to consider in our study. The first aspect is regarding data availability. All data used for model estimation and validation is digitized data, thus the true data points might be slightly different from the ones we used. In many cases there was a lot of covariate data missing (such as BMI, age, and HOMA-IR), which excluded several studies from the analysis done (4,9). Furthermore, a delimitation of the study was not to include plasma TAG data, which further excluded a lot of available datasets. By including plasma TAG, the model complexity would also need to increase, for inclusion of various lipoproteins, i.e. VLDL. We have previously developed models for NEFA metabolism (27), lipolysis (28), glucose-insulin regulation (29,30), and multi-level and multi-timescale models for insulin resistance progression (31,32). An integration of all these models into a singular framework could be capable of describing mixed meal postprandial metabolism and therefore connecting more complex mechanisms of obesity progression and dieting. Herein, we present a step towards such a more complete mathematical model of human meal responses.
In conclusion, we have created a model for plasma CM appearance capable of describing dynamics related to the SME. The mechanism behind CM, such as person-specific variability and form of consumed type of dietary fats (oil, emulsion, or spread), is still not fully understood, and there seems to be large inter-patient, and inter-cohort variations in CM dynamics that will need further work to be elucidated. Nevertheless, our presented model is usable and could be integrated into any mathematical model for meal simulations, needing a description of CM appearance and postprandial lipid dynamics.
Method
The code to reproduce all reported observations can be found in our GitHub repository (https://github.com/chrsi30/CMDT).
Modelling approach and software
The model was formulated using ordinary differential equations (ODEs). The mathematical analysis, model simulation, model formulation and the model parameter estimation were all performed in MATLAB (2022a; The MathWorks, Natick, MA), using the systems biology toolbox (IQM) (33). The parameter estimation was performed using the extended scatter search (ESS) optimization algorithm from the MEIGO toolbox (34). The quality of the resulting optimal simulations was evaluated using a χ2-test (35). Lastly, the parameter and model uncertainty analysis was done using Markov-Chain Monte Carlo sampling implemented in the PESTO toolbox (36).
Parameter estimation
Parameter estimation was done by quantifying the model performance, using the model output ŷ to calculate a weighted least squares cost function (Eqs. 1-2). where θ is the model parameters; yi(tj) is the measured data from a study i, at time point j and from one type of measure k; ŷi(tj,q) is the simulation value for a given experiment setup i and time point j, and SEM is the standard error of the mean, which is the sample standard deviation,σi(tj) divided with the square root of the number of repeats, ni(tj) at each time point. The value of the cost function, V(θ), is then minimized by tuning the values of the parameters.
To evaluate the model, a χ2-test was performed for the size of the residuals, with the null hypothesis that the experimental data have been generated by the model. The test is performed under the assumption that the experimental noise is additive and normally distributed (35). In practice, the cost function value was compared to a χ2-test statistic, T°, with a cumulative density function (Eq. 3).
where is the inverse density function; and α is the significance level (α = 0.05 was used) and v is the degrees of freedom, equal to the number of data points in the training dataset. In practice, the model is rejected if the model cost is larger than the χ2-threshold (T°).
Model description – the chylomicron model
The ODE model was developed by the iterative process of physiologically based data-driven hypothesis testing. The capability of several different candidate model structures was tested using the Barrows et al. dataset. The candidates were evaluated based on their ability to accurately describe the experimental data using a χ2-test, with special emphasis on the ability to reproduce the sharp second meal peak. Furthermore, the candidate structure complexity was also evaluated based on the Akaike information criterion (AIC). The results of this evaluation can be seen in the supplementary material (Table S2). The chosen model (‘M9’) had the lowest χ2-test statistic, and it could best describe the sharp second meal peak observed in data. This section will detail the model equations.
The model has a total of seven model states and five model parameters (Fig. 1A, right). The first state represents TAG that has entered the stomach and is moving towards the intestine and has notation LS (Eq. 4). where; ufat is the model input, the mass of triglycerides from the meal, expressed as mmol (converted from grams using a TAG molar mass of 845mol/g), and kL1 (min-1) is the rate parameter for the absorption and transfer of lipids into and within the enterocytes. TAG is then moved through the Enterocytes and prepared for storage or plasma release. This part is represented by the three model states: L1 (Eq. 5), L2 (Eq. 6), and L3 (Eq. 7). Where kL1 again is the rate of transfer. where kL1 is reused from (Eq. 5) to describe the transfer rate moving into the next enterocyte state. In (Eq. 6) enterocyte lipid storage is introduced, which is the movement of TAG into lipid droplets represented by the state LIP. The transfer into and from the LIP state is represented by the parameters kLIP1 (min-1) and kLIP2 (min-1) respectively. Here we also introduce the phenomenological effect of the second meal, denoted sng, which affects the efflux from the LIP state. The sng effect is a model state which has the following derivative (Eq. 8): where ucarb and uprotein are the carbohydrate and protein portions of the meal. The last enterocyte state, L3, has the following derivative (Eq. 9): where kL1 is reused from (Eqs. 5 and 6) as the rate of TAG-chylomicron excretion into plasma. The plasma TAG-chylomicron state is denoted Chp and has the following derivative (Eq. 10): Where, kBreakdown (min-1) is the breakdown and usage of TAG-chylomicrons in the body. To obtain the concertation in plasma the Chp state is divided into with plasma volume, calculated via the Nadler formula, based on sex, height and weight (37).
Data curation
Data was gathered from a total of seven published studies (4,5,9,12,13,16–18). The estimation data includes six of these studies and was comprised of a total of 10 different cohorts (Table S1), and various interventions. The interventions that were investigated in each study are described in the supplementary material (Table S3). Inclusion criteria for data were that TAG chylomicrons concentration should be measured in connection to a meal, and detailed information about the meal, and any prior meal should be detailed. All data was digitized using the WebPlotDigitzer.
Funding information
GC acknowledges support from the Swedish Research Council (2023-03186, 2023-05460), the Horizon Europe project STRATIF-AI (101080875), VINNOVA (VisualSweden), ALF (RÖ-1001928), and the Exploring Inflammation in Health and Disease (X-HiDE) Consortium - a strategic research profile at Örebro University funded by the Knowledge Foundation (20200017). EN acknowledges support from Zenith and the Swedish Fund for Research without Animal Experiments.
Conflict of interest
Gunnar Cedersund is the owner of SUND sound medical decisions AB (unrelated to the work presented herein). Oscar Silfvergren and Kajsa Tunedal are employees of SUND sound medical decisions AB.