Advertisement
Research ArticleArticle

Highlighting Vitamin D Receptor–Targeted Activities of 1α,25-Dihydroxyvitamin D3 in Mice via Physiologically Based Pharmacokinetic-Pharmacodynamic Modeling

Qi Joy Yang, Paola Bukuroshi, Holly P. Quach, Edwin C. Y. Chow and K. Sandy Pang
Drug Metabolism and Disposition January 2018, 46 (1) 75-87; DOI: https://doi.org/10.1124/dmd.117.077271

Abstract

We expanded our published physiologically based pharmacokinetic model (PBPK) on 1α,25-dihydroxyvitamin D3 [1,25(OH)2D3], ligand of the vitamin D receptor (VDR), to appraise VDR-mediated pharmacodynamics in mice. Since 1,25(OH)2D3 kinetics was best described by a segregated-flow intestinal model (SFM) that described a low/partial intestinal (blood/plasma) flow to enterocytes, with feedback regulation of its synthesis (Cyp27b1) and degradation (Cyp24a1) enzymes, this PBPK(SFM) model was expanded to describe the VDR-mediated changes (altered/basal mRNA expression) of target genes/responses with the indirect response model. We examined data on 1) renal Trpv5 (transient receptor potential cation channel, subfamily V member 5) and Trpv6 and intestinal Trpv6 (calcium channels) for calcium absorption; 2) liver 3-hydroxy-3-methyl-glutaryl-coenzyme A reductase (Hmgcr) and cytochrome 7α-hydroxylase (Cyp7a1) for cholesterol synthesis and degradation, respectively; and 3) renal and brain Mdr1 (multidrug-resistance protein that encodes the P-glycoprotein) for digoxin disposition after repetitive intraperitoneal doses of 120 pmol 1,25(OH)2D3. Fitting, performed with modeling software, yielded reasonable prediction of a dominant role of intestinal Trpv6 in calcium absorption, circadian rhythm that is characterized by simple cosine models for Hmgcr and Cyp7a1 on liver cholesterol, and brain and renal Mdr1 on tissue efflux of digoxin. Fitted parameters on the Emax, EC50, and turnover rate constants of VDR-target genes [zero-order production (kin) and first-order degradation (kout) rate constants] showed low coefficients of variation and acceptable median prediction errors (4.5%–40.6%). Sensitivity analyses showed that the Emax and EC50 values are key parameters that could influence the pharmacodynamic responses. In conclusion, the PBPK(SFM)-pharmacodynamic model successfully characterized VDR gene activation and serves as a useful tool to predict the therapeutic effects of 1,25(OH)2D3.

Introduction

Vitamin D is an essential hormone for health and diseases, and is regarded as important for the maintenance of bone and cellular homeostasis, longevity, and antiproliferation. Vitamin D deficiency is associated with bone disease, hyperparathyroidism, cardiovascular disease, hypertension, inflammation, diabetes, and cancer (Valdivielso et al., 2009; Sarkinen, 2011). Concentrations of the active vitamin D receptor (VDR) ligand 1α,25-dihydroxyvitamin D3 [1,25(OH)2D3] are tightly regulated by the rate-limiting, synthetic enzyme CYP27B1 or 1α-hydroxylase in kidney for activation of 25-hydroxyvitamin D3, the first and relatively inert intermediate that is formed from vitamin D in the liver (Shinki et al., 1992; Lemay et al., 1995; Henry, 2001). High concentrations of 1,25(OH)2D3 lead to induction of the degradation enzyme CYP24A1 or 24-hydroxylase, which can inactivate both 25-hydroxyvitamin D3 and 1,25(OH)2D3 (Makin et al., 1989; Jones et al., 1998). These feedback pathways prevent the accumulation of 1,25(OH)2D3 that would otherwise result in hypercalcemia due to excess reabsorption of calcium via the intestinal TRPV6 (transient receptor potential cation channel subfamily V member 6) (den Dekker et al., 2003) and renal TRPV5 and TRPV6. High calcium levels would result in indirect inhibition of the parathyroid hormone for feedback inhibitory control of CYP27B1 to curtail the synthesis of 1,25(OH)2D3 (Brenza and DeLuca, 2000; Turunen et al., 2007).

In addition to calcium and bone homeostasis, two new VDR targets have been identified in our laboratory. First, 1,25(OH)2D3 exerts a direct role on inhibition of the small heterodimer partner (Shp) that normally represses the transcription of liver cytochrome 7α-hydroxylase (Cyp7a1), the rate-limiting enzyme in cholesterol metabolism in hypercholesterolemic mice, when these were fed a high-fat/high-cholesterol diet (Chow et al., 2014). The mRNA/protein expression levels of the rate-limiting enzyme for cholesterol synthesis, 3-hydroxy-3-methyl-glutaryl-coenzyme A (Hmgcr), and Cyp7a1 exhibit circadian rhythm (Mayer, 1976; Noshiro et al., 1990). In rodents, these circadian oscillations result in maximum rates of cholesterol synthesis around midnight (Edwards et al., 1972; Ho, 1979) due to maximal Hmgcr and Cyp7a1 expression levels occurring at around 9 p.m. to 12 a.m. in C57BL/6 mice (Chow et al., 2014). Second, the expression of brain and renal Mdr1 (multidrug resistance protein 1) transcription was induced by the Vdr, and the resultant, elevated P-glycoprotein (P-gp) levels hastened the renal excretion and brain efflux of [3H]digoxin in mice when treated with repeated i.p. doses of 120 pmol or 2.5 μg·kg−1 1,25(OH)2D3 (Chow et al., 2011). Treatment of mice with this 1,25(OH)2D3 dose raised plasma calcium levels slightly due to elevated expression of the Trpv6 in the intestine and Trpv6 and Trpv5 the kidney (Chow et al., 2013). A mechanistically based pharmacokinetic-pharmacodynamic (PK-PD)–linked model that can quantitatively describe these events, however, is missing.

Compartmental PK-PD (Quach et al., 2015) and physiologically based pharmacokinetic (PBPK) (Ramakrishnan et al., 2016) models that have incorporated the induction of Cyp24a1 and inhibition of Cyp27b1 in the description of 1,25(OH)2D3 kinetic profiles after repeated i.v. (Quach et al., 2015) and i.p. (Ramakrishnan et al., 2016) doses in mice, respectively, have been compared. The PBPK model, comprising the kidney, intestine, liver, and brain, was found superior over the compartmental model. Improvement in goodness of fit was obtained, especially when the PBPK model contained the intestine nested as two tissular regions: the enterocyte region, which receives a low, partial intestinal blood flow wherein the enzymes and transporters are located, and an inert, serosal region. This intestinal model, known as the segregated flow model (SFM) (Cong et al., 2000; Ramakrishnan et al., 2016; Yang et al., 2016), was developed to better describe route-dependent intestinal metabolism. These findings greatly support the need to account for Vdr-mediated feedback regulation of the enzymes, Cyp24a1 and Cyp27b1, within multiple tissues to fully describe the kinetics of 1,25(OH)2D3.

In this study, we expanded the published PBPK(SFM) model (Ramakrishnan et al., 2016), which has been parameterized with physiologic constants (such as volume and blood flow), tissue partitioning coefficients (KT), and key tissue compartments comprising the intestine, liver, kidney, and brain, to describe the various Vdr-related dynamic activities in C57BL/6 mice toward 1) Trpv6 and Trpv5 for calcium absorption (Chow et al., 2013); 2) induction of Cyp7a1, which decreases liver cholesterol (Chow et al., 2013, 2014); and 3) increased renal and brain Mdr1/P-gp expression in the excretion/efflux of [3H]digoxin (Chow et al., 2011). The circadian rhythmic control of liver Hmgcr, Cyp7a1, and cholesterol levels was accounted for in this integrated pharmacodynamic model. The PBPK(SFM)-PD model was reasonable in describing the complex regulation of Vdr-related gene expression contributing to the various therapeutic outcomes.

Materials and Methods

The PBPK(SFM)-PD Model.

Previously, we have shown that the full PBPK(SFM) model was superior over the traditional PBPK intestinal model in describing the dose- and route-dependent kinetics of 1,25(OH)2D3 (Ramakrishnan et al., 2016). This published model incorporated Cyp24a1 (induction) and Cyp27b1 (inhibition) in their feedback regulation by the 1,25(OH)2D3-liganded Vdr, with use of indirect response equations that embellished the sigmoidal Emax, EC50, and Hill coefficients (see Fig. 1 for model scheme). The intestine, composed of two tissular regions (the serosa or inert region and the enterocyte region, where transporters and enzymes reside), are nested within the PBPK(SFM) model. With this model, 100% of the orally administered drug traverses the enterocyte region, but about 10%–20% of the intestinal flow brings blood-borne 1,25(OH)2D3 to the enterocyte region for intestinal processing. The model explains the lesser extent of intestinal removal of the systemically delivered drug versus drug delivered into the intestinal lumen, a phenomenon known as route-dependent intestinal drug metabolism (Cong et al., 2000).

Fig. 1.
Fig. 1.

Extension of the PBPK(SFM)-PD model of Ramakrishnan et al. (2016) to describe 1,25(OH)2D3-mediated induction of intestinal Trpv6 (A) and renal Trpv5 and Trpv6 (B), which mediate the absorption of calcium; induction of liver Cyp7a1 transcription (C), leading to the lowering of liver cholesterol; and renal (D) and brain (E) Mdr1, which facilitate the secretion of digoxin, in mice. The model also incorporates upregulation of its own degradation enzyme, Cyp24a1, in intestine, kidney, brain, and liver, as well as inhibition of its degradation enzyme, Cyp27b1, in kidney. Induction and inhibition are represented by the white and black boxes, respectively; kin and kout are zero-order synthesis and first-order degradation rate constants of enzyme, respectively; ka and kdeg are first-order absorption and luminal degradation rate constants for 1,25(OH)2D3, respectively, in the intestine compartment; the rhythmic changes of liver Hmgcr, Cyp7a1, and cholesterol are incorporated in the liver compartment; CLint,met is the metabolic intrinsic clearance for 1,25(OH)2D3; Embedded Imageis the secretory intrinsic clearance for digoxin; Q denotes plasma flow; Rsyn is the zero-order synthesis rate for 1,25(OH)2D3; τ is the time-delay term for liver Cyp7a1 transcription. Subscripts I, K, Br, and L represent intestine, kidney, brain, and liver tissue, respectively.

Data.

Data sets (Chow et al., 2011, 2013, 2014) that were used for pharmacodynamic modeling originated in studies whereby male C57BL/6 mice were fed a normal diet and treated with corn oil (vehicle for baseline levels, control group) or 2.5 μg·kg−1 1,25(OH)2D3 [(1R,3S,5Z)-4-methylene-5-[(2E)-2-[(1R,3aS,7aS)-octahydro-1-[(1R)-5-hydroxy-1,5-dimethylhexyl]-7a-methyl-4H-inden-4-ylidene]ethylidene]-1,3-cyclohexanediol (The Merck Index Online: https://www.rsc.org/merck-index) in corn oil–treated group] i.p. at 9 a.m. every other day for 8 days (four doses). Pharmacodynamic effects on 1) plasma calcium levels, assayed by inductively coupled plasma atomic emission spectroscopy, and the mRNA expression levels of (intestine and kidney) Trpv6, determined with quantitative polymerase chain reaction (qPCR) (Chow et al., 2013), and 2) the mRNA expression levels of liver Cyp7a1 (baseline and treated) and liver cholesterol levels, assayed by the extraction method of Folch et al. (1957), were obtained from the studies by Chow et al. (2013, 2014). Last, for Mdr1, data from a study involving [3H]digoxin levels and excretion (normalized by the injected, i.v. dose) in murine kidney and brain were used (Chow et al., 2011). These data defined the mRNA expression levels of Mdr1 (kidney and brain). [3H]Digoxin [(3β,5β,12β)-3-[(O-2,6-dideoxy-β-d-ribo-hexopyranosyl-(1→4)-O-2,6-dideoxy-β-d-ribo-hexopyranosyl-(1→4)-2,6-dideoxy-β-d-ribo-hexopyranosyl)oxy]-12,14-dihydroxycard-20(22)-enolide (The Merk Index Online: https://www.rsc.org/merck-index)], which is not metabolized in mice, was given as an intravenous bolus dose (0.1 mg·kg−1) at 24 hours after the cessation of 1,25(OH)2D3 treatment.

Parameters and 1,25(OH)2D3 levels for modeling were identical to those used previously for PBPK(SFM) modeling (Ramakrishnan et al., 2016). The assigned and fitted constants used are summarized in Table 1. Since calcium absorption also occurs with renal Trpv5, in addition to Trpv6 in the intestine and kidney, we took the control (for baseline or untreated levels) and 1,25(OH)2D3-treated samples from the study by Chow et al. (2013) and assayed for the relative mRNA expression levels of renal Trpv5. Although the Vdr-mediated cholesterol-lowering effect is usually examined in our laboratory among mice fed the high-fat/high-cholesterol diet (Chow et al., 2014), comparable temporal data for mice under the Western diet were not available. Hence, we restricted our pharmacodynamic analysis to mice that were fed the normal diet only. We also determined, from the liver samples of Chow et al. (2013), the temporal profiles of liver cholesterol and relative mRNA expression levels of Hmgcr by qPCR in control and treated mice. Additional samples that were collected at 6 hours after the last i.p. injection (after fourth dose) from another treatment study [same vehicle and same i.p. 1,25(OH)2D3 doses] were assayed for liver Hmgcr, Cyp7a1 mRNA expression, and cholesterol levels (Bukuroshi and Pang, unpublished study) and added to the data pool.

View this table:
TABLE 1

Physiologic volumes, blood flows, and pharmacokinetic parameters used for continued, PD modeling and simulation of mice data [see Ramakrishnan et al. (2016) for details]

Liver Cholesterol Assay.

Cholesterol levels of C57BL/6 mice given corn oil (vehicle) or 2.5 μg·kg−1 1,25(OH)2D3 [samples from Chow et al, (2013)] were measured according to published methods (Folch et al., 1957; Cho, 1983), as described earlier (Chow et al., 2014). Around 0.1–0.2 g of liver tissue was homogenized in 4 ml of chloroform:methanol (2:1, v/v) mixture to extract the lipids. One milliliter of 50 mM NaCl was added to the homogenate and then centrifuged, and the organic phase was removed and washed with 1 ml of 0.36 M CaCl2/methanol. With repeated centrifugation and removal of the organic phase, the washing was repeated with another 1 ml of 0.36 M CaCl2/methanol. The final, organic phase was placed in a volumetric flask and made up to a total volume of 5 ml with chloroform. A 100μl aliquot of this chloroform solution was removed onto a glass tube, and 10 μl of chloroform:Triton X-100 (1:1, v/v) was added and air dried overnight; the same was repeated for 10 μl of the standards. A colorimetric enzymatic assay for liver cholesterol was then performed using commercial reagents based on the manufacturer’s protocol (Infinity, Ca#TR13421; Thermo Scientific, Mississauga, ON, Canada).

Real-Time qPCR.

The relative mRNA expression of renal Trpv5 and liver Hmgcr mRNA over the course of 8 days of treatment with corn oil or 1,25(OH)2D3 from the study by Chow et al. (2013) was determined. The primer sequence for Trpv5 was as follows: forward, 5′ CATGATGGGCGACACTCACT 3′ and reverse, 5′ GGTGGTGTTCAACCCGTAAGA 3′. The primer sequence for Hmgcr was as follows: forward, 5′ CAAGGAGCATGCAAAGACAA 3′ and reverse, 5′ GCCATCACAGTGCCACATAC 3′ (Chow et al., 2014). Liver Hmgcr and Cyp7a1 mRNA expression and cholesterol levels at 6 hours after the last dose of the same 1,25(OH)2D3 dosing regimen (Bukuroshi and Pang, unpublished study) were also assayed. The total mRNA was extracted from liver tissue using the TRIzol extraction method according to the manufacturer’s protocol (Sigma-Aldrich, St. Louis, MO) with modification (Chow et al., 2011, 2013, 2014). The mRNA data were normalized to that of cyclophilin (Chow et al., 2013).

Data Fitting and Simulations: Application of PBPK(SFM)-PD Model to Describe 1,25(OH)2D3 Pharmacodynamics.

The published PBPK(SFM) model (Ramakrishnan et al., 2016) was used as a template to expand upon for pharmacodynamic modeling. As described in earlier PBPK modeling, the intrinsic clearance estimates (for influx or efflux or metabolism) were expressed as the product of the unbound fraction in plasma or tissue (fP or fT) × intrinsic clearance. The plasma protein binding of 1,25(OH)2D3 to the vitamin D binding protein was expected to be linear, inasmuch as an excess of vitamin D binding protein (5.1–9.7 μM) present in serum (Faict et al., 1986); the tissue unbound fractions of 1,25(OH)2D3, however, remained largely undetermined. Changes in the relative expression of Vdr-target genes were expressed as fold-change (FC) in mRNA expressions, obtained upon normalization of the treated by the control or nontreated data (Ramakrishnan et al., 2016). The initial (at day 0) mRNA FC was set as unity. The fitted parameters from the published model were then related to the temporal plasma and tissue 1,25(OH)2D3 profiles, as well as Cyp24a1 and Cyp27b1 relative mRNA expression changes in C57BL/6 mice (see Table 1). These estimates were adopted and fixed as constants in this extended PBPK(SFM)-PD model.

The indirect response model (Sharma and Jusko, 1998) was applied to fit these pharmacodynamics responses (R) using the naïve pooled data maximum likelihood method in ADAPT5 (Biomedical Simulations Resource, University of South California, Los Angeles, CA) (see Fig. 1 for the model scheme, and detailed events among the individual tissues in subsequent figures). The model contains the kin (zero-order production rate constant) and kout (first-order degradation rate constant) to describe the response (R) time profile of Vdr-target genes. Here, the rate of change of the response is given as the synthesis rate minus the degradation rateEmbedded Image(1)where the stimulatory function [E(t)]Embedded Image(2)and inhibitory function [I(t)]Embedded Image(3)are incorporated to denote the zero-order synthesis (kin) and first-order degradation (kout) rate constants of the Vdr target that would elicit the pharmacodynamic response (see detailed mass balance equations in Appendix). Notably, kin = kout·R0, where R0 is the initial response or average value of control samples; kin has the same value as kout when R0 = 1 in describing the mRNA fold-change. The units, however, would differ. For plasma calcium and liver cholesterol levels, values of R0 equal the average basal level (see Appendix for details). The Emax (or Imax), EC50 (or IC50), and γ denote the maximum stimulatory (inhibitory) effect, the tissue 1,25(OH)2D3 concentrations whereby 50% Emax (or Imax) is achieved, and the Hill coefficient, respectively. For simplicity, γ was set to 1.

Fitting.

For the modeling of Trpv5 and Trpv6 FC, fitting was performed with FC values from all Trpv5 and Trpv6 mRNA expression levels to determine the extent of calcium absorption. The average of the calcium concentrations for the control samples was set as R0, which is fixed at 10.7 mg·dl−1. We then expressed kout as kin/R0 to reduce the number of parameters (see Appendix for additional details). Then simulation was performed to show the relative contribution of Trpv5 or Trpv6 in calcium absorption.

The stimulatory approach was applied to examine renal and brain efflux of digoxin due to Mdr1 FC upon induction by the Vdr. For describing renal and brain digoxin PD effects, we transformed the digoxin tissue concentration data in kidney and brain data after treatment and took the difference (treated value − control value) divided by the control value (R0) × 100% as the percent reduction (%reduced). Hence, the data of Chow et al. (2011) on digoxin, with/without 1,25(OH)2D3 treatment, in the kidney and brain were recalculated in this fashion to reflect the effect of 1,25(OH)2D3 treatment on reducing digoxin accumulation in the brain and kidney.

For the modeling of cholesterol lowering, we applied a simple cosine model (D’Argenio et al., 2009) for calculating the mean liver Hmgcr and Cyp7a1 endogenous synthesis rates that display circadian rhythm (Edwards et al., 1972; Shefer et al., 1972; Mayer, 1976; Kai et al., 1995; Panda et al., 2002; Aoyama et al., 2010). The equations used for describing the circadian baseline are not harmonic (or multiple cosine) equations, as described by D’Argenio et al. (2009) (please see page 283 of the ADAPT5 User’s Guide (D’Argenio et al., 2009), with contributions from Dr. Wojciech Krzyzanski). The kin values for Hmgcr and Cyp7a1 were expressed by the simple cosine model (see eqs. A10 and A13 in Appendix), as reported by Chakraborty et al. (1999); however, the Rmean value in our model [equivalent to the Rm of Chakraborty et al. (1999)] for the mean synthesis rate was re-parameterized in our ADAPT5 fit, according to Krzyzanski (see Appendix, eqs. A9 and A12 on how we defined Rmean for Hmgcr and Cyp7a1 using the ADAPT5 User’s Guide). Similarly, other researchers have also applied these Rmean equations, as we did, to describe circadian rhythmic changes of plasma mevalonic acid following rosuvastatin treatment (Aoyama et al., 2010). We described diurnal variation for both Hmgcr and Cyp7a1 (Fig. 1, see Appendix), and a series of transit compartments and time-delay functions were added to improve the fit of Cyp7a1 mRNA FC to explain the Cyp7a1-mediated induction of cholesterol metabolism. Model complexity can be readily appreciated even when one input parameter was modified with circadian rhythm. Therefore, the circadian variation in other genes that regulate cholesterol turnover (such as farnesoid X receptor (Fxr) and small heterodimer partner (Shp)) for the mouse was not considered in the extended PBPK(SFM)-PD model for sake of simplicity and avoidance of over-parameterization, although these genes would influence the kin for cholesterol. For Hmgcr mRNA expression levels, which were found unchanged with 1,25(OH)2D3 treatment (Chow et al., 2014), the modeling of Hmgcr on cholesterol synthesis was simplified, and Embedded Image was assumed to be influenced only by Hmgcr FC values. For Cyp7a1, input and output rate constants (Embedded Imageand Embedded Image) were used to describe Cyp7a1 mRNA turnover, which would increase the degradation of cholesterol Embedded Image.

The variance model was modeled withEmbedded Image(4)where Y(t) is the model output; var(t) is the variance function associated with the output; and σinter and σslope denote the two variance parameters describing a linear relationship between the S.D. of the model output [Y(t)]. The prediction capacity was validated by calculating the median value of absolute percentage prediction error, or %PE:Embedded Image(5)For examination of the validity of the extended PBPK(SFM)-PD model, a visual predictive check (Cox et al., 1999; Duffull and Aarons, 2000; Yano et al., 2001; Westerhout et al., 2012) was generated after simulating the data 1000 times with the final fitted parameter estimates (Table 2). The 5th and 95th percentile of the predicted concentrations were calculated to obtain the 90% prediction interval, which reflects the precision of parameter estimates.

View this table:
TABLE 2 

Fitted parameters [estimate and (coefficient of variation; CV%)] obtained from the extended PBPK(SFM)-PD model to include the pharmacodynamic responses

Sensitivity Analysis.

The sensitivity analysis of model outputs (pharmacodynamic responses) was evaluated by calculating the percentage change of the area under the curve (AUC) of model outputs after increasing or decreasing the model parameters by 2-fold:Embedded Image(6)where AUCsim is the AUC of the pharmacological response calculated using fitted parameters, and AUC±2-fold was obtained after altering the fitted parameters by 2-fold (Emond et al., 2006; Urva et al., 2010).

Results

PBPK(SFM)-PD Model to Describe Vdr-Regulated Pharmacodynamic Responses

The PBPK(SFM) model (Ramakrishnan et al., 2016) was expanded upon, incorporating additional data on the mRNA relative expression of Vdr-target genes as well as other changes in plasma calcium and liver cholesterol levels. The assigned (Table 1) and optimized (Table 2) parameters for all of the pharmacodynamic effects from fitting are summarized here. The coefficients of variation for the estimated parameters were mostly within the accepted range. The fitted values of EC50 for Vdr-target gene expressions (such as intestinal Trpv6 and renal Trpv5 and Trpv6, renal and brain Mdr1, and hepatic Cyp7a1) are around 11.4–1,135 pM (Table 2), whereas the average plasma is 11,360 pM (Chow et al., 2013). All of the estimated EC50 values are lower than the average plasma and the relevant tissue 1,25(OH)2D3 concentrations, suggesting these enzymes are likely saturated by the high concentrations attained with the dosing regimen (120 pmol i.p. repeated every other day for four doses). Moreover, the EC50 estimates for Trpv6I, Trpv5K, and Trpv6K are similar (136, 163, and 104 pmol·kg−1 tissue), whereas the Emax is the highest for Trpv6I > Trpv5K > Trpv6K [361, 11.7, and 2.9 (FC)], a ranking that is consistent with observations on the relative abundance of Trpv6 protein expression (Chow et al., 2013).

The model overpredicted the induction of intestinal Trpv6 transcription for the first two doses, but underpredicted the FC changes that were dramatically higher for the third and fourth doses, failing to describe the systematic changes of the intestinal Trpv6 (Fig. 2A, linear scale). We further plotted the intestinal Trpv6 data on a semilog scale, and found that the plot overemphasized deviations for lower values and de-emphasized the underpredictions in the third and fourth doses (data not shown). We suspect that the systematic trend is due to receptor sensitization with the hysteresis loop plots [FC vs. tissue concentration of 1,25(OH)2D3]. Similarly, the model successfully predicted the induction of renal Trpv5 and Trpv6 mRNA FCs for the first two doses, and again slightly underpredicted the greater changes for the latter two doses. The median values of the prediction error remained, however, acceptable (<50%) (Sager et al., 2015), ranging from 4.5% to 40.6% (Table 3) amidst the trends. The observed and predicted data from the PBPK(SFM)-PD model and the fitted responses over time following multiple i.p. doses of 1,25(OH)2D3 are summarized in Figs. 2–4.

Fig. 2.
Fig. 2.

The intestine and kidney compartments were used to describe intestinal Trpv6 mRNA (A) and renal Trpv5 and Trpv6 mRNA fold changes (B), and increased plasma calcium levels (mg·dl−1) (C) after repeated 1,25(OH)2D3 i.p. administration. The fits to the fold changes in mRNA expression of intestinal Trpv6 (A), renal Trpv6 and Trpv5 (B), and plasma calcium concentrations (C) are shown. We further simulated the profiles for calcium absorption (calcium concentration) due to intestinal Trpv6 (solid red line) only, renal Trpv5 (dotted red line) only, and renal Trpv6 (dashed red line) only when the other calcium channels were absent (C). The circles represent the data of Chow et al. (2013); renal Trpv5 mRNA expression was assayed by qPCR. In all the figures, the solid line represents the fitted line, and the dashed lines represents the 90% confidence interval of prediction generated based on 1000 simulations with parameters summarized in Table 2.

Fig. 3.
Fig. 3.

The liver compartment (A) describing the basal, temporal liver Hmgcr and Cyp7a1 mRNA expression (expressed as FCs, normalized by the 9 a.m. value) and liver cholesterol concentrations in mice fed a normal diet without treatment (B), and Hmgcr and Cyp7a1 mRNA and cholesterol levels (C) after repeated 1,25(OH)2D3 i.p. administration. The insets in (B) show the observed FCs between 9 a.m. and 9 p.m., with solid lines and dots representing fitted profiles and observations of Hmgcr and Cyp7a1 mRNA and cholesterol (without treatment), showing that both exhibit circadian rhythms with a 24-hour span. The fitted parameters based on the first 24 hours were used to simulate profiles of the next 8 days, and were presented against observations. (C) Treatment with 1,25(OH)2D3 showed relatively no change in (measured) Hmgcr but higher Cyp7a1 mRNA expression levels. The highlighted black solid circles represent the additional samples collected at 6 hours after the fourth dose of a treatment study (Bukuroshi and Pang, unpublished study). In (B) and (C), the dashed line represents the 90% confidence interval of prediction generated based on 1000 simulations with parameters summarized in Table 2. The solid, purple lines are simulated profiles of Hmgcr, Cyp7a1, and cholesterol in liver when there was no circadian rhythm, either for basal levels (B) or those after 1,25(OH)2D3 treatment (C).

Fig. 4.
Fig. 4.

(A) The expanded kidney compartment for describing temporal renal Mdr1 mRNA FCs and %reduction of [3H]digoxin accumulation in kidney tissue. (B) The expanded brain compartment for describing temporal brain Mdr1 mRNA fold changes and %reduction of digoxin accumulation in brain tissue (%reduction = (treated - control)/control·100%) after repeated 1,25(OH)2D3 i.p. administration. The solid lines and dots represent fitted and observed temporal data. The dashed lines represent the 90% confidence interval of prediction generated based on 1000 simulations with parameters summarized in Table 2 [data of Chow et al. (2011)].

View this table:
TABLE 3

Calculated median %PE for the PBPK(SFM)-PD model

Trpv5 and Trpv6 for Calcium Absorption.

When individual components (the intestinal Trpv6, kidney Trpv5, and kidney Trpv6 mRNA expression) were simulated for their contributions to the absorbed calcium concentration, intestinal Trpv6, and not renal Trpv5 and Trpv6 (see Appendix, eqs. A1–A4) (Sharma and Jusko, 1998), was found responsible for the bulk of the absorbed calcium concentration (Fig. 2C, see red lines), which was elevated 10%–40% compared with the basal level (∼10 mg·dl−1) (Fig. 2C). This is reasonable since the Emax for Trpv6I mRNA is the highest and >> Trpv5K > Trpv6K mRNA (Table 2). Results from the visual predictive check suggest that the extended PBPK(SFM)-PD model described the majority of the observed data reasonably well within the 90% prediction interval, which reflects the precision of parameter estimation (Fig. 2, dashed lines and shaded area).

The best fit for calcium absorption occurred when all intestinal Trpv6 and renal Trpv5/6 were involved in calcium absorption. When the PBPK(SFM)-PD model was further tested on the relative importance of intestinal Trpv6 versus renal Trpv6 or Trpv5 as the only calcium channel contributing to calcium absorption (see Appendix, eq. A7), the fit based on intestinal Trpv6 alone was better statistically compared with those for renal Trpv6 or renal Trpv5 alone (fit not shown), as revealed by the prediction error (%PE) and weighted residual sum of squares (Table 4). This was substantiated by the relative abundance of protein expression level of intestinal Trpv6 versus that for renal Trpv6 (Chow et al., 2013).

View this table:
TABLE 4

Summary of fitted results of the PBPK(SFM)-PD model when only intestinal Trpv6, only renal Trpv6, or only renal Trpv5 is present versus when all Trpv5 and Trpv6 channels of intestine and kidney are present

Cyp7a1 and Liver Cholesterol.

In similar fashion, the model was extended to describe the upregulation of Cyp7a1 that lowered cholesterol due to the inhibition of Shp in the liver. Since Shp is unstable (Miao et al., 2009), Cyp7a1 and not Shp mRNA FC was used to illustrate the pharmacodynamic changes. For the modeling of Vdr-mediated cholesterol-lowering effects, the description of cholesterol synthesis was simplified. We used Embedded Image to incorporate the influence of Hmgcr. The expression of Hmgcr was previously found not to be significantly changed upon 1,25(OH)2D3 treatment, and therefore, stimulatory function (with Emax and EC50) was not applied to describe its turnover (Chow et al., 2014). By contrast, Cyp7a1 was affected by 1,25(OH)2D3 treatment, and the parameterEmbedded Image would be modified due to Cyp7a1 upregulation with 1,25(OH)2D3 repeated i.p. dosing (Fig. 3) (Chow et al., 2014) (see Appendix).

The liver compartment (Fig. 3A) was used to model Vdr-mediated induction of Cyp7a1 expression, which, in turn, would lead to changes in liver cholesterol in mice in vivo. Over a 24-hour span, basal levels of liver Hmgcr and Cyp7a1 mRNA expression, which display circadian rhythm, showed that peak levels were achieved at around 9 p.m., and basal levels of cholesterol rose and fell in parallel (Fig. 3B, see insets). The fits over the first 24 hours described circadian rhythm for the basal, Hmgcr and Cyp7a1 mRNA expression levels (Fig. 3B, see insets) and the resultant parameters (Table 2) were used to simulate the time courses of Hmgcr and Cyp7a1 over the 8 days for the control mice (Appendix, eqs. A8–13); these simulated profiles were used for comparison with the fitted profiles on the treated mice. When compared with controls, treated mice showed unchanged Hmgcr expression (Fig. 3B vs. 3C) but higher Cyp7a1 mRNA expression, with inductive effects being more prominent during Days 2–8 rather than on the first day (Fig. 3C). The fitted, liver cholesterol temporal profile was well described, although there was high variability of the observed data. Although high intersubject variability was observed for all of the doses, the fit was well captured and mostly fell within the 90% prediction interval (Fig. 3C). When circadian rhythm was absent, the basal levels of Hmgcr, Cyp7a1, and cholesterol remained time-invariant (Fig. 3B, purple lines). Statistically, the fitted data lacking circadian rhythm contained higher weighted residual sum of squares and PE% and were not as well described by the model that considered circadian rhythm (Fig. 3, B and C; Table 5).

View this table:
TABLE 5

Summary of fitted results for liver Cyp7a1 and cholesterol when circadian rhythm is present or absent

Mdr1/P-gp and Digoxin.

The brain compartment was applied to examine the upregulation of Mdr1 expression for increased clearance of digoxin, a P-gp substrate (Chow et al., 2011), in the kidney (Fig. 4A) and brain (Fig. 4B). Elevated Mdr1 expression was observed only in the kidney and brain. Although the Vdr mRNA expression is high in the intestine (Chow et al., 2013), we failed to observe Mdr1 induction in the intestine (Chow et al., 2011, 2013), and hence restricted our analysis to the kidney and brain only. The model slightly overpredicted renal Mdr1 expression for the first two doses but was able to describe upregulation of Mdr1 mRNA in the third and fourth doses. The model also overpredicted the brain Mdr1 mRNA expression with a median prediction error of 11.5% (Table 3). The %reduction of digoxin accumulation in kidney and brain tissue fell within the 90% prediction interval (denoted by dashed lines and shaded areas), and model prediction aptly captured the observations (Fig. 4).

Sensitivity Analysis

For determination of which pharmacodynamic parameter (such as Emax, EC50, and turnover rate constants) is relatively more sensitive, a sensitivity analysis was conducted. The results (Fig. 5) show that the Emax and EC50 values for intestinal Trpv6, renal Trpv5 and Trpv6, and liver Cyp7a1 were key parameters influencing the extent of induction of their transcription after repeated i.p. doses. The FCs of Trpv5 and Trpv6 were relatively insensitive to the corresponding turnover rate constants (kin and kout). On the other hand, renal and brain Mdr1 gene levels were found to be sensitive to changes in Emax and turnover rate constants, but less sensitive to the changes in EC50. Additionally, the FCs of liver Hmgcr and Cyp7a1 were relatively sensitive to the peak mRNA synthesis rate (Ramp). The increase and decrease in Tpeak values (or the time to reach the peak mRNA synthesis rate) were associated with reduced AUC of Cyp7a1 and Hmgcr mRNA turnover when considering circadian rhythm. Values of the liver Hmgcr, cholesterol expression levels, and plasma calcium concentrations were less sensitive to their corresponding turnover rate constants. Overall, the pharmacodynamic responses were relatively insensitive to the changes of the Hill coefficient associated with upregulation of Trpv5/6 and Mdr1 transcription (data not shown), and therefore, the Hill coefficient was set to 1 for simulation of the mRNA expression.

Fig. 5.
Fig. 5.

Sensitivity analyses. The plot represents the percentage change (%change) in pharmacodynamic changes (represented by AUC from data of Figs. 2–4) following a 2-fold increase or decrease of the pharmacodynamic parameters (Table 2). Small changes in the model output suggest that the corresponding pharmacodynamic parameter is less sensitive to the adjusted parameters. Turnover rate constants (kin and kout) are represented by a single term (k) for simplicity since the values of kin and kout are identical, although the units differed (kin = kout·FCbaseline, where FCbaseline = 1).

Discussion

The development of pharmacokinetic-pharmacodynamic–linked models has gained more and more utility in the determination of drug utilization during drug development, especially in relating to changes in relative mRNA or protein expression (Jin and Jusko, 2007). PBPK modeling can provide mechanistic insight as it interrelates physiologically meaningful parameters and dynamic responses (e.g., Emax and EC50). With our extensive database on 1,25(OH)2D3, especially on temporal changes with 1,25(OH)2D3 treatment (Chow et al., 2013), we first established a sound model that relates tissue 1,25(OH)2D3 concentrations to the synthetic and degradative enzymes, Cyp27b1 and Cyp24a1, which tightly regulate levels of 1,25(OH)2D3, wherein sigmoidal Emax/Imax equations of the indirect response model were used to account for the concentration-dependent pharmacodynamic behavior of 1,25(OH)2D3 (Ramakrishnan et al., 2016). Since we had also obtained mRNA expression of enzymes and transporters of interest within the rich tissue data sets and information about the corresponding Vdr-mediated pharmacological effects (Chow et al., 2011, 2013, 2014), the present work was undertaken to extend our PBPK(SFM)-PD model for describing turnover steps, time-delay factors due to the intermediate signaling transduction processes, and circadian rhythm characteristics.

The model-fitted profiles agreed reasonably well with our measurements. Calcium absorption (Embedded Image) was modified by factors related to increased relative mRNA expression of intestinal and renal Trpv6 mRNA (Trpv6I,FC and Trpv6K,FC) as well as renal Trpv5 mRNA (Trpv5K,FC) (Fig. 2). Although the elevated plasma calcium concentration can in turn activate the calcium-sensing receptor, which inhibits calcium transport (Blankenship et al., 2001), and decreases the production of the parathyroid hormone that normally increases the expression of Cyp27b1 (Jones et al., 1998; Abraham et al., 2009, 2011), our simplified model is able to provide a reasonably good prediction of the temporal plasma calcium changes. Accordingly, plasma calcium levels are attributed more to intestinal Trpv6 and its level of induction than to renal Trpv6 or Trpv5 (Fig. 2C). However, our model failed to capture the greater increase in the intestinal Trpv6 expression in the latter two doses. To understand whether these patterns might be due to receptor sensitization, we examined hysteresis loops of the effect-concentration curve (data not shown). Counterclockwise hysteresis loops were found for the Trpv6 FC, and a greater extent of induction occurred with the third and fourth doses compared with the first two doses. The threshold and magnitude of response accompanying receptor activation could also be increased when sensitization develops (Gold and Gebhart, 2010). But 1,25(OH)2D3 treatment was found to upregulate the mRNA expression of Vdr only by 1- to 2-fold (Chow et al., 2013), which in turn can further increase Trpv6. Vdr activation may not be a driving force since the intestine tissue is Vdr-rich, and the extent of Vdr upregulation (1- to 2-fold) may not be significant enough to elicit the ∼600-fold change in intestinal Trpv6. Other possibilities are the delayed pharmacological response or the formation of other plausible active metabolite of 1,25(OH)2D3 that can bind to the Vdr and elicit the gene changes (Louizos et al., 2014). However, incorporation of turnover delay functions for Trpv6 for delayed effects (described by multiple transient compartments) failed to characterize the sudden increase in gene expression for the latter two doses, and no major active metabolite of 1,25(OH)2D3 was described in the literature. As our PBPK(SFM)-PD model is already quite complicated, the Vdr and its occupancy by 1,25(OH)2D3 for the upregulation processes have not been considered. The mechanism for this systematic change for the latter doses is currently unknown.

Previous studies have shown diurnal variations in liver cholesterol biosynthesis and metabolism (Back et al., 1969; Kandutsch and Saucier, 1969; Edwards et al., 1972; Mayer, 1976). For the modeling of cholesterol, the kinetic profiles of Hmgcr and Cyp7a1, the rate-limiting enzymes for cholesterol synthesis (Shapiro and Rodwell, 1969; Edwards et al., 1972) and metabolism (Gielen et al., 1975; Noshiro et al., 1990), respectively, need to be considered. These key hepatic, circadian genes are regulated by the endogenous glucocorticoids which are under control of the hypothalamus-pituitary-adrenal axis (Van Cantfort and Gielen, 1979; Oishi et al., 2005). Several PK-PD models have successfully considered circadian rhythm for the input rate with complex harmonic functions using an indirect response model (Yao et al., 2006; Hazra et al., 2007; Aoyama et al., 2010; Scheff et al., 2010), namely, for endogenous corticosterone production and the transcription of the hepatic glucocorticoid receptor (Yao et al., 2006; Hazra et al., 2007) and cytokines and cyclic production of hormones (cortisol and melatonin). Model complexity can be readily appreciated when more than one input parameter modify the circadian rhythm, especially with other genes that regulate cholesterol turnover [such as farnesoid X receptor, SHP, and FGF19 (Bookout et al., 2006; Lundåsen et al., 2006; Yang et al., 2006)]. Although the “two rates” model has been used to represent a high zero-order production of cortisol in the morning and low production rate for the rest of the day (Scheff et al., 2010), we used the simple cosine model (D’Argenio et al., 2009) for Hmgcr and Cyp7a1, whose circadian patterns influence cholesterol levels with a peak around 9 p.m. to 12 a.m. and a nadir at 9 a.m. (Nakano et al., 1990) (Fig. 3B). Moreover, it was shown that the circadian oscillation of Cyp7a1 enzymatic activity is well correlated with the mRNA and paralleled the protein expression levels in rodents, with enzymatic activities being higher during nighttime and relatively lower during the daytime (Noshiro et al., 1990; Sundseth and Waxman, 1990; Kai et al., 1995; Chow et al., 2014). Our model was able to predict Hmgcr and Cyp7a1 mRNA expression levels and cholesterol levels and the cyclic variations well (Fig. 3). Although we have also studied cholesterol lowering in mice fed a high-fat/high-cholesterol diet to elevate both plasma and liver cholesterol and changes after treatment with 1,25(OH)2D3 (Chow et al., 2014), we only obtained end data at 48 hours after the last 1,25(OH)2D3 dose on cholesterol levels and altered gene (Shp and Cyp7a1) expression and not the temporal data. We did not fit these end data, and we recognized that the diet may modify the baseline turnover and metabolism of cholesterol or other feedback controls (Jones et al., 1996).

In addition, we used the PBPK(SFM)-PD model to describe the brain compartment. Our model readily predicted the induction (fold-change) of Mdr1 (Mdr1FC) mRNA expression over time in the kidney and the brain (Fig. 4), events that led to reduced renal and brain accumulation of [3H]digoxin, a P-gp substrate (Chow et al., 2011), due to the increased first-order excretion rate constant or efflux clearance, consequences of elevated Mdr1 transcription in the kidney and brain (or Mdr1FC; see Appendix). Similar to observations on intestinal Trpv6, renal Mdr1 expression also showed a greater extent of induction for the third and fourth doses. The reason is again unknown, except that there was a higher Vdr mRNA induction with subsequent doses (Chow et al., 2013). Additionally, brain Mdr1 mRNA expression is low (<1.5-fold) and tended to be overpredicted. The dynamic range for the increase in brain Mdr1 occurred over a narrow range, although the change in Mdr1 expression in kidney was considerably larger (Chow et al., 2011; Chow et al., 2013; Durk et al., 2014). However, the median prediction error for the brain Mdr1 fit is 11.5%, which is within the acceptable range (<50%), suggesting reasonably good model performance. The increase of brain Mdr1, although <1.5-fold, was significant (Durk et al., 2014).

In conclusion, the previously established PBPK(SFM)-PD model (Ramakrishnan et al., 2016) was extended for characterizing multiple data sets from the literature on 1,25(OH)2D3 pharmacodynamics (Chow et al., 2011, 2013, 2014). With the developed model, we demonstrate good capability of describing induction of Trpv5/6 for calcium absorption, Hmgcr and Cyp7a1 mRNA expression for cholesterol synthesis/degradation, and Mdr1 in reducing the brain accumulation of the P-gp substrate digoxin. Sensitivity analyses (Fig. 5) suggest that the Emax and EC50 values are important parameters governing pharmacodynamic responses. This is not unexpected, as Emax and EC50 reflect the maximum response and drug concentration producing 50% of maximum response, respectively (Sharma and Jusko, 1998). The circadian rhythmic variations of liver Cyp7a1 mRNA expression as well as cholesterol turnover were found important in cholesterol homeostasis. The model was relatively less sensitive to mRNA turnover rate constants and Hill coefficients. The developed model is quite robust and may find application in predicting biomarker (such as VDR-target genes) profiles after 1,25(OH)2D3 treatment in preclinical or clinical studies for the exploration of VDR therapeutic targets and treatment of diseases, especially cancer, and for in vivo–in vitro extrapolation or interspecies scaling of 1,25(OH)2D3 data to humans. We have used the PBPK(SFM) model (Ramakrishnan et al., 2016), based on mouse 1,25(OH)2D3 data, and scaled the model to predict clinical data in cancer patients after high i.v. and oral doses (Yang and Pang Abstract #W4117, presented at AAPS Annual Meeting and Exposition, November 12-15 2017, San Diego, CA), and will apply the PD component to the model in the near future.

Acknowledgments

We thank Dr. Keumhan Noh, Leslie Dan Faculty of Pharmacy, University of Toronto, for designing the primers for murine Trpv5, and Dr. Wojciech Krzyzanski and Dr. Rob Bies, Department of Pharmaceutical Science, University of Buffalo, and Dr. David D’Argenio, University of Southern California, for discussion.

Appendix

In the following equations, kin and kout denote zero-order production and first-order degradation rate constants, respectively, as described in the indirect response model (Sharma and Jusko, 1998; Mager et al., 2003; Ramakrishnan et al., 2016); SM1 and SM2 are the scaling factor on the plasma calcium and liver cholesterol turnover, respectively. Emax is the maximum inductive FC; EC50 is the 1,25(OH)2D3 tissue concentration (CT) that results in 50% of Emax; CL, CI, CK, and CBr are the 1,25(OH)2D3 concentrations in liver, intestine, kidney, and brain, respectively. Initial estimates (kin, kout, Emax, and EC50) were obtained from plotting Vdr-target gene FC versus tissue 1,25(OH)2D3 (Ramakrishnan et al., 2016).

The extended PBPK(SFM)-PD model incorporates renal Trpv5 (eq. A1, Trpv5FC,K) as well as renal and intestinal (eqs. A2 and A3) Trpv6 mRNA FC (Trpv6FC,K and Trpv6FC,I, respectively):

Embedded Image(A1)Embedded Image(A2)Embedded Image(A3)

For describing plasma calcium concentration ([Ca2+]) due to multiple calcium channels,

Embedded Image(A4)

where Embedded Image= Embedded Image=Embedded Image (where R0 is the initial response = Embedded Image·VP; Embedded Image is 10.7 mg·dl−1 in control samples; VP is the plasma volume). Here, the definition of kin for calcium is not a rate constant for the synthetic enzyme but for intake or reabsorption of calcium by the calcium channels in the intestine and kidney; and kout for calcium is the degradation rate constant associated with the elimination or utilization of calcium.

Upon dividing eq. A4 by ([Embedded Image] VP), the equation becomes:

Embedded Image(A5)Embedded Image

We further simplified the above equation, with Embedded Image,Embedded Image(A6) The output plasma calcium is expressed as [Ca2+] = (Embedded Image·Embedded Image), unit of Embedded Image. When only intestinal Trpv6, renal Trpv5, or Trpv6 is considered for calcium absorption, then the respective Trpv6FC,I, Trpv5FC,K, or Trpv6FC,K term was included for describing plasma calcium:

Embedded Image(A7)

The extended PBPK(SFM)-PD model was modified to incorporate circadian rhythmic changes of liver Hmgcr (HmgcrFC,L in eqs. A8–A10) and Cyp7a1 mRNA levels (Cyp7a1FC,L, eqs. A11–A15) and liver cholesterol (eq. A16). Circadian rhythm effects, embellished onto the input rate (D'Argenio et al., 2009; Aoyama et al., 2010), were applied to describe the 24-hour diurnal variation of liver Hmgcr and Cyp7a1 mRNA expression (eqs. A8–A10 and A11–13, respectively): Embedded Image(A8)Embedded Image(A9)Embedded Image(A10)Embedded Image(A11)Embedded Image(A12)Embedded Image(A13)where Embedded Image and Embedded Image are the zero-order synthesis and first-order degradation rate constants of liver Hmgcr mRNA; Embedded Image and Embedded Image are the zero-order synthesis and first-order degradation rate constants of liver Cyp7a1 mRNA, respectively; Embedded Image and Embedded Image are the peak and the mean values and of the mRNA synthesis rate, respectively; R0 is the initial condition of the mRNA expression, Embedded Image = Hmgcrbaseline,L and Embedded Image = Cyp7a1baseline,L; t is the decimal clock time; and Embedded Image is the time to reach peak mRNA synthesis rate. A series of transit compartments and a time-delay function were added to improve fitting of Hmgcr and Cyp7a1 mRNA FC and to explain the Hmgcr- and Cyp7a1-mediated control of cholesterol synthesis and metabolism. In the absence of circadian rhythm, liver Hmgcr and Cyp7a1 mRNA are readily described by eqs. A8 and A11, respectively, and eqs. A9–10 and A12–A13 are irrelevant.

We further used time-delay functions (eqs. A14 and A15) to provide an improved fit for Cyp7a1 expression, where τ is the time-delay term in hours, and Atransit1 and Atransit2 are the amounts in transit compartments 1 and 2, respectively. The transit compartments improved the fit to CYp7a1 only:

Embedded Image(A14)Embedded Image(A15)

Subsequently, Cyp7a1 expression induced degradation of liver cholesterol (Embedded Image) with the endogenous synthesis of cholesterol (Embedded Image) being controlled by Hmgcr. The amount of cholesterol was modeled in a similar fashion as plasma calcium (eqs. A4–A6):Embedded Image(A16)where the basal concentration of cholesterol [Cholesterolbaseline] is 2.35 mg·g−1 liver at time = 0 and is fixed; the value was estimated as the average value of measured liver cholesterol concentration at baseline without treatment, and CholesterolFC = Embedded Image

The output of cholesterol is expressed as [Cholesterol] = (Embedded Image·[Cholesterolbaseline]), unit of mg·g−l liver.

The extended PBPK(SFM)-PD model that incorporates renal and brain Mdr1 (Mdr1FC,K and Mdr1FC,Br, respectively; eqs. A17 and A18) and digoxin accumulation extents (expressed as %reduced over control value) in kidney (Digoxin%reduced,K; eq. A19) and brain (Digoxin%reduced,Br; eq. A20) are described. We transformed the digoxin tissue concentration data in the kidney and brain after treatment by taking the difference (treated value − control value), normalized to the control value (R0) x 100% as (%reduced). Hence, the data of Chow et al. (2011) on digoxin, with/without 1,25(OH)2D3 treatment, in kidney and brain were recalculated to reflect the effect of 1,25(OH)2D3 treatment on reducing digoxin accumulation in tissue:

Embedded Image(A17)Embedded Image(A18)Embedded Image(A19)Embedded Image(A20)

Authorship Contributions

Participated in research design: Pang, Yang, Bukuroshi.

Provided samples and data for fitting: Chow.

Conducted experiments: Yang, Bukuroshi.

Performed data analysis: Yang, Quach, Bukuroshi.

Wrote or contributed to the writing of the manuscript: Yang, Bukuroshi, Pang.

Footnotes

    • Received June 27, 2017.
    • Accepted October 26, 2017.

  • This work was supported by the Canadian Institutes of Health Research (K.S.P., P.B.) the Centre for Collaborative Drug Research (K.S.P.), the National Sciences and Engineering Research Council of Canada (E.C.Y.C., H.P.Q.), and the Ontario Graduate Scholarship Program (Q.J.Y.).

  • Authors declare no conflict.

  • https://doi.org/10.1124/dmd.117.077271.

Abbreviations

AUC
area under the curve
f
unbound fraction
FC
fold-change of mRNA expression
Hmgcr
rodent 3-hydroxy-3-methyl-glutaryl-coenzyme A reductase
MDR1/Mdr1
human/rodent multidrug-resistance protein 1
1,25(OH)2D3
1α,25-dihydroxyvitamin D3
PBPK
physiologically based pharmacokinetic model
PE
prediction error
P-gp
P-glycoprotein
PK-PD
pharmacokinetic-pharmacodynamic
qPCR
quantitative polymerase chain reaction
SFM
segregated flow model for intestine
SHP
small heterodimer partner
TRPV/Trpv
human/rodent transient receptor potential cation channel, subfamily V
VDR/Vdr
human/rodent vitamin D receptor

References

  1. ↵
    1. Abraham AK,
    2. Mager DE,
    3. Gao X,
    4. Li M,
    5. Healy DR, and
    6. Maurer TS
    (2009) Mechanism-based pharmacokinetic/pharmacodynamic model of parathyroid hormone-calcium homeostasis in rats and humans. J Pharmacol Exp Ther 330:169–178.
    OpenUrlAbstract/FREE Full Text
  2. ↵
    1. Abraham AK,
    2. Maurer TS,
    3. Kalgutkar AS,
    4. Gao X,
    5. Li M,
    6. Healy DR,
    7. Petersen DN,
    8. Griffith DA, and
    9. Mager DE
    (2011) Pharmacodynamic model of parathyroid hormone modulation by a negative allosteric modulator of the calcium-sensing receptor. AAPS J 13:265–273.
    OpenUrlCrossRefPubMed
  3. ↵
    1. Aoyama T,
    2. Omori T,
    3. Watabe S,
    4. Shioya A,
    5. Ueno T,
    6. Fukuda N, and
    7. Matsumoto Y
    (2010) Pharmacokinetic/pharmacodynamic modeling and simulation of rosuvastatin using an extension of the indirect response model by incorporating a circadian rhythm. Biol Pharm Bull 33:1082–1087.
    OpenUrlCrossRefPubMed
  4. ↵
    1. Back P,
    2. Hamprecht B, and
    3. Lynen F
    (1969) Regulation of cholesterol biosynthesis in rat liver: diurnal changes of activity and influence of bile acids. Arch Biochem Biophys 133:11–21.
    OpenUrlCrossRefPubMed
  5. ↵
    1. Blankenship KA,
    2. Williams JJ,
    3. Lawrence MS,
    4. McLeish KR,
    5. Dean WL, and
    6. Arthur JM
    (2001) The calcium-sensing receptor regulates calcium absorption in MDCK cells by inhibition of PMCA. Am J Physiol Renal Physiol 280:F815–F822.
    OpenUrlAbstract/FREE Full Text
  6. ↵
    1. Bookout AL,
    2. Jeong Y,
    3. Downes M,
    4. Yu RT,
    5. Evans RM, and
    6. Mangelsdorf DJ
    (2006) Anatomical profiling of nuclear receptor expression reveals a hierarchical transcriptional network. Cell 126:789–799.
    OpenUrlCrossRefPubMed
  7. ↵
    1. Brenza HL and
    2. DeLuca HF
    (2000) Regulation of 25-hydroxyvitamin D3 1alpha-hydroxylase gene expression by parathyroid hormone and 1,25-dihydroxyvitamin D3. Arch Biochem Biophys 381:143–152.
    OpenUrlCrossRefPubMed
  8. ↵
    1. Chakraborty A,
    2. Krzyzanski W, and
    3. Jusko WJ
    (1999) Mathematical modeling of circadian cortisol concentrations using indirect response models: comparison of several methods. J Pharmacokinet Biopharm 27:23–43.
    OpenUrlCrossRefPubMed
  9. ↵
    1. Cho BH
    (1983) Improved enzymatic determination of total cholesterol in tissues. Clin Chem 29:166–168.
    OpenUrlAbstract/FREE Full Text
  10. ↵
    1. Chow EC,
    2. Durk MR,
    3. Cummins CL, and
    4. Pang KS
    (2011) 1α,25-dihydroxyvitamin D3 up-regulates P-glycoprotein via the vitamin D receptor and not farnesoid X receptor in both fxr(-/-) and fxr(+/+) mice and increased renal and brain efflux of digoxin in mice in vivo. J Pharmacol Exp Ther 337:846–859.
    OpenUrlAbstract/FREE Full Text
  11. ↵
    1. Chow EC,
    2. Magomedova L,
    3. Quach HP,
    4. Patel R,
    5. Durk MR,
    6. Fan J,
    7. Maeng HJ,
    8. Irondi K,
    9. Anakk S,
    10. Moore DD, et al.
    (2014) Vitamin D receptor activation down-regulates the small heterodimer partner and increases CYP7A1 to lower cholesterol. Gastroenterology 146:1048–1059.
    OpenUrlCrossRefPubMed
  12. ↵
    1. Chow EC,
    2. Quach HP,
    3. Vieth R, and
    4. Pang KS
    (2013) Temporal changes in tissue 1α,25-dihydroxyvitamin D3, vitamin D receptor target genes, and calcium and PTH levels after 1,25(OH)2D3 treatment in mice. Am J Physiol Endocrinol Metab 304:E977–E989.
    OpenUrlAbstract/FREE Full Text
  13. ↵
    1. Cong D,
    2. Doherty M, and
    3. Pang KS
    (2000) A new physiologically based, segregated-flow model to explain route-dependent intestinal metabolism. Drug Metab Dispos 28:224–235.
    OpenUrlAbstract/FREE Full Text
  14. ↵
    1. Cox EH,
    2. Veyrat-Follet C,
    3. Beal SL,
    4. Fuseau E,
    5. Kenkare S, and
    6. Sheiner LB
    (1999) A population pharmacokinetic-pharmacodynamic analysis of repeated measures time-to-event pharmacodynamic responses: the antiemetic effect of ondansetron. J Pharmacokinet Biopharm 27:625–644.
    OpenUrlCrossRefPubMed
  15. ↵
    1. D’Argenio DZ,
    2. Schumitzky A, and
    3. Wang X
    (2009) ADAPT 5 User’s Guide: Pharmacokinetic/Pharmacodynamic Systems Analysis Software, Biomedical Simulations Resource, Los Angeles, CA.
  16. ↵
    1. den Dekker E,
    2. Hoenderop JG,
    3. Nilius B, and
    4. Bindels RJ
    (2003) The epithelial calcium channels, TRPV5 & TRPV6: from identification towards regulation. Cell Calcium 33:497–507.
    OpenUrlCrossRefPubMed
  17. ↵
    1. Duffull SB and
    2. Aarons L
    (2000) Development of a sequential linked pharmacokinetic and pharmacodynamic simulation model for ivabradine in healthy volunteers. Eur J Pharm Sci 10:275–284.
    OpenUrlCrossRefPubMed
  18. ↵
    1. Durk MR,
    2. Han K,
    3. Chow EC,
    4. Ahrens R,
    5. Henderson JT,
    6. Fraser PE, and
    7. Pang KS
    (2014) 1α,25-Dihydroxyvitamin D3 reduces cerebral amyloid-β accumulation and improves cognition in mouse models of Alzheimer’s disease. J Neurosci 34:7091–7101.
    OpenUrlAbstract/FREE Full Text
  19. ↵
    1. Edwards PA,
    2. Muroya H, and
    3. Gould RG
    (1972) In vivo demonstration of the circadian thythm of cholesterol biosynthesis in the liver and intestine of the rat. J Lipid Res 13:396–401.
    OpenUrlAbstract
  20. ↵
    1. Emond C,
    2. Birnbaum LS, and
    3. DeVito MJ
    (2006) Use of a physiologically based pharmacokinetic model for rats to study the influence of body fat mass and induction of CYP1A2 on the pharmacokinetics of TCDD. Environ Health Perspect 114:1394–1400.
    OpenUrlCrossRefPubMed
  21. ↵
    1. Faict D,
    2. De Moor P,
    3. Bouillon R,
    4. Heyns W,
    5. Heiniger HJ,
    6. Corrow D, and
    7. Lesaffre E
    (1986) Transcortin and vitamin D-binding protein levels in mouse serum. J Endocrinol 109:141–147.
    OpenUrlAbstract/FREE Full Text
  22. ↵
    1. Folch J,
    2. Lees M, and
    3. Sloane Stanley GH
    (1957) A simple method for the isolation and purification of total lipides from animal tissues. J Biol Chem 226:497–509.
    OpenUrlFREE Full Text
  23. ↵
    1. Gielen J,
    2. Van Cantfort J,
    3. Robaye B, and
    4. Renson J
    (1975) Rat-liver cholesterol 7alpha-hydroxylase. 3. New results about its circadian rhythm. Eur J Biochem 55:41–48.
    OpenUrlPubMed
  24. ↵
    1. Gold MS and
    2. Gebhart GF
    (2010) Nociceptor sensitization in pain pathogenesis. Nat Med 16:1248–1257.
    OpenUrlCrossRefPubMed
  25. ↵
    1. Hazra A,
    2. Pyszczynski N,
    3. DuBois DC,
    4. Almon RR, and
    5. Jusko WJ
    (2007) Modeling receptor/gene-mediated effects of corticosteroids on hepatic tyrosine aminotransferase dynamics in rats: dual regulation by endogenous and exogenous corticosteroids. J Pharmacokinet Pharmacodyn 34:643–667.
    OpenUrlCrossRefPubMed
  26. ↵
    1. Henry HL
    (2001) The 25(OH)D(3)/1alpha,25(OH)(2)D(3)-24R-hydroxylase: a catabolic or biosynthetic enzyme? Steroids 66:391–398.
    OpenUrlCrossRefPubMed
  27. ↵
    1. Ho KJ
    (1979) Circadian rhythm of cholesterol biosynthesis: dietary regulation in the liver and small intestine of hamsters. Int J Chronobiol 6:39–50.
    OpenUrlPubMed
  28. ↵
    1. Ette EI and
    2. Williams PJ
    1. Jin YJ and
    2. Jusko WJ
    (2007) Pharmacogenomics and pharmacoinetic/pharmacodynamic modeling, in Pharmacometrics (Ette EI and Williams PJ eds) pp 509–528, John Wiley & Sons, Hoboken, NJ.
  29. ↵
    1. Jones G,
    2. Strugnell SA, and
    3. DeLuca HF
    (1998) Current understanding of the molecular actions of vitamin D. Physiol Rev 78:1193–1231.
    OpenUrlAbstract/FREE Full Text
  30. ↵
    1. Jones PJ,
    2. Pappu AS,
    3. Hatcher L,
    4. Li ZC,
    5. Illingworth DR, and
    6. Connor WE
    (1996) Dietary cholesterol feeding suppresses human cholesterol synthesis measured by deuterium incorporation and urinary mevalonic acid levels. Arterioscler Thromb Vasc Biol 16:1222–1228.
    OpenUrlAbstract/FREE Full Text
  31. ↵
    1. Kai M,
    2. Eto T,
    3. Kondo K,
    4. Setoguchi Y,
    5. Higashi S,
    6. Maeda Y, and
    7. Setoguchi T
    (1995) Synchronous circadian rhythms of mRNA levels and activities of cholesterol 7 alpha-hydroxylase in the rabbit and rat. J Lipid Res 36:367–374.
    OpenUrlAbstract
  32. ↵
    1. Kandutsch AA and
    2. Saucier SE
    (1969) Prevention of cyclic and triton-induced increases in hydroxymethylglutaryl coenzyme A reductase and sterol synthesis by puromycin. J Biol Chem 244:2299–2305.
    OpenUrlAbstract/FREE Full Text
  33. ↵
    1. Lemay J,
    2. Demers C,
    3. Hendy GN,
    4. Delvin EE, and
    5. Gascon-Barré M
    (1995) Expression of the 1,25-dihydroxyvitamin D3-24-hydroxylase gene in rat intestine: response to calcium, vitamin D3 and calcitriol administration in vivo. J Bone Miner Res 10:1148–1157.
    OpenUrlPubMed
  34. ↵
    1. Louizos C,
    2. Yáñez JA,
    3. Forrest ML, and
    4. Davies NM
    (2014) Understanding the hysteresis loop conundrum in pharmacokinetic/pharmacodynamic relationships. J Pharm Pharm Sci 17:34–91.
    OpenUrlCrossRefPubMed
  35. ↵
    1. Lundåsen T,
    2. Gälman C,
    3. Angelin B, and
    4. Rudling M
    (2006) Circulating intestinal fibroblast growth factor 19 has a pronounced diurnal variation and modulates hepatic bile acid synthesis in man. J Intern Med 260:530–536.
    OpenUrlCrossRefPubMed
  36. ↵
    1. Mager DE,
    2. Wyska E, and
    3. Jusko WJ
    (2003) Diversity of mechanism-based pharmacodynamic models. Drug Metab Dispos 31:510–518.
    OpenUrlAbstract/FREE Full Text
  37. ↵
    1. Makin G,
    2. Lohnes D,
    3. Byford V,
    4. Ray R, and
    5. Jones G
    (1989) Target cell metabolism of 1,25-dihydroxyvitamin D3 to calcitroic acid. Evidence for a pathway in kidney and bone involving 24-oxidation. Biochem J 262:173–180.
    OpenUrlAbstract/FREE Full Text
  38. ↵
    1. Mayer D
    (1976) The circadian rhythm of synthesis and catabolism of cholesterol. Arch Toxicol 36:267–276.
    OpenUrlCrossRefPubMed
  39. ↵
    1. Miao J,
    2. Xiao Z,
    3. Kanamaluru D,
    4. Min G,
    5. Yau PM,
    6. Veenstra TD,
    7. Ellis E,
    8. Strom S,
    9. Suino-Powell K,
    10. Xu HE, et al.
    (2009) Bile acid signaling pathways increase stability of small heterodimer partner (SHP) by inhibiting ubiquitin-proteasomal degradation. Genes Dev 23:986–996.
    OpenUrlAbstract/FREE Full Text
  40. ↵
    1. Nakano A,
    2. Tietz PS, and
    3. LaRusso NF
    (1990) Circadian rhythms of biliary protein and lipid excretion in rats. Am J Physiol 258:G653–G659.
    OpenUrlAbstract/FREE Full Text
  41. ↵
    1. Noshiro M,
    2. Nishimoto M, and
    3. Okuda K
    (1990) Rat liver cholesterol 7 alpha-hydroxylase. Pretranslational regulation for circadian rhythm. J Biol Chem 265:10036–10041.
    OpenUrlAbstract/FREE Full Text
  42. ↵
    1. Oishi K,
    2. Amagai N,
    3. Shirai H,
    4. Kadota K,
    5. Ohkura N, and
    6. Ishida N
    (2005) Genome-wide expression analysis reveals 100 adrenal gland-dependent circadian genes in the mouse liver. DNA Res 12:191–202.
    OpenUrlCrossRefPubMed
  43. ↵
    1. Panda S,
    2. Antoch MP,
    3. Miller BH,
    4. Su AI,
    5. Schook AB,
    6. Straume M,
    7. Schultz PG,
    8. Kay SA,
    9. Takahashi JS, and
    10. Hogenesch JB
    (2002) Coordinated transcription of key pathways in the mouse by the circadian clock. Cell 109:307–320.
    OpenUrlCrossRefPubMed
  44. ↵
    1. Quach HP,
    2. Yang QJ,
    3. Chow EC,
    4. Mager DE,
    5. Hoi SY, and
    6. Pang KS
    (2015) PKPD modelling to predict altered disposition of 1α,25-dihydroxyvitamin D3 in mice due to dose-dependent regulation of CYP27B1 on synthesis and CYP24A1 on degradation. Br J Pharmacol 172:3611–3626.
    OpenUrlCrossRefPubMed
  45. ↵
    1. Ramakrishnan V,
    2. Yang QJ,
    3. Quach HP,
    4. Cao Y,
    5. Chow EC,
    6. Mager DE, and
    7. Pang KS
    (2016) Physiologically-based pharmacokinetic-pharmacodynamic modeling of 1α,25-dihydroxyvitamin D3 in mice. Drug Metab Dispos 44:189–208.
    OpenUrlAbstract/FREE Full Text
  46. ↵
    1. Sager JE,
    2. Yu J,
    3. Ragueneau-Majlessi I, and
    4. Isoherranen N
    (2015) Physiologically based pharmacokinetic (PBPK) modeling and simulation approaches: a systematic review of published models, applications, and model verification. Drug Metab Dispos 43:1823–1837.
    OpenUrlAbstract/FREE Full Text
  47. ↵
    1. Sarkinen BE
    (2011) Vitamin D deficiency & cardiovascular disease. Nurse Pract 36:46–53.
    OpenUrlPubMed
  48. ↵
    1. Scheff JD,
    2. Calvano SE,
    3. Lowry SF, and
    4. Androulakis IP
    (2010) Modeling the influence of circadian rhythms on the acute inflammatory response. J Theor Biol 264:1068–1076.
    OpenUrlCrossRefPubMed
  49. ↵
    1. Shapiro DJ and
    2. Rodwell VW
    (1969) Diurnal variation and cholesterol regulation of hepatic HMG-CoA reductase activity. Biochem Biophys Res Commun 37:867–872.
    OpenUrlCrossRefPubMed
  50. ↵
    1. Sharma A and
    2. Jusko WJ
    (1998) Characteristics of indirect pharmacodynamic models and applications to clinical drug responses. Br J Clin Pharmacol 45:229–239.
    OpenUrlCrossRefPubMed
  51. ↵
    1. Shefer S,
    2. Hauser S,
    3. Lapar V, and
    4. Mosbach EH
    (1972) Diurnal variation of HMG CoA reductase activity in rat intestine. J Lipid Res 13:571–573.
    OpenUrlAbstract
  52. ↵
    1. Shinki T,
    2. Jin CH,
    3. Nishimura A,
    4. Nagai Y,
    5. Ohyama Y,
    6. Noshiro M,
    7. Okuda K, and
    8. Suda T
    (1992) Parathyroid hormone inhibits 25-hydroxyvitamin D3-24-hydroxylase mRNA expression stimulated by 1 alpha,25-dihydroxyvitamin D3 in rat kidney but not in intestine. J Biol Chem 267:13757–13762.
    OpenUrlAbstract/FREE Full Text
  53. ↵
    1. Sundseth SS and
    2. Waxman DJ
    (1990) Hepatic P-450 cholesterol 7 alpha-hydroxylase. Regulation in vivo at the protein and mRNA level in response to mevalonate, diurnal rhythm, and bile acid feedback. J Biol Chem 265:15090–15095.
    OpenUrlAbstract/FREE Full Text
  54. ↵
    1. Turunen MM,
    2. Dunlop TW,
    3. Carlberg C, and
    4. Väisänen S
    (2007) Selective use of multiple vitamin D response elements underlies the 1 alpha,25-dihydroxyvitamin D3-mediated negative regulation of the human CYP27B1 gene. Nucleic Acids Res 35:2734–2747.
    OpenUrlCrossRefPubMed
  55. ↵
    1. Urva SR,
    2. Yang VC, and
    3. Balthasar JP
    (2010) Physiologically based pharmacokinetic model for T84.66: a monoclonal anti-CEA antibody. J Pharm Sci 99:1582–1600.
    OpenUrlCrossRefPubMed
  56. ↵
    1. Valdivielso JM,
    2. Cannata-Andía J,
    3. Coll B, and
    4. Fernández E
    (2009) A new role for vitamin D receptor activation in chronic kidney disease. Am J Physiol Renal Physiol 297:F1502–F1509.
    OpenUrlAbstract/FREE Full Text
  57. ↵
    1. Van Cantfort J and
    2. Gielen JE
    (1979) Comparison of rat and mouse circadian rhythm of cholesterol-7 alpha-hydroxylase activity. J Steroid Biochem 10:647–651.
    OpenUrlCrossRefPubMed
  58. ↵
    1. Westerhout J,
    2. Ploeger B,
    3. Smeets J,
    4. Danhof M, and
    5. de Lange EC
    (2012) Physiologically based pharmacokinetic modeling to investigate regional brain distribution kinetics in rats. AAPS J 14:543–553.
    OpenUrlCrossRefPubMed
  59. ↵
    1. Yang QJ,
    2. Fan J,
    3. Chen S,
    4. Liu L,
    5. Sun H, and
    6. Pang KS
    (2016) Metabolite kinetics: the segregated flow model for intestinal and whole body physiologically based pharmacokinetic modeling to describe intestinal and hepatic glucuronidation of morphine in rats in vivo. Drug Metab Dispos 44:1123–1138.
    OpenUrlAbstract/FREE Full Text
  60. ↵
    1. Yang X,
    2. Downes M,
    3. Yu RT,
    4. Bookout AL,
    5. He W,
    6. Straume M,
    7. Mangelsdorf DJ, and
    8. Evans RM
    (2006) Nuclear receptor expression links the circadian clock to metabolism. Cell 126:801–810.
    OpenUrlCrossRefPubMed
  61. ↵
    1. Yano Y,
    2. Beal SL, and
    3. Sheiner LB
    (2001) Evaluating pharmacokinetic/pharmacodynamic models using the posterior predictive check. J Pharmacokinet Pharmacodyn 28:171–192.
    OpenUrlCrossRefPubMed
  62. ↵
    1. Yao Z,
    2. DuBois DC,
    3. Almon RR, and
    4. Jusko WJ
    (2006) Modeling circadian rhythms of glucocorticoid receptor and glutamine synthetase expression in rat skeletal muscle. Pharm Res 23:670–679.
    OpenUrlCrossRefPubMed
Back to top

In this issue

Download PDF
Article Alerts
Email Article
Citation Tools

Research ArticleArticle

PBPK-PD on VDR targets in mice

Qi Joy Yang, Paola Bukuroshi, Holly P. Quach, Edwin C. Y. Chow and K. Sandy Pang
Drug Metabolism and Disposition January 1, 2018, 46 (1) 75-87; DOI: https://doi.org/10.1124/dmd.117.077271
del.icio.us logo Digg logo Reddit logo Twitter logo Facebook logo Google logo Mendeley logo

Jump to section

Advertisement