Valproic acid physiological pharmacokinetics simulation for epilepsy patients via a refined seven-compartmental model: A pilot study

Abstract

This study simulated valproic acid (VPA) physiological pharmacokinetics for epilepsy patients via a refined seven-compartmental model. The seven compartments were defined as oral, GI Tract, liver, whole body (WB), cerebrospinal fluid (CSF), kidney, and bladder. A first-order system of seven differential equations was defined according to this specific model and solved by a self-developed program run in MATLAB to evaluate the VPA time-dependent change among compartments. Each compartment's preset dissolving half-life (in hours) was as follows: oral 0.05, GI Tract 0.10, liver 0.5, WB 2.2, CSF 0.8, kidney 1.2, and bladder 0.5, respectively. The derived results were reorganized according to various presets of dissolving half-lives of the liver, WB, CSF, or kidney to analyze the VPA time-dependent changes under various scenarios. Either the liver or WB was the dominant compartment in this model, which mainly controlled the VPA changes. In contrast, others compartments followed the principle of secular equilibrium in the chain decay of radioactive nuclides. Accordingly, the VPA degradation changes in WB (“mother” compartment) mainly affected the VPA changes in other (“daughter”) compartments. Thus, the predicted VPA degradation in CSF, kidney, or bladder was always longer than in WB despite its dissolving half-life changes. The predicted results of VPA changes in various compartments were also compared with other studies, and a reasonable agreement was reached on whether the dissolving half-life of the liver or WB ranged from the original 0.2/2.2 to 1.4/0.9 h. The programable capability of this self-developed program allows one to easily modify the primary preset to comply with findings from other studies and has great potential in similar applications.

Keywords

Valproic acid pharmacokinetics model interpretation simulation self-developed program

1. Introduction

Valproic acid (VPA) physiological pharmacokinetics for epilepsy patients was simulated via a refined seven-compartmental model in this study. VPA, also known as 2-propylvaleric acid, in the form of sodium valproate, is a medication used to treat epilepsy and bipolar disorder and to prevent migraines. It prevents the absence of seizures, partial seizures, and generalized seizures. VPA can be administered intravenously or orally as well. It is included in the World Health Organization's Standard List of Essential Medicines and is one of the essential medicines for the basic public health system. The anticonvulsant effect of VPA is generally attributed to the blockade of sodium channels and the increase of gamma-aminobutyric acid (GABA) levels in the brain. However, the mechanism of sodium valproate is not fully understood. Sodium valproate increases the levels of the inhibitory neurotransmitter GABA in the brain and cerebellum in the animal test, most likely through the inhibition of GABA-degrading enzymes, such as GABA aminotransferase, succinate semialdehyde dehydrogenase and through the inhibition of GABA reuptake by neuronal cells.^1–6

Many researchers have tried to construct the pharmacokinetic model of VPA inside the human body from theoretical simulations or practical surveys. Ibarra et al.,⁷ Ogungbenro et al.,⁸ and Simeon et al.⁹ adopted the pharmacokinetic model with 4–10 compartments to define the mechanism of VPA inside the human body and solve the simultaneous differential equation group to evaluate the time-dependent VPA concentration in every compartment theoretically. In doing so, many essential parameters in solving the differential equation group were estimated or surmised from other published articles with similar topics. In contrast, the crucial challenge for verifying the derived pharmacokinetic model needs a more practical survey of the real human metabolic mechanism. Thus, the time-dependent VPA concentration in internal organs is too difficult to quantify. Sopasakis et al.¹⁰ and Azizi¹¹ developed a theory of fractional kinetics in the pharmacokinetics field and tried to find the numerical solution using the Riemann-Liouville integral with discrete-time finite memory approximations. In contrast, Jones and Rowland-Yeo¹² proposed a concept of Physiologically Based Pharmacokinetic (PBPK) modeling in their tutorial study, claiming that a clearance rate [h⁻¹] of essential organs, liver, and kidney dominated the time-dependent concentration of drugs inside the human body.

Accordingly, a concise multi-compartmental model was proposed in this study. The refined seven-compartmental model was built up to compose a close loop of the pharmacokinetics of VPA undertaken by patients. The corresponding system of 1^st-order differential equations was defined according to the compartmental model, while the essential parameters or correlated constants to derive the mathematical solution were mostly quoted from similar pharmacokinetic studies. The theoretical simulation complied with some preliminary results that were mostly verified by surveys of practical VPAs administered to patients. The defined system of differential equations was solved by MATLAB default program, where each parameter or constant could be adjusted to reflect various scenarios (discussed in the next sections). Furthermore, the model elaboration, definitions, theoretical presumption of parameters, and the potential application of this model were discussed in detail.

2. Theoretical simulation

2.1 Presumption of the program

Once the radioactivated solution is administered into the patient's body, it can be easily quantified by gamma camera from a radiological viewpoint. In contrast, the VPA time-dependent change in the patient's body can only be quantified by blood testing from various organs. Most researchers proceeded only with the theoretical simulation, whereas Ogungbenro et al.,⁸ Soars et al.,¹³ and Cloyd et al.^14,15 accomplished both theoretical simulation and practical blood tests in their studies. Thus, most essential parameters as preset in this simulated program were quoted from the above sources. Yet, the pharmacokinetic model was also revised based on their suggested prototype. However, as defined herein, the revised model was constructed layer by layer according to the preset principle of a system of differential equations for simplifying the correlation among compartments.^16–18

2.2 The structure of the multi-compartmental model

Assuming the human body is an all-inclusive vessel, internal organs can be categorized into several compartments, implying the importance of operating or transferring the VPA inside the human body. If so, whole body fluid (WB) can be defined as similar to the transient compartment to store or deliver the VPA after it is administered orally or injected. Accordingly, as implied in Figure 1, if the VPA is administrated orally, the first compartment is oral; the next is the gastrointestinal tract (GI Tract), where it is fully dissolved. The epithelial cells of the villi of the small intestine are absorbed into the microvasculature and enter the liver via the hepatic portal vein, then from the hepatic vein into the inferior vena cava, and finally through the heart to the cells of the whole-body fluid (WB), then split into three paths as to cerebrospinal fluid (CSF) (40%), kidney (20%) or feedback to GI Tract (40%) again. The kidney is transferred to the bladder and then excreted outside the body. Notably, CSF or kidneys have a feedback path back to WB to fulfill the human metabolic mechanism (CSF and kidney feedback to WB is preset as 50% and 20%, respectively). The three feedback paths (WB to GI Tract and CSF/kidney to WB) make the numerical analysis too complicated to calculate the VPA time-dependent change among compartments; however, with the assistance of a self-developed program run in MATLAB,¹⁹ it can easily derive the solution and fit the practical requirements in clinical observation.

Figure 1.

A seven-compartmental model adopted in this study. The first compartment is the oral, followed by the gastrointestinal tract (GI Tract). The epithelial cells of the villi of the small intestine are absorbed into the microvasculature and enter the liver via the hepatic portal vein, then from the hepatic vein into the inferior vena cava, and finally through the heart to the cells of the whole-body fluid (WB), then split into three paths as to cerebrospinal fluid (CSF), kidney, or feedback to GI Tract again. The kidney is transferred to the bladder and then excreted outside the body.

2.3 The definition of simultaneous differential equations group

A system of seven first-order simultaneous differential equations was proposed in this study according to the theoretical scenario described above. Defining decay constant (λ [h⁻¹]) is inversely proportional to a medicine or compartment dissolving half-life; thus, λ_R = Ln 2/T_1/2(radiological) and T_1/2(radiological) is the radiological half-life of VPA dissolving. The value is preset as 200 h in this study because it provided a neglectable contribution to affect the simulation. In contrast, if it is radiolabel medicine to fulfill special pharmaceutical needs in reality, the value can be manually changed to go with the radiological half-life of any specific isotope as I-131 or Tc-99 m.²⁰ Decay constant λ_x represents VPA dissolving in the specific compartment (x) in the theoretical scenario, and the accompanied half-life of VPA dissolving is as follows: 1. Oral (T_1/2 = 0.05 h), 2. GI tract (T_1/2 = 0.6 h), 3. Liver (T_1/2 = 0.5 h), 4. Whole body (WB) (T_1/2 = 2.2 h), 5. cerebrospinal fluid (CSF) (T_1/2 = 0.8 h), 6. Kidney (T_1/2 = 1.2 h), and 7. Bladder (T_1/2 = 0.5 h). Furthermore, i_xy implies the transmitted path from compartment x to y, and the suggested values of i₁₂, i₂₃, i₃₄, i₄₅, i₄₆, i₄₂, i₅₄, i₆₇, i₆₄, and i₇₈ are 1.0, 1.0, 1.0, 0.4, 0.2, 0.4, 0.5, 0.8, 0.2, and 1.0, respectively. Specifically, the branching ratio is preset, assuming only 40% (i₄₅ = 0.4) of dissolved VPA is transmitted from WB to CSF. Note that either CSF or kidney is also assigned the feedback path (i₅₄ = 0.5) or (i₆₄ = 0.2) to WB, whereas the remaining VPA inside the kidney totally transfers to the bladder (i₇₈ = 0.1) (cf. Figure 1). This scenario presumes that VPA is mainly absorbed by the CSF and partly delivered to the kidneys and bladder, followed by excretion. The liver is a relay organ that metabolizes and transmits the VPA to WB. The WB can give feedback to the GI Tract or transmit to CSF and kidney only. In addition, all the percentages can be adjusted to fulfill the practical examination from previous surveys or theoretical calculations in other studies.

\begin{aligned} \frac{{dq}_{1}}{dt} & = - (λ_{R} + λ_{12}) \cdot q_{1} \end{aligned}

(1)

\begin{aligned} \frac{{dq}_{2}}{dt} & = λ_{12} \cdot q_{1} + λ_{42} \cdot q_{4} - (λ_{R} + λ_{23}) \cdot q_{2} \end{aligned}

(2)

\begin{aligned} \frac{{dq}_{3}}{dt} & = λ_{23} \cdot q_{2} - (λ_{R} + λ_{34}) \cdot q_{3} \end{aligned}

(3)

\begin{aligned} \frac{{dq}_{4}}{dt} & = λ_{34} \cdot q_{3} + λ_{54} \cdot q_{5} + λ_{64} \cdot q_{6} - (λ_{R} + λ_{45} + λ_{46} + λ_{42}) \cdot q_{4} \end{aligned}

(4)

\begin{aligned} \frac{{dq}_{5}}{dt} & = λ_{45} \cdot q_{4} - (λ_{R} + λ_{54}) \cdot q_{5} \end{aligned}

(5)

\begin{aligned} \frac{{dq}_{6}}{dt} & = λ_{46} \cdot q_{4} - (λ_{R} + λ_{67} + λ_{64}) \cdot q_{6} \end{aligned}

(6)

\begin{aligned} \frac{{dq}_{7}}{dt} & = λ_{67} \cdot q_{6} - (λ_{R} + λ_{78}) \cdot q_{7} \end{aligned}

(7)

2.4 Running MATLAB program

Accordingly, Eqs (1)–(7) were reorganized following the input format of MATLAB to compute the numerical solution. In reality, the initial condition was preset at 0.0, implying that there was no residual dose before the patient undertook the first VPA dosage. The standard input dataset was integrated via Eq. (8). As implied, a system of differential equations with a dataset matrix [8 × 1] = [8 × 8]·[8 × 1]. In addition, all input parameters were incorporated as variables, not fixed values, to simplify the user input in later debugging or revised processes.²⁰

\begin{aligned} [\begin{array}{l} d q_{1} / d t \\ d q_{2} / d t \\ d q_{3} / d t \\ d q_{4} / d t \\ d q_{5} / d t \\ d q_{6} / d t \\ d q_{7} / d t \end{array}] \\ = [\begin{array}{ccccccc} - (λ_{R} + λ_{12}) & 0 & 0 & 0 & 0 & 0 & 0 \\ λ_{12} & - (λ_{R} + λ_{23}) & 0 & λ_{42} & 0 & 0 & 0 \\ 0 & λ_{23} & - (λ_{R} + λ_{34}) & 0 & 0 & 0 & 0 \\ 0 & 0 & λ_{34} & - (λ_{R} + λ_{45} + λ_{46} + λ_{42}) & λ_{54} & λ_{64} & 0 \\ 0 & 0 & 0 & λ_{45} & - (λ_{R} + λ_{54}) & 0 & 0 \\ 0 & 0 & 0 & λ_{46} & 0 & - (λ_{R} + λ_{64} + λ_{67}) & 0 \\ 0 & 0 & 0 & 0 & 0 & λ_{67} & - (λ_{R} + λ_{78}) \end{array}] \\ \times [\begin{array}{l} q_{1} \\ q_{2} \\ q_{3} \\ q_{4} \\ q_{5} \\ q_{6} \\ q_{7} \end{array}] \end{aligned}

(8)

Table 1.

The program's essential parameters simulate various scenarios, such as standard, liver, kidney, whole-body disorder, changing path ratio into compartments, or transmitting to VPA concentration.

	Scenario
Parameters $(L n 2 / T_{1 / 2} (x)) [h r^{- 1}];$ X:	I: Standard case	II: liver disorder (A)	III: liver disorder (B)	IV: WB disorder (A)	V: WB disorder (B)	VI: CSF disorder (A)	VII CSF disorder (B)	VIII: kidney disorder (A)	IX: kidney disorder (B)	X: various paths compose (A)	XI: various paths compose (B)	XII: various paths compose (C)
Radiological	200	200	200	200	200	200	200	200	200	200	200	200
Oral	0.05	0.05	0.05	0.05	0.05	0.05	0.05	0.05	0.05	0.05	0.05	0.05
GI Tract	0.10	0.10	0.10	0.10	0.10	0.10	0.10	0.10	0.10	0.10	0.10	0.10
Liver	0.5	1.0	3.0	3.0	3.0	0.5	0.5	3.0	3.0	0.5	0.5	0.5
Whole body (WB)	2.2	2.2	2.2	1.1	4.5	2.2	2.2	2.2	2.2	2.2	2.2	2.2
CSF	0.8	0.8	0.8	0.8	0.8	0.6	1.2	0.8	0.8	0.8	0.8	0.8
Kidney	1.2	1.2	1.2	1.2	1.2	1.2	1.2	1.8	3.0	1.2	1.2	1.2
Bladder	0.5	0.5	0.5	0.5	0.5	0.5	0.5	0.5	0.5	0.5	0.5	0.5
Path ratio
i12, oral to GI Tract	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0
i23, GI Tract to Liver	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0
i34, Liver to WB	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0
i45, WB to CSF	0.4	0.4	0.4	0.4	0.4	0.4	0.4	0.4	0.4	0.3	0.5	0.4
i46, WB to kidney	0.2	0.2	0.2	0.2	0.2	0.2	0.2	0.2	0.2	0.1	0.2	0.3
i42; WB GI Tract	0.4	0.4	0.4	0.4	0.4	0.4	0.4	0.4	0.4	0.6	0.3	0.3
i54, CSF to WB	0.5	0.5	0.5	0.5	0.5	0.5	0.5	0.5	0.5	0.5	0.5	0.5
i67, Kidney to bladder	0.8	0.8	0.8	0.8	0.8	0.8	0.8	0.8	0.8	0.8	0.8	0.8
i64, Kidney to WB	0.2	0.2	0.2	0.2	0.2	0.2	0.2	0.2	0.2	0.2	0.2	0.2
i78, Bladder to excretion	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0

Bold represents the significance difference in terms of numbers in the same two groups with the same disorders.

3. Simulation results

Table 1 implies all essential parameters as adopted in the program to simulate various scenarios as standard case (I), liver (II, III), WB (IV, V), CSF (VI, VII), or kidney (VIII, IX) disorder, as well as changing path ratio from WB into three compartments (CSF, Kidney or GI Tract) (X, XI, XII). As shown in the first column, every scenario included seven dissolving half-lives of specific compartments in the program. The preset value was determined from the previous studies or theoretical assumptions. The computational result of the standard scenario was automatically plotted by MATLAB, as depicted in Figure 2. As clearly illustrated, the VPA time-dependent change of WB attributed to the largest part and other compartments also have their unique characteristics, as shown in this model. In contrast, the eight disorder scenarios were preset as either small or large dissolving half-life of the specific compartment to imply a fast or slow metabolic mechanism. For instance, the half-life of the liver changed from 0.5 to 2 or 3 h. to represent a slow metabolic mechanism of the liver (scenario I, II, or III in this study). The computation process also defined a similar presumption for other essential organs/tissues. Scenarios X, XI, or XII represented alternative combinations of WB distributed to CSF, kidney, or GI Tract in different branching ratios. Specifically, it changed from the original [0.4, 0.2, 0.4] (standard case) to [0.3 0.1 0.6] (X), [0.5 0.2 0.3] (XI), or [0.4 0.3 0.3] (XII), respectively.

Figure 2.

The computational result of the standard scenario case was automatically plotted by MATLAB. The VPA time-dependent change of WB attributed to the largest part and other compartments also have their unique characteristics, as shown in this model.

Figure 3 demonstrates the correlated VPA time-dependent changes among compartments. As depicted, the eleven VPA time-dependent changes among compartments show various combinations of parameters in different scenarios. The MATLAB default function automatically plotted Figures 2 and 3. In doing so, the self-developed program was compiled with the subroutine “plot” to derive specific results in this study.¹⁹

Figure 3.

The eleven VPA time-dependent changes among compartments show different combinations of parameters under different presumptions.

4. Discussion

4.1 Interpretation of the VPA time-dependent changes among compartments

4.1.1 Various preset of half-life in the liver (scenario I, II, and III)

The liver is undoubtedly the most complex and dominant compartment in this model. When VPA is absorbed in the GI Tract, it is first absorbed by the epithelial cells of the small intestine's villi. It enters into microvessels, passes through the hepatic portal vein to the liver, and then through the hepatic vein to the inferior vena cava. Thus, VPA is finally distributed to WB by the heart and then to other compartments. Even so, part of the VPA is still fed back toward the liver and transported again to the WB, so VPA time-dependent changes cannot be dealt with as a unidirectional amount of time elapsed; it must be calculated by a group of differential equations as proposed in this study.

Accordingly, the program computation result is very surprising. The VPA in liver absorption causes a high peak, then degrades quickly. The larger half-time of the liver means a longer holding time of VPA in the liver, and the more the liver is absorbed. Thus, the high peak center of liver absorption is also shifted backward by the large half-life itself. As implied in Tab 1 and 2, once the half-life of the liver changes from 0.5 to 1.0 or 3.0, then the corresponding ratio of maximal peak high, Ratio_max, and time of peak reaching maximal, T_max (h), also changes from 0.64, 0.45; 0.76, 0.45 to 0.88, 0.60, respectively. In contrast, the derived Ratio_max, T_max [h] and T_1/2 (=ln2/m) [h] where “m” is inversely proportional to the order of regressed exponential function (i.e., exp(-mx), T_1/2 = ln2/m) for the VPA degradation of WB changes apparently from original (Ratio_max,) 0.70, 0.59, to 0.40; (T_max [h]) 2.10, 3.3. to 7.20; and (T_1/2 [h]) 17.20, 18.73 to 26.65. All the parameters of WB change increasingly with the long dissolving half-life of the liver because the long dissolving half-life can slow down the liver's metabolic mechanism and increase the deposited quantity into WB (cf. Figure 4A). Restated that the long dissolving half-life of the liver can delay the VPA into WB to reach its maximal peak high.

Figure 4.

Four kinds of various VPA dissolving half-lives in specific compartments: (A) liver, (B) WB, (C) CSF, (D) kidney, and (E) four different combinations of branching ratios of WB to CSF, kidney or feedback to GI tract.

4.1.2 Various preset of half-life in WB (scenario I, IV, and V)

WB is another dominant compartment in this model. It attributes the major contribution from itself to others. The model initiates from the oral, GI Tract, liver, then WB; however, the liver is mainly manipulated by the feedback from WB; in addition, the first two compartments (oral and GI Tract) have a minor contribution to the VPA time-dependent change due to their small half-lives. Thus, the WB is practically the dominant compartment. A similar interpretation of the correlation can be analogous to the secular equilibrium in the chain decay of radioactive nuclide.²¹ The WB compartment is equal to the mother nuclide in the series of chain decay; thus, it manipulates the VPA degradation into other compartments. Accordingly, the precise data of every derived VPA time-dependent change among compartments in various scenarios are listed in Table 2. The various VPA changes according to different dissolving half-lives of WB are illustrated in Figure 4B.

Table 2.
The derived numerical data for three dominant compartments: whole body, liver, and CSF in various scenarios.

scenario

Numerical data for the dominant compartment I: Standard case II: liver disorder (A) III: liver disorder (B) IV: WB disorder (A) V: WB disorder (B) VI: CSF disorder (A) VII CSF disorder (B) VIII: kidney disorder (A) IV: kidney disorder (B) X: various paths compose (A) XI: various paths compose (B) XII: various paths compose (C)

Whole body

Ratio_max 0.70 0.59 0.40 0.56 0.81 0.71 0.69 0.70 0.70 0.72 0.69 0.69

T_max (hr) 2.10 3.30 7.20 1.95 2.55 2.10 1.95 2.10 2.10 2.25 1.95 1.95

T_1/2 (hr) 17.20 18.73 26.65 11.74 27.72 16.90 17.36 17.24 17.21 30.13 17.46 11.95

Liver

Ratio_max 0.64 0.76 0.88 0.64 0.64 0.64 0.64 0.64 0.64 0.64 0.64 0.64

T_max (hr) 0.45 0.45 0.60 0.45 0.45 0.45 0.45 0.45 0.45 0.45 0.45 0.45

CSF

Ratio_max 0.16 0.15 0.11 0.23 0.10 0.13 0.21 0.16 0.16 0.13 0.20 0.15

T_max (hr) 6.15 7.35 10.80 5.10 7.50 5.40 7.50 6.15 6.00 7.35 6.00 5.55

T_1/2 (hr) 18.72 20.38 28.88 12.83 31.50 17.77 21.66 19.25 19.25 34.13 19.25 13.08

	scenario
Whole body
Ratio_max	0.70	0.59	0.40	0.56	0.81	0.71	0.69	0.70	0.70	0.72	0.69	0.69
T_max (hr)	2.10	3.30	7.20	1.95	2.55	2.10	1.95	2.10	2.10	2.25	1.95	1.95
T_1/2 (hr)	17.20	18.73	26.65	11.74	27.72	16.90	17.36	17.24	17.21	30.13	17.46	11.95
Liver
Ratio_max	0.64	0.76	0.88	0.64	0.64	0.64	0.64	0.64	0.64	0.64	0.64	0.64
T_max (hr)	0.45	0.45	0.60	0.45	0.45	0.45	0.45	0.45	0.45	0.45	0.45	0.45
CSF
Ratio_max	0.16	0.15	0.11	0.23	0.10	0.13	0.21	0.16	0.16	0.13	0.20	0.15
T_max (hr)	6.15	7.35	10.80	5.10	7.50	5.40	7.50	6.15	6.00	7.35	6.00	5.55
T_1/2 (hr)	18.72	20.38	28.88	12.83	31.50	17.77	21.66	19.25	19.25	34.13	19.25	13.08

Bold represents the “standard case” column represents the original values in normal conditions, comparing to those of abberant numbers in different experiment groups.

As implied in Table 2, once the dissolving half-life of WB changes to small or large (cf. Tab. 1, II, 1.1, or III, 4.5), then the derivation of VPA time-dependent change in WB is Ratio_max: 0.7, 0.56, 0.81; T_max: 2.10, 1.95, 2.55 h, and T_1/2: 17.20, 11.74, 27.72 h, respectively (cf. Tab. 2, scenario I, IV, and V). To further inspect the VPA time-dependent change of WB (cf. Figure 4B), the long dissolving half-life of WB delays the metabolic mechanism of the human body. Therefore, the VPA can deposit more quantity and increase the concentration in WB, although the liver's VPA concentration has barely changed.

4.1.3 Various preset of half-life in CSF (scenario I, VI, and VII)

The correlation between WB and CSF is comparatively simple, among others. It can be interpreted as mother and daughter nuclides in the chained decay of radioactive nuclides.²¹ The CSF is just like a daughter-to-mother compartment (WB). It reaches secular equilibrium with the mother compartment after time elapses, so the VPA change in CSF follows the same trend as WB. As implied in Table 2, only changing the dissolving half-life of CSF barely causes any obvious difference in VPA change in either WB or CSF. In reality, the VPA degradation of CSF is always slower than that of WB. The derived VPA degradation half-life of CSF to WB was changed from 17.77 to 16.90, 18.72 to 17.20, and 21.66 to 17.36 h, referred to the CSF own dissolving half-life as 0.6, 0.8, and 1.2 h, respectively (cf. Tab. 2 scenario VI, VII; and Figure 4C). These interesting phenomena also fulfill the principle of secular equilibrium in chain decay that the decay of the daughter nuclide follows not it's own but the half-life of the mother nuclide. Therefore, the dissolution of VPA in CSF in various scenarios only slightly reflects on the Ratio_max and T_max of either CSF or WB. However, the dissolving half-life of VPA in CSF changes dramatically from 0.6, 0.8 to 1.2 h (cf. Tab. 1, scenario I, VI, and VII, Figure 4C).

4.1.4 Various preset of half-life in the kidney (scenario I, VIII, and IX)

The VPA distribution in the kidney is similar to that in CSF; even though CSF only has one feedback path to WB, the kidney still has another direct path to the bladder. Thus, the VPA in the kidney degrades more rapidly than in CSF. However, the mother compartment (i.e., WB) still dominates the VPA time-dependent change in the daughter compartment (kidney). In addition, the VPA remaining inside the kidney only occupies a minor part. Therefore, any changes in the dissolving half-life of the kidney cause insignificant changes in WB, liver, or CSF at all (cf. Tab. 2, scenario I, VIII, and IX; Figure 4D)

4.1.5 Other combinations of paths of WB (scenario I, X, XI, and XII)

The branching ratio from WB to CSF, kidney, or feedback to GI Tract still needs to be discovered, and no solid information can be quoted from correlated research. Therefore, the combination can be arbitrarily preset as [0.4, 0.2, 0.4] in the initial assumption and also changed to [0.3,0.1, 0.6], [0.5, 0.2, 0.3], or [0.4, 0.3, 0.3], as shown in Table 1 (Scenarios X, XI, and XII). The computational results from the MATLAB program are listed in Figure 3 (Scenarios I, X, XI, and XII). The VPA time-dependent change among various compartments truthfully reflects the different combinations of branching ratios. However, the WB compartment always dominates the whole VPA time-dependent changes as depicted in Figure 4E; the different preset of branching ratios causes the various trends of VPA in WB, specifically, combination of [0.3 0.1 0.6] has the longest VPA degradation T_1/2 of WB among all, (17.20, 30.13, 17.46, 11.95) to Scenarios I, X, XI and XII, respectively, even though the dissolving half-life of WB is still 2.2 h. In addition, either Ratio_max or T_max maintained almost the same and fluctuated slightly between 0.69–0.72 or 1.95–2.10. Thus, the large feedback ratio (cf. Table 1, Scenario X, i₄₂ = 0.6, WB to GI Tract) in Scenario X can effectively slow down the VPA degradation of WB (cf. Figure 4E).

4.2 Other models of VPA time-dependent changes among compartments

A four-compartmental model was proposed by Ibarra et al.,⁷ and the four compartments were defined as the GI Tract, central, peripheral, and bladder. A group of two differential equations was defined and calculated according to the first-order conditional estimation (FOCE) algorithm to be compiled with practical surveyed data for parameter estimation. The evaluated results of VPA time-dependent changes were split into two trends according to seven females and males, respectively, and did not agree with the standard scenario in this study. However, once the dissolving half-lives of the liver and WB were adjusted to 1.4 and 0.9 h, respectively, then the VPA time-dependent changes in WB perfectly fitted the empirical data of Ibarra's finding in the central compartment (i.e., WB in this study) as depicted in Figure 5. Notably, the practical survey in Ibarra's study had a 1-h lag after VPA was given to volunteers. Thus, the two different computational results should have more coincidence once shifted 1 h forward to eliminate the lag timing.

Figure 5.

Other research results on VPA in compartments are also included for comparison. These data were mainly acquired from practical surveys using various clinical techniques.

Another similar research was introduced by Ogungbenro et al..⁸ They defined a 6-compartmental model to analyze the VPA time-dependent changes as systemic blood, liver, kidney, brain, muscle, and the rest of the body, and the computational results were implied according to various ages. They observed the systemic blood and simulated the VPA time-dependent change according to the 2.5^−th, 50^−th, and 97.5^−th percentile of VPA after an oral dose in adults. The corresponding trend of VPA time-dependent changes in systemic blood (i.e., WB in this study) is also shown in Figure 5, where the VPA time-dependent change is plotted according to 50% oral dose in their study.

Cloyd¹³ and his research team explored the VPA time-dependent change in patients with epilepsy, and the results were demonstrated in both total and unbound status. As also attached in Figure 5, the VPA time-dependent change was shown according to unbound status at a maximum of 15%. The semi-empirical analysis was performed using a simplified theoretical basis of loading dose, body mass, and the protocol-specified volume of distribution. They surveyed the first 6 h after dose administration to 12 female and male patients.

4.3 The multi-compartmental models and their solution

It is quite common to define a multi-compartmental model in solving the time-dependent change of solution inside the human body, and it is even more popular in analyzing the patients’ metabolic mechanism in nuclear medicine examination because the traceable feature of radioactive solution can be easily collected from the gamma camera imaging scan of patient's body. The corresponding semi-empirical solution ws often expressed via a linear equation,²² polynomial form,²³ or exponential terms^24–26 to illustrate the complicated status of time-dependent concentration changes. Moreover, for a good regression, some researchers even defined two variables in exponential function or two linear equations in separate sections to have a convincible fit. However, the simple mathematical expression cannot interpret the complex correlations of time-dependent concentration among various compartments. Therefore, we adopt the numerical solution of time-dependent concentrations in this model with the assistance of the MATLAB program and obtain a perfect computational result. This is superior to other researchers using Laplace transform or simple linear regression to analyze the simultaneous differential equations group.^27–29

Furthermore, the simple programmable setup in applying MATLAB makes it easy for the researchers to revise the parameters to match the empirical data. This is essential because the simultaneous differential equations group has to be solved altogether and cannot adjust any single one to barely fit limited empirical data. Note that the correlation between the liver and WB needs precise investigation in future studies because these two compartments have connected feedback paths. Thus, it has to be solved consistently to fulfill the requirements in numerical analysis.

4.4 Future perspective of applying this model

The most important compartments are the liver and WB, whereas other compartments are minor in this model. To solidify and integrate the theoretical estimation in this study and future clinical results from others, we suggest first putting lagging needles into the patient's wrist and hepatic veins, then drawing blood every 30 min during the first four hours after the VPA oral administration. Eventually, the practical VPA time-dependent changes among compartments, especially in the liver or WB, can be recorded and analyzed; then, secondly, by precisely adjusting parameters in the seven-compartmental model, the theoretical distribution of VPA in these two compartments can be experimentally verified. Similar results are demonstrated in Figure 5, covering other researchers’ findings.

5. Conclusions

In this study, VPA physiological pharmacokinetics for epilepsy patients was simulated via a refined seven-compartmental model. The self-developed program run in MATLAB was specifically defined to solve the system of differential equations of VPA time-dependent changes among compartments. The computational results showed quite unexpected VPA concentration changes in the human body. Either the liver or WB was the dominant compartment in this model, and the dissolving half-life was 0.5 and 2.2 h, respectively. In contrast, the other compartments could be dealt with as minor ones that follow the secular equilibrium principle as chain decay of radioactive nuclide. The WB could be defined as similar to the mother compartment, which dominates the daughter compartments in this model. Accordingly, the derived VPA degradation of minor compartments always went with the VPA degradation of WB with a large half-life. All parameters, including the half-life or branching ratio of every compartment, were preset and changeable to satisfy the practical survey in future prospective studies.

Footnotes

Acknowledgements

None to report.

Ethical approval

This study did not involve human or animal subjects directly. All simulations were based on hypothetical or previously published pharmacokinetic parameters. Therefore, ethical approval was not required.

Author contributions

Tsung-Han Lin and Ya-Hui Lin conceived and designed the study. Keng-Yi Wu and Tzu-Hsien Chao performed the investigation and data collection. Bing-Ru Peng provided resources. Lung-Kwang Pan conducted the formal analysis and visualization. The manuscript was drafted by Ya-Hui Lin and all authors contributed to reviewing, editing, and final approval. Ya-Hui Lin supervised the project and acquired funding.

Funding

The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the Taichung Armed Forces General Hospital in Taiwan, (grant number TCAFGH-D-113034).

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Data availability

The MATLAB code used for simulation in this study is available from the corresponding author upon reasonable request.

Copyright and originality

The authors declare that this manuscript is original, has not been published previously, and is not under consideration for publication elsewhere. All authors have read and approved the final version of the manuscript.

ORCID iD

Ya-Hui Lin

References

Singh

Gupta

Verma

, et al. Hidden pharmacological activities of valproic acid: a new insight. Biomed Pharmacother. 2021; 142: 112021.

Romoli

Mazzocchetti

D’Alonzo

, et al. Valproic acid and epilepsy: from molecular mechanisms to clinical evidences. Curr Neuropharmacol. 2018; 17: 926–946.

Zhu

Shi

, et al. The pharmacogenomics of valproic acid. J Hum Genet. 2017; 62: 1009–1014.

Bertelsen

Landau

Vase

, et al. Acute in vivo effect of valproic acid on the GABAergic system in rat brain: a [11C]Ro15-4513 microPET study. Brain Res. 2018; 1680: 110–114.

Kay

Greene

Kang

, et al. M-current preservation contributes to anticonvulsant effects of valproic acid. J Clin Invest. 2015; 125: 3904–3914.

Chang

Walkeere

Williams

RSB

. Seizure-induced reduction in PIP3 levels contributes to seizure-activity and is rescued by valproic acid. Neurobiol Dis. 2014; 62: 296–306.

Ibarra

Vazquez

Fagiolino

, et al. Sex-related differences in valproic acid pharmacokinetics after oral single dose. J Pharmacokinet Pharmacodyn 2013; 40: 479–486.

Ogungbenro

Aarons

; the CRESim & Epi-CRESim project Groups. A physiologically based pharmacokinetic model for Valproic acid in adults and children. Eur J Pharm Sci. 2014; 63: 45–52.

Simeon

Beaudouin

Brotzmann

, et al. Multistate models of developmental toxicity: application to valproic acid-induced malformations in the zebrafish embryo. Toxicol Appl Pharmacol. 2021; 414: 115424.

10.

Sopasakis

Sarimveis

Macheras

, et al. Fractional calculus in pharmacokinetics. J. Pharmacokinet Pharmacodyn 2018; 45: 107–125.

11.

Azizi

. Study on the application of the fractional calculus in pharmacokinetic modeling Commun Biomath. Technol Innov Eng Res 2022; 5: 52–70.

12.

Jones

Rowland-Yeo

. Basic concepts in physiologically based pharmacokinetic modeling in drug discovery and development. Pharmacometrics & Systems Pharmacology 2013; 2: e63.

13.

Soars

Burchell

Riley

. In vitro analysis of human drug glucuronidation and prediction of in vivo metabolic clearance. J. Pharmacology and Experimental Therapeutics 2002; 301: 382–390.

14.

Cloyd

Fischer

Kriel

, et al. Valproic acid pharmacokinetics in children IV. Effects of age and antiepileptic drugs on protein binding and intrinsic clearance. Clin Pharmacol 1993; 53: 22–29.

15.

Cloyd

Dutta

Cao

, et al. The depacon study group. Valproate unbound fraction and distribution volume following rapid infusions in patients with epilepsy. Epilepsy Res. 2003; 53: 19–27.

16.

Huang

Lin

Kittipayak

, et al. Biokinetic model of radioiodine I-131 in nine thyroid cancer patients subjected to in-vivo gamma camera scanning: a simplified five compartmental model. PLOS ONE 2020; 15: e0232480.

17.

Yeh

Chen

Tang

, et al. A quantitative evaluation of multiple biokinetic models using an assembled water phantom: a feasibility study. PLOS One 2017; 12: e0189244.

18.

Chiang

Chung

, et al. Quantitative analysis of multiple biokinetic models using a dynamic water phantom: a feasibility study. Bioengineered 2016; 7: 304–313.

19.

MATLAB, V7.0.1.24704 (R14). Matrix Laboratory Developed by MathWorks: Natick, MA, USA, 2004. 2004.

20.

Pan

Chiang

Huang

, et al. Thyroid biokinetics for radioactive I-131 in twelve thyroid cancer patients via the refined nine-compartmental model. Appl Sci. 2022; 12: 5538.

21.

Cember

Jojnson

. Introduction to Health Physics 4th Ed. Ch. 4. New York, NY: McGraw-Hill NY, 2008, ISBN: 978-0-07-164323-8.

22.

Bucci

Rebelos

Oikonen

, et al. Kinetic modeling of brain [18-F] FDG positron emission tomography time activity curves with input function recovery (IR) method. Metabolites 2024; 114: 14020114.

23.

GaoHua

Abduljalil

Jamei

, et al. Pregnancy physiologically based pharmacokinetic (p-PBPK) model for disposition of drugs metabolized by CYP1A2, CYP2D6 and CYP3A4. British J Clinical Pharmacology 2012; 74: 873–885.

24.

Donnet

Samson

. A review on estimation of stochastic differential equations for pharmacokinetic/pharmacodynamic models. Adv Drug Delivery Rev. 2013; 65: 1.

25.

Beni

Watabe

. Compvision: an open-source five-compartmental software for biokinetic simulations. Open Physics 2021; 19: 454–459.

26.

Shuryak

Dadachova

. Usefulness of continuous probability distributions of rates for modelling radionuclide biokinetics in humans and animals. Sci Rep. 2019; 9: 1218.

27.

Voicu

Jiquidi

Mircioiu

, et al. Experimental evaluation of 65Zn decorporation kinetics following rapid and delayed Zn-DTPA interventions in rats: biphasic compartmental and square-root law mathematical modeling. Pharmaceutics 2021; 13: 1830.

28.

Vieru

Fetecau

Ahmed

, et al. A generalized kinetic model of the advection-dispersion process in a sorbing medium. Mathematical Modelling of Nature Phenomena 2021; 16: 39.

29.

Vinnett

Waters

. Representation of kinetics models in batch flotation as distributed first-order reactions. Minerals 2020; 10: 913.

	scenario
Numerical data for the dominant compartment	I: Standard case	II: liver disorder (A)	III: liver disorder (B)	IV: WB disorder (A)	V: WB disorder (B)	VI: CSF disorder (A)	VII CSF disorder (B)	VIII: kidney disorder (A)	IV: kidney disorder (B)	X: various paths compose (A)	XI: various paths compose (B)	XII: various paths compose (C)
Whole body
Ratio_max	0.70	0.59	0.40	0.56	0.81	0.71	0.69	0.70	0.70	0.72	0.69	0.69
T_max (hr)	2.10	3.30	7.20	1.95	2.55	2.10	1.95	2.10	2.10	2.25	1.95	1.95
T_1/2 (hr)	17.20	18.73	26.65	11.74	27.72	16.90	17.36	17.24	17.21	30.13	17.46	11.95
Liver
Ratio_max	0.64	0.76	0.88	0.64	0.64	0.64	0.64	0.64	0.64	0.64	0.64	0.64
T_max (hr)	0.45	0.45	0.60	0.45	0.45	0.45	0.45	0.45	0.45	0.45	0.45	0.45
CSF
Ratio_max	0.16	0.15	0.11	0.23	0.10	0.13	0.21	0.16	0.16	0.13	0.20	0.15
T_max (hr)	6.15	7.35	10.80	5.10	7.50	5.40	7.50	6.15	6.00	7.35	6.00	5.55
T_1/2 (hr)	18.72	20.38	28.88	12.83	31.50	17.77	21.66	19.25	19.25	34.13	19.25	13.08