Development of a Compartment Model to Study the Pharmacokinetics of Medical THC after Oral Administration

1. Introduction

Delta9-tetrahydrocannabinol (THC), one of the cannabinoids found in the cannabis plant, possesses several clinically useful pharmacological characteristics [1]. As a result, THC is extracted and developed into a medical product, serving as an alternative medicine to treat patients who do not respond to modern medicine. This includes patients experiencing nausea and vomiting from chemotherapy, medication-resistant epilepsy, anorexia in AIDS patients, and those seeking to improve the quality of life in palliative care. Several studies show that medical THC has positive effects in treating various diseases [2,3]. On the other hand, THC is a psychoactive chemical that can be addictive if taken in inappropriate amounts. Therefore, the administration of medical THC to patients must be approached with caution because the precise and appropriate dosage for various patient groups in the administration of different medical THC products has not yet been established. The general advice is to start with a small dosage, for example, administering it once a day, and gradually increase the dosage when the patient does not experience side effects. Therefore, having a tool to help describe the pharmacokinetics of medical THC will increase our understanding of medical THC pharmacokinetics and can also be used to guide the determination of the dosage and administration duration of medical THC products.

Currently, various forms of medical THC products are available, catering to different treatment needs that demand rapid bloodstream absorption, ranging from oral, inhalation, and smoking to injection. Among these, oral dosage forms like capsules or tablets are widely preferred due to their convenience and straightforward dosage determination. While previous research has delved into THC pharmacokinetics post-smoking [4], this study will focus on investigating the pharmacokinetics of clinical THC administered orally. After oral administration, THC is absorbed through the gastrointestinal tract (stomach and small intestine) with a bioavailability of about 0.1-0.2% [5]. Subsequently, it enters the liver via the hepatic portal vein for first-pass metabolism. In the liver, THC is metabolized by the enzymes CYP2C and CYP3A into its main metabolites, THC-OH and THCCOOH. THC-OH is the primary psychoactive metabolite, while THCCOOH is a secondary metabolite and is not psychoactive [6]. The remaining THC and both of its metabolites pass through the heart before entering the bloodstream. Both THC and THC-OH simultaneously enter the brain. Approximately 20% of THC and other cannabinoids are excreted in urine, with more than 65% eliminated in feces within 5 days [7].

Compartment modeling is another method that can be employed to study pharmacokinetics, encompassing processes such as absorption, distribution, metabolism, and excretion. It is a mathematical technique used to represent the entire body by dividing it into compartments. This modeling approach helps avoid challenges and reduces costs associated with clinical studies. Previous research studies have proposed compartment models to describe the complex pharmacokinetics of THC. For instance, in 2015, Jules and their group [8] provided a four-compartment pharmacokinetic model for THC after oral, intravenous, and pulmonary dosing. In 2017, Rakesh [9] conducted a study on the pharmacokinetics and pharmacodynamics of THC and its active metabolite using effect compartment modeling-based approaches. In 2020, Cristina and their team [10] constructed a population multi-compartment pharmacokinetic model for THC and its active and inactive metabolites after controlled smoked cannabis delivery. Additionally, in 2023, we created a compartment model to study the pharmacokinetics of THC and its metabolites after smoking [4]. Usually, we can solve analytical solutions when building a compartment model with one or two compartments. However, because the system of differential equations for a three- or multi-compartment model is large, it becomes difficult to derive an analytical solution. Therefore, numerical methods were introduced to help reduce the challenge of solving analytical solutions for the multi-compartment model. One of the numerical methods used to solve the ODE system of the compartment pharmacokinetic model is the nonstandard finite difference method (NSFD). This method evolved from the standard finite difference method (SFD) developed by Mickens [11,12] to construct numerical solution schemes. Previous studies have shown that the NSFD method provides more stable numerical results compared to the SFD method, regardless of the chosen step size [13,14,15,16].

In this study, we developed a compartmental model based on our previous research, which focused on the compartmental pharmacokinetic model of medical THC following smoking [4]. The objective of this work was to develop a compartmental model to comprehensively investigate the pharmacokinetics of medical THC and its metabolites in the body after oral administration. The model was converted into ordinary differential equations (ODEs) to depict the rate of concentration change of THC and its metabolites in bodily compartments, utilizing the law of mass action. For the simulation, we will find numerical solutions by applying the NSFD method to create numerical solution schemes. These schemes were then implemented in the MATLAB program, along with the estimated pharmacokinetic rate constants of THC and its metabolites, to generate simulation curves.

2. Methodology

In this section, we develop a multi-compartment model to describe the pharmacokinetics of THC after oral administration. The compartment model is converted into ordinary differential equations (ODEs) to describe the rate of change of THC and its metabolite concentrations in each compartment, using the law of mass action. Afterward, we apply the nonstandard finite difference method (NSFD) to transform the ODEs into numerical solution schemes. Additionally, we provide the estimated pharmacokinetic rate constants used to simulate the results, including examples of actual plasma concentrations of THC and THC-OH that are compared with our simulation results.

2.1. Compartmental Modeling

Here, we considered dividing the compartment model representing the body into four phases: the first phase encompassed the absorption compartment, representing the gastrointestinal tract, to elucidate THC absorption. The second phase dealt with the pharmacokinetics of THC, while the third phase focused on the pharmacokinetics of THC-OH. The fourth phase examined the pharmacokinetics of THCCOOH in plasma. In this model, we considered the pharmacokinetics of THC and THC-OH in the blood and tissues by dividing the body into three physiologically significant compartments: the central compartment, representing the blood (plasma) or systemic circulation; the rapidly equilibrating tissue compartment, representing organs like the liver, lungs, and kidneys; and the slowly equilibrating tissue compartment, representing tissues such as muscle, fat, and bone. Additionally, we created effect compartments for THC and THC-OH to describe their concentrations at the sites of effect.

We defined rate constants for the absorption, distribution, metabolism, and excretion of THC and its metabolites in each compartment. The concentrations were represented as follows: $C_a\left(t\right)$ for THC in the absorption compartment (gastrointestinal tract), $C_1\left(t\right)$ for THC in the central compartment (plasma), $C_2\left(t\right)$ for THC in the rapidly equilibrating tissue compartment, $C_3\left(t\right)$ for THC in the slowly equilibrating tissue compartment, and $C_{e1}\left(t\right)$ for THC in the effect compartment. $C_4\left(t\right)$ for THC-OH in the central compartment, $C_5\left(t\right)$ for THC-OH in the rapidly equilibrating tissue compartment, $C_6\left(t\right)$ for THC-OH in the slowly equilibrating tissue compartment, $C_{e2}\left(t\right)$ for THC-OH in the effect compartment, and $C_7\left(t\right)$ for THCCOOH in the central compartment. We let $k_a$ be the absorption rate constant of THC. The distribution rate constants of THC in each compartment were denoted as $k_{12}$, $k_{21}$, $k_{13}$, $k_{31}$, and $k_{1e}$. The metabolic rate constants of THC to THC-OH and THCCOOH were represented by $k_{14}$ and $k_{17}$, respectively. For THC-OH, the distribution rate constants in each compartment were $k_{45}$, $k_{54}$, $k_{46}$, $k_{64}$, and $k_{4e}$, while $k_{47}$ was the metabolic rate constant of THC-OH to THCCOH. The elimination rate constants of THC-OH and THCCOOH from the central compartment were $k_{40}$ and $k_{70}$, respectively. Additionally, $k_{e10}$ and $k_{e20}$ denoted the elimination rate constants of THC and THC-OH from the effect compartment. This representation was illustrated in Figure 1.

By applying the law of mass action to the compartment model illustrated in Figure 1, we derived a set of ordinary differential equations representing the rate of change of THC and its metabolite concentrations in each compartment, as follows:

$$\left.\begin{array}{ccc}\hfill \frac{d{C}_a\left(t\right)}{dt}=& -k_a{C}_a\left(t\right);\hfill & \hfill C_a\left(0\right)=C_0\\ \hfill \frac{d{C}_1\left(t\right)}{dt}=& k_a{C}_a\left(t\right)-(k_{12}+k_{13}+k_{14}+k_{17}+k_{1e})C_1\left(t\right)\hfill \\ \hfill & +k_{21}{C}_2\left(t\right)+k_{31}{C}_3\left(t\right);\hfill & \hfill C_1\left(0\right)=0\\ \hfill \frac{d{C}_2\left(t\right)}{dt}=& k_{12}{C}_1\left(t\right)-k_{21}{C}_2\left(t\right);\hfill & \hfill C_2\left(0\right)=0\\ \hfill \frac{d{C}_3\left(t\right)}{dt}=& k_{13}{C}_1\left(t\right)-k_{31}{C}_3\left(t\right);\hfill & \hfill C_3\left(0\right)=0\\ \hfill \frac{d{C}_4\left(t\right)}{dt}=& k_{14}{C}_1\left(t\right)-(k_{40}+k_{45}+k_{46}+k_{47}+k_{4e})C_4\left(t\right)\hfill \\ \hfill & +k_{54}{C}_5\left(t\right)+k_{64}{C}_6\left(t\right);\hfill & \hfill C_4\left(0\right)=0\\ \hfill \frac{d{C}_5\left(t\right)}{dt}=& k_{45}{C}_4\left(t\right)-k_{54}{C}_5\left(t\right);\hfill & \hfill C_5\left(0\right)=0\\ \hfill \frac{d{C}_6\left(t\right)}{dt}=& k_{46}{C}_4\left(t\right)-k_{64}{C}_6\left(t\right);\hfill & \hfill C_6\left(0\right)=0\\ \hfill \frac{d{C}_7\left(t\right)}{dt}=& k_{17}{C}_1\left(t\right)+k_{47}{C}_4\left(t\right)-k_{70}{C}_7\left(t\right);\hfill & \hfill C_7\left(0\right)=0\\ \hfill \frac{d{C}_{e1}\left(t\right)}{dt}=& k_{1e}{C}_1\left(t\right)-k_{e10}{C}_{e1}\left(t\right);\hfill & \hfill C_{e1}\left(0\right)=0\\ \hfill \frac{d{C}_{e2}\left(t\right)}{dt}=& k_{4e}{C}_4\left(t\right)-k_{e20}{C}_{e2}\left(t\right);\hfill & \hfill C_{e2}\left(0\right)=0\end{array}\right\}$$

(1)

where $C_0\mathrm{is}\mathrm{the}\mathrm{initial}\mathrm{concentration}\mathrm{of}\mathrm{THC}(\mathrm{ng}/\mathrm{ml}).$

Next, we will apply the NSFD method to System 1 to construct numerical solution schemes, thereby reducing the challenge of solving analytical solutions for large systems of equations.

2.2. NSFD Schemes

The nonstandard finite difference method (NSFD), a numerical approach, has been used in continuous models. It was developed from the standard finite difference method (SFD) by Mickens [11,12]. It serves as an alternative method to ensure the stability of numerical solutions, irrespective of the chosen step size. The NSFD procedures are based on the following rules:

1) The following represents the discrete first derivative:

$$\frac{dC}{dt}\to \frac{C_{k+1}-\psi C_k}{\varphi }$$

where $\psi$ and $\varphi$ depends on step-size $△t=h$ and satisfy the conditions:

$$\psi =1+\mathcal{O}\left(h\right),\varphi =h+\mathcal{O}\left(h^2\right).$$

$\psi$ and $\varphi$ can have different functions from one another. While specific forms for a given equation can be readily found, there are currently no general guidelines for choosing the functions $\psi \left(h\right)$ and $\varphi \left(h\right)$. Frequently employed functional representations for $\psi \left(h\right)$ and $\varphi \left(h\right)$ are

$$\varphi \left(h\right)=\frac{e^{\lambda h}-1}{\lambda },\psi \left(h\right)=1,$$

where $\lambda$ is some parameter appearing in the differential equations.

2) Both linear and nonlinear terms may require a nonlocal representation on the discrete computational; for example:

$$\begin{array}{cc}\hfill & C⟹C_{k+1}.\hfill \end{array}$$

Preliminary rules exist for constructing denominator functions for a system of coupled first-order ordinary differential equations:

1) Form an initial, finite difference model by replacing all first-derivatives by discrete forward-Euler terms.

2) For a particular discrete equation, in general, its dependent variable will occur linearly in its evaluation at the $(k+1)$-th time step. Solve for this dependent variable at the $(k+1)$-th time step in terms of all other dependent variables evaluated at the k-th time step.

3) The denominator function can be chosen as follows

$$\varphi (h,\lambda )=\frac{e^{\lambda h}-1}{\lambda }.$$

If a factor in a given discrete equation has an expression of the form $(1+\lambda h)$, where $\lambda$ is made up of one or more parameters that occur in the original differential equations.

4) In the discrete finite-difference schemes constructed in (1), replace h by the appropriate $\varphi (h,\lambda )$.

5) For the case where $\lambda =0$, use the denominator function $\varphi \left(h\right)=h$.

Therefore, from the ordinary differential equations in Equation (1), by applying the NSFD method to convert them into numerical solution schemes, we obtain

$$\left.\begin{array}{ccc}\hfill \frac{C_{a,k+1}-C_{a,k}}{{\varphi }_a}=& -k_a{C}_{a,k+1};\hfill & \hfill C_{a,0}=C_0\\ \hfill \frac{C_{1,k+1}-C_{1,k}}{{\varphi }_1}=& k_a{C}_{a,k+1}-(k_{12}+k_{13}+k_{14}+k_{17}+k_{1e})C_{1,k+1}\hfill \\ \hfill & +k_{21}{C}_{2,k}+k_{31}{C}_{3,k};\hfill & \hfill C_{1,0}=0\\ \hfill \frac{C_{2,k+1}-C_{2,k}}{{\varphi }_2}=& k_{12}{C}_{1,k+1}-k_{21}{C}_{2,k+1};\hfill & \hfill C_{2,0}=0\\ \hfill \frac{C_{3,k+1}-C_{3,k}}{{\varphi }_3}=& k_{13}{C}_{1,k+1}-k_{31}{C}_{3,k+1};\hfill & \hfill C_{3,0}=0\\ \hfill \frac{C_{4,k+1}-C_{4,k}}{{\varphi }_4}=& k_{14}{C}_{1,k+1}-(k_{40}+k_{45}\hfill \\ \hfill & +k_{46}+k_{47}+k_{4e})C_{4,k+1}+k_{54}{C}_{5,k}+k_{64}{C}_{6,k};\hfill & \hfill C_{4,0}=0\\ \hfill \frac{C_{5,k+1}-C_{5,k}}{{\varphi }_5}=& k_{45}{C}_{4,k+1}-k_{54}{C}_{5,k+1};\hfill & \hfill C_{5,0}=0\\ \hfill \frac{C_{6,k+1}-C_{6,k}}{{\varphi }_6}=& k_{46}{C}_{4,k+1}-k_{64}{C}_{6,k+1};\hfill & \hfill C_{6,0}=0\\ \hfill \frac{C_{7,k+1}-C_{7,k}}{{\varphi }_7}=& k_{17}{C}_{1,k+1}+k_{47}{C}_{4,k+1}-k_{70}{C}_{7,k+1};\hfill & \hfill C_{7,0}=0\\ \hfill \frac{C_{e1,k+1}-C_{e1,k}}{{\varphi }_{e1}}=& k_{1e}{C}_{1,k+1}-k_{e10}{C}_{e1,k+1};\hfill & \hfill C_{e1,0}=0\\ \hfill \frac{C_{e2,k+1}-C_{e2,k}}{{\varphi }_{e2}}=& k_{4e}{C}_{4,k+1}-k_{e20}{C}_{e2,k+1};\hfill & \hfill C_{e2,0}=0\end{array}\right\}$$

(2)

By rearranging System (2), we obtain the solution schemes as follows:

$$\left.\begin{array}{ccc}\hfill & C_{a,k+1}=\frac{C_{a,k}}{1+k_a{\varphi }_a};\hfill & \hfill C_{a,0}=C_0\\ \hfill & C_{1,k+1}=\frac{C_{1,k}+k_a{\varphi }_1{C}_{a,k+1}+k_{21}{\varphi }_1{C}_{2,k}+k_{31}{\varphi }_1{C}_{3,k}}{1+(k_{12}+k_{13}+k_{14}+k_{17}+k_{1e}){\varphi }_1};\hfill & \hfill C_{1,0}=0\\ \hfill & C_{2,k+1}=\frac{C_{2,k}+k_{12}{\varphi }_2{C}_{1,k+1}}{1+k_{21}{\varphi }_2};\hfill & \hfill C_{2,0}=0\\ \hfill & C_{3,k+1}=\frac{C_{3,k}+k_{13}{\varphi }_3{C}_{1,k+1}}{1+k_{31}{\varphi }_3};\hfill & \hfill C_{3,0}=0\\ \hfill & C_{4,k+1}=\frac{C_{4,k}+k_{14}{\varphi }_4{C}_{1,k+1}+k_{54}{\varphi }_4{C}_{5,k}+k_{64}{\varphi }_4{C}_{6,k}}{1+(k_{40}+k_{45}+k_{46}+k_{47}+k_{4e}){\varphi }_4};\hfill & \hfill C_{4,0}=0\\ \hfill & C_{5,k+1}=\frac{C_{5,k}+k_{45}{\varphi }_5{C}_{4,k+1}}{1+k_{54}{\varphi }_5};\hfill & \hfill C_{5,0}=0\\ \hfill & C_{6,k+1}=\frac{C_{6,k}+k_{46}{\varphi }_6{C}_{4,k+1}}{1+k_{64}{\varphi }_6};\hfill & \hfill C_{6,0}=0\\ \hfill & C_{7,k+1}=\frac{C_{7,k}+k_{17}{\varphi }_7{C}_{1,k+1}+k_{47}{\varphi }_7{C}_{4,k+1}}{1+k_{70}{\varphi }_7};\hfill & \hfill C_{7,0}=0\\ \hfill & C_{e1,k+1}=\frac{C_{e1,k}+k_{1e}{\varphi }_{e1}{C}_{1,k+1}}{1+k_{e10}{\varphi }_{e1}};\hfill & \hfill C_{e1,0}=0\\ \hfill & C_{e2,k+1}=\frac{C_{e2,k}+k_{4e}{\varphi }_{e2}{C}_{4,k+1}}{1+k_{e20}{\varphi }_{e2}};\hfill & \hfill C_{e2,0}=0\end{array}\right\}$$

(3)

where the various $\varphi$ depend on step-size (h) and are written in the functional form as follows:

$$\begin{array}{cc}\hfill & {\varphi }_a={\displaystyle \frac{e^{k_{a}h}-1}{k_a}},{\varphi }_1={\displaystyle \frac{e^{(k_{12}+k_{13}+k_{14}+k_{17}+k_{1e})h}-1}{k_{12}+k_{13}+k_{14}+k_{17}+k_{1e}}},\hfill \\ \hfill & {\varphi }_2={\displaystyle \frac{e^{k_{21}h}-1}{k_{21}}},{\varphi }_3={\displaystyle \frac{e^{k_{31}h}-1}{k_{31}}},\hfill \\ \hfill & {\varphi }_4={\displaystyle \frac{e^{(k_{40}+k_{45}+k_{46}+k_{47}+k_{4e})h}-1}{k_{40}+k_{45}+k_{46}+k_{47}+k_{4e}}},{\varphi }_5={\displaystyle \frac{e^{k_{54}h}-1}{k_{54}}},\hfill \\ \hfill & {\varphi }_6={\displaystyle \frac{e^{k_{64}h}-1}{k_{64}}},{\varphi }_7={\displaystyle \frac{e^{k_{70}h}-1}{k_{70}}},\hfill \\ \hfill & {\varphi }_{e1}={\displaystyle \frac{e^{k_{1e}h}-1}{k_{1e}}},{\varphi }_{e2}={\displaystyle \frac{e^{k_{4e}h}-1}{k_{4e}}}.\hfill \end{array}$$

2.3. Pharmacokinetic Parameters

To simulate the results, we needed to know the pharmacokinetic rate constants of THC and its metabolites. In this study, we derived these rate constants from our previous investigation on the pharmacokinetics of THC and its metabolites after smoking [4], adjusting them to obtain simulation results as close to the actual data as possible. Therefore, we obtained the estimated pharmacokinetic rate constants for THC and THC-OH, as shown in Table 1.

2.4. Actual Data

The actual data sample we used for comparison with the simulation results in this study consisted of the plasma concentrations of THC and THC-OH after the oral administration of 10 mg of THC, based on a study conducted by Guy and Robson [17]. Experimental data from the clinical study are shown in Table 2.

From the actual data sample as shown in Table (Table 2), we will compare it with the simulated results of the concentration of THC and THC-OH to investigate any errors, which we will discuss in the conclusion and discussion section.

3. Simulation Results

For the simulations, we utilized the NSFD scheme specified in System (3) and implemented it in the MATLAB program. We employed the estimated pharmacokinetic rate constants from Table 2 and selected a step size (h) of 0.05. Consequently, we obtained simulated concentrations of THC and THC-OH in plasma and other tissues following the oral administration of 10 mg of THC, replicating the experimental conditions of the clinical studies conducted by Guy and Robson [17]. These results are presented in Figure 2 and Figure 3 for THC and THC-OH, respectively. Additionally, the simulated concentrations of THC and THC-OH are illustrated in Figure 4 and Figure 5, respectively.

4. Conclusion and Discussion

In this study, we developed a compartment model to describe the pharmacokinetics of delta9-tetrahydrocannabinol (THC) and its metabolites following oral administration. The compartmental model was converted into ordinary differential equations (ODEs), which represent the rate of change in the concentration of THC and its metabolites in each compartment, using the law of mass action. Subsequently, we applied the nonstandard finite difference method (NSFD) to create numerical solution schemes. These schemes were then implemented in a MATLAB program, along with the estimated pharmacokinetic rate constants, to generate the simulation results. The rate constants used in the simulation were adjusted from our previous study on the pharmacokinetics of THC and its metabolites following smoking.

In the simulation results of THC and THC-OH following the oral administration of 10 mg THC under the conditions of the clinical study by Guy and Robson [17], using the NSFD scheme as depicted in System (3) with a chosen h of 0.05, we observed the following trends:

For THC, the concentration in plasma, denoted as $C_1\left(t\right)$, rapidly increased, reaching a maximum concentration of 4.52 ng/ml between 50 and 55 minutes. Subsequently, it decreased rapidly within the 55-360 minute interval before gradually declining over time. When comparing the simulation results to the actual data from the study by Guy and Robson, we found that the simulated plasma THC concentration curves closely resembled the actual data, with an R-squared value of 0.8167 and an RMSE value of 0.6626. In the absorption compartment, denoted as $C_a\left(t\right)$, the THC concentration rapidly decreased from the initial concentration of 24.85 ng/ml within the 0-200 minute interval and continued to decrease gradually until reaching 0 ng/ml after 200 minutes. The THC concentration in other tissues gradually increased, depending on the distribution rate constants. The maximum concentrations of THC in the rapidly equilibrating tissue, the slowly equilibrating tissue, and the effect site were 3.36 ng/ml at 50-55 minutes, 1.68 ng/ml at 50-55 minutes, and 0.043 ng/ml at 160-235 minutes, respectively. After reaching their maximum concentrations, it gradually decreased to 0 ng/ml over time, as illustrated in Figure 2. The simulated concentration of THC after oral administration of 10 mg in each tissue is depicted in Figure 4.

For THC-OH, we observed a rapid increase in the THC-OH concentration in plasma ($C_4\left(t\right)$), reaching a maximum concentration of 6.76 ng/ml at 105 minutes. When comparing the simulation results to the actual data from the study by Guy and Robson, we found that the simulated plasma THC-OH concentration curve closely resembled the actual data, albeit with a maximum time discrepancy. The R-squared value was 1.0189, and the RMSE value was 0.7848. This discrepancy may be related to the estimated pharmacokinetic parameters. Additionally, we noted an increase in the THC-OH concentration in other tissues to their respective maximum points, followed by a gradual decrease over time. The maximum concentrations of THC-OH in the rapidly equilibrating tissue, the slowly equilibrating tissue, and the effect site were 0.042 ng/ml at 105 minutes, 0.022 ng/ml at 105 minutes, and 1.60 ng/ml at 250-260 minutes, respectively, as depicted in Figure 3. The simulated concentration of THC-OH after oral administration of 10 mg THC in each tissue is shown in Figure 5.

Furthermore, when comparing the simulation results with another study that administered the same initial dosage of 10 mg oral medical THC, such as the previously published study by Nadulski and their team [18], we noted differences in the mean maximum concentrations of THC and THC-OH in plasma ($C_1$ and $C_4$, respectively). In their study, the mean maximum concentration of THC was 3.20 ng/ml at a mean maximum time of 63 minutes, while the mean maximum concentration of THC-OH was 4.48 ng/ml at a mean maximum time of 90 minutes. These slight discrepancies between our simulation results and Nadulski’s study may be attributed to variations in the clinical study conditions. Based on our simulation results and the findings from Nadulski’s study, we can infer that following the oral administration of 10 mg of THC, the maximum concentration of THC will range from approximately 3.2 to 4.52 ng/ml, occurring approximately between 50 and 63 minutes. Similarly, the maximum concentration of THC-OH will range from approximately 4.88 to 6.76 ng/ml, occurring approximately between 90 and 105 minutes. Therefore, in subsequent medical THC administrations, it is advisable to wait for the concentrations of THC and THC-OH in plasma to decrease, thereby ensuring that they do not reach excessively high levels.

Based on the simulation results, our model can produce pharmacokinetic curves of THC and its metabolites in plasma after oral administration that are close to the actual data from the study conducted by Guy and Robson. Furthermore, our model can predict the concentration of THC and its metabolites in other bodily compartments, providing an advantage over most clinical studies that only consider or investigate plasma concentrations of THC and its metabolites. Consequently, our model can serve as a tool to advance our understanding of the pharmacokinetics of medical THC, including assisting in determining the dosage and administration duration for oral medical THC products. Moreover, the model shows potential for further refinement and development in the future.