Two-compartment model (PKPD)

Simulations example for a simple 2-compartment PKPD model¶

This notebook demonstrates a pharmacokinetic/pharmacodynamic (PKPD) simulation using a 2-compartment model with a downstream PD effect. The model equations govern the absorption and distribution of Drug X across two compartments, as well as its pharmacodynamic response. The simulation is carried out using the proprietary SUND simulation toolbox.

In [1]:

import numpy as np
import matplotlib.pyplot as plt
import sund

import numpy as np import matplotlib.pyplot as plt import sund

Defining the PKPD Model¶

The model used in this example is defined in the code cell below. The system consists of:

• An absorption compartment: A,
• A central compartment: C1,
• A peripheral compartment: C2,
• A pharmacodynamic (PD) response variable: R.

The model includes first-order absorption and elimination, inter-compartmental exchange between the central and peripheral compartments, and a sigmoidal inhibitory $I_{max}$ model to represent the pharmacodynamic response.

This code block uses the SUND-toobox to install the model, and then load a model object.

In [2]:

# Define the model
two_compartment_PKPD_model = """\
    ########## NAME
    TwoComp_PKPDmodel
    ########## METADATA
    time_unit = h

    ########## MACROS
    ########## STATES
    ddt_A = -ka*A
    ddt_C1 = ka*(A/V1)-(CL/V1)*C1 - k1*C1 + k2*C2*(V2/V1)
    ddt_C2 = k1*C1*(V1/V2) - k2*(C2)
    ddt_R = kin*(1-Imax*(C2^n/(C2^n+IC50^n)))-kout*R

    A(0) = 0
    C1(0) = 0
    C2(0) = 0
    R(0)=kin/kout
    ########## PARAMETERS
    kin = 1
    kout = 1
    Imax = 1
    IC50 = 1
    n = 1
    k1 = 2.9
    k2 = 0.85

    ########## VARIABLES
    ka = 1.3095 //220 // [1/week] // parameter values recalculated to be per hour rather than per week
    CL = 0.0008333 // [L/week/kg]
    V1 = 0.2 // [L/kg]
    V2 = 0.1 // [L/kg]

    ########## FUNCTIONS
    ########## EVENTS
    ########## OUTPUTS
    ########## INPUTS
    ########## FEATURES
    PK_conc = C2
    PD_Effect = R
    """

# Install the model using the SUND-toolbox.
sund.install_model(two_compartment_PKPD_model)
# Load the installed model
model_object = sund.load_model('TwoComp_PKPDmodel')

# Create a new simulation object for the model using the constructor 'Simulation(model)' from the sund package, and store the object in a variable
simulation_object = sund.Simulation(models = model_object,  time_unit = 'h')

# Define the model two_compartment_PKPD_model = """\ ########## NAME TwoComp_PKPDmodel ########## METADATA time_unit = h ########## MACROS ########## STATES ddt_A = -ka*A ddt_C1 = ka*(A/V1)-(CL/V1)*C1 - k1*C1 + k2*C2*(V2/V1) ddt_C2 = k1*C1*(V1/V2) - k2*(C2) ddt_R = kin*(1-Imax*(C2^n/(C2^n+IC50^n)))-kout*R A(0) = 0 C1(0) = 0 C2(0) = 0 R(0)=kin/kout ########## PARAMETERS kin = 1 kout = 1 Imax = 1 IC50 = 1 n = 1 k1 = 2.9 k2 = 0.85 ########## VARIABLES ka = 1.3095 //220 // [1/week] // parameter values recalculated to be per hour rather than per week CL = 0.0008333 // [L/week/kg] V1 = 0.2 // [L/kg] V2 = 0.1 // [L/kg] ########## FUNCTIONS ########## EVENTS ########## OUTPUTS ########## INPUTS ########## FEATURES PK_conc = C2 PD_Effect = R """ # Install the model using the SUND-toolbox. sund.install_model(two_compartment_PKPD_model) # Load the installed model model_object = sund.load_model('TwoComp_PKPDmodel') # Create a new simulation object for the model using the constructor 'Simulation(model)' from the sund package, and store the object in a variable simulation_object = sund.Simulation(models = model_object, time_unit = 'h')

Model 'TwoComp_PKPDmodel' succesfully installed.

Simulating the PKPD¶

In this example, the drug X is administered once per week for four weeks, mimicking a common clinical trial dosing scenario. The two-compartment model is used to simulate the tissue concentrations, i.e. the concentration in the peripheral compartment (C2), of the drug for different dose levels: 1 mg, 3 mg, 10 mg, and 30 mg. To ensure the clearance of the drug is captured, the simulation covers 10 weeks to capture both the dosing period and the post-dosing observation period.

The following code block defines a function that runs the model simulation for the above described dosing schedule for all defined doses. The inputs to the function are the model object and the parameter values. Outputs from the function include the both the simulated drug concentrations and the PD responses over the 10 week period.

In [3]:

def simulate_multipleDoses(parameterValues: list[float], simulation_object: sund.Simulation) -> tuple[dict, dict, dict]:
    """Simulate multiple doses of the TwoComp_PKPD model.
    This function simulates the model for multiple doses and returns the concentration and PD effect over time.
    The function takes in the parameter values and the simulation object as inputs.
    The function simulates the model for different doses (0, 1, 3, 10, 30 mg) and for different time points (3 days, 4 days, 1 week, and 6 weeks).
    The function returns the concentration and PD effect over time for each dose.
    
    Args:
        parameterValues (list[float]): Model parameter values
        simulation_object (sund.Simulation): The simulation object

    Returns:
        tuple[dict, dict, dict]: The concentration and PD effect over time for each dose
    """
    
    timepoints ={}
    timepoints['3Days'] = np.arange(0, 24*3, 0.1) # timepoints for simulation
    timepoints['4Days'] = np.arange(0, 24*4, 0.1) # timepoints for simulation
    timepoints['1Week'] = np.arange(0, 24*7, 0.1) # timepoints for simulation
    timepoints['6Week'] = np.arange(0, 24*7*6, 0.1) # timepoints for simulation

    conc_sims = {} # Preallocate an empty dictionary for the output
    PD_effect = {} # Create an empty dictionary
    simTimeVector ={}
    
    for doses_size in [0, 1, 3, 10, 30]: #  Loop over the different doses
        conc_sims[str(doses_size)] =[]
        PD_effect[str(doses_size)] = []
        simTimeVector[str(doses_size)] =[0]
        
        simulation_object.reset_states() # Reset the state variables to their initial values
        simulation_object.state_values[[1,2,3]] = [0, 0, (parameterValues[0]/parameterValues[1])] # Set the initial values of the states based on the parameter values. C1(0) = 0, C2(0) = 0, and R(0) = parameterValues[0]/parameterValues[1]
        for simuationTime in ['1Week', '1Week', '1Week', '1Week', '6Week']: # Loop over the different time points
            simulation_object.state_values[0] = doses_size # Set the initial values of the Input state to the current dose size.

            # Simulate the model for the current dose size and time point. 
            simulation_object.simulate(
                time_vector=timepoints[simuationTime], 
                parameter_values=parameterValues, 
                reset=False # The reset=False argument indicates that the state variables should not be reset to their initial values after each simulation.
                ) 

            #Store the resulting simulation in the output variables.    
            conc_sims[str(doses_size)] = np.concatenate((conc_sims[str(doses_size)], simulation_object.feature_data[:,0])) # concatenate the two simulations
            PD_effect[str(doses_size)] = np.concatenate((PD_effect[str(doses_size)], (simulation_object.feature_data[:,1]/(parameterValues[0]/parameterValues[1]))*100)) 
            simTimeVector[str(doses_size)] = np.concatenate((simTimeVector[str(doses_size)], simulation_object.time_vector + simTimeVector[str(doses_size)][-1])) # concatenate the two simulations
        
        simTimeVector[str(doses_size)] = simTimeVector[str(doses_size)][1:] # remove the first element of the time vector
    
    return conc_sims, PD_effect, simTimeVector

def simulate_multipleDoses(parameterValues: list[float], simulation_object: sund.Simulation) -> tuple[dict, dict, dict]: """Simulate multiple doses of the TwoComp_PKPD model. This function simulates the model for multiple doses and returns the concentration and PD effect over time. The function takes in the parameter values and the simulation object as inputs. The function simulates the model for different doses (0, 1, 3, 10, 30 mg) and for different time points (3 days, 4 days, 1 week, and 6 weeks). The function returns the concentration and PD effect over time for each dose. Args: parameterValues (list[float]): Model parameter values simulation_object (sund.Simulation): The simulation object Returns: tuple[dict, dict, dict]: The concentration and PD effect over time for each dose """ timepoints ={} timepoints['3Days'] = np.arange(0, 24*3, 0.1) # timepoints for simulation timepoints['4Days'] = np.arange(0, 24*4, 0.1) # timepoints for simulation timepoints['1Week'] = np.arange(0, 24*7, 0.1) # timepoints for simulation timepoints['6Week'] = np.arange(0, 24*7*6, 0.1) # timepoints for simulation conc_sims = {} # Preallocate an empty dictionary for the output PD_effect = {} # Create an empty dictionary simTimeVector ={} for doses_size in [0, 1, 3, 10, 30]: # Loop over the different doses conc_sims[str(doses_size)] =[] PD_effect[str(doses_size)] = [] simTimeVector[str(doses_size)] =[0] simulation_object.reset_states() # Reset the state variables to their initial values simulation_object.state_values[[1,2,3]] = [0, 0, (parameterValues[0]/parameterValues[1])] # Set the initial values of the states based on the parameter values. C1(0) = 0, C2(0) = 0, and R(0) = parameterValues[0]/parameterValues[1] for simuationTime in ['1Week', '1Week', '1Week', '1Week', '6Week']: # Loop over the different time points simulation_object.state_values[0] = doses_size # Set the initial values of the Input state to the current dose size. # Simulate the model for the current dose size and time point. simulation_object.simulate( time_vector=timepoints[simuationTime], parameter_values=parameterValues, reset=False # The reset=False argument indicates that the state variables should not be reset to their initial values after each simulation. ) #Store the resulting simulation in the output variables. conc_sims[str(doses_size)] = np.concatenate((conc_sims[str(doses_size)], simulation_object.feature_data[:,0])) # concatenate the two simulations PD_effect[str(doses_size)] = np.concatenate((PD_effect[str(doses_size)], (simulation_object.feature_data[:,1]/(parameterValues[0]/parameterValues[1]))*100)) simTimeVector[str(doses_size)] = np.concatenate((simTimeVector[str(doses_size)], simulation_object.time_vector + simTimeVector[str(doses_size)][-1])) # concatenate the two simulations simTimeVector[str(doses_size)] = simTimeVector[str(doses_size)][1:] # remove the first element of the time vector return conc_sims, PD_effect, simTimeVector

The following code block runs the above-defined simulation function with predefined parameter values.

Parameter values used for this simulation example:

Name	Value	Unit
k1	1	[1/h]
k2	3	[1/h]
V1	0.2	[L]
V2	0.1	[L]
CL	0.0008	[L/h]
kin	0.0011	[1/h]
kout	0.0046	[1/h]
ICMax	0.859	[n\a]
IC50	8.91	[mg/L]
n	1.24	[n\a]

In [4]:

parameterValues = [0.00108, 0.00459, 0.859, 8.913, 1.241, 1, 3] # kin, kout, ICMax, IC50, n, k1, k2

# Call the function to simulate the model for dosing protocol and for multiple doses
conc_sims, PD_effect, simTimeVector = simulate_multipleDoses(parameterValues, simulation_object) 

parameterValues = [0.00108, 0.00459, 0.859, 8.913, 1.241, 1, 3] # kin, kout, ICMax, IC50, n, k1, k2 # Call the function to simulate the model for dosing protocol and for multiple doses conc_sims, PD_effect, simTimeVector = simulate_multipleDoses(parameterValues, simulation_object)

Plotting Tissue Concentration Over Time¶

This section visualizes the concentration of Drug X in the peripheral (tissue) compartment (C2) for each dose. These plots help assess how tissue exposure scales with increasing doses.

In [5]:

plot_colour = [('#0000ff'),
               ('#002db3'),
               ('#ff9900'),
               ('#b30000'),
               ('#b36b00')] # Define the colours for the plot


doses_sizes = [1, 3, 10, 30] # Define the dose sizes for the plot

plt.figure()

for i in [0,1,2,3,4]: #plot a vertical line for each dosing time.
    plt.axvline(x=i, color=(0.8, 0, 0.4, 0.6), linestyle='-') # plot a vertical line for each dosing time. Note the the time vector is divided by 24*7 to convert the time from hours to weeks.
plt.text(0.5, 0.06, 'Dosing', rotation=0, ha='center', va='center', fontsize=12, color=(0.8, 0, 0.4, 1), bbox=dict(facecolor='w', alpha=1, edgecolor = 'none')) # add text to the plot

for i in doses_sizes: # plot the tissue concentration for each dose size
    plt.plot(simTimeVector[str(i)]/(24*7),
              conc_sims[str(i)],
              label=f'{i} mg',
              color=plot_colour[np.where(np.array(doses_sizes) == i)[0][0]],
              linewidth = 3) # plot the tissue concentration for each dose size. Note the the time vector is divided by 24*7 to convert the time from hours to weeks.

plt.xlabel('Time (weeks)')
plt.ylabel('Concentration (ug/g tissue)')
plt.yscale('log')
plt.grid(which='major', color='0.8', linestyle='-')
plt.grid(which='minor', color='0.8', linestyle='--')
plt.legend(title='Dose Size (mg)', loc='upper right')
plt.title('Concentration in tissue over time')

plot_colour = [('#0000ff'), ('#002db3'), ('#ff9900'), ('#b30000'), ('#b36b00')] # Define the colours for the plot doses_sizes = [1, 3, 10, 30] # Define the dose sizes for the plot plt.figure() for i in [0,1,2,3,4]: #plot a vertical line for each dosing time. plt.axvline(x=i, color=(0.8, 0, 0.4, 0.6), linestyle='-') # plot a vertical line for each dosing time. Note the the time vector is divided by 24*7 to convert the time from hours to weeks. plt.text(0.5, 0.06, 'Dosing', rotation=0, ha='center', va='center', fontsize=12, color=(0.8, 0, 0.4, 1), bbox=dict(facecolor='w', alpha=1, edgecolor = 'none')) # add text to the plot for i in doses_sizes: # plot the tissue concentration for each dose size plt.plot(simTimeVector[str(i)]/(24*7), conc_sims[str(i)], label=f'{i} mg', color=plot_colour[np.where(np.array(doses_sizes) == i)[0][0]], linewidth = 3) # plot the tissue concentration for each dose size. Note the the time vector is divided by 24*7 to convert the time from hours to weeks. plt.xlabel('Time (weeks)') plt.ylabel('Concentration (ug/g tissue)') plt.yscale('log') plt.grid(which='major', color='0.8', linestyle='-') plt.grid(which='minor', color='0.8', linestyle='--') plt.legend(title='Dose Size (mg)', loc='upper right') plt.title('Concentration in tissue over time')

Out[5]:

Text(0.5, 1.0, 'Concentration in tissue over time')

Plotting the PD Effect Over Time¶

The PD response (R) is plotted over the 10-week simulation for each dose level. These plots demonstrate the temporal dynamics of the drug's effect and the impact of dose escalation on PD inhibition.

In [6]:

plt.figure()
plt.axhline(y=100, color='k', linestyle='--', label='Control')

for i in [0,1,2,3,4]: #plot a vertical line for each dosing time.
    plt.axvline(x=i, color=(0.8, 0, 0.4, 0.6), linestyle='-') # plot a vertical line for each dosing time. Note the the time vector is divided by 24*7 to convert the time from hours to weeks.
plt.text(0.5, 6, 'Dosing', rotation=0, ha='center', va='center', fontsize=12, color=(0.8, 0, 0.4, 1), bbox=dict(facecolor='w', alpha=1, edgecolor = 'none')) # add text to the plot

for i in doses_sizes:
    plt.plot(simTimeVector[str(i)]/(24*7),
              PD_effect[str(i)], 
              label=f'{i} mg',
              color=plot_colour[np.where(np.array(doses_sizes) == i)[0][0]],
              linewidth = 3) # plot the PD effect for each dose size. Note the the time vector is divided by 24*7 to convert the time from hours to weeks.


plt.xlabel('Time (weeks)')
plt.ylabel('PD  Effect ("%" of controll)')
plt.grid(which='major', color='0.8', linestyle='-')
plt.legend(title='Drug effect', loc='upper right')
plt.title('PD Effect over time')
plt.ylim(0, 110)

plt.figure() plt.axhline(y=100, color='k', linestyle='--', label='Control') for i in [0,1,2,3,4]: #plot a vertical line for each dosing time. plt.axvline(x=i, color=(0.8, 0, 0.4, 0.6), linestyle='-') # plot a vertical line for each dosing time. Note the the time vector is divided by 24*7 to convert the time from hours to weeks. plt.text(0.5, 6, 'Dosing', rotation=0, ha='center', va='center', fontsize=12, color=(0.8, 0, 0.4, 1), bbox=dict(facecolor='w', alpha=1, edgecolor = 'none')) # add text to the plot for i in doses_sizes: plt.plot(simTimeVector[str(i)]/(24*7), PD_effect[str(i)], label=f'{i} mg', color=plot_colour[np.where(np.array(doses_sizes) == i)[0][0]], linewidth = 3) # plot the PD effect for each dose size. Note the the time vector is divided by 24*7 to convert the time from hours to weeks. plt.xlabel('Time (weeks)') plt.ylabel('PD Effect ("%" of controll)') plt.grid(which='major', color='0.8', linestyle='-') plt.legend(title='Drug effect', loc='upper right') plt.title('PD Effect over time') plt.ylim(0, 110)

Out[6]:

(0.0, 110.0)

Tissue Concentration vs PD Effect¶

This plot illustrates the relationship between the concentration of Drug X in the tissue compartment (C2) and its corresponding PD effect (R). The characteristic sigmoidal shape of the curve reflects the inhibitory Imax model, where higher tissue concentrations lead to a greater suppression of the response, up to a maximum inhibition (Imax). This visualization helps link drug exposure directly to pharmacodynamic outcome, providing insight into the effective concentration range and supporting dose optimization.

In [7]:

plt.figure()

for i in doses_sizes:
    plt.plot(conc_sims[str(i)][-(24*7*5*10):-1],
             PD_effect[str(i)][-(24*7*5*10):-1],
             label=f'{i} mg',
             color=plot_colour[np.where(np.array(doses_sizes) == i)[0][0]],
             linewidth = 3)

plt.ylabel('PD  Effect ("%" of controll)')
plt.xlabel('Concentration (ug/g tissue)')
plt.xscale('log')
plt.grid(which='major', color='0.8', linestyle='-')
plt.grid(which='minor', color='0.8', linestyle='--')
plt.title('Tissue Concentration vs PD Effect')
plt.legend(title='Dose Size (mg)', loc='upper right')

plt.show()

plt.figure() for i in doses_sizes: plt.plot(conc_sims[str(i)][-(24*7*5*10):-1], PD_effect[str(i)][-(24*7*5*10):-1], label=f'{i} mg', color=plot_colour[np.where(np.array(doses_sizes) == i)[0][0]], linewidth = 3) plt.ylabel('PD Effect ("%" of controll)') plt.xlabel('Concentration (ug/g tissue)') plt.xscale('log') plt.grid(which='major', color='0.8', linestyle='-') plt.grid(which='minor', color='0.8', linestyle='--') plt.title('Tissue Concentration vs PD Effect') plt.legend(title='Dose Size (mg)', loc='upper right') plt.show()