Physiologically-Based Pharmacokinetic Modelling of Serum 25-Hydroxyvitamin D Concentrations in Schoolchildren Receiving Weekly Oral Vitamin D3 Supplementation

1. Introduction

Vitamin D plays crucial roles in bone health and development in children, as well as in maintaining calcium, phosphate, and magnesium homeostasis in the body [1,2]. Vitamin D deficiency is defined as serum 25-hydroxyvitamin D₃ (25(OH)D₃) concentrations less than 50 nmol/L by the Institute of Medicine (IOM) [3]. Following this definition, a worldwide study spanning 2000 to 2022 revealed 48.5% (90% CI: [42.5-54.5%]) of individuals <18 years of age are deficient [4]. In particular, regional studies found deficiency in 39.1% of under 18-year-olds in the United Arab Emirates [5], 64% of the 9-to 13-year-olds in Poland in March [6], and 65% children aged 6 to 12 years in Pakistan [7]. Deficiency is partly due to insufficient sun exposure during winter and spring [5,6,7,8]. By even the most conservative serum 25(OH)D₃ threshold of 25 nmol/L, vitamin D deficiency is still frequent in both low- and high-income countries [8].

A meta-regression of 1429 subjects across 11 trials with vitamin D-fortified foods, conducted by Cashman and colleagues [9], revealed that the dietary requirement for 97.5% of children aged 2-17 years to maintain serum 25(OH)D₃ over 50 nmol/L is 24.4 µg/day (95% CI: 16.8-38.2 µg/day) [9]. This model was adjusted for baseline levels, age, and BMI from the population living between latitudes 40° and 63°N [9].

Compared with regression modelling, pharmacokinetic (PK) modelling is often used to predict longitudinal changes in circulating concentrations of metabolites following administration of a drug or nutritional supplement, and has been adopted to inform dosing guidelines for paediatric clinical care [10,11]. Previously, a one-compartment model of 25(OH)D₃ PK in children with chronic kidney disease suggested that a weight-based dosing strategy is needed [12]. Unfortunately, population predictions using this model systematically overestimated serum 25(OH)D₃ concentrations, indicating that the model structure needed changing in order to make reliable predictions. Taken together, these factors motivated us to develop a PK model and use it to examine the 24.4 µg/day dose recommended by Cashman and colleagues [9].

For healthy adults, our previous physiologically-based pharmacokinetic (PBPK) modelling accurately predicted the mean serum 25(OH)D₃ concentrations across a wide range of vitamin D₃ doses in 155 trial arms, from 12.5 µg/day to 1250 µg/day for daily dosing, and up to 50,000 µg for bolus dosing [13]. This naïve average model highlighted the non-linear dose-exposure relationship for 25(OH)D₃. We further developed the model into a nonlinear mixed effects model and identified weight as a predictive covariate for PK in healthy adults [14]. Here, we modified this adult model to fit data from a randomised controlled trial of vitamin D supplementation, which was conducted in Cape Town, South Africa to test the prevention of tuberculosis infection and disease [15], and compared model predictions using data from a sister trial in Mongolia [16].

2. Materials and Methods

2.1. Study Participants and Clinical Data Collection

Children in Grades 1-4 (aged 6 to 11 years) at participating schools in the Klipfontein district of Cape Town, South Africa, were recruited to the ViDiKids trial [15]. Over a 3-year period, each child randomised to the intervention arm received a weekly oral capsule containing 250µg (10,000 IU) vitamin D₃. Capsule administration was recorded by the trial team during term-times and was supervised by parents/guardians during school holidays and COVID-19 lockdowns. Dosing records for each participant were used for model fitting and simulation. All received doses were assumed to be taken. Blood samples were collected at baseline and 36 months for all participants; and at 6, 12 and 24 months in a subset. Serum 25(OH)D₂ and 25(OH)D₃ concentrations were quantified with liquid chromatography tandem mass spectrometry [17]. The concentration of serum 25(OH)D₂ was comparable to approximately 1% that of the serum 25(OH)D₃, and serum 25(OH)D₂ was not consistently detected. The fitting and testing of the model only used the 25(OH)D₃ data. For simplicity, vitamin D₃ and 25(OH)D₃ are referred to as vitamin D and 25(OH)D onwards, respectively. In addition, the 6-month time point was confounded by seasonal changes and was excluded from modelling. The model was fitted to 77 children and was tested on data from 463 other children.

2.2. Physiological Parameters

We averaged the cardiac output for boys (4.8 L/min) and girls (4.3 L/min) aged 7 to 9 [18] and took the mean 4.55 L/min (273 L/h) for our model. The average blood flow for children in this model was estimated assuming the fraction of blood flow going into each compartment was the same as a 70 kg man [18]. For simplicity, we assumed the proportion of each compartment’s volume was the same for children and adults, with a 1 kg/L apparent density. All physiological parameters are presented in Table S3.

2.3. Model Structure

We developed a series of 11 nested models with different covariate structures to select the best PBPK model.

The base model (Model 1) was based on our published adult model [13], which considered the arterial blood (V_art), venous blood (V_ven), the liver (V_l), and the rest-of-the-body compartment (V_rb) which lumped together all non-elimination organs. These four compartments were applied to models 1-8 and 11 (Figure S1). In models 9 and 10, non-elimination organs were divided into fat mass (V_fm) and lean mass (V_lm) (Figure S2). All model equations are listed in the supplementary materials.

Partition coefficients K_p were introduced to describe the ratio between specific compartments and the venous blood. K_p>1 indicates the concentration is higher in a compartment than in the venous blood. K_prb and K_p25rb are for vitamin D and 25(OH)D in the rest of the body compartment, respectively; K_pl and K_p25l are for vitamin D and 25(OH)D in the liver, respectively; K_pfm is for 25(OH)D in fat mass (models 9-10); and K_plm is for 25(OH)D in lean mass (models 9-10).

F_m refers to the fraction of vitamin D metabolised into 25(OH)D, which was assumed to be 0.33 [13]. CL_H is the hepatic clearance rate (L/h).

The clearance of 25(OH)D was described by a sigmoidal function, to describe the nonlinear relationship between dose and 25(OH)D serum concentration:

$${\displaystyle \frac{{\mathrm{CL}}_{\max }\text{×}{\mathrm{C}}_{25\left(\mathrm{OH}\right)\mathrm{D}}^{\mathsf{\gamma }}}{{{\mathrm{C}}_{50}}^{\mathsf{\gamma }}+{\mathrm{C}}_{25\left(\mathrm{OH}\right)\mathrm{D}}^{\mathsf{\gamma }}}}\text{×}{\mathrm{C}}_{25\left(\mathrm{OH}\right)\mathrm{D}}$$

To speed up simulation, we assumed vitamin D level in venous blood was kept in a steady state with arterial blood and the rest of the body. Thus, vitamin D concentrations of these two compartments were explicit functions of the venous blood concentration (equations 1 and 2). Simulated 25(OH)D was not affected by this assumption and the computing time using the “daspk” ordinary differential equations (ODE) solver, an implicit method for stiff ODEs and differential-algebraic equations (DAEs) implemented in R version 4.2.2, was reduced by ~90%. See supplementary section 1 for details.

$${\mathrm{A}}_{\mathrm{art}}={\displaystyle \frac{{\mathrm{A}}_{\mathrm{ven}}}{{\mathrm{V}}_{\mathrm{ven}}}}\text{×}{\displaystyle \frac{{\mathrm{Q}}_{\mathrm{co}}}{({\mathrm{Q}}_{\mathrm{l}}+{\mathrm{Q}}_{\mathrm{rb}})\text{×}{\mathrm{V}}_{\mathrm{art}}}}$$

(1)

$${\mathrm{A}}_{\mathrm{rb}}={\mathrm{V}}_{\mathrm{rb}}\text{×}{\displaystyle \frac{{\mathrm{K}}_{\mathrm{prb}}}{{\mathrm{V}}_{\mathrm{art}}\text{×}{\mathrm{A}}_{\mathrm{art}}}}$$

(2)

A_art and A_ven are the amount of vitamin D in the arterial and venous blood, respectively. Q_co , Q_l , and Q_rb are the cardiac output, the liver blood flow and the blood flow to the rest of the body, respectively. V_art and V_ven are the volume of distribution for the arterial and venous blood, respectively.

All models assumed combined additive and proportional residual errors. The fixed and random effects are tabulated (Tables 1 and S5).

The base model (Model 1) did not include any covariate to estimate parameters. All other models assumed that the volume of distribution of each compartment was directly proportional to weight (WT), and the maximum clearance rate for serum 25(OH)D (CL_MAX) was proportional to WT^0.75. Models 2-4 and 11 fixed the weight to its baseline value for each child. Models 5 and 6 assumed linear weight gain:

$$\mathrm{WT}\left(\mathrm{t}\right) = \mathrm{WT}\left(0\right) + 3.08 \mathrm{kg}/\mathrm{year} \text{×} \mathrm{time} \left(\mathrm{year}\right) \mathrm{for} \mathrm{boys}$$

(3)

$$\mathrm{WT}\left(\mathrm{t}\right)=\mathrm{WT}\left(0\right)+3.23 \mathrm{kg}/\mathrm{year} \text{×} \mathrm{time} \left(\mathrm{year}\right) \mathrm{for} \mathrm{girls}$$

(4)

Models 7-10 interpolated weight from baseline to 3 years linearly for each child:

$$\mathrm{WT}\left(\mathrm{t}\right) = \mathrm{WT}\left(0\right) + \left[\mathrm{WT}\right(3\mathrm{y}) - \mathrm{WT}(0\left)\right] \mathrm{kg} \text{×} {\displaystyle \frac{\mathrm{time} \left(\mathrm{year}\right)}{3 \mathrm{years}}}$$

(5)

In Model 6, ${\mathrm{K}}_{\mathrm{p}25\mathrm{rb}}$ was related to the BMI-for-age Z-score (ZBMI) via ${\mathrm{K}}_{\mathrm{p}25\mathrm{rb}}^{\mathrm{Z}\mathrm{B}\mathrm{M}\mathrm{I}}$:

$${\mathrm{K}}_{\mathrm{p}25\mathrm{rb}}={\mathrm{e}}^{{\mathrm{TK}}_{\mathrm{p}25\mathrm{rb}}+\mathrm{ZBMI}\text{×}{\mathrm{K}}_{\mathrm{p}25\mathrm{rb}}^{\mathrm{Z}\mathrm{B}\mathrm{M}\mathrm{I}}+{\mathsf{\eta }}_{{\mathrm{K}}_{\mathrm{p}25\mathrm{rb}}}}$$

(6)

In Models 9 and 10, the fraction of fat mass was reported to be related to ZBMI [19]:

$${\mathrm{f}}_{\mathrm{FM}}=\left\{\begin{array}{c}\mathrm{m}\mathrm{a}\mathrm{x}({\displaystyle \frac{28.61+7.82\text{×}\mathrm{ZBMI}-0.91\text{×}{\mathrm{ZBMI}}^2+0.03\text{×}{\mathrm{ZBMI}}^3}{100}},0.05), if ZBMI\le 6.19\\ 0.49, if ZBMI>6.19\end{array}\right.$$

(7)

2.4. Model Assessment and Selection

All nonlinear mixed-effects (NLME) models were fitted using the stochastic approximation expectation-maximization (SAEM) routine in nlmixr2 R package. To assess the fitting, we visualised the goodness-of-fit, distribution of conditional weighted residuals (CWRES), normalised prediction distribution errors (NPDEs) and performed a visual predictive check (VPC). We deemed η shrinkage <30% to be acceptable. Akaike information criterion (AIC) and Bayesian information criterion (BIC) were assessed for model selection. As conclusions based on the AIC and BIC were consistent, we did not prioritise one over the other.

For VPC, 1,000 samples of random effects (${\eta }_{{\mathrm{Kp}}_{25\mathrm{rb}}}$ and ${\eta }_{{\mathrm{CL}}_{\mathrm{MAX}}}$) and residual errors were drawn to simulate the 5th, 50th and 95th percentiles of serum 25(OH)D, together with 90% prediction intervals (PI) around each percentile.

2.5. Simulations

To test the model independently, we simulated 9 children in the treatment arm of the safety sub-study (baseline, 1, and 2 years) and another 454 children in the treatment group (baseline and 3 years). For each child, we generated 30 virtual studies to estimate the 95% prediction intervals of the simulated 2.5th, 50th and 97.5th percentiles. In each study, a realisation of the fixed effect uncertainty (${\mathrm{CL}}_{\max }$ and ${\mathrm{K}}_{\mathrm{p}25\mathrm{fm}}$) was made, and 2,500 samples of the IIV (${\mathsf{\eta }}_{{\mathrm{CL}}_{\max }}$) and the additive residual error were generated to simulate the PK according to the dosing record of that child. The proportional residual error was not sampled to ensure simulations started with the observed baseline values.

To simulate the European children, a virtual population was generated for 1,250 boys and 1,250 girls for the 6-year-olds and 11-year-olds, respectively. Baseline serum 25(OH)D was sampled as per the reported 5th, 25th, 50th, 75th and 95th percentiles for European children [20]. For simplicity, points were sampled uniformly between each pair of neighbouring percentiles. To cover the whole range, we assumed the lower and upper bounds of 10 and 90 nmol/L, respectively. The WHO provides weight-for-age information from birth to 10 years. For the 6-year group, weight was sampled as per the 1st, 3rd, 5th, 15th, 25th, 50th, 75th, 85th, 95th, 97th and 99th percentiles reported by the WHO [21]. Similarly, uniform distribution was assumed between two percentiles. We used the weight-for-age Z-score as a proxy for the BMI Z-score. For the 11-year group, the WHO does not provide weight distribution, but instead reports BMI and height distribution. Hence, we initially sampled the body mass index-for-age [21] and height-for-age [22] percentiles (both at 1st, 3rd, 5th, 15th, 25th, 50th, 75th, 85th, 95th, 97th and 99th ) assuming independent distribution, and calculated the weight for each virtual child. As was done for the simulation of each child, 30 virtual studies were generated to sample the fixed effect uncertainty, random effects and the additive errors to estimate the 95% prediction intervals of the simulated 2.5th, 50th and 97.5th percentiles.

2.6. Comparing Trial Results

To compare the results from the Cape Town trial with its sister Mongolian trial, we compared the mean weight and serum 25(OH)D at baseline, weekly dose (250µg Cape Town v 350µg Mongolia), compliance (i.e. received doses / number of weeks ×100%), serum 25(OH)D increase. We graphed the 3-year changes in serum 25(OH)D versus with compliance. Baseline weight and age were fitted to a linear regression model to infer the mean weight increase each year (boy: 3.0 kg/year; girl: 3.4 kg/year). Similar to the European children described in section 2.5, we simulated the Mongolian children where each child started with their baseline weight and baseline serum 25(OH)D and was assumed to conform to the mean weight increase per year. Apart from this, simulation was also carried out at $4\times {CL}_{max}$.

2.7. Computing Software and Environment

Baseline characteristics were analysed descriptively using IBM SPSS software, version 25 (IBM, NY, US). Categorical variables were summarised as frequencies and percentages, while continuous variables were presented as means with standard deviations. ZBMI for each child was computed in the WHO 2007 Z-Score calculator shiny app [23]. For NLME fitting, we used the nlmixr2 package version 2.0.9 in R (version 4.2.2) on a Linux high-performance cluster with 60 cores. Model fitting was performed using the SAEM routine implemented in nlmixr2. VPC was carried out using vpcSim function from the nlmixr2est R package version 2.1.5. Simulations were carried out with the rxode2 R package version 2.0.13 [24]. Graphical diagnostics were conducted using RStudio version 4.3.1.

3. Results

Out of the 1,682 participants in the ViDiKids trial, the first 200 participants enrolled on the study were also enrolled on a safety study. Among the 94 children in the treatment arm of this safety study, 77 children had their serum 25(OH)D consistently assessed at baseline, 1, 2 and 3 years, and were used to fit the model; 9 were assessed at baseline, 1 and 2 years, and were used to test the model; 8 had missing values and were excluded (Figure 1) [17]. Another 454 children were assessed at baseline and 3 years. They were also used to test the model (Figure 1). The baseline age and serum 25(OH)D of the 77 children for fitting were 8.9±1.3 years and 64.7±14.9 nmol/L, respectively. The baseline data were collected between March and September in 2017 (Table S1).

In the 94 children in the intervention arm of the safety study, mean serum 25(OH)D showed statistically significant increases (Figure S3). Baseline 25(OH)D was somewhat anti-correlated with ZBMI among the 77 children for model fitting (slope = -2.644 nmol/L; p-value = 0.10) (Figure S4).

3.1. The Final Model

The final model encompasses 5 compartments for vitamin D, including the GI tract, arterial blood, venous blood, the liver, and the rest of the body which lumps together all non-eliminating organs. For 25(OH)D, non-eliminating organs were divided into fat mass and lean mass (Figure 2). Predictive covariates include weight (to estimate the volume of distribution) and ZBMI (to estimate the fat mass and lean mass) [25]. Similar to our previous healthy adult PBPK model [13], a combined basal rate of vitamin D cutaneous synthesis and intake in each child was calculated so that the model was in balance with the baseline 25(OH)D prior to dosing (see supplementary materials Section 1 for details).

Model fitting made good parametric inference, with η shrinkage <30% (Table 2). Additive error (0.00249 nmol/L) was negligible, and proportional residual error was 0.109 which was small (Table 2). The random effect ${\eta }_{{CL}_{max}}$ was independent from covariates (Figure S5), and the model was fitted at the optimal value of Kp_25lm = 4 (Figure S6).

* IIV CV% is $\sqrt{e^{\omega }-1}\times 100\%$, where $\omega$ is the variance of random effects. Hence, ${\omega }_{{CL}_{max}}=0.332$.

The final model successfully repeated the observed serum 25(OH)D in 77 children for 3 years. Population predictions are in good agreement with observed serum 25(OH)D (Figure 3A). The concordance with data is higher for individual predictions (R2=0.70 in Figure 3B). The CWRES and NPDEs did not exhibit systemic errors, suggesting the model predictions were unbiased (Figure 2C-2F). Population and individual predictions for each child were plotted together with doses received during the 3-year span (Figure S7).

3.2. Visual Predictive Check

VPC verified model predictions. The observed 5th and 95th percentiles (red dashed lines in Figure 4) were within the 90% prediction intervals for the predicted 5th and 95th percentiles (shaded areas in Figure 4). The observed medians (black solid line in Figure 4) were also within 90% prediction intervals for the predicted medians. This confirms that model predictions were in good agreement with the fitting data.

3.3. Simulation of 9 Cape Town Children

To test the model, we first simulated 9 children from the same trial arm, who were not fitted. Here, we found observations for 5 out of 9 children were entirely within the 2.5th - 97.5th percentile range (Figure 5A, 5D, 5F, 5H, 5I). For the other 4 children, a 1-year observation (Figure 5C) and a 2-year observation (Figure 5E) were within the 2.5th - 97.5th predicted percentiles, and observations outside the 2.5th - 97.5th predicted percentiles were within ±20 nmol/L of the range (Figure 5B, 5C, 5E, 5G). This provided confidence in model predictions for up to 2 years of dosing.

In addition, we simulated further 454 children from the same trial, and found 76.9% of observed serum 25(OH)D were between the predicted 2.5th-97.5th percentiles (Figure S8). The mean difference between observation and the predicted 50th percentile is 5.4 nmol/L (90% CI: [-26.2, 39.5] nmol/L, Figure S9). The range of the predicted 2.5th to 97.5th percentiles changed over time. The range was smaller at the end of year 1, compared with the end of years 2 and 3 (Figure S10). Such range exhibited strong anti-correlation with baseline weight (R²=0.56) and was somewhat anti-correlated with baseline ZBMI (R²=0.18) (Figure S11).

3.4. Simulate Cashman’s 24.4µg Daily Dose

We simulated the proposed 24.4µg daily dose in European children of 6 and 11 years of age (Figure 6). Baseline 25(OH)D were sampled according to a reported European study [20]. Weight was based on WHO reports [21]. At 100% compliance, over 97.5% of the 6-year-old children were predicted to reach the 50 nmol/L threshold within 21 weeks (Figure 6A). For the 11-year-old children, 97.5% were predicted to reach the same threshold around 39 weeks (Figure 6B). This is intuitive, as older children are heavier. Their larger volume of distribution reduces the increase in serum 25(OH)D. In short, model simulation results were consistent with the expectation that 24.4µg daily vitamin D supplementation would be sufficient to reach the 50 nmol/L threshold among children aged 6-11 years.

3.5. Comparing Trials in Cape Town and Ulaanbaatar

Finally, we compared the Cape Town trial with the sister Mongolian trial (Figure S12). The mean baseline weight was similar for both studies (Figure S12C and Figure S12H), but baseline serum 25(OH)D was higher among the Cape Town children (mean = 64 nmol/L) than Mongolian children (mean = 38 nmol/L). Despite the lower compliance (Figures S12D and S12I) and the lower weekly dose (250µg Cape Town v 350µg Mongolia) for the Cape Town trial, serum 25(OH)D increase was higher among Cape Town children on average (mean = 32 nmol/L, Figure S12G) than Mongolian children (mean = 25 nmol/L, Figure S12B). In fact, many Mongolian children reported a 3-year serum 25(OH) that is lower than baseline (Figure S12B). In addition, the 3-year changes in serum 25(OH)D were not positively correlated with compliance (Figures S12E and S12J). Consistent with these findings, model simulations highlighted the difference between the two trials (Figure S13). The model overestimated the increase in the serum 25(OH)D at the final time point for Mongolian children (Figure S13A). To test whether this might be attributed to a higher clearance rate among the Mongolian children alone, we increased ${CL}_{max}$ by 4 folds, but the simulation results were not in perfect agreement with the Mongolian data (Figure S13B). This suggests other factors warrant considerations.

4. Discussion

Currently, there is no vitamin D PK model for healthy children. Here, we reported a PBPK model for serum 25(OH)D based on 3-year observations among 77 Cape Town children. The model gave reasonable predictions for another 463 children among the same trial, for up to 3 years, with the 95% CI covering 77.3% of the observed values.

The model is concerned with the long-term changes in serum 25(OH)D. We identified both weight and ZMBI as covariates to predict the volume of distribution. In addition, the range of the 2.5th to 97.5th predicted percentiles at 3 years was strongly inversely correlated with weight (Figure S11A). The range was also somewhat anti-correlated with ZBMI (Figure S11B), probably due to the fact that weight and ZBMI were associated. This is intuitive: the narrower range for a heavy child might be due to the lipophilic nature of vitamin D, which tends to accumulate in adipose tissue, thereby reducing variability in serum levels.

Our simulation of the European children who received 24.4µg daily vitamin D was consistent with the finding by Cashman that this dose would be sufficient to bring the serum 25(OH)D above the 50 nmol/L threshold in the majority of European children. Here are some additional considerations for model parameter and potential limitations.

4.1. 25(OH)D Disposition in Children

Previously, population PK modelling of children with chronic kidney disease reported the volume of distribution (V/F) for 25(OH)D was inferred to be 292 L (95% CI: [199,572] L) for a 24 kg child [12]. They used a 1-compartment model to link vitamin D₃ dose directly to serum 25(OH)D, assuming each molecule of vitamin D₃ produces one molecule of 25(OH)D. Our model assumes that every 3 molecules of vitamin D produce 1 molecule of 25(OH)D, which is more accurate. Had Wan et al made the same assumption, their V/F would have been 97 L for a 24 kg child, or 122 L for a 30 kg child. These values would agree with our parameter values given the reasoning below.

In our PBPK model, volume of distribution V_d = Σ(V_organ × K_p). For a 30kg child with ZBMI=0, 28.61% of weight is expected to be fat mass (i.e. 8.58 kg) [19]. Assuming 0.6 L arterial blood, 1.80 L venous blood and a 0.77 L liver (with K_p=1) (Table S3 and S4), the lean mass should be 30 x (100%–28.61%)–1.80–0.77–0.6 = 18.25 kg. Hence, V_d = Σ(V_organ × K_p)= 1.8 + 0.6 + 0.77 × 1 + 8.58 × 4 + 18.25 × 4.71 = 123 L for 25(OH)D. This is consistent with what was reported for children with chronic kidney disease. We speculate the high K_p for 25(OH)D in lean mass (K_p25lm = 4) may indicate 25(OH)D efficiently binds with vitamin D receptors in children to meet the demand of rapid growth. This leads to the high V_d in children. In contrast, we previously reported the volume of distribution for vitamin D was between 10.0-16.6 L for a 70 kg man from a non-compartmental analysis [14]. In that model, the volume of distribution for the rest of the body (including all non-eliminating organs) V_rb = 62.2 L. Because V_d = Σ(V_organ × K_p), we reasoned the partition coefficient for vitamin D K_prb must be lower than 1 [14]. Since 25(OH)D has similar physicochemical properties with vitamin D, K_p25rb<1 for adults was also expected.

4.2. Maximum Serum 25(OH)D Clearance Rate Constant

Here, we reported the maximum 25(OH)D clearance rate constant CLmax for a 30 kg child is 0.0119 h^-1 (95% CI: [0.0102, 0.0139] h^-1) (Table 1). For a 70 kg child, this extrapolates to 0.0225 h^-1 (95% CI: [0.0193, 0.0262] h^-1), and is smaller than a 70 kg adult we reported earlier (0.0463 h^-1, 95% CI: [0.0359, 0.0597] h^-1) [14]. This suggests there might be genuine differences in clearance between children and adults.

4.3. Limitation of the Data Used for Fitting and Testing

The data we fitted were sampled once a year for 3 years. The sampling was sparse, and we did not have observations for PK changes in shorter terms such as once every week. Hence, although the predictions made in years were reliable, the model does not necessarily capture short-term kinetics.

Secondly, the recorded dispensed doses were used to fit the model. Following long gaps in dosing (school holidays, COVID-19 lockdowns), serum 25(OH)D was predicted to diminish. Graphing the dosing information together with model predictions helped make sense of the fitting data (Figure S7). However, there were cases we could not explain, such as the lack of increase in the 1-year observations in children 41 and 92, and the 2-year observations in children 36 and 164 (Figure S7). These anomalous cases warrant further investigation.

Thirdly, the model does not consider seasonal variations. We fitted the serum 25(OH)D in the control group of the ViDiKids safety study (n=105) to a sinusoidal curve and inferred an amplitude of 11.5 nmol/L (Figure S14). Therefore, significant fluctuation is expected for Cape Town children within a year.

5. Conclusions

Here, we developed a population PBPK model to successfully fit the sparsely sampled 3-year time course of serum 25(OH)D concentrations in healthy Cape Town children started between 6-11 years. This model accurately predicted serum 25(OH)D for up to 3 years in 463 children from the same study who were not used to fit the model. Model simulation results agreed with Cashman’s finding that 24.4µg daily vitamin D is needed to achieve the 50 nmol/L sufficiency threshold in over 97.5% European children. The inferred volume of distribution was in agreement with another study on children with chronic kidney disease [12]. This model may be used to generate long-term forecasts in Cape Town children. However, Ulaanbaatar children exhibited significant differences in 25(OH)D dose-PK relationship which warrants further investigation.