Chiang Mai Journal of Science

Design and Implementation of Individualized Drug Dosage Regimen Based on Bayesian Hierarchical Model — A Case Study of Vancomycin

1. INTRODUCTION

     VANCOMYCIN has been recognized as a first-line therapeutic agent against severe Gram-positive infections for some time, due to its potent activity against methicillin-resistant Staphylococcus aureus (MRSA). However, recent clinical investigations have reported a progressive increase in the minimum inhibitory concentration (MIC) of Staphylococcus aureus to vancomycin, alongside a rising prevalence of heteroresistent strains such as heterogeneous vancomycin-intermediate Staphylococcus aureus (hVISA) and Vancomycin-intermediate Staphylococcus aureus (VISA). These findings have highlighted the urgent need for more precise and individualised dosing strategies, as they have led to increased risks of therapeutic failure and diminished efficacy.
     As a time-dependent antibiotic, the pharmacodynamic effect of vancomycin is closely linked to its exposure profile in vivo. The drug is primarily eliminated via renal excretion, and its pharmacokinetic (PK) parameters are significantly influenced by renal function, particularly in vulnerable populations such as the elderly, patients with renal insufficiency, obese individuals, ICU patients, and those with malignancies. Errors in dose adjustment can result in treatment failure or the escalation in the incidence of nephrotoxicity and ototoxicity. Conventional empirical dosing or trough-guided therapeutic drug monitoring (TDM) strategies are inadequate for high-risk populations.
     To address these challenges, the adoption of individualized dosing guided by population pharmacokinetics (PopPK) models has become widespread. PopPK frameworks consist of structural, covariate, and residual error models, which are based on the nonlinear mixed effects model (NONMEM). These models quantitatively describe the drug–physiology relationships. Common estimation methods include first-order (FO) and first-order conditional estimation with interaction (FOCE-I), coupled with validation techniques such as bootstrap resampling, visual predictive checks (VPC), and normalized prediction distribution errors (NPDE) to ensure robustness and predictive accuracy.
     Significant factors affecting vancomycin disposition have been identified as including body weight, age, sex, creatinine clearance (CLCr), hepatic biomarkers, albumin levels, concomitant medications, and MIC values. Rodvold et al. (2009) proposed that, particularly in critically ill patients with fluctuating renal function, dosing should be guided by the pharmacodynamic target of AUC0–24/MIC rather than trough levels. Subsequently, guidelines published by the Infectious Diseases Society of America (IDSA) in 2009 and updated in 2020 endorsed Bayesian-based AUC-guided dosing as a preferred strategy.
     The basis for contemporary PopPK modelling was established by Sheiner et al. [1] with the development of the NONMEM framework utilising FORTRAN. Subsequent to this seminal work, a plethora of estimation approaches were developed, encompassing semi-parametric, nonparametric, and two-stage methods.
     Bayesian inference has been incorporated into PopPK modelling in a progressive manner [2]. Gradually, studies have conducted a comprehensive assessment of the Bayesian methods in PopPK, thereby laying the foundation for uncertainty quantification and individualized inference. Subsequently, Frazee et al. [3] proposed a generalized linear model framework that assumed gamma-distributed PK data, thus relaxing the conventional normality assumptions.
     However, the majority of contemporary models continue to presuppose normality in residual structures, thereby constraining their adaptability to non-Gaussian or outlier-prone datasets. To address this, Rybak et al. [4] introduced heavy-tailed multivariate distributions to enhance model robustness. Subsequently, Tobler et al. [5] employed t-distributed error structures in panel data analysis using Bayesian inference for both location and scale parameters, thereby enriching the toolkit for modelling complex clinical data.
     Internationally, Yamamoto [6], Broeker [7], and others developed PopPK models for adults and ICU populations, identifying strong correlations between clearance (CL), age, and CLCr. In contrast, Edginton [8] established a one-compartment model for paediatric populations, emphasizing the need to account for developmental PK variability in extrapolations.
     Domestically, Li Y et al. [9-10] developed a one compartment PopPK model based on TDM data from 316 elderly patients, identifying age, body weight, and CLCr as major predictors of clearance and constructing corresponding dosage nomograms. In a further development, Li j et al. [11] enhanced the robustness of the model by incorporating t-distributed error structures in neonatal populations. Concurrently, Yang’s team developed a Bayesian posterior inference system using WinBUGS for individualized dosing in patients with renal impairment.
     Furthermore, Bayesian inference and Markov chain Monte Carlo (MCMC) methods are increasingly applied to estimate individual PK parameters via maximum a posteriori (MAP) or fully Bayesian frameworks, particularly in sparse sampling or limited TDM scenarios. The integration of PK/PD models using AUC/MIC as dosing targets has emerged as a critical decision-making tool for clinical optimisation.
     Research on vancomycin individualised dosing has evolved from empirical adjustment to model-driven precision. The integration of population pharmacokinetics (PopPK) and Bayesian methodologies provides a robust theoretical and computational foundation for dose optimisation. Nevertheless, there remains a paucity of modelling efforts for special populations, such as elderly Chinese patients and individuals with renal impairment, thus underscoring the necessity for enhanced model generalizability and external validation, a need which this study is designed to address. Vancomycin, a first-line antibiotic for the treatment of methicillin-resistant Staphylococcus aureus  (MRSA) infections, exemplifies the challenges inherent in individualised therapy. Maintaining therapeutic concentrations (10–20 mg/L ) through therapeutic drug monitoring (TDM) is critical to ensure clinical efficacy while minimizing toxicity. However, traditional pharmacokinetic models frequently lack precision due to sparse sampling and high interpatient variability. Recent studies have demonstrated that one-compartment PopPK models, when coupled with Bayesian estimation techniques, are particularly effective in addressing these limitations under sparse TDM conditions.
     In order to surmount the aforementioned challenges, the present study proposes a Bayesian hierarchical modelling framework that integrates Hamiltonian Monte Carlo (HMC) sampling with adaptive prior regularisation. This innovative approach dynamically adjusts prior intensity, thereby enhancing the accuracy of posterior estimations and enabling precise individualised dosing recommendations. The validity of this approach is assessed using simulated clinical data from patients with heart failure, with a view to determining the effectiveness of this novel method in improving therapeutic precision and patient safety.

2. MATERIALS AND METHODS

2.1 Materials
2.1.1 Data source
     The dataset is derived from the anonymised clinical dataset of heart failure patients, originally published by Federspiel (2021) and hosted on a mirrored site of the Hugging Face platform. All procedures were conducted in strict accordance with the “Ethical Review Measures for Biomedical Research Involving Human Subjects” to ensure compliance with ethical standards. Personally identifiable information was removed, and all experimental protocols adhered to clinical safety regulations, ensuring participant safety.
     The dataset was initially published on Kaggle under a Creative Commons license, comprising detailed clinical data for 293 patients, encompassing age, sex, serum creatinine concentrations, and survival status.

2.2 Methods
2.2.1 Pharmacokinetics
     Pharmacokinetics studies the movement of drugs within the human body, starting from the processes of absorption, distribution, metabolism, and excretion (abbreviated as ADME). Pharmacokinetic models provide a rigorous interpretation of the relationship between drug concentrations in plasma and pharmacological responses. The simplest pharmacokinetic model is the one-compartment model.
     The one-compartment pharmacokinetic model assumes that the drug distributes instantaneously throughout the body’s fluid space, achieving dynamic equilibrium with the site of action immediately after administration. This model is typically suitable for drugs that mix rapidly and distribute uniformly within the body. It simplifies the description of the drug’s kinetic behavior, enabling an initial understanding and calculation of fundamental pharmacokinetic parameters such as absorption, distribution, metabolism, and excretion. The mathematical representation of the one-compartment pharmacokinetic model [12] is often based on the following equation:

\(\dfrac{dc(t)}{dt} = -k \cdot c(t),\)      (1)

where \(c(t)\) is the concentration of the plasma drug at time \(t\), \(k\) represents the rate constant for elimination, and its definition formula is \(k = \dfrac{CL}{V}\). \(CL\) is clearance and \(V\) is the volume of distribution.

2.2.2 Bayesian hierarchical model
     In order to achieve the goal of vancomycin individualized [13–15], in this study, Bayesian multi-level random effects model was adopted, which was divided into four layers: observation layer, individual parameter layer, population distribution layer and prior distribution [15–16] layer. The following details the model building methods and formulas for each layer.
     In the hierarchical model of Lunn et al. [17], which is based on the Bayesian framework of Hastings [18], posterior distribution of the model parameters is given by:

\(\pi(\theta_i \mid c_{ij}) \propto \pi(c_{ij} \mid \theta_i) \cdot \pi(\theta_i),\)      (2)

where \(c_{ij}\) represents observed drug concentration of person \(i\) at time \(t_j\) [19].

1) Level 1: Observation Layer
     Assume that the random variations among individuals follow an exponential normal distribution [9][11], the relationship between the observed concentration and the predicted concentration was modeled as:

\(c_{ij,obs} = c_{ij,pre} \cdot e^{\varepsilon_{ij}}, \quad \varepsilon_{ij} \sim \mathcal{N}(0,\sigma^2),\) (3)

where \(c_{ij,obs}\) represents observed drug concentration of person \(i\) at time \(t_j\), \(c_{ij,pre}\) represents prediction of drug concentration of person \(i\) at time \(t_j\), which come from the one-compartment pharmacokinetic model. After taking the logarithm of formula (3):

\(\ln(c_{ij,obs}) = \ln(c_{ij,pre}) + \varepsilon_{ij},\) (4)

\(\ln(c_{ij,obs}) \sim \mathcal{N}(\ln(c_{ij,pre}), \sigma^2).\) (5)

2) Level 2: Individual Parameter Layer
     Since the i-th individual can be regarded as a sample drawn from the overall population, and the sample and the resident are identically distributed, based on the first layer of assumptions, it is easy to conclude:

\((CL_i, V_i)^T \sim \mathcal{N}_2(\mu_i, \Sigma).\) (6)

3) Level 3: Population Distribution Layer
     The modeling of individual variability driven by covariates was derived from the literature [13]. Given that the age range of the population in this study's dataset is consistent with the age range described in the literature, this formula is also applicable to the dataset of this study, specific form is as follows:

\(\mu_i = (\mu_{CL_i}, \mu_{V_i}),\) (7)

\(\mu_{CL_i} = \left[1.71 \cdot e^{\eta_{1i}} + 8.31 \cdot \left(1 - e^{-0.0113 \cdot \ln(2) \cdot CLcr_i}\right)\right] \cdot \left(\dfrac{AGE_i}{72}\right)^{-0.475},\) (8)

\(\mu_{V_i} = 26.2 \cdot 2.09 \cdot \left(\dfrac{AGE_i}{72}\right) \cdot e^{\eta_{3i}},\) (9)

where \(\eta_{ki} \sim \mathcal{N}(0,\omega_k^2)\), for \(k = 1, 2, 3\), are inter-individual variability in pharmacokinetic parameters, \(CLcr\) are calculated by the Cockcroft–Gault formula with age, weight and SCR. Notice that \(CLcr\) is a function of weight since there are only age and SCR is not weight in the dataset.

4) Level 4: Prior Distribution Layer
     The prior distribution layer is predicated on the prior hypothesis of parameters, which is established in accordance with statistical data and physiological significance.
It is a postulate of this study, based on Bayesian theory, that random effects follow standard normal distributions:

\(\eta_k \sim \mathcal{N}(0,\omega_k^2), \quad k = 1,2,3.\)       (10)

and the specific variability values are:

\(\omega_1 = \sqrt{Var(\eta_1)} = 0.415,\quad \omega_2 = \sqrt{Var(\eta_2)} = 0.766,\quad \omega_3 = \sqrt{Var(\eta_3)} = 0.381,\)

where the value of \(\omega_k\), \(k = 1, 2, 3\), come from the references [13]. These parameters were selected based on the population pharmacokinetics (PopPK) model described in reference [13]. These parameters are derived from Table 1 in the Section 2 of that reference and are specifically tailored to the biological characteristics of the adult patient population. Given that the observed data in this study mainly pertains to adults, the selected parameter values are in line with the pharmacokinetic characteristics of this population, thereby ensuring the applicability and accuracy of the model.
     According to the theory of conjugate priors in statistics, it is established that the covariance matrix \(\Sigma\) follows an inverse-Wishart distribution [20–21]:

\(\Sigma \sim \mathcal{W}_2^{-1}(m, W),\) (11)

where \(W\) is the scale matrix, and \(m \ge p = 2\) represents the degrees of freedom. The prior distribution of the residual variance, denoted by \(\sigma^2\), is modelled using inverse gamma distribution,

\(\sigma^2 \sim \mathrm{Inv\text{-}Gamma}(\alpha, \beta),\) (12)

where \(\alpha\) and \(\beta\) are the hyperparameters of the prior distribution.

5) Level 5: Posterior Distribution
     According to formulas (2)–(12), and based on Bayesian prior theory, we arrive as:

\(\pi(\theta_i \mid c_{i}) \propto \pi(c_{i} \mid \theta_i)\cdot \pi(\Sigma)\cdot \pi(CL_i, V_i \mid \eta_{1}, \eta_{2}, \eta_{3}) \cdot \pi(\eta_1)\cdot \pi(\eta_2)\cdot \pi(\eta_3)\)
\(= \mathcal{N}(\ln c_{pre}, \sigma^2)\mathcal{W}_m^{-1}(m+n, A+\Omega)\mathcal{N}_2(\mu, \Sigma)\prod_k \mathcal{N}(0,\omega_k^2),\) (13)

with \(\theta_i = (CL_i, V_i, \eta_{1i}, \eta_{2i}, \eta_{3i}, \Sigma)\) and

\(A = \sum_i (x_i - \mu)(x_i - \mu)^T.\) (14)

This hierarchical Bayesian formulation facilitates the comprehensive integration of inter- and intra-individual variability into the posterior inference process.

3. RESULTS AND DISCUSSION

     In recent years, researchers have enhanced the robustness of models by using adaptive prior enhancement methods. For instance, Chen [21] proposed dynamically adjusting the sampling strategy in MCMC, and Turner et al. [22] explored adaptive adjustments based on sample size. However, these methods based on traditional MCMC have limitations in sampling efficiency and handling high-dimensional parameter spaces. Most methods rely on predefined heuristic rules or require post-sampling adjustments, lacking a mechanism to dynamically and real-time adjust the prior influence based on the sample distribution density during the sampling process. This is inefficient and lacks a real-time dynamic adjustment mechanism for the prior in high-dimensional spaces.
     In this context, this paper proposes a Hamiltonian Monte Carlo method for dynamically adjusting the prior strength. We base it on the classical Hamiltonian Monte Carlo sampling as the basic framework, and focus on modifying its core potential energy function, and introducing a parameter \(\lambda\) that can be dynamically adjusted. This improvement draws on the ideas from Chen [21]'s research on the dynamic adjustment of traditional MCMC sampling, combines the advantages of Hamiltonian Monte Carlo in exploring high-dimensional parameter spaces, and meets the demand for more stable algorithms in the precise vancomycin administration field [7].

3.1 Prediction of Population Parameters Using Hamiltonian Monte Carlo (HMC) Sampling with Dynamically Adjusted Prior Strength
3.1.1 Classical Hamiltonian Monte Carlo sampling
     Hamiltonian Monte Carlo (HMC) sampling is a probabilistic approach that employs physical dynamics to simulate posterior distributions, according to Hamiltonian dynamics:

\(H(\theta, p) = U(\theta) + K(p),\) (15)

where \(U(\theta) = -\log p(\theta)\) represents potential energy, \(K(p) = \frac{1}{2} p^T M^{-1} p\) denotes kinetic energy. Parameters and momenta are updated via the Leapfrog integration method:

\(p_{t+\frac{1}{2}} = p_t - \frac{\delta}{2} \nabla U(\theta_t).\) (16)

\(\theta_{t+\delta} = \theta_t + \delta\, p_{t+\frac{\delta}{2}},\) (17)

\(p_{t+\delta} = p_{t+\frac{\delta}{2}} - \frac{\delta}{2}\nabla U(\theta_{t+\delta}).\) (18)

3.1.2 With the dynamic adjustment \(\lambda\) Hamiltonian Monte Carlo sampling
     To mitigate conflicts between sparse data and prior information (or to strengthen the robustness of the model in sparse data scenarios), we introduce a dynamic adjustment to the prior strength \(\lambda\) during the Hamiltonian Monte Carlo (HMC) sampling process.

The potential energy is modified to:

\(H(\theta, p) = U_\lambda(\theta) + K(p),\) (19)

where \(U(\theta) = -\lambda \log p(\theta)\) represents potential energy, \(K(p) = \frac{1}{2} p^T M^{-1} p\) denotes kinetic energy. Parameters and momenta are updated via the Leapfrog integrator method:

\(p_{t+\frac{\delta}{2}} = p_t - \frac{\delta}{2}\nabla U_\lambda(\theta_t),\) (20)

\(\theta_{t+\delta} = \theta_t + \delta\, p_{t+\frac{\delta}{2}},\) (21)

\(p_{t+\delta} = p_{t+\frac{\delta}{2}} - \frac{\delta}{2}\nabla U_\lambda(\theta_{t+\delta}).\) (22)

3.1.3 Dynamic Adjustment of Prior Strength in Hamiltonian Monte Carlo (HMC) sampling
     In order to enhance the efficiency of the sampling process, the prior intensity was dynamically adjusted within the framework of the classical Hamiltonian Monte Carlo (HMC) sampling procedure. The adjustment procedure comprises the following steps:

Step 1: Initialization of Joint Distribution: \(\pi(\theta, p)\).
     The joint distribution combines the posterior distribution \(\pi(\theta)\) with the auxiliary momentum distribution:

\(\pi(\theta, p) \propto \pi(\theta)\cdot \mathcal{N}(p \mid 0, M),\)

where \(M\) denotes the mass matrix, and initial values \(\theta^{(0)}, p^{(0)}\) are randomly assigned.

Step 2: Dynamic Adjustment of Prior Intensity.
     The prior distribution is dynamically adjusted via a scaling parameter \(\lambda\):

\(\pi_{\text{posterior}}(\theta; \lambda) \propto \exp(-\lambda \cdot U(\theta)),\)

where \(\lambda\) controls the prior intensity. Higher \(\lambda\) increases prior influence, whereas lower values emphasize data-driven inference.

Step 3: Simulation of Hamiltonian Dynamics.
     Parameters and momentum evolve according to Hamiltonian dynamics:

\(\dfrac{d\theta}{dt} = \dfrac{\partial H}{\partial p},\)

\(\dfrac{dp}{dt} = -\dfrac{\partial H}{\partial \theta}.\)

with Hamiltonian defined as:

\(H(\theta, p) = U(\theta) + \frac{1}{2} p^T M^{-1} p.\) (23)

Step 4: Sampling and Metropolis Acceptance.
     The Leapfrog integrator solves the Hamiltonian equations to propose new states, accepted with probability:

\(\alpha = \min\left\{1, \frac{\pi(\theta^{(k+1)}, p^{(k+1)})}{\pi(\theta^{(k)}, p^{(k)})}\right\}.\) (24)

Step 5: Updating Prior Strength \(\lambda\).
     The \(\lambda\) value adjusts dynamically—increased if samples concentrate in high-density regions or reduced when samples are scattered—thus enhancing efficiency.

Step 6: Iterative Sampling until Convergence.
     Sampling continues iteratively until sufficient \(\theta\) samples are collected, ensuring accurate posterior representation for statistical inference.
     The dynamic adjustment scheme adapts the prior strength \(\lambda\) in real-time based on the density of the sampled data, optimizing the influence of prior information during the model training process. This adjustment reduces reliance on sparse priors, allowing the model to more accurately reflect individual variability, thereby improving the precision of drug concentration predictions.

3.2 Prediction Model
     Using the posterior estimates of individual clearance (\(CL_i\)) and volume of distribution (\(V_i\)), the time-concentration relationship was described by a one-compartment model with first-order elimination kinetics:

\(c_i(t) = \frac{D_i}{V_i} \cdot e^{-\frac{CL_i}{V_i} t}.\) (25)

To determine the required dose to \(D_i\) reach a specified target concentration \(c_{\text{target}}\) at time \(t_{\text{target}}\), the equation is algebraically rearranged as:

\(D_i = c_{\text{target}} \cdot V_i \cdot e^{\frac{CL_i}{V_i} t_{\text{target}}}.\) (26)

For sequential doses administered at intervals \(T_n\), the plasma concentration dynamics integrate residual drug from prior administrations. The concentration after the n-th dose is modeled by a piecewise function:

\(c_n = \begin{cases} c_{n-1} e^{-\frac{CL_n}{V_n} t} + \dfrac{D_n}{T_n V_n}\left(1 - e^{-\frac{CL_n}{V_n} t}\right), & t \le T_n, \\ \left[c_{n-1} e^{-\frac{CL_n}{V_n} T_n} + \dfrac{D_n}{T_n V_n}\left(1 - e^{-\frac{CL_n}{V_n} T_n}\right)\right] e^{-\frac{CL_n}{V_n}(t-T_n)}, & t > T_n, \end{cases}\) (27)

     Here, \(c_{n-1}\) denotes the residual concentration (or trough concentration) from the preceding dose, while \(CL_n\) and \(V_n\) represent the clearance and volume of distribution specific to the n-th dosing interval. This model accounts for both infusion periods (\(t \le T_n\)) and post-infusion elimination phases (\(t > T_n\)), enabling precise dose adaptation for complex therapeutic schedules. The value of \(T_n\) is typically determined based on the drug’s infusion rate, desired therapeutic target, and clinical guidelines. For example, in the case of intravenous vancomycin administration, \(T_n\) could be set based on the recommended infusion duration (e.g., 1 to 2 hours), with adjustments made depending on patient-specific factors such as renal function and age.
     This individualized dosing algorithm provides a flexible and data-driven framework for personalized vancomycin therapy, enabling real-time dose optimization based on Bayesian pharmacokinetic inference. By continuously integrating prior population knowledge with individual patient data through hierarchical Bayesian updating, the algorithm dynamically refines pharmacokinetic parameter estimates such as clearance and volume of distribution.

4. MODEL SIMULATION ANALYSIS

4.1 Simulation Data
     In this study, core indicators including body weight, age, height, and serum creatinine (SCR) were extracted from 293 patients in the dataset. Firstly, it is easy to calculate the creatinine clearance (ClCr) for each patient by the Cockcroft-Gault formula. Subsequently, with the ClCr value as a key input parameter, the drug clearance (CL) and apparent volume of distribution (V) of each patient were further derived by combining Equations (8) and (9) in the Bayesian Hierarchical Model. Furthermore, the population-level Clearance (μ_CL, L/h) and Volume (μ_V, L), along with their 95% confidence intervals, were statistically derived. Detailed statistical results are presented in Table 1:

     The simulation of therapeutic drug monitoring (TDM) data for vancomycin was conducted using a one-compartment pharmacokinetic model with first-order elimination kinetics. The simulation incorporated a half-life range of 6–8 hours. Intravenous infusion doses ranged from 250 mg to 1000 mg, with sparse blood concentration sampling at 0.5 hours to 48 hours post-administration, in order to mimic real-world clinical scenarios (1–3 timepoints per patient). the final model estimated the clearance (CL) to be 3.39 L/h (range: 2.59 L/h – 4.19 L/h), and the volume of distribution (V) to be 66.9 L (range: 44.87 L – 76.69 L). Vancomycin mimic therapeutic drug monitoring (TDM) data are generated using an intravenous infusion dosing regimen, as outlined by the pharmacokinetic algorithm, as follows:

\(c_n = \begin{cases} c_{n-1} e^{-\frac{CL_n}{V_n} t} + \dfrac{D_n}{T_n V_n}\left(1 - e^{-\frac{CL_n}{V_n} t}\right), & t \le T_n, \\ \left[c_{n-1} e^{-\frac{CL_n}{V_n} T_n} + \dfrac{D_n}{T_n V_n}\left(1 - e^{-\frac{CL_n}{V_n} T_n}\right)\right] e^{-\frac{CL_n}{V_n}(t-T_n)}, & t > T_n, \end{cases}\)

     In this equation, \(c_n\) represents the plasma drug concentration at the n-th dose timepoint, \(c_{n-1}\) is the plasma drug concentration at the previous dose, \(CL_n\) is the clearance of the drug at the n-th dose, \(V_n\) is the volume of distribution at the n-th dose, \(D_n\) is the dose of the drug administered at the n-th dose, \(T_n\) is the infusion duration of the n-th dose, and \(t\) is the time post-dose in hour.

4.2 Predicted Concentration
     Based on the Bayesian hierarchical model and the HMC sampling method that dynamically adjusts the prior strength, combined with the simulation data simulated in Section 4.1, we further derived and analyzed the individual pharmacokinetic parameters \(\theta_i\). Based on experience, it is found that after discarding (i.e., burning-in) the first 1000 samples in the collected data, the Markov chain tends to stabilize. That is, the target samples obtained in the posterior distribution are effective. It is worth noting that the effective sample size is determined based on the reference minimum required threshold of 500 [22–24]. In this study the effective sample size \(n = 550\) is adopted to ensure adequate estimation precision. Therefore, the mean of the effective sample values serves as an estimate of the individual pharmacokinetic parameters \(\theta_i\), i.e., \(\theta_i = \frac{1}{550}\sum_{j=1}^{550}\theta_i^{(j)}\), where \(\theta_i^{(j)}\) is effective sample. So we obtain the \(CL_i\) and \(V_i\) belonging to the i-th individual, because \(CL_i\) and \(V_i\) are part of \(\theta_i\). Furthermore, we use the individual \(CL_i\) and \(V_i\) of the i-th subject, along with Equation (27), to derive the concentration at a specific observation time for the i-th individual.

4.3 Model Accuracy Analysis
     In order to evaluate the impact of dynamically adjusted prior strength on the model’s performance, a comparison was made of its predictive accuracy before and after adjustment using standard error metrics, including mean prediction error (MPE) and root mean square error (RMSE) [25]. These metrics quantify deviations between predicted concentrations \(\hat{c}_{pred}\) and observed values \(c_{obs}\) as defined below:

\( \mathrm{MPE} = \frac{1}{NM}\sum_{i=1}^{N}\sum_{j=1}^{M}\left(c_{\text{pred},ij} - c_{\text{obs},ij}\right) \) (29)

Where \(c_{\text{obs},ij}\) represents the observed drug concentration for person \(i\) at time \(t_j\), and \(c_{\text{pred},ij}\) represents the predicted drug concentration for person \(i\) at time \(t_j\), derived from the compartment pharmacokinetic model. The calculation is performed across all individuals and time points, where \(N\) is the total number of individuals and \(M\) represents the number of observation time points for each individual.

\(\mathrm{RMSE} = \sqrt{\frac{1}{N}\sum_{i=1}^{N}\sum_{j=1}^{M}\left(c_{\text{pred},ij} - c_{\text{obs},ij}\right)^2}\) (30)

Where \(c_{\text{obs},ij}\) represents the observed drug concentration for person \(i\) at time \(t_j\), and \(c_{\text{pred},ij}\) represents the predicted drug concentration for person \(i\) at time \(t_j\), derived from the compartment pharmacokinetic model. The RMSE is calculated across all individuals and time points, where \(N\) is the total number of individuals and \(M\) represents the number of observation time points for each individual.
     The dynamic adjustment of the model resulted in enhanced robustness and improved accuracy in conditions of sparse sampling. Prior to adjustment, the mean percentage error (MPE) and root mean square error (RMSE) were 2.4959 mg/L and 2.7822 mg/L, respectively. Following the adjust–ment process, these metrics underwent a substantial reduction, reaching 1.3452 mg/L and 1.2673 mg/L, as outlined in Table 2.

     The incorporation of posterior sampling and dynamic prior adjustment has been demonstrated to enhance the precision of predicted drug concentrations. Integration of TDM data with HMC sampling facilitates more accurate estimation of pharmacokinetic parameters. The initial 100 sample results are illustrated in the subsequent figure:

     This figure compares the predicted vancomycin concentrations before and after dynamic prior adjustment using the Hamiltonian Monte Carlo (HMC) method. The black interpolation line represents the connecting line of the first 100 sample observations in the sparse treatment drug monitoring (TDM) database. It is noteworthy that the performance of 100 sample points is sufficient to reflect the accuracy of the dynamically adjusted Hamiltonian Monte Carlo (HMC) sampling method. The red crosses represent the predicted values before the adjustment, and the blue crosses represent the predicted values after the adjustment. As illustrated in the figure, the blue data points on the ordinate, representing the adjusted predicted values, exhibit a closer proximity to the actual observed values than the red data points do. This indicates that dynamically adjusted Hamiltonian Monte Carlo (HMC) sampling improves the estimation of patient-specific pharmacokinetic parameters, thereby enhancing the predictive performance under sparse data conditions.
     The integration of therapeutic drug monitoring (TDM) data with Hamiltonian Monte Carlo (HMC) sampling has demonstrated substantial improvements in individualized pharmacokinetic parameter estimation, such as clearance and volume. By dynamically adjusting the prior distribution, the model enhances its ability to account for inter-patient variability, leading to more accurate predictions of drug concentrations. The results show that the dynamic adjustments significantly improve the accuracy of dose recommendations, reducing prediction errors compared to non-adjusted models. This approach offers a robust framework for real-time personalized therapy, with potential for optimizing therapeutic outcomes and improving patient safety.

4.4 Personalized Drug Delivery System Design
     Implemented using R 4.2.1 with R stan integration, the developed decision support platform employs three core components: the Bayesian hierarchical model for pharmacokinetic parameter estimation; Hamiltonian Monte Carlo (HMC) sampling with dynamically adjustable prior strength for posterior sampling efficiency; and the individualized dosing algorithm. The computational core of the system lies in its dual calibration design, synchronizing prior distribution optimization with posterior sampling efficiency, which effectively handles non-linear relationships between creatinine clearance and drug elimination rates. The platform also incorporates real-time feedback mechanisms to adjust the model’s predictions as new patient data becomes available. This dynamic adaptability allows the system to remain accurate and responsive to changes in patient conditions. The pseudo-code design based on the model algorithm is as Algorithm 1:

     The individualized dosing system design consists of three main components: the basic patient information interface, the prior information input interface, and the individualized dosing recommendation interface. The basic patient information interface collects patient-specific data, such as name, age, body weight, and serum creatinine concentration. The prior information input interface includes the input of population pharmacokinetic (PK) values for clearance (CL), volume of distribution (V), dosing intervals, and initial dosing time.
     To demonstrate the system, the following simulation case is provided. Suppose the patient’s basic information is shown in the table below:

     The table shows the individualized dosing recommendations for a representative patient at four specific time points: 8:00, 14:00, 20:00, and 02:00 (next day). The recommended doses vary across time points, reflecting the adjustments made by the system based on the dynamic pharmacokinetic changes and therapeutic goals.

5. CONCLUSIONS

     Although sampling methods and prediction models have been applied in the medical field from various perspectives [26-28], MCMC (Markov Chain Monte Carlo) sampling provides a new perspective for individualized drug administration [29]. This study, based on the traditional MCMC, introduces the HMC (Hamiltonian Monte Carlo) sampling method with adjustable prior strength. By dynamically adjusting the prior strength, we can optimize the prior distribution based on the observed data, reduce prior bias, and improve model accuracy. Compared with previous methods, the improved HMC sampling has significantly enhanced accuracy and reliability. This method not only strengthens the individualized drug administration model for vancomycin but also provides a more effective strategy for future individualized drug therapy.

ACKNOWLEDGEMENTS

     The project was supported by the Provincial-level Undergraduate Innovation and Entrepreneurship Training Program Project at Guilin University of Electronic Technology, entitled "Design and Implementation of Individualized Drug Dosage Regimen Based on Bayesian Hierarchical Model — A Case Study of Vancomycin" (Project No. S202310595317). This project is funded by the Science and Technology Project of Guangxi (Grant Number Guike AD23023002). This project is funded by the Guangxi Higher Education Reform Project (Grant Number 2024JGA182). This research is funded by Course Construction Project of GUET Graduate Education (Grant Number YKC202502). This project is funded by Innovation Project of Guangxi Graduate Education (Grant Number YCSW2025371).

AUTHOR CONTRIBUTIONS

Daxin Gui: Conceptualization, Methodology, Software, Formal analysis, Writing - Original draft preparation. Pengcheng Ma: Conceptualization, Methodology, Resources, Supervision, Funding acquisition, Writing - Original draft preparation, Writing - Reviewing and Editing. Jing Chao: Investigation, Validation, Supervision, Writing - Reviewing and Editing. Haobang Huang: Data curation, Visualization. Yuyong Tan: Investigation, Validation. Xintai Guo: Software, Formal analysis. Furong Yang: Visualization, Project administration.

CONFLICT OF INTEREST STATEMENT

     The authors declare that they hold no competing interests.

FUNDING

This research was financially supported by the Provincial-level Undergraduate Innovation and Entrepreneurship Training Program Project at Guilin University of Electronic Technology (Project No. S202310595317).
This project is funded by the Science and Technology Project of Guangxi (Grant Number Guike AD23023002).
This project is funded by the Guangxi Higher Education Reform Project (Grant Number 2024JGA182).
This research is funded by Course Construction Project of GUET Graduate Education (Grant Number YKC202502).
This project is funded by Innovation Project of Guangxi Graduate Education (Grant Number YCSW2025371).

REFERENCES

[1] Sheiner L.B., Rosenberg B. and Marathe V.V., Estimation of population characteristics of pharmacokinetic parameters from routine clinical data. Journal of Pharmacokinetics and Biopharmaceutics, 1977; 5(5): 445–479. DOI 10.1007/BF01061728.

[2] Zhu W.B., Yang L.L. and He Z.K., HPLC determination of vancomycin in human plasma. Chinese Journal of Pharmaceutical Analysis, 2006; 26(8): 1106–1108. DOI 10.16155/j.0254-1793.2006.08.025.

[3] Frazee E.N., Kuper P.J., Schramm G.E., Larson S.L., Kashani K.B., Osmond D.R., et al., Effect of continuous venovenous hemodiafiltration dose on achievement of adequate vancomycin concentrations. Antimicrobial Agents and Chemotherapy, 2012; 56(12): 6181–6185. DOI 10.1128/AAC.00459-12.

[4] Rybak M.J., Lomaestro B., Rotschafer J.C., Moellering J.R., Craig W.A., Billeter M., et al., Therapeutic monitoring of vancomycin in adult patients: A consensus review of the American Society of Health-System Pharmacists, the Infectious Diseases Society of America, and the Society of Infectious Diseases Pharmacists. American Journal of Health-System Pharmacy, 2009; 66(1): 82–98. DOI 10.2146/ajhp080434.

[5] Tobler A. and Mühlebach S., Intravenous phenytoin: A retrospective analysis of Bayesian forecasting versus conventional dosing in patients. International Journal of Clinical Pharmacy, 2013; 35(5): 790–797. DOI 10.1007/s11096-013-9809-5.

[6] Yamamoto M., Kuzuya T., Baba H., Yamada K. and Nabeshima T., Population pharmacokinetic analysis of vancomycin in patients with gram-positive infections and the influence of infectious disease type. Journal of Clinical Pharmacy and Therapeutics, 2009; 34(4): 473–483. DOI 10.1111/j.1365-2710.2008.01016.x.

[7] Broeker A., Nardecchia M., Klinker K., Derendorf H., Day R., Marriott D., et al., Towards precision dosing of vancomycin: A systematic evaluation of pharmacometric models for Bayesian forecasting. Clinical Microbiology and Infection, 2019; 25(10): e1–e7. DOI 10.1016/j.cmi.2019.02.029.

[8] Edginton A.N., Schmitt W., Voith B. and Willmann S., A mechanistic approach for the scaling of clearance in children. Clinical Pharmacokinetics, 2006; 45(7): 683–704. DOI 10.2165/00003088-200645070-00004.

[9] Li Y., Jiao X., Sun G., Wang F., Wu X., Dong W., et al., Population pharmacokinetics and dosing optimization of norvancomycin for Chinese patients with community-acquired pneumonia. Infection and Drug Resistance, 2024; 17: 5881–5893. DOI 10.2147/IDR.S496776.

[10] Khangtragool A., Santidherakul S. and Leesawat P., Stability of vancomycin 31 mg/mL in extemporaneous eye drops determined with capillary electrophoresis. Chiang Mai Journal of Science, 2011; 38(4): 533–540.

[11] Li J.J., Liu Y.X., Tang L., Weng X.H., Wang S.N., Jiao Z., et al., Population pharmacokinetics of vancomycin in Chinese neonates. Chinese Pharmaceutical Journal, 2017; 52(16): 1434–1441. DOI 10.11669/cpj.2017.16.012.

[12] Lindsey J. and Jones B., Modeling pharmacokinetic data using heavy-tailed multivariate distributions. Journal of Biopharmaceutical Statistics, 2000; 10(3): 369–381. DOI 10.1081/BIP-100102500.

[13] Gao Y.C., Jiao Z., Huang H., Xie C., Gao J.J., Zhang L., et al., Development of decision system for individualization of vancomycin dosage. Acta Pharmaceutica Sinica, 2018; 53(1): 104–110. DOI 10.16438/j.0513-4870.2017-0673.

[14] Karlsson M. and Savic R., Diagnosing model diagnostics. Clinical Pharmacology & Therapeutics, 2007; 82(1): 17–20. DOI 10.1038/sj.clpt.6100241.

[15] Thammakornk P., Srion A., Chokvivat W., Hemstapat R., Morales N.P. and Suwannateeb J., Development of polycaprolactone infiltrated anti-tuberculosis drug-loaded 3D-printed hydroxyapatite for localized and sustained drug release in bone and joint tuberculosis treatment. Chiang Mai Journal of Science, 2022; 49(1): 105–121. DOI 10.12982/CMJS.2022.009.

[16] Mohammad R.E.A., Gubartallah E.A., Elbashir A.A., Saad B., Yahaya N. and N.N.M., et al., Simultaneous determination of nonsteroidal anti-inflammatory drugs using vortex-assisted liquid-liquid microextraction with back extraction coupled with HPLC-UV. Chiang Mai Journal of Science, 2021; 48(5): 1363–1373.

[17] Lunn D.J., Best N., Thomas A., Wakefield J. and Spiegelhalter D., Bayesian analysis of population PK/PD models: General concepts and software. Journal of Pharmacokinetics and Pharmacodynamics, 2002; 29(3): 271–307. DOI 10.1023/A:1020206907668.

[18] Wang F., Hu Y. and Qin Y., Spatio-temporal prediction of air quality using spatio-temporal clustering and hierarchical Bayesian model. Chiang Mai Journal of Science, 2024; 51(5): e2024083. DOI 10.12982/CMJS.2024.083.

[19] Hastings W.K., Monte Carlo sampling methods using Markov chains and their applications. Biometrika, 1970; 57(1): 97–109. DOI 10.1093/biomet/57.1.97.

[20] Durieux P., Trinquart L., Colombet I., Nies J., Walton R., Rajeswaran A., et al., Computerized advice on drug dosage to improve prescribing practice. Cochrane Database of Systematic Reviews, 2008; 3: CD002894. DOI 10.1002/14651858.CD002894.pub2.

[21] Chen S.W., Sampling Method Design Based on PPK Model: Applied to Individualized Administration of Vancomycin, PhD Thesis, Guilin University of Electronic Technology, China, 2023.

[22] Turner B.M. and Van Zandt T., A tutorial on approximate Bayesian computation. Journal of Mathematical Psychology, 2012; 56(2): 69–85. DOI 10.1016/j.jmp.2012.02.005.

[23] Tahir M., Aslam M., Hussain Z. and Khan A.A., On the 3-component mixture of exponential, Rayleigh and Burr type-XII distributions: A simulation study in Bayesian framework. Chiang Mai Journal of Science, 2018; 45(2): 1161–1180.

[24] Wei L.S. and Zhang W.P., Bayesian Analysis, University of Science and Technology of China Press, Hefei, 2013.

[25] Ding J.J., Jiao Z. and Wang Y., Estimation of individual pharmacokinetic parameters using maximum a posteriori Bayesian method with D-optimal sampling strategy. Acta Pharmaceutica Sinica, 2011; 46(12): 1493–1500. DOI 10.16438/j.0513-4870.2011.12.015.

[26] Pachthong C., Supyen D., Buddhasukh D. and Jatisatienr A., Isolation and characterization of brassinolide and castasterone in the pollen of pumpkin. Chiang Mai Journal of Science, 2006; 33(1): 95–101.

[27] Liu J.S., Liang F. and Wong W.H., The multiple-try method and local optimization in Metropolis sampling. Journal of the American Statistical Association, 2000; 95(449): 1211–1234. DOI 10.1080/01621459.2000.10473908.

[28] Ratchanapun P., Kumsuk N., Thipo K., Warchalapupan P., Prediction models for shelf life of pumpkin crackers in different packages based on its moisture content. Chiang Mai Journal of Science, 2010; 37(3): 410–420.

[29] Deng C.H., Ji S.M., Wang S.Y., Li L., Zhou T.Y. and Lu W., Development of therapeutic drug monitoring software of vancomycin for adult patients. Chinese Pharmaceutical Journal, 2014; 49(10): 881–885. DOI 10.11669/cpj.2014.10.020.

Related Articles