Pharmacometrics modeling and simulation to assist phenobarbital dose optimization in dogs

Introduction: Phenobarbital (PB) remains the first-line treatment for canine epilepsy due to its efficacy, affordability, and favorable pharmacokinetics. However, its narrow therapeutic index and substantial interindividual variability necessitate therapeutic drug monitoring (TDM) and individualized dosing. This study aimed to develop and validate a population pharmacokinetic (popPK) model of PB in dogs to support model-informed precision dosing (MIPD) in clinical practice.

Methods: A total of 121 serum samples from 100 dogs receiving PB monotherapy at steady state were used to build the model. An external dataset comprising 53 samples from 50 dogs was used for validation. Modeling was performed using nonlinear mixed-effects (NLME) techniques in MonolixSuite 2023R1. Covariate analysis included age, sex, and body weight (WT). Model performance was assessed using goodness-of-fit plots, prediction-corrected visual predictive checks (pcVPC), and calculation of mean error (ME), mean relative error (MRE), and root mean square error (RMSE). Monte Carlo simulations were conducted to evaluate the probability of target attainment (PTA) under different dosing regimens.

Results: A one-compartment model with autoinductive clearance (CL) best described PB pharmacokinetics. WT and age were significant covariates on apparent clearance (CL/F). The model accurately predicted PB concentrations in the external dataset (ME = −0.08 mg/L, MRE = 0.04%, RMSE = 2.04%). Simulations identified optimal dosing regimens stratified by age and WT, including recommendations for loading doses to accelerate attainment of therapeutic concentrations.

Discussion: The validated popPK model enables individualized PB dosing in dogs, accounting for variability in WT and age. This approach supports the implementation of MIPD in veterinary practice and may improve therapeutic outcomes while minimizing toxicity.

1 Introduction

Pharmacometrics, also referred to as “the science of quantitative pharmacology” (1), is a multidisciplinary field that develops and applies mathematical and statistical models to characterize variability in drug exposure, response, and disease progression. It supports data-driven decision-making in both drug development and clinical pharmacotherapy. This discipline integrates principles of pharmacokinetics (PK), pharmacodynamics (PD), pharmaceutical sciences, therapeutics, and pathophysiology, using mathematical and computational tools to describe and predict drug–organism interactions and clinical outcomes (2). Such models provide a rational framework to anticipate variability in drug exposure and response across different dose levels and populations. The use of computational tools enhances the value of data obtained from clinical trials, therapeutic drug monitoring, and pharmacovigilance (20).

Phenobarbital (PB) is the first-line treatment for seizures in canine epilepsy due to its efficacy, low cost, and favorable pharmacokinetic properties, which allow dosing every 12–24 h, unlike other classic anticonvulsants such as phenytoin or carbamazepine (3). Following oral administration in dogs, PB is rapidly and almost completely absorbed (bioavailability >85%), reaching peak plasma concentrations between 4 and 8 h post-dose. Approximately 50% of PB binds to plasma proteins. Its primary elimination pathway is hepatic enzymatic biotransformation, with about 33% excreted unchanged in the urine (4). PB is also an auto-inducer of the microsomal CYP450 enzyme complex, which leads to a progressive reduction in its elimination half-life with chronic use (5). According to Thurman et al. (21), the mean elimination half-life is 46.3 ± 11.3 h after a single dose and 29.3 ± 4.6 h at steady state, typically reached after 3 weeks of daily dosing.

Due to its narrow therapeutic index (15–45 mg/L) and significant pharmacokinetic interindividual variability observed in dogs, therapeutic drug monitoring (TDM) of PB is strongly recommended to optimize dosing, maximizing the probability of therapeutic efficacy while minimizing the risks of toxicity including sedation, lethargy, ataxia and hepatotoxicity (6, 7).

Traditionally, therapeutic drug monitoring relies on trough concentration, which provides limited information for individual pharmacotherapeutic interpretation and clinical outcome prediction. Nonlinear mixed effects (NLME) models implemented in pharmacometrics analyses allow the use of sparse data in the development of tools based on the population approach. In this context, population pharmacokinetic parameters can be estimated from sparse observations performed in patients belonging to the same population, and sources of interindividual and interoccasion variability can be recognized and described in a quantitative fashion. Furthermore, population pharmacokinetic models can be employed in the clinical setting to obtain maximum a posteriori estimates of individual parameters using sparse observations from the patient and the a priori distribution provided by the model. Bayesian forecasting is increasingly implemented within the context of model-informed precision dosing (MIPD) (8, 9).

Previous population pharmacokinetic studies in humans have identified several factors that significantly influence phenobarbital C, including total body weight, body surface area, and concomitant use of other antiepileptic drugs (10). To date, however, such models have not been applied to support dose optimization of phenobarbital in dogs. In this study, we aimed to develop a population pharmacokinetic model for phenobarbital in canine patients, evaluate current dosage regimens, and propose individualized dosing strategies based on patient-specific characteristics. Our goal was to provide a practical tool to support evidence-based dosing adjustments in clinical veterinary practice.

2 Materials and methods

2.1 Patients and data collection

A total of 121 serum samples were obtained from 100 client-owned dogs receiving phenobarbital monotherapy at steady state, defined as at least 20 days after treatment initiation. These data were used to develop the popPK model.

An independent validation dataset, consisting of 53 serum samples from 50 dogs, was also available. Dogs included in both datasets varied in age, body weight, sex, and neuter status; no breed restrictions were applied, and breed was not considered as a covariate in the analysis. All samples were collected as part of routine therapeutic drug monitoring with informed owner consent, under a protocol approved by the institutional Animal Use Ethics Committee.

Dogs received commercial oral phenobarbital tablets (100, 40 or 15 mg), administered according to the prescribing veterinarian’s instructions. Dosing was typically based on whole or halved tablets, reflecting standard clinical practice.

To assess demographic comparability between the training and validation populations, sex distribution was evaluated using the chi-square test, and continuous variables (age and BW) were compared using unpaired t-tests. No statistically significant differences were observed (p > 0.05), indicating that the two populations were demographically comparable. A summary of baseline characteristics is presented in Table 1.

Table 1

Table 1. Demographic characteristics of the canine population undergoing monotherapy with phenobarbital at steady state.

2.2 Phenobarbital determination

Phenobarbital serum concentrations were determined using the ARCHITECT® iPhenobarbital assay (Abbott), a chemiluminescent microparticle immunoassay (CMIA) performed on the ARCHITECT i1000SR analyzer. The assay’s validated analytical range is 1.1–80.0 μg/mL, with a lower limit of quantification (LLOQ) of 1.1 μg/mL. Intra- and inter-assay coefficients of variation (CV) were below 4%, and accuracy was within ±10% across the measurement range. Analyses were conducted at the Therapeutic Drug Monitoring Service, Department of Pharmaceutical Sciences, University Hospital, Universidad de la República (Uruguay), in accordance with the manufacturer’s validation standards.

2.3 Pharmacokinetic modeling

2.3.1 Structural model

Phenobarbital serum concentrations were analyzed using NLME implemented in MonolixSuite 2023R1 (Lixoft, a Simulations Plus Inc., company) to estimate the typical CL/F of PB and to perform the covariate analysis. Due to the sparse nature of the dataset, a one-compartment disposition model was selected. The absorption rate constant (ka) and the apparent volume of distribution (Vd/F) were fixed at 0.61 h⁻¹ and 20.0 L, respectively, based on previously published data (3, 5, 11). Both linear and autoinductive CL models were evaluated.

Interindividual variability (IIV) in CL/F was described using a log-normal distribution, as shown in Equation 1:

$\begin{array}{ll}\mathit{CL}/F_i=\mathit{CL}/F_{\mathit{pop}}\ast e^{{\eta }_i}& (1)\end{array}$

Where $\mathit{CL}/F_i$ and $\mathit{CL}/F_{\mathit{pop}}$ represent the ith individual estimate and the typical (population) parameter, respectively. The random effects ${\eta }_i$ describe the ith individual’s deviation from $\mathit{CL}/F_{\mathit{pop}}$ with a normal distribution: ${\eta }_i~N(0,{\omega }^2),$ being ${\omega }^2$ the IIV variance.

Residual unexplained variability was modeled using an additive error model, as described in Equation 2:

$\begin{array}{ll}{C}_{i,j}=C_{\mathit{pre}{d}_{i,j}}+\mathrm{a}\ast {\epsilon }_{i,j}& (2)\end{array}$

In this equation, C_ij represents the observed concentration for the ith subject at time j, Cpred_ij is the corresponding model-predicted concentration, a is the estimated additive residual error, and ε_ij is a standard normal random variable with mean 0 and variance 1.

2.3.2 Model evaluation

The model-building process was guided by diagnostic metrics and graphical diagnostics. The corrected Bayesian Information Criterion (BICc) computed from the estimated log-likelihood was used to optimize model parsimony, assessing the tradeoff between data fit and model complexity. Uncertainty in parameter estimation, calculated with the Fisher Information Matrix in Monolix® via linearization and expressed as relative standard error (RSE) was also considered as an essential diagnostic, along with plausibility of estimates. Basic goodness-of-fit plots included observed concentrations versus individual and population predictions, residuals plotted against time and predicted concentrations, and the distributions of both residuals and random effects. In addition, simulation-based diagnostics were performed, including prediction-corrected visual predictive checks (pcVPC) and normalized prediction distribution errors (NPDE).

2.3.2.1 Covariates analysis

The impact of age (years), sex, and body weight (WT), in kg on CL/F (L/h) was evaluated once the base structural model was developed.

Continuous covariates were incorporated using a power model with median centering, as described in Equation 3:

$\begin{array}{ll}\mathit{CL}/F_i=\mathit{CL}/F_{\mathit{pop}}\ast {\left(\frac{{\mathit{COV}}_i}{\mathit{CO}{V}_{\mathit{pop}}}\right)}^{{\beta }_{\mathrm{cov}}}{e}^{{\eta }_{{\theta }_i}}& (3)\end{array}$

Where ${\mathit{COV}}_i$ is the individual value of the continuous covariate for the ith animal, $\mathit{CO}{V}_{\mathit{pop}}$ the population median value (use for centering), ${\beta }_{\mathrm{cov}}$ represents the estimated effect of the covariate on CL/F ηᵢ the random effect representing unexplained IIV.

Categorical covariates such as sex a categorical covariate, were evaluated according an exponencial model, as shown in Equation 4:

$\begin{array}{ll}\mathit{CL}/F_i=\mathit{CL}/F_{\mathit{pop}}\ast e^{{\beta }_{\mathit{cat}}\ast {\mathit{sex}}_i}{e}^{{\eta }_{{\theta }_i}}& (4)\end{array}$

Where ${\beta }_{\mathit{sex}}$ SEX_i is coded as 1 for males and 0 for females, ${\beta }_{\mathit{sex}}$ s the estimated effect associated with sex, and η_i is the random effect accounting for IIV not explained by the covariate.

A forward selection method guided by graphical and statistical diagnostics was used to build the final model. The likelihood ratio test (LRT) and the Wald test were implemented to assess the statistical significance of each evaluated covariate relation. Covariate effects were retained in the model during the forward selection if its inclusion resulted in a reduction of the objective function value (OFV: −2*log-likelihood) by at least 3.84 units (p < 0.05) for the likelihood ratio test or LRT with one degree of freedom (12), and if the estimated parameter was significantly different than 0 (p < 0.05 in the Wald test). A backward deletion was then implemented to refine the full model based on the LRT, retaining significant effects with p < 0.01 (6.63-unit change in the objective function).

A stepwise approach was used to build the final model, beginning with forward selection and followed by backward elimination. The inclusion of a covariate was retained during forward selection if it reduced the objective function value (OFV = −2 × log-likelihood) by at least 3.84 units (p < 0.05) for the Likelihood Ratio Test (LRT) with one degree of freedom, and if the estimated parameter was significantly different from zero (p < 0.05 in the Wald test). In the backward elimination step, covariates were retained if their removal increased the OFV by at least 6.63 units (p < 0.01 for the LRT). Graphical and statistical diagnostics supported the selection process throughout.

2.3.3 External model validation

The validation dataset was used to assess the predictive performance of the final model. Both the a priori prediction (i.e., model predictions obtained for each individual using covariate data only) and the prediction obtained after Bayesian updating of pharmacokinetic parameters (i.e., after estimating individual parameters using individual characteristics and observations) were evaluated with goodness of fit plots and metrics. Mean error (ME, Equation 4), mean relative error (MRE, Equation 2), and relative root mean squared error (RMSE, Equation 3) were calculated to assess prediction bias and precision, as shown in Equations 5–7:

$\begin{array}{ll}\mathit{ME}\left(\frac{\mathit{mg}}{L}\right)=\frac{\sum (\mathit{\text{Cpred}}-\mathit{\text{Cobs}})}{n}& (5)\end{array}$

$\begin{array}{ll}\mathit{MRE}(\%)=\frac{\sum \left(\frac{\mathit{\text{Cpred}}-\mathit{\text{Cobs}}}{\mathit{\text{Cobs}}}\right)}{n}.100& (6)\end{array}$

$\begin{array}{ll}\mathit{\text{RMSE}}(\%)=\sqrt{\frac{\sum {\left(\frac{\mathit{\text{Cpred}}-\mathit{\text{Cobs}}}{\mathit{\text{Cobs}}}\right)}^2}{n}.100}& (7)\end{array}$

Where Cpred is the concentration predicted by the model, Cobs is the observed concentration and n is the number of observations.

2.4 Simulation assessment for different dose regimens

Monte Carlo simulations were performed using the final model in Simulx 2023R1 (Lixoft SAS, a Simulations Plus company) to evaluate PB pharmacokinetics under various dosing scenarios. A total of 1,000 individual profiles were simulated to explore different dosage regimens and estimate the probability of target attainment (PTA), defined as the proportion of individuals achieving steady-state concentrations within the therapeutic range. Therapeutic concentrations were defined as those between 15 and 45 μg/mL, based on published clinical guidelines (4). A PTA of 90% was considered an acceptable threshold for regimen selection. Additionally, the inclusion of a loading dose was evaluated to accelerate the achievement of therapeutic concentrations.

3 Results

A one-compartment model with an elimination process that is enhanced over time due to autoinduction was found to best describe PB concentration – time data. Autoinduction models have been proposed by Smythe et al. (13), and in Hassan et al. (14) to describe drugs that stimulate its own metabolism. In this model, CL/F increases over time as PB induces the expression of metabolic enzymes responsible for its own elimination.

The system is described by the following differential equations:

$\begin{array}{ll}\frac{\mathit{dAc}}{\mathit{dt}}=\mathit{ka}-\frac{\mathrm{C}L}{V}\times \mathit{Enz}.\mathit{Ac}& (8)\end{array}$

$\begin{array}{ll}\mathit{\text{dEnz}}=\mathit{ken}{z}_{\mathit{\text{in}}}-\mathit{ken}{z}_{\mathit{\text{out}}}\times (1-\frac{\mathit{Cc}}{\mathit{Cc}+I{C}_{50}})& (9)\end{array}$

The time-course of the system is described by two differential equations. Equation 8 models the amount of drug in the central compartment (Ac), which decreases over time due to CL modulated by the level of metabolizing enzymes (Enz). Equation 9 describes the dynamics of the enzyme pool. Although Kenz_in and Kenz _out are distinguished to reflect the mechanistic interpretation of enzyme synthesis and degradation, the model assumes that both parameters are equal (Kenz_in = Kenz_out = Kenz_pop), implying steady-state enzyme levels in the absence of drug.

This simplification allows the autoinductive effect of phenobarbital (PB) to be represented as an inhibition of enzyme degradation, governed by a sigmoidal function in which IC₅₀ represents the PB concentration that reduces enzyme degradation by 50%. As PB concentrations increase, enzyme degradation slows down, resulting in progressive enzyme accumulation over time. Since these enzymes mediate PB metabolism, their accumulation leads to a time-dependent increase in apparent CL/F. This dynamic behavior, captured by the model, aligns with clinical observations of reduced PB exposure over time despite constant dosing (Figure 1).

Figure 1

Figure 1. Schematic representation of the autoinduction model with enzyme turn over phenobarbital CL model. Ka = absorption rate constant, CL/V = elimination rate constant, Kenz,_in = enzyme production constant rate, Kenz,_out = enzyme degradation constant rate.

The apparent volume of distribution (Vd/F) was fixed to a literature value and scaled linearly with WT according to Equation 10.

$\begin{array}{ll}V=\mathit{\text{Vpop}}.\left(\frac{\mathit{WT}}{20}\right)& (10)\end{array}$

Where 20 kg represents the weighted mean WT in the study population.

Both age and WT significantly influenced CL/F, while sex was not a significant covariate (Figure 2, generated using R version 4.5.0 (22)). Incorporating WT into the model reduced the objective function value (OFV) by 62 units. Adding age led to a further decrease of 8 units. CL was significantly influenced by both body weight and age, while sex showed no significant effect. Inclusion of body weight using allometric scaling with a fixed exponent of 0.75 reduced the objective function value (OFV) by 62 units. Adding age as a covariate resulted in an additional reduction of 8 units.

Figure 2

Figure 2. Variations of CL (L/h) with age (years) and weight (kg) in different therapeutic range based on age and weight for the canine population undergoing monotherapy with phenobarbital at steady state.

Equation 11 represents the final model equation for CL/F:

$\begin{array}{ll}\mathit{CL}/F=\mathit{\text{CLpop}}{\left(\frac{\mathit{WT}}{20}\right)}^{3/4}.{\left(\frac{\mathit{AGE}}{5}\right)}^{{\beta }_{\mathit{AGE}}}& (11)\end{array}$

Where 20 kg and 5 years correspond to the population median values used for centering.

The final population pharmacokinetic model characterized the disposition of PB in dogs receiving monotherapy at steady state. The estimated model parameters are summarized in Table 2 and illustrated in Figure 3.

Table 2

Table 2. Estimated population parameters for the canine population undergoing monotherapy with phenobarbital at steady state (final model).

Figure 3

Figure 3. Goodness-of-fit and residual diagnostic plots for the final model applied to the training dataset (left column) and to the external validation dataset (right column) in dogs receiving phenobarbital monotherapy at steady state. (A) Top left: Observed concentrations versus individual predictions (OBS vs. PRED) for the training dataset. (B) Top right: OBS vs. PRED for the validation dataset. (C) Bottom left: Normalized prediction distribution errors (NPDE) versus population predictions for the training dataset. (D) Bottom right: NPDE versus population predictions for the validation dataset. In the OBS vs. PRED plots (top row), the black dashed line represents the line of identity, while the dotted lines indicate the 90% prediction interval. The bottom row shows residual diagnostics, with orange loess lines indicating trends in residuals, and the horizontal gray line representing the zero residual level.

Parameter precision was assessed using the Fisher Information Matrix, and the results indicated reliable estimates.

External validation demonstrated that the model adequately predicted observed concentrations, with a mean error (ME) of −2.10 mg/L, a mean relative error (MRE) of −1.74, and a relative root mean squared error (RMSE) of 24.80% (Figure 3).

The predictive performance of the final model was explored using a simplified visual predictive check (VPC) adapted in R version 4.5.0 (22), focusing on predose concentrations. Figure 4 presents observed serum concentrations at 12 h post-dose.

Figure 4

Figure 4. Simplified visual predictive check (VPC) for phenobarbital predose concentrations.

Dosage regimens derived from simulation-based analysis using the final model are summarized in Table 3 and illustrated in Figure 5, highlighting the influence of age and weight on phenobarbital exposure and supporting the rationale for stratified dosing. The resulting recommendations aim to guide individualized phenobarbital therapy in dogs to rapidly and consistently achieve therapeutic concentrations.

Table 3

Table 3. Dosage regimen recommendations to attain a therapeutic range of phenobarbital based on the age and weight of the individual dogs.

Figure 5

Figure 5. Simulated phenobarbital concentration-time profiles (n = 1,000) across nine age and weight strata in dogs receiving optimized dosing regimens. Each panel represents a subgroup based on weight (kg) and age (years), as defined in Table 3. The solid blue line denotes the median predicted serum concentration, and the shaded area represents the 90% prediction interval. Red horizontal lines indicate thresholds at 15, 30, and 45 μg/mL.

4 Discussion

Variability in drug disposition, both between and within individuals, poses a major challenge in the treatment of chronic conditions such as epilepsy. In this study, a popPK model was developed to characterize PB concentrations in dogs receiving monotherapy, based on real TDM data. The model demonstrated adequate predictive performance, as confirmed by external validation, and identified relevant covariates that explain a significant portion of the observed variability, supporting its clinical applicability.

The selected structural model was a one-compartment model with an elimination process subject to autoinduction, in which CL/F increases progressively over time. This type of model, previously proposed by Smythe et al. (13) and Hassan et al. (14), simulates the upregulation of metabolizing enzymes as a result of their degradation being inhibited by the drug. The IC₅₀ parameter represents the PB concentration that reduces enzyme degradation rate by 50% and plays a central role in capturing the system’s dynamics and adequately simulating progressive enzyme induction. This mechanism helps explain the plasma concentration profiles observed after treatment initiation or dose adjustments with PB. Initially, concentrations rises due to drug accumulation, but then progressively decline as PB induces the expression of metabolizing enzymes. This increase in CL stabilizes over time, such that at steady state, plasma concentrations remain at levels lower than those reached during the first days of treatment. This behavior contrasts with linear pharmacokinetics, where concentrations gradually rise until stabilizing at a level proportional to the administered dose, typically within 5–7 elimination half-lives. The autoinductive structural model captured this temporal pattern accurately, whereas linear models, though statistically acceptable, failed to represent the observed dynamic or provide a mechanistic explanation. Moreover, the autoinductive model yielded better fit metrics, supporting its selection as a more realistic predictive tool. Importantly, the use of a one-compartment model is consistent with both biological plausibility and previous evidence. A systematic review by Methaneethorn and Leelakanok (15), which evaluated 18 pharmacokinetic studies of phenobarbital in humans, found that all employed one-compartment models, even in the presence of rich sampling designs. This reinforces the empirical and translational robustness of our structural model selection.

In addition, chronic administration of antiepileptic drugs—including PB, phenytoin, and carbamazepine—has been associated with increased expression of efflux transporters such as P-glycoprotein (Pgp) in the brains of rats after 60 days of treatment (16). Efflux transporters are membrane proteins widely distributed in tissues, acting in concert with metabolic enzymes to regulate intracellular drug concentrations by actively exporting their substrates. Pgp is constitutively expressed in the liver, bile canaliculi, intestines, kidneys, blood–brain barrier, and other tissues (17) and significantly influences PB pharmacokinetics by limiting absorption and enhancing biliary and renal elimination. Alvariza et al. (18) further demonstrated that Pgp overexpression occurs in a concentration-dependent manner across multiple tissues, contributing to increased CL/F through a combined effect of enzyme induction and enhanced excretory CL.

The final model also incorporated significant effects of WT and age on CL/F. WT was accounted for using allometric scaling with a fixed exponent of 3/4, in line with well-established physiological principles. Age showed a negative effect on CL, consistent with age-related decline in renal function. In this regard, Laroute et al. (19) reported a 34% decrease in glomerular filtration rate between young and elderly Beagle dogs—a relevant finding given that approximately 25% of PB elimination occurs via renal excretion.

To assess the model’s clinical utility, Monte Carlo simulations were performed in Simulx® under different dosing schemes. A PTA ≥ 90% was considered an acceptable threshold for regimen selection, defined as the proportion of individuals achieving steady-state concentrations within the recommended therapeutic range (15–45 μg/mL). The simulations also included evaluation of loading doses to accelerate attainment of effective concentrations. According to the model, PB concentrations tend to stabilize after approximately 15 days of continuous treatment.

The results showed that dogs weighing less than 10 kg—and particularly those under 1 year of age—require higher doses to reach therapeutic levels, whereas older dogs may require lower doses. In clinical practice, PB is typically prescribed at doses ranging from 2.5 to 5.5 mg/kg, with body weight being the primary factor guiding dose selection Podell et al. (4). However, our findings suggest that age also plays a relevant role in drug disposition and should be considered when individualizing therapy, especially in very young or elderly dogs.

Finally, although other factors such as breed or reproductive status may be clinically relevant, these were not included as covariates due to limitations in the available data and the objective of developing a broadly applicable model for diverse canine populations.

In summary, this study demonstrates the value of pharmacometric tools in optimizing veterinary pharmacotherapy. The proposed model enables individualized dosing strategies based on patient-specific characteristics such as weight and age, which may improve the safety and effectiveness of antiepileptic therapy with phenobarbital.

5 Conclusion

This study presents a population pharmacokinetic model that successfully characterizes the disposition of phenobarbital (PB) in dogs receiving monotherapy, incorporating key covariates such as body weight and age. The model not only demonstrated good predictive performance but also mechanistically captured the time-dependent changes in PB CL due to autoinduction. Based on this model, individualized dosing regimens were proposed, allowing for more precise and clinically relevant dose adjustments from the outset and throughout therapy. While current clinical practice primarily considers body weight when selecting PB doses, our findings highlight the importance of also accounting for age, especially in puppies and senior dogs—to optimize exposure. Additionally, the inclusion of loading dose strategies supports a faster achievement of therapeutic levels, which may be particularly beneficial in the management of acute epileptic conditions. Overall, this model provides a robust framework to support evidence-based, patient-tailored dosing of PB, with the potential to enhance both treatment efficacy and safety in canine epilepsy management.

Data availability statement

The data used in this study are not publicly available due to confidentiality restrictions. However, they may be provided by the corresponding authors upon reasonable request and subject to confidentiality agreements.

Ethics statement

The animal studies were approved by Ethics Committee on Animal Use – Faculty of Veterinary Medicine, Universidad de la República (UdelaR). The studies were conducted in accordance with the local legislation and institutional requirements. Written informed consent was obtained from the owners for the participation of their animals in this study.

Author contributions

SA: Writing – review & editing, Supervision, Conceptualization, Software, Formal analysis, Investigation, Methodology, Writing – original draft, Funding acquisition, Data curation, Project administration, Validation, Visualization. MI: Visualization, Data curation, Software, Writing – original draft, Methodology, Formal analysis, Validation, Writing – review & editing. NG: Resources, Investigation, Writing – review & editing, Supervision. CM: Supervision, Resources, Writing – review & editing. MV: Project administration, Writing – review & editing, Supervision, Resources. GF: Resources, Investigation, Writing – review & editing. GS: Project administration, Software, Writing – review & editing, Writing – original draft, Visualization.

Funding

The author(s) declare that financial support was received for the research and/or publication of this article. This study was funded by the Pharmacology and Therapeutics Unit, Faculty of Veterinary Medicine, and the Therapeutic Drug Monitoring Service, Pharmaceutical Sciences Department, Faculty of Chemistry, University Hospital, Montevideo, Uruguay. No external funding was received for this study.

Acknowledgments

The authors thank SimulationPlus for providing an academic license for MonolixSuite® software.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The authors declare that Gen AI was used in the creation of this manuscript. This manuscript benefited from language editing and structural feedback provided with the assistance of ChatGPT (GPT-4, OpenAI), under the supervision and full responsibility of the corresponding author.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. Ette, EI, and Williams, PJ. Pharmacometrics: the science of quantitative pharmacology. New York, NY: Wiley-Interscience (2007).

Google Scholar

2. Ibarra, M, Lorier, M, and Trocóniz, IF. Pharmacometrics: definition and history In: A Talevi, editor. The ADME encyclopedia. Cham: Springer (2021)

Google Scholar

3. Al-Tahan, F, and Frey, HH. Absorption kinetics and bioavailability of phenobarbital after oral administration to dogs. J Vet Pharmacol Ther. (1985) 8:205–7. doi: 10.1111/j.1365-2885.1985.tb00946.x

PubMed Abstract | Crossref Full Text | Google Scholar

4. Podell, M, Volk, HA, Berendt, M, Löscher, W, Muñana, KR, Patterson, EE, et al. 2015 ACVIM small animal consensus statement on seizure Management in Dogs. J Vet Intern Med. (2016) 30:477–90. doi: 10.1111/jvim.13841

PubMed Abstract | Crossref Full Text | Google Scholar

5. Hojo, T, Ohno, R, Shimoda, M, and Kokue, E. Enzyme and plasma protein induction by multiple oral administrations of phenobarbital at a therapeutic dosage regimen in dogs. J Vet Pharmacol Ther. (2002) 25:121–7. doi: 10.1046/j.1365-2885.2002.00385.x

PubMed Abstract | Crossref Full Text | Google Scholar

6. Boothe, DM, Perkins, J, and Carpenter, J. Pharmacokinetics of phenobarbital in dogs with idiopathic epilepsy. J Vet Pharmacol Ther. (1998) 21:85–91.

Google Scholar

7. Muñana, KR. Update: seizure management in small animal practice. Vet Clin North Am Small Anim Pract. (2013) 43:1127–47. doi: 10.1016/j.cvsm.2013.04.008

Crossref Full Text | Google Scholar

8. Darwich, AS, Polasek, TM, Aronson, JK, Ogungbenro, K, Wright, DFB, Achour, B, et al. Model-informed precision dosing: background, requirements, validation, implementation, and forward trajectory of individualizing drug therapy. Annu Rev Pharmacol Toxicol. (2021) 61:225–245. doi: 10.1146/annurev-pharmtox-033020-113257

Crossref Full Text | Google Scholar

9. Keizer, RJ, Ter Heine, R, Frymoyer, A, Lesko, LJ, Mangat, R, and Goswami, S. Model-informed precision dosing at the bedside: scientific challenges and opportunities. CPT Pharmacometrics Syst Pharmacol. (2018) 7:785–7. doi: 10.1002/psp4.12353

Crossref Full Text | Google Scholar

10. Teixeira- da-Silva, P, Santos-Buelga, D, Otero, M, and Garcia, M. Population pharmacokinetics of phenobarbital in Caucasian patients with epilepsy. Eur J Pharm Sci (2020); 153:105484. doi: 10.1016/j.ejps.2020.105484

Crossref Full Text | Google Scholar

11. Ravis, WR, Nachreiner, RF, Pedersoli, WM, and Houghton, NS. Pharmacokinetics of phenobarbital in dogs after multiple oral administrations. Am J Vet Res. (1984) 45:1283–6.

Google Scholar

12. Ooi, QX, Wright, DFB, Isbister, GK, and Duffull, SB. Evaluation of assumptions underpinning pharmacometric models. AAPS J. (2019) 21:97. doi: 10.1208/s12248-019-0366-2

PubMed Abstract | Crossref Full Text | Google Scholar

13. Smythe, W, Khandelwal, A, Merle, C, Rustomjee, R, Gninafon, M, Lo, MB, et al. A semimechanistic pharmacokinetic-enzyme turnover model for rifampin autoinduction in adult tuberculosis patients. Antimicrob Agents Chemother. (2012) 56:2091–8. doi: 10.1128/AAC.05792-11

PubMed Abstract | Crossref Full Text | Google Scholar

14. Hassan, M, Svensson, USH, Ljungman, P, Björkstrand, B, Olsson, H, Bielenstein, M, et al. A mechanism-based pharmacokinetic-enzyme model for cyclophosphamide autoinduction in breast cancer patients. J Clin Pharmacol. (1999) 48:669–77.

Google Scholar

15. Methaneethorn, J, and Leelakanok, N. Pharmacokinetic variability of phenobarbital: a systematic review of population pharmacokinetic analysis. Eur J Clin Pharmacol. (2021) 77:291–309. doi: 10.1007/s00228-020-03011-x

PubMed Abstract | Crossref Full Text | Google Scholar

16. Yang, H, Liu, H, Liu, X, Zhang, D, Liu, Y, Liu, X, et al. Increased P-glycoprotein function and level after long-term exposure of four antiepileptic drugs to rat brain microvascular endothelial cells in vitro. Neurosci Lett. (2008) 434:299–303. doi: 10.1016/j.neulet.2008.01.071

PubMed Abstract | Crossref Full Text | Google Scholar

17. Staud, F, Ceckova, M, Micuda, S, and Pavek, P. Expression and function of P-glycoprotein in Normal tissues: effect on pharmacokinetics. Methods Mol Biol. (2010) 596:133–72. doi: 10.1007/978-1-60761-416-6_10

Crossref Full Text | Google Scholar

18. Alvariza, S, Fagiolino, P, Vázquez, M, Feria-Romero, I, and Orozco-Suárez, S. Chronic administration of phenytoin induces efflux transporter overexpression in rats. Pharmacol Rep. (2014) 66:946–51. doi: 10.1016/j.pharep.2014.06.007

PubMed Abstract | Crossref Full Text | Google Scholar

19. Laroute, V, Chetboul, V, Roche, L, Maurey, C, Costes, G, Pouchelon, JL, et al. Quantitative evaluation of renal function in healthy beagle puppies and mature dogs. Res Vet Sci. (2005) 79:161–7. doi: 10.1016/j.rvsc.2004.11.011

PubMed Abstract | Crossref Full Text | Google Scholar

20. Martín-Jiménez, T, and Riviere, JE. Population pharmacokinetics in veterinary medicine: potential use for therapeutic drug monitoring and prediction of tissue residues. J Vet Pharmacol Ther. (1998) 21:167–89. doi: 10.1046/j.1365-2885.1998.00121.x

PubMed Abstract | Crossref Full Text | Google Scholar

21. Thurman, GD, McFadyen, ML, and Miller, R. The pharmacokinetics of phenobarbitone in fasting and non-fasting dogs. J S Afr Vet Assoc. (1990) 61:86.

Google Scholar

22. R Core Team. (2025). R: A language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing.

Google Scholar