Abstract
A physiologically based pharmacokinetic (PBPK) model of vitamin D3 and metabolites [25(OH)D3, 1,25(OH)2D3, and 24,25(OH)2D3] is presented. In this study, patients with 25(OH)D3 plasma concentrations below 30 ng/ml were studied after a single dose of 5000 I.U. (125 µg) cholecalciferol, provided with 5000 I.U. daily cholecalciferol supplementation until vitamin D replete [25(OH)D3 plasma concentrations above 30 ng/ml], and had serial plasma samples were collected at each phase for 14 days. Total concentrations of vitamin D3 and metabolites were measured by ultra-high performance liquid chromatography tandem mass spectrometry. A nine-compartment PBPK model was built using MATLAB to represent the triphasic study nature (insufficient, replenishing, and sufficient). The stimulatory and inhibitory effect of 1,25(OH)2D3 were incorporated by fold-changes in the primary metabolic enzymes CYP27B1 and CYP24A1, respectively. Incorporation of dynamic adipose partition coefficients for vitamin D3 and 25(OH)D3 and variable enzymatic reactions aided in model fitting. Measures of model predictions agreed well with data from metabolites, with 97%, 88%, and 98% of the data for 25(OH)D3, 24,25(OH)2D3, and 1,25(OH)2D3, respectively, within twofold of unity (fold error values between 0.5 and 2.0). Bootstrapping was performed and optimized parameters were reported with 95% confidence intervals. This PBPK model could be a useful tool for understanding the connections between vitamin D and its metabolites under a variety of clinical situations.
SIGNIFICANCE STATEMENT This study developed a physiologically based pharmacokinetic (PBPK) model of vitamin D3 and metabolites for patients moving from an insufficient to a repleted state over a period of 16 weeks.
Introduction
Vitamin D is a fat-soluble, prohormone that plays an essential role in regulating calcium and phosphorus to maintain musculoskeletal health. The main source of vitamin D is through endogenous production of cholecalciferol (VitD3) in the skin upon UV B exposure from the sun. The active form of VitD3, 1,25(OH)2D3 (calcitriol, 1,25D3) has been found to exhibit many pleiotropic actions beyond calcium and phosphorus homeostasis (e.g., musculoskeletal health), that are often referred to as “nonclassical” actions. Additionally, 1,25D3 controls over 200 genes, including those responsible for the regulation of hormone secretion, immune function, and cell proliferation and differentiation (Holick, 2006).
VitD insufficiency is a global health issue, and it has been estimated that up to 80% of men and women in the United States, Canada, and Europe meet this classification (Ganji et al., 2012; van Schoor and Lips, 2018). VitD status is based on concentrations of 25-hydroxyvitamin D, 25(OH)D (calcifediol, 25D), given that it is the major circulating form and it is readily assayed in hospital laboratories. There is a lack of consensus on the serum concentrations of 25D considered to be adequate, but it is commonly agreed that 25D concentrations should be above 20 ng/ml (1 ng/ml = 2.5 nmol/l) (Ganji et al., 2012, van Schoor and Lips, 2018). The Institute of Medicine (IOM) defines deficiency as 25D < 12 ng/ml and recommends a target serum concentration of 20 ng/ml (Ross, 2011). The Endocrine Society (ENDO) defines VitD deficiency as 25D < 20 ng/ml and recommends a target concentration of 30 ng/ml (Holick et al., 2011). However, other leading experts define VitD insufficiency as 25D concentrations between 20 and 30 ng/ml and deficiency as concentrations <20 ng/ml (Dawson-Hughes et al., 2005; Holick, 2007; Holick et al., 2011; Cianferotti and Marcocci, 2012). Although many people with VitD insufficiency take supplements, treatment recommendations are not consistent. The IOM recommended a daily intake of 600 I.U. VitD for children, adolescents, and adults, and 800 I.U. for adults over the age of 70 to maintain 25D concentrations of 20 ng/ml (Jernigan and Andress, 2003). ENDO recommended 1500–2000 I.U./d to prevent or treat VitD insufficiency and achieve 25D concentrations above 30 ng/ml, with a preferred range of 40–60 ng/ml (Holick et al., 2011). Whereas these are the common guidelines in the United States, recommendations in other regions of the world differ based on many population-specific factors, such as sunlight exposure, skin pigmentation, clothing, and dietary practices (Pérez-López et al., 2012; Society, 2012; Płudowski et al., 2013; Rizzoli et al., 2013; Munns et al., 2016; Haq et al., 2018).
Activation of VitD occurs through two sequential hydroxylation steps to generate the metabolically active metabolite 1,25D. Vitamin D is hydroxylated in the liver primarily by the cytochrome P450 (P450) enzyme CYP2R1 but also through pathways with CYP27A1 (Sakaki et al., 2005) and CYP3A4 (Battault et al., 2013; Wang et al., 2013), to form 25D, the major circulating metabolite. 25D can then be hydroxylated in the kidney by either CYP27B1 to the active metabolite, 1,25D, or by CYP24A1 to the inactive metabolite 24,25-dihydroxyvitamin D (24,25D). Phase 2 metabolism pathways have also been described for VitD but appear to be of lesser importance (Wang et al., 2014; Gao et al., 2017; Wong et al., 2018). 1,25D is tightly regulated in a feedback loop with CYP27B1 and CYP24A1; high concentrations of 1,25D will respectively suppress and stimulate expressions of CYP27B1 and CYP24A1.
Despite VitD being used clinically for decades, there are a limited number of published studies that include assessments of the parent VitD compound and metabolites using pharmacokinetic or physiologically based pharmacokinetic (PBPK) approaches (Kimura et al., 1991; Levine and Song, 1996; Bailie and Johnson, 2002; Armas et al., 2004; Ilahi et al., 2008; Roth et al., 2012; Benaboud et al., 2013; Jetter et al., 2014; Meekins et al., 2014; Ocampo-Pelland et al., 2016, 2017; Ramakrishnan et al., 2016; Fassio et al., 2020; Hsu et al., 2021). The current study sought to develop a PBPK model of VitD3 and metabolites in VitD insufficient patients in the United States who were treated with moderate (5000 I.U.) daily doses of VitD3 for up to 16 weeks to achieve replacement as defined by 25D ≥ 30 ng/ml. Through establishing a model that might align well with a triphasic clinical experience of a patient–an initial visit for establishing VitD concentrations, time to achieve repletion on a dosing regimen, and a follow up visit–this study aimed to better understand the longitudinal effects of daily supplementation of moderate doses of VitD3.
Materials and Methods
Design Overview
The parent study from which the data were obtained was conducted in accordance with the Declaration of Helsinki. All subjects provided informed consent and the protocols were approved by the Institutional Review Boards at the University of Colorado and the University of Pittsburgh. The parent study (NCT02360644) was registered under ClinicalTrials.gov and can be referenced for additional details.
Study Participants
Healthy subjects from the University of Colorado and University of Pittsburgh with VitD insufficiency defined in the study as total 25D below 30 ng/ml and not receiving VitD replacement therapy were recruited. Other eligibility criteria included age 18 to 75 years, predicted compliance with study visits, not pregnant or lactating, no changes in medications within 4 weeks, no predisposition to hypercalcemia, and hemoglobin ≥10 g/dl. Exclusion criteria were active autoimmune disease, active or recent infections requiring antimicrobial treatment, and hepatic insufficiency. Note that while the parent clinical study focused on patients with chronic kidney disease, this study used data from a control group of healthy study participants from the parent clinical study. Baseline characteristics for study participants (n = 11) can be found in Table 1.
Baseline characteristics of study participants (n = 11)
Clinical Study Design
Vitamin D insufficient subjects were admitted to the Clinical and Translational Research Centers (CTRC) at the University of Colorado or University of Pittsburgh for a 12 hour stay, followed by visits at 24, 48, 168, and 336 hours. Participants were required to be fasting at the beginning of the study. Any prescribed medications were withheld for the first two hours of the study. Subjects were given a single 5000 I.U. (125 µg) oral dose of VitD3 (Jarrow Formulas, Los Angeles, CA). Serial blood samples (7.5 ml) were collected at baseline and at 0.5, 1, 2, 4, 8, 12, 24, 48, 168, and 336 hours into heparinized vacutainers. Blood samples were centrifuged immediately following collection for 10 minutes at 3000 ×G at 4°C and plasma samples were stored at −80◦C until analysis. After the blood draw at 336 hours, participants were given up to 16 weeks of daily supplementation of 5000 I.U. VitD3 to attain replete concentrations (≥30 ng/ml). At the repleted study phase, participants were given a final dose of 5000 I.U. VitD3 followed by collection of blood samples as previously described for up to 336 hours. The total study reflected three phases–the insufficient phase, the replenishing phase, and the repleted phase. Data are available for subjects in the insufficient and repleted phases.
Analytical Assay
Total concentrations of VitD3 and metabolites [25D3, 1,25D3, and 24,25D3] were determined by a novel ultra-high performance liquid chromatography-tandem mass spectrometry (UHPLC-MS/MS) assay capable of detecting all four analytes simultaneously as previously described (Stubbs et al., 2014) with minor modifications. UHPLC was performed with a Waters Acquity UPLC I-class (Waters, Milford, MA, USA), which includes a sample manager and a binary solvent manager. Briefly, 500 µl samples were precipitated with acetonitrile, extracted with methyl tert butyl ether, then derivatized with 4-phenyl-1,2,4-triazoline-3,5dione. Separation of derivatized VitD analytes was achieved using a Waters Acquity BEH C18 column (150 mm × 2.1 mm, 1.7 µm particles) with a gradient elution of water with 0.1% formic acid and acetonitrile. The flow rate was 500 µl per min and the total run time was 8 minutes. Detection of analytes was achieved using positive atmospheric pressure chemical ionization and selected reaction monitoring on a TSQ Quantum Ultra triple quadrupole mass spectrometer (Thermo Scientific, San Jose, CA). Standards and quality control samples were constructed using blank human serum. Standard curve ranges were 0.10–15 ng/ml for VitD3 and 24,25D3, 0.0100–0.500 ng/ml for 1,25D3, and 1.0–100 ng/ml for 25D3. Mean correlation coefficients were ≥0.994 for all calibration curves. The within-run and between-run accuracy and precision percentage coefficient of variation were <10.6% for all analytes.
Physiologically Based Pharmacokinetic Model Development
A nine-compartment PBPK model (adipose tissue, brain, heart, kidney, liver, rapidly perfused tissue, slowly perfused tissue, and plasma) was developed for VitD3 and its metabolites using MATLAB (version R2020a, The Mathworks Inc, Natick, MA). A schematic of the overall PBPK model for each compound is shown in Fig. 1 with the connections between each metabolite PBPK model shown in Fig. 2.
General physiologically based pharmacokinetic model diagram for VitD3 and metabolites. Symbols are defined as the following: Q is the plasma flow rate, Cap is the arterial concentration of the compound; Cvp:T = AT/(VT ∗ PT) is the apparent concentration of the compound in tissue T defined as the amount of compound divided by the volume of the tissue times the tissue:plasma partition coefficient; ka is the absorption rate from the gut.
Network diagram for the VitD3 metabolic cascade model. The network contains compounds with measured concentrations and black arrows define conversion steps with kinetic equations as defined in Table 2. Note that concentrations of 1,24,25(OH)3D3 were not measured in this study.
Physiologic parameters were obtained from a previous publication (Davies and Morris, 1993). The slowly perfused compartment was composed of skin and muscle. The rapidly perfused compartment served as a mass balance compartment and comprised all tissue not previously named. Fractional blood flow rate for the slowly perfused compartment was equal to the sum of the rates of its components. The fractional blood flow rate of the rapidly perfused compartment was equal to the remaining fractional blood flow rates. Following the methods of Ramakrishnan et al. (2016), this study assumed no binding to red blood cells and scaled the blood flow to tissues by multiplying the blood flow QB and the blood volume VB by BP, the blood:plasma ratio, to model the rate of flow of plasma, as shown in eqs. 1 and 2:
The PBPK model was comprised of conventional mass balance equations, assuming well-stirred distribution into each model compartment with detailed equations for the PBPK model found in the Appendix (eqs. A1 to A14). Calculations for relevant partition coefficients and physiologic values used in this study are found in the supplemental materials (Supplemental Tables 1 and 2).
Metabolic Routes
The metabolic cascade for VitD3 considered in this study consists of three P450 enzymes: CYP2R1, CYP27B1, and CYP24A1. Although the P450 enzyme CYP2R1, a microsomal enzyme found mainly in the liver, is primarily responsible for the 25-hydroxylation of VitD3 to 25D3 following Michaelis-Menten kinetics (Cheng et al., 2004; Shinkyo et al., 2004), there are several other enzymatic pathways that may be involved in the 25-hydroxylation pathway, including CYP27A1 (Sakaki et al., 2005), CYP2J2, and CYP3A4 (Battault et al., 2013; Wang et al., 2013). This study chose to focus on the actions of CYP2R1, with the assumption that Vmax may vary depending on levels of sufficiency (Abramson, 1986). The variable rate of Vmax is given by eq. 3:
where V' is the maximum value of Vmax, [D] is the total liver concentration of VitD3, D50 is the concentration at which 50% of the inhibition occurs, and hV is the Hill coefficient for this function. The conversion of VitD3 to 25D3 by CYP2R1 is given by standard Michalis-Menten kinetics (eq. 4) (Ocampo-Pelland et al., 2016):
Endogenous production of VitD3 through incidental exposure to sunlight and additional dietary sources are incorporated into the model using eq. 4. This follows the assumption that the input rate for VitD3 is equal to the output rate of VitD3 to 25D3 when there is no external supplementation at low baseline concentrations of VitD3 (Ocampo-Pelland et al., 2016). Following the approach of Ramakrishnan et al. (2016), the baseline concentration of VitD3 and its metabolites in tissue were calculated using the relationship between measured total plasma concentrations and partition coefficients as CT,base = PT:p ·Cplasma,base, where CT,base is the baseline concentration in tissue T, PT:p is the calculated tissue:plasma partition coefficient, and Cplasma,base is the baseline concentration of the compound in the plasma compartment.
The metabolite 25D3 is further converted to either 1,25D3 by the 1α-hydroxylase CYP27B1, or deactivated to the metabolite 24,25D3 by the 24-hydroxylase CYP24A1. Both CYP27B1 and CYP24A1 are renal mitochondrial enzymes which have the capability to demonstrate Michaelis-Menten kinetics (Inouye and Sakaki, 2001; Sakaki et al., 2005). However, since observed concentrations of 25D3 range roughly 20–50 times smaller than literature values of Km for CYP27B1 (Inouye and Sakaki, 2001) and 200–400 times smaller than literature values of Km for CYP24A1 (Sakaki et al., 2005), this study assumes linear first-order kinetics for CYP27B1 and CYP24A1 acting on 25D3.
Supplementaldoses of VitD3 during the replenishment period were incorporated into the model through pulsing the initial condition for Agut, the amount of oral VitD3 effectively introduced in the intestine en route to the liver. During the initial and final portions of the study where no daily supplements were given, the only source of input to the model is through the endogenous production of VitD3 throughout sunlight and diet. For the purposes of continuity of the model, all participants were assumed to have 16 weeks of supplementation to achieve sufficient concentrations of 25D3.
Dynamic Adipose Partition Coefficients
Previous publications have suggested the possibility that partition coefficients in tissues, in particular the kidney, may change with levels of sufficiency (Quach et al., 2015) and the adipose tissue (Sawyer et al., 2017). We chose to incorporate dynamic partitioning for the adipose tissue compartment for VitD3 and 25D3 similar to (Sawyer et al., 2017) using the following equations (eq. 5 and 6):
Here, AdiPC[](t) are the adipose partition coefficients at any point in t related directly to the total plasma concentration of either VitD3 or 25D3 at time t (indicated by [D] and [25D], respectively), Adi[],50 is the concentration of VitD3 or 25D3 at which 50% of the inhibition occurs, and h[] is the Hill coefficient for these functions.
Stimulatory and Inhibitory Effects of 1,25D3
Circulating and tissue concentrations of 1,25D3 are tightly regulated through its synthesis by CYP27B1 and degradation by CYP24A1. CYP27B1 expression is tightly regulated by parathyroid hormone and 1,25D3, where high levels of 1,25D3 suppress expression of CYP27B1 (Schuster, 2011). Conversely, the expression of CYP24A1 is enhanced by the presence of 1,25D3; upregulation of CYP24A1 by 1,25D3 serves as feedback control to reduce concentrations of 1,25D3 (Ramakrishnan et al., 2016). To incorporate the regulatory effects of 1,25D3 plasma concentrations, this study follows the approach previously outlined (Ramakrishnan et al., 2016) and adjusts the fold change (FC; ratio of changed or basal mRNA levels) of CYP27B1 and CYP24A1 in proportion to total 1,25D3 plasma concentrations calculated from the partition coefficient in the relevant tissues. The rate of CYP27B1 fold change was expressed as shown in eq. 7:
with Imax the maximum inhibitory effect, IC50 the total plasma concentration of 1,25D3 in the kidney to achieve 50% of Imax, γ2 is the Hill coefficient for CYP27B1 in the kidney, and [1,25DK] is the total concentration of 1,25D3 in the kidney (Ramakrishnan et al., 2016). For CYP24A1, eq. 8 gives the fold change in tissue T, for T liver, kidney, intestine, or brain:
where Emax,T is the maximum stimulatory effect in tissue T, EC50,T is the total plasma concentration of 1,25D3 in the tissue to achieve 50% of Emax, γ1,T is the Hill coefficient for CYP24A1 in the tissue, and [1,25DT] is the total concentration of 1,25D3 in the tissue (Ramakrishnan et al., 2016). Following the methods of Noh et al. (2020) and Ramakrishnan et al. (2016), kin,[] and kout,[] are presumed to have the same value but different (appropriate) units for the differential equation.
Model Simulations and Parameter Estimations
This study follows the approach outlined by McNally et al. (2011) to explore the global sensitivity of the kinetic parameters for this PBPK model. The model equations were coded in MATLAB (version R2020a, The Mathworks Inc, Natick, MA) and sensitivity analysis was conducted using the SAFE Toolbox (Pianosi et al., 2015). For the initial analysis, the Morris’ Method was used to distinguish between the set of influential and noninfluential parameters through ranking (Morris, 1991). Parameters in this study were considered influential if the normalized mean sensitivity measure obtained using Morris’ Method was greater than 0.1 (Hsieh et al., 2018). Using this method, the set of influential parameters were chosen for optimization and all other parameters were held at available literature values. If a parameter did not have an available literature value, it was also included in the set chosen for optimization.
This set of parameters was optimized using nonlinear least squares where the cost of fit was calculated using the least squares difference between the observed and predicted model data. To avoid unfairly skewing the model, each compound’s data were transformed by subtracting the minimum observed value and dividing by the interquartile range of the observed data for that compound. Data were further weighted by the count of the number of available points for optimization in each phase; this ensured that data below the analytical level of quantification were not considered in and did not negatively affect the optimization routine. This was present primarily in the insufficient VitD3 data and comprised approximately 15% of the overall data points (3% of points without the insufficient VitD3 data) used for optimization. The lack of data caused by patient VitD3 concentrations below the limit of quantification poses some issue for initialization of the model for this compartment; however, the set of influential parameters span the entire timespan of the model, so we used visual inspection of the curve of the model prediction to assess fit in this area for this compartment. The model was initialized to the mean of the first two time points (t = 0 and t = 0.5 hours).
To generate 95% confidence intervals for the optimized parameters, a bootstrapping method was used after initial optimization of the model. For this method, parameters identified as candidates for optimization were allowed to vary twofold from their optimized values, which were based on available literature values, or the best result from multiple iterative guesses. Bootstrapping was performed by systematically running the optimization solver over randomly sampled parameters within these constraints, holding nonoptimized parameters constant at literature values. The resulting parameter output range for each optimized parameter was then subject to standard bootstrapping techniques and the 95% confidence intervals were found for each parameter. To generate a 95% confidence interval around the optimized model, parameters were randomly sampled in a normal distribution from a twofold range of their optimized value for 50 samples, the model run over these parameter sets, and then the 95th percentile of the total model outputs at every half-hour mark was used to generate the total band for the duration of the time-course data.
Assessment of Prediction Accuracy
The accuracy of the PBPK predictions was evaluated on predicted plasma concentration fit to observations for the three measured compounds. The predicted plasma concentration was formulated using the set of optimized parameters developed from the bootstrapping results. The goodness of prediction for each compound was based on the average fold error (AFE), root mean square error (RMSE), and the normalized root mean square error (NRMSE) (Sheiner and Beal, 1981). The fold error (eq. 9), AFE (eq. 10), RMSE, (eq. 11), and NRMSE (eq. 12) were calculated as follows:
where Pred is the predicted value, Obs is the observed value, IQR(Obs) is the interquartile range of the observed values, and n is the number of observed samples.
Results
A basic PBPK approach was applied to VitD3 and metabolites using intrinsic clearance parameters determined from the published literature and physiologic distribution parameters (Fig. 1). Data from eleven subjects (demographics given in Table 1) were used to fit the model; two subjects withdrew from the study before the conclusion of the experiment.
Because of limited VitD3 time course data, the nearly instantaneous conversion of VitD3 to 25D3 in VitD deficient individuals (Heaney et al., 2008), and the adipose tissue as a storage compartment that may contribute to circulating levels of VitD3 (Best et al., 2020), our model was insufficient at predicting available VitD3 time course data (results not shown). However, clinically, VitD3 levels are not measured and not likely to be important since therapy decisions are based on 25D3 levels. All optimized and literature model parameters are reported in Table 2. Bootstrapping was performed and optimized parameters are reported with 95% confidence intervals.
Kinetic parameters for the PBPK model
Optimized parameters are given with bootstrapped 95% confidence intervals and all other parameters are held to available literature values or estimates as notated. Note the 95% CIs are not necessarily symmetric about the optimized parameter due to the method of sampling. See the Model simulations and parameter estimations section for details.
The observed and final model predicted concentrations for 25D3, 24,25D3 and 1,25D3 during the insufficient and repleted periods are shown in Fig. 3, A–C, respectively, with a 95% confidence interval band around the optimized parameter model as described in section 7.5. Detailed model predictions for the insufficient and repleted periods are shown as insets in Fig. 3. The wavy line in the center of the plots indicate a break in the supplementation period for better clarity on the insufficient and repleted period.
Predicted PBPK plasma concentrations for VitD3 metabolites (25D3, 24,25D3, and 1,25D3). Box plots indicate the middle 50% of the data with the median data indicated by a line, and the shading indicates the 95% confidence interval of the model predictions. Insets indicate enhanced view of data for predicted PBPK plasma concentrations for VitD3 metabolites (25D3, 24,25D3, and 1,25D3), clipping out the supplementation period from days 20 through 92 of the experiment, where the wavy line in the center indicates the clipped portion of the model during the supplementation period where no data were observed. (A) 25D3 serum levels (nmol l−1). (B) 24,25D3 serum levels (nmol l−1). (C) 1,25D3 serum levels (nmol l−1).
Predictive performance of the PBPK model for the insufficient and repleted periods are shown in Fig. 4. The average fold error (AFE) for the VitD3 metabolites ranges between 1.01 and 1.05 and the normalized root mean square error (NRMSE) ranges between 0.26 and 0.38. AFE values close to one and NRMSE close to zero are indicative of better model fits to data. The correlation of determination (R2) for each metabolite is 0.77, 0.67, and 0.63, respectively, and 0.92 when considering all data points together. When considering the fold error values for the data, the model performed well with 97%, 88%, and 98% of the data for 25D3, 24,25D3, and 1,25D3, respectively, within twofold of unity (fold error values between 0.5 and 2.0). When considering a tighter range of fold error (1.5-fold of unity, with fold error values between 0.67 and 1.5), the model performs adequately with nearly 80% of 25D3 and 1,25D3 captured within a 1.5-fold error of unity. 24,25D3 data were captured less well, with a 1.5-fold error value of 63%.
Observed plasma concentrations versus predicted PBPK plasma concentrations in vitamin and three metabolites. Red x-marks indicate observations from the deficient period, blue circles indicate observations from the replete period, the dotted line represents the identity line, and the green lines represent the twofold interval.
Discussion
Vitamin D insufficiency is highly prevalent in the community, afflicting up to 80% of men and women in the United States, Canada, and Europe (Ganji et al., 2012; van Schoor and Lips, 2018). While the definitions of VitD insufficiency (and deficiency) are not standardized across medical organizations, targeted concentrations are generally in the 20 ng/ml to 40 ng/ml range for 25D, the primary metabolite used for classification (Holick et al., 2011; Ross et al., 2011; Ganji et al., 2012; van Schoor and Lips, 2018). Whereas various treatment recommendations have also been proposed according to the IOM (Jernigan and Andress, 2003) and ENDO (Holick et al., 2011), there is variability in success of achieving target concentrations of 25D3 in patients. There is currently not a clinically established approach to enable prediction of plasma 25D concentrations that might result from a given treatment regimen for a given patient. There is further complication since the plasma concentration of 25D will be impacted by the function of numerous P450 enzymes through activation and deactivation pathways, as well as through sunlight exposure, skin pigment, diet, and clothing. PBPK modeling has a potential to predict expected plasma concentrations of 25D and subsequent metabolites after administration of oral VitD therapy. The current study reports the development of a PBPK model of VitD3 and metabolites in VitD insufficient subjects in the United States who were treated with moderate doses (5000 I.U.) of daily VitD3 for up to 16 weeks to achieve replacement as defined by 25D ≥ 30 ng/ml. The comprehensive PBPK model for VitD3 incorporated dynamic adipose:plasma partition coefficients, atypical kinetics for CYP2R1 using a variable Vmax based on the total liver concentration of VitD3, and inhibitory and stimulatory effects of 1,25D3 for CYP27B1 and CYP24A1. The resultant model sufficiently predicted the concentrations of VitD3 metabolites well in healthy subjects throughout a multiphase study that incorporated an insufficient phase and replete phase. The model informs an understanding of the disposition of VitD3 metabolites and could be used to predict plasma concentrations that might result from a given dosing regimen of VitD3 in human patients.
The developed PBPK model for VitD3 and its metabolites consisted of nine compartments (adipose tissue, brain, heart, intestines, kidney, liver, rapidly perfused tissue, slowly perfused tissue, and plasma). Partition coefficients were either selected from the literature or calculated as described in the Supplemental Material. Although a previous publication in mice suggested that partition coefficients may differ in states of vitamin D deficiency (Quach et al., 2015), our simulations showed less than a 2% change in liver partition coefficients (as defined by the ratio of liver concentration to plasma concentration) between insufficient and sufficient states for VitD3 and metabolites, and less than 1% for kidney partition coefficients (data not shown). However, we did find a difference in model outcomes with the inclusion of variable partitioning in the adipose tissue compartment. As subjects in our study were initially at levels of insufficiency, there is evidence to suggest that normal daily inputs and endogenous production of vitamin D is not sufficient for accumulation in tissues (Heaney et al., 2009); however, at higher levels of sufficiency, the body is able to store VitD3 and 25D3 in adipose tissue for use at a later time (Mawer et al., 1972; Heaney et al., 2009; Abbas, 2017). The addition of a variable adipose partition coefficient for VitD and 25D, as described in section 7.2.2 and discussed in Sawyer et al. (2017), led to enhanced fits for the model. In addition, our model agrees well with previously published human studies without Kp values defined (Holick et al., 2008, Fig. 2B).
As VitD3 is acquired endogenously through the skin from the effects of UV-B light on 7-dehydrocholesterol (7-DHC) and dietary sources including dairy and fish, these contribute to baseline concentrations of VitD3 in the plasma. The PBPK model incorporated endogenous VitD3 levels with an assumption of output to 25D3. After oral ingestion and absorption of supplemental VitD3, the free fraction is taken up in the liver where it is hydroxylated by 25-hydroxylase to form calcifediol or 25D3. Subsequently, 25D3 is hydroxylated by the kidneys to form calcitriol or 1,25D3, which is considered the most active metabolite of VitD3. The current PBPK model used linear first order kinetics for CYP27B1 and CYP24A1 based on the significantly lower 25D3 concentrations versus published Km values for the respective enzymes. Supplemental dosing of cholecalciferol was incorporated into the model through pulsing the amount introduced into the intestines.
Circulating and tissue concentrations of 1,25D3 are tightly regulated through the synthesis by CYP27B1 and degradation by CYP24A1. Although CYP27B1 expression is tightly regulated by the binding of 1,25D3 and parathyroid hormone (Schuster, 2011), the expression of CYP24A1 is enhanced by the presence of 1,25D3; upregulation of CYP24A1 by 1,25D3 serves as a feedback control to reduce concentrations of 1,25D3 (Ramakrishnan et al., 2016). The current PBPK model incorporated the regulatory effects of 1,25D3 plasma concentrations by adjusting the fold change of CYP27B1 and CYP24A1 in proportion to 1,25D3 plasma concentrations in the liver, kidney, intestine, and brain, as previously identified (Ramakrishnan et al., 2016). Whereas the fold-change data were simulated and although the dosing regimen and the supplementation in Ramakrishnan et al. (2016) is different than our dosing regimen and supplementation, we observed the same general behavior for the fold change of each enzyme (results not shown), namely that as the level of 1,25D3 increases, the fold change of CYP27B1 decreases and the levels of CYP24A1 increase in each compartment. When the supplementation was removed at the repleted portion of the study, we observed an expected increase and decrease, respectively, in the fold changes of the CYP24A1 and CYP27B1 enzymes comparative to the concentration of 1,25D3 in the appropriate compartments. This behavior follows the expected directions for CYP27B1 and CYP24A1 given the critical concentration of 1,25D3 governing these physiologic effects on these enzymes.
Recent publications have discussed the advancement of atypical Michaelis-Menten kinetics for a variety of P450 enzymes, including several enzymes involved in the secondary metabolism cascade for vitamin D. This includes CYP3A4 (Arendse et al., 2013) and CYP2J2 (Leow and Chan, 2019; Leow et al., 2021). In particular, CYP2J2 shares 72.5% sequence similarity with CYP2R1, so it is not unreasonable to assume that CYP2R1 may also exhibit atypical Michaelis-Menten kinetics at varying sufficiency levels. To model this, we chose to use a modification on substrate inhibition kinetics, as shown in eq. 3. This led to strong improvement of model predictions, particularly in latter time points (data not shown).
Since the metabolites versus the parent VitD3 compound are measured clinically to assist clinicians in determining whether patients are VitD3 insufficient and/or whether they require changes to dosing regimens, we performed some simulations of predicted concentrations with common regimens. A simulation using the current model was performed inputting a dose of 1000 I.U./daily, which is consistent with dosing recommendations, according to the IOM. Our model predicted 25D3 plasma concentrations well (Fig. 5) and is consistent with a previous study with this dosing regimen (Holick et al., 2008, Fig. 2B). Note subjects in (Holick et al., 2008) began the study with sufficient levels of VitD; we hypothesis this explains the slight underprediction of our model to participant data at later time points as the ratio of VitD3 and metabolites may not remain the same across levels of sufficiency for initialization of the model.
Simulation of daily 1000 I.U. dosing regimen. Solid line indicates model predictions, data taken from (Holick et al., 2008) are mean ± SEM.
Overall, the developed PBPK model of VitD3 that incorporated data from insufficient participants into a repleted phase after daily administration of 5000 I.U. VitD3 for 12–16 weeks led to acceptable and satisfactory predictions in healthy human subjects. When considering the fold error values for the data, the model performed well with 97%, 88%, and 98% of the data for 25D3, 24,25D3, and 1,25D3, respectively, within twofold of unity (fold error values between 0.5 and 2.0). Because of limited VitD3 time course data, the nearly instantaneous conversion of VitD3 to 25D3 in VitD deficient individuals (Heaney et al., 2008), and the adipose tissue as a storage compartment that may contribute to circulating levels of VitD3 (Best et al., 2020), our model was insufficient at predicting available VitD3 time course data (results not shown). However, this may not be overly relevant since VitD3 levels are not measured or used clinically. An additional limitation to this model is the lack of available experimental data for comparison of atypical kinetic parameters for CYP2R1; however, the inclusion of the atypical kinetics greatly improved our model behavior and is within the realm of biologic possibility. The model serves to inform an understanding of VitD3 and metabolite disposition and could be used to predict 25D3 plasma concentrations that might result from a given dosing regimen in human patients. A beneficial future addition to this study would include incorporation of patient-specific data throughout the supplementation phase to better inform the model predictions. This model has applications in studying the effects of repletion schemes for populations of patients with impaired enzymatic abilities, such as chronic kidney disease, an important diseased population we are currently evaluating.
Acknowledgments
The authors thank University of Colorado and the University of Pittsburgh for their help in the conduct of the clinical studies. Jarrow Formulas, Inc. (Los Angeles, CA) is gratefully acknowledged for providing study drug (cholecalciferol 5000 I.U.).
Appendix
PBPK Model Equations
The parameters in the PBPK model equations refer to tissue (T), arterial blood flow (Q), compound amount (A), concentration (C), volume (V), and tissue:plasma partition coefficients 
. The concentration in each tissue is given by eq. A1, where the tissue:plasma partition coefficients are given in Supplemental Table 2.
For compartments that have the same equation form across all compounds, general mass balance differential equations are given in eqs. A2 and A3 based on Fig. 1 in the main text. Definitions and values for parameter values for eqs. A2 through A14 are given in Table 2 and Supplemental Tables 1 and 2.
Noneliminating tissue in the set Z (adipose, heart, rapidly perfused, slowly perfused)
Plasma:
Equations A4 through A7 describe the mass balance equations for the remaining four compartments (brain [Br], intestines [I], kidney [K], and liver [L]) for VitD3 and each metabolite, where additional components are described in eqs. A8 through A14.VitD3 equations:
25D3 equations:
24,25D3 equations:
1,25D3 equations:
Additional differential equations for the PBPK model are given in eqs. A8 to A14, where concentrations in tissue T are indicated by [MT], where M is VitD3 or a metabolite:Absorption of VitD3:
Endogenous production of VitD3:
Conversion of VitD3 to 25D3 in the liver (L):
Conversion of 25D3 to 1,25D3 in kidney by CYP27B1:
Conversion of 25D3 to 24,25D3 in tissue T by CYP24A1:
Clearance of 1,25D3 in tissue T by CYP24A1:
Clearance of 24,25D3 in the kidney:
The initial conditions for the model were generated by taking the mean of the data points to generate an initial estimate for the median plasma concentration.
Authorship Contributions
Participated in research design: Nolin, Joy.
Conducted experiments: Tuey, West, Nolin, Joy.
Performed data analysis: Sawyer.
Wrote or contributed to the writing of the manuscript: Sawyer, Tuey, Nolin, Joy.
Footnotes
- Received August 12, 2021.
- Accepted June 13, 2022.
This research was supported by the National Institutes of Health National Institute of General Medical Sciences [Grant R01-GM107122] (to M.S.J., T.D.N.).
No author has an actual or perceived conflict of interest with the contents of this article.
Primary laboratory of origin: M.S.J. and T.D.N.
↵
This article has supplemental material available at dmd.aspetjournals.org.
Abbreviations
- 25D3
- 25(OH)D3
- 1
- 25D3
- 1
- 25(OH)2D3
- 24
- 25D3
- 24
- 25(OH)2D3
- P450
- cytochrome P450
- VitD3
- vitamin D3
- Copyright © 2022 by The American Society for Pharmacology and Experimental Therapeutics
