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 a 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.
The AUC refers to the area under the curve. This is a metric used to calculate the total dosage delivered to the patient. It is just the integral
\[ AUC=\int_0^T C(t)dt. \]
\(C_{max}\) refers to the maximum concentration experienced by the patient.
In a real world setting constraints might be placed on AUC and \(C_{max}\) (e.g. we can play around with the timing of delivery but the max concentration must never exceed some critical value).
#| 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.sidebar(
ui.input_slider(id="k",label="Drug half life (h)",min=1,max=12,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=18.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.output_plot("plot"),
),
)
def server(input, output, session):
@render.plot
def plot():
fig, ax = plt.subplots()
#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.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.set_xlabel('$t$ (h)')
ax.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)1.2 A specific modelling problem
Here is a modelling challenge?
You work for a drug company who are designing a new drug. The drug has to be administered twice a day, at 12 hour intervals, and the administered drug concentration is fixed at a value of 4 (for reasons outside your control).
Can you use Figure 2 to answer the following questions:
- What must the half-life of the drug be such that the concentration does not fall below a critical level 2?
- How rigidly ought the patient stick the 12 hour intervals (e.g. will the minimal drug concentration be reached if the drugs are taken with a gap of 10 hours)?
- What happens to the AUC as the half-life of the drug increases?
#| 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.sidebar(
ui.input_slider(id="k",label="Drug half life (h)",min=1,max=12,value=6,step=1),
ui.input_slider(id="t1",label="First time (h)",min=0.0,max=10.0,value=6.0,step=1.0),
ui.input_slider(id="t2",label="Second time (h)",min=12.0,max=24.0,value=16.0,step=1.0),
),
ui.output_plot("plot"),
),
)
def server(input, output, session):
@render.plot
def plot():
fig, ax = plt.subplots()
#ax.set_ylim([-2, 2])
# Filter fata
drug_half_life=float(input.k())
t1=float(input.t1())
t2=float(input.t2())
c1=4.0
c2=4.0
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,0]
t_sort_day=[0.0,t1,t2]
#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.plot(t, N,'b')
AUC=AUC+simpson(N, x=t)
Cmax=np.max([Cmax,np.max(N)])
#ax[1].set_xlabel('$t$ (h)')
ax.set_xlabel('$t$ (h)')
ax.set_ylabel('$C(t)$')
ax.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)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.