Modeling Drug Interaction Using the Pharmacokinetic One-Compartment Model and Differential Equations

• Pharmacokinetics => branch of Pharmacology, specifically the study of how drugs move throughout the body.
• We will be using first order ODEs to model drug interaction using the One Compartment Model(OCM). The One Compartment Model is the most simple form of interaction of a drug with a patient's body. Using the OCM, as well as in other models, the elimination of a drug causes the amount inside the compartment to decrease exponentially. While modeling drug interaction, the ODEs being used will be IVPs since we will be given a value at time t=0 in the problem statement. A visual of the OCM is shown on the next slide:

The One-Compartment Model (OCM) is the basic idea of how a concentration of a drug moves throughout the body. This happens in three stages: The administration to the patient, entrance to the 'compartment' which is the drug taking effect in the patient, and the excretion (elimination) from the compartment. In mathematical terms these three stages are:

• The Initial Amount
• The Current Concentration in the body
• The rate at which the drug leaves the body (Often referred to as half-life) Since the OCM is a first order model, and decreases exponentially, we begin with the exponential decay function:

\[ \begin{aligned} \frac{dP}{dt}=-Pk \end{aligned} \] where P is the current amount inside the compartment, -k is the decay rate, and \( \frac{dP}{dt} \) is the overall rate of change. We can see that the form of our solution equations will be exponentially decaying functions by finding the general solution of the exponential decay function. Using the method of separation of variables, we arrive at: \[ \begin{aligned} P=P_{0}e^{-kt} \end{aligned} \] This shows that as we have done previously, we can assume the solution to our IVPs is in the form \( y=e^{-kt} \). Now we can pick a method suitable to solve our ODEs. The Runge-Kutta Predictor/Corrector Method will be what we use to numerically solve using python.

Problem Statement: The amount of drug in the blood of a patient administered via intravenous line is governed by the IVP \( y'(t)=-0.2y+3 \) where \( y(0)=0 \) (The initial concentration is 0mg). The amount of drug is measured in milligrams and t is measured in hours. Find and graph the solution of the drug concentration over a 24 hour interval.

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import norm
import pandas as pd

def rk4(N):
    f = lambda x,y: -0.2*y+3    #Given derivative function from IVP
    g = lambda x:  np.exp(-x/5)*(15)+15 #Analytically derived exact solution

    a = 0    #Left endpoint of [a,b]
    b = 24    #Right endpoint of [a,b]
    h = (b-a)/N  #Initialize step size
    x0 = a   #Initial x value given in IVP
    y0 = 0    #Initial y value given in IVP
    x = np.zeros(N+1)    #Initialize x values for numerical solution

    y = np.zeros(N+1)   #Initialize y values for numerical solution
    gn = np.zeros(N+1)  #Initialize gn values for exact solution
    en = np.zeros(N+1)  #Initialize en values for error vector
    x[0] = x0       #Initialize x value
    y[0] = y0     #Initialize y value

    gn[0] = g(x0)   #Initialize g value
    en[0] = g(x0) - y0  #Initialize e value

    #RK4 Loop
    for n in range(0,N):
        x[n+1] = x[n] + h   #Next x value
        k1 = h*f(x[n],y[n])  #Predictor 1
        k2 = h*f(x[n]+0.5*h, y[n]+0.5*k1)    #Predictor 2
        k3 = h*f(x[n]+0.5*h, y[n]+0.5*k2)   #Predictor 3
        k4 = h*f(x[n]+h, y[n]+k3)     #Predictor 4
        y[n+1] = y[n] + (k1+2*k2+2*k3+k4)/6   #Next y value
        gn[n+1] = g(x[n+1])    #Next g value
        en[n+1] = g(x[n+1]) - y[n+1]   #Next e value

    #Calculate and Display PRE using error vector
    NormE = norm(en)
    print('Norm of Error Vector = ', NormE)
    PRE = NormE/norm(gn)*100   #PRE Value
    print('Percent Relative Error = ', PRE)

    #Plot the numerical and analytical solutions
    L = N*h   #[x0,x0+L] = x-axis interval for solution
    xvalues = np.linspace(x0,x0+L,1000)   
    plt.figure(figsize=(12,6))
    plt.plot(xvalues, g(xvalues), 'r', x, y, 'bo')
    plt.xlabel('X')
    plt.legend(('Exact', 'RK4'), loc=0)
    plt.show()

rk4(20)

For some reason the code is not running the way it is supposed to. With a 56% relative error there needs to be some debugging.

1 - Academic Journal: https://core.ac.uk/download/pdf/234815388.pdf

2 - Online Textbook “Applied Pharmaceutics and Pharmacokinetics” https://accesspharmacy.mhmedical.com/content.aspx?bookid=513&sectionid=41488021#:~:text=The%20one%2Dcompartment%20open%20model,like%20a%20single%2C%20uniform%20compartment.

3 - Academic Journal https://www.sciencedirect.com/science/article/pii/S2090506816300136

4 - Pharmacokinetic Information https://accesspharmacy.mhmedical.com/content.aspx?bookid=513&sectionid=41488021#:~:text=The%20one%2Dcompartment%20open%20model,like%20a%20single%2C%20uniform%20compartment.

5 - Academic Journal http://math.colorado.edu/~stade/CLS/Section_3_1.pdf

6 - Academic Journal https://math.dartmouth.edu/~m3cod/ionescuslides/Lecture18.pdf

7 - Image https://www.nuventra.com/resources/blog/pk-compartmental-modeling/

8 - Ch21.1 and 21.2 Notes and Textbook