1 Pharmaco kinetic modelling

We can use calculus to study how to optimise drug delivery.

Suppose you are working as part of a team who are designing clinical trials for a drug. It is proposed that patients will orally ingest the drug at a number of times every day. The delivery protocol needs to be optimised so that the drug concentration is as smooth as possible whilst maintaining the total dose constant.

A first proposal is that the patient will take the drug 12 units of the drug once per day.

1.1 Model development

Let \(t\) represent time and \(C(t)\) represent the drug concentration in the blood stream.

Suppose that the drug has a half-life of 6 hours and that at some set of times, \(\{t_i\}\), a concentration of drug \(c_i\) is delivered to the patient.

Consider the model \[ \frac{dC}{dt}=-kC+\sum_i c_i \delta(t_i), \quad N(0)=0. \]

The first term describes the linear degradation of the drug. The second term represent instantaneous delivery of drug at a prescribed set of times.

In Figure 1 you can explore the effect of different treatment protocols.

#| standalone: true
#| components: [viewer]
#| viewerHeight: 500

from shiny import App, Inputs, Outputs, Session, render, ui
from shiny import reactive

import numpy as np
from pathlib import Path
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from scipy.integrate import simpson

app_ui = ui.page_fluid(
    ui.layout_sidebar(
        ui.panel_sidebar(
    ui.input_slider(id="k",label="Drug half life (h)",min=1,max=8,value=6,step=1),
    ui.input_slider(id="t1",label="First time (h)",min=0.0,max=10.0,value=7.0,step=1.0),
    ui.input_slider(id="t2",label="Second time (h)",min=10.0,max=17.0,value=12.0,step=1.0),
    ui.input_slider(id="t3",label="Third time (h)",min=17.0,max=22.0,value=18.0,step=1.0),
    ui.input_slider(id="c1",label="First conc.",min=0.0,max=12.0,value=4.0,step=1.0),
    ui.input_slider(id="c2",label="Second conc.",min=0.0,max=12.0,value=4.0,step=1.0),
    ui.input_slider(id="c3",label="Third conc.",min=0.0,max=12.0,value=4.0,step=1.0),
    
                 
          
            ),

        ui.panel_main(ui.output_plot("plot"),),
    ),
)

def server(input, output, session):
    
    @render.plot
    def plot():
        fig, ax = plt.subplots(2,1)
        #ax.set_ylim([-2, 2])
        # Filter fata
        
        
        drug_half_life=float(input.k())
        t1=float(input.t1())
        t2=float(input.t2())
        t3=float(input.t3())
        c1=float(input.c1())
        c2=float(input.c2())
        c3=float(input.c3())

        k=np.log(2)/drug_half_life

        
        # Define rhs of LV ODEs
        def rhs_pop_model(x,t,k,r):
          rhs=np.zeros_like(x,dtype=float)
          N=x[0]
          dN_dt=-k*N
          rhs[0]=dN_dt
          return rhs

        c1_default=12
        c2_default=0.0
        c3_default=0.0

        T_final=96
        num_days=int(T_final/24)

        drug_conc_default_day=[c1_default,c2_default,c3_default,0]
        drug_conc_day=[c1,c2,c3,0]
        t_sort_day=[0.0,t1,t2,t3]

        drug_conc_default=[]
        drug_conc=[]
        t_sort=[]
        Cmax=0.0
        for i in range(num_days):
            drug_conc_default=drug_conc_default+drug_conc_default_day
            drug_conc=drug_conc+drug_conc_day

            t_sort_day_i = [x + i*24.0 for x in t_sort_day]


            t_sort=t_sort+t_sort_day_i
        
        t_sort.append(T_final)        #t_sort=[0.0,t1,t2,t3,24.0,t1+24.0,t2+24.0,t3+24,48.0]
        N_0=0.0
        AUC=0.0

        for i in range(len(t_sort)-1):
        # Define discretised t domain
            t = np.linspace(t_sort[i], t_sort[i+1], 1000)
            print(t)

            # define initial conditions
            init_cond=[N_0]
        
            # Compute numerical solution of ODEs
            sol1 = odeint(rhs_pop_model, init_cond,t,args=(k,1))

            # Plot results
            N=sol1[:,0]

            N_0=N[-1]
            N_0=N_0+drug_conc[i]
            ax[0].plot(t, N,'b')
            AUC=AUC+simpson(N, x=t)
            Cmax=np.max([Cmax,np.max(N)])

        
        N_0=0.0
        AUC_t=0.0
        cmax_t=0.0
        for i in range(len(t_sort)-1):
        # Define discretised t domain
            t = np.linspace(t_sort[i], t_sort[i+1], 1000)

            # define initial conditions
            init_cond=[N_0]
        
            # Compute numerical solution of ODEs
            sol1 = odeint(rhs_pop_model, init_cond,t,args=(k,1))

            # Plot results
            N=sol1[:,0]

            N_0=N[-1]
            N_0=N_0+drug_conc_default[i]
        
        
            ax[1].plot(t, N,'r')
            AUC_t=AUC_t+simpson(N, x=t)
            cmax_t=np.max([cmax_t,np.max(N)])


        ax[1].set_xlabel('$t$ (h)')
        ax[0].set_xlabel('$t$ (h)')

        ax[0].set_title('AUC = ' + str(int(AUC)) +', Cmax = ' + str(int(Cmax)))
        ax[1].set_title('AUC = ' + str(int(AUC_t)) + ', Cmax = ' + str(int(cmax_t)))

        #ax.set_ylim([0,max_inf*1.4])

        plt.grid()
        plt.show()
    
app = App(app_ui, server)

Figure 1

Explore the effect of the drug half life on hte concentration profile.
Is taking the drug three times per day ‘better’ than taking it twice?

Note

At Dundee, core concepts from calculus (e.g. differential equations) are studied in the modules Maths 1A and Maths 1B and developed further in the modules Maths 2A and Maths 2B.

At Level 2 in the modules Computer algebra and dynamical systems and Introduction to Programming you would be introduced to techniques that are used to compute numerical solutions to differential equations.

At Level 3 in the module Mathematical Biology you would learn how to formulate and study mathematical models of biological systems.

You can find out more about these modules here.