PK/PD Modeling and Simulation to Guide Dosing Strategy for Antibiotics

This example shows how to perform a Monte Carlo simulation of a pharmacokinetic/pharmacodynamic (PK/PD) model for an antibacterial agent. This example is adapted from Katsube et al. [1] This example also shows how to use the SimBiology exported model to perform parameter scans in parallel.

This example requires Statistics Toolbox™. The performance can be improved if you have the Parallel Computing Toolbox™ software.

Background

Katsube et al. [1] used a PK/PD modeling and simulation approach to determine the most effective dosage regimens for doripenem, a carbapenem antibiotic. The objectives of their study were:

  • Develop a PK/PD model to describe the antibacterial effect of doripenem against several Pseudomonas aeruginosa strains

  • Use Monte Carlo simulations to compare the efficacy of four common antibiotic dosage regimes, and to determine the most effective dosing strategy

  • Investigate the effect of renal function on the antibacterial efficacy of the treatments

In this example, we will implement the antibacterial PK/PD model developed by Katsube et al. [1] in SimBiology®, and replicate the results of the Monte Carlo simulation described in their work.

References

[1] T. Katsube, Y. Yano, T. Wajima, Y. Yamano and M. Takano. Pharmacokinetic/pharmacodynamic modeling and simulation to determine effective dosage regimens for doripenem. Journal of Pharmaceutical Sciences (2010) 99(5), 2483-91.

PK/PD Model

Katsube et al. assumed a two-compartment infusion model with linear elimination from the central compartment to describe the pharmacokinetics of the doripenem. For the bacterial growth model, they assumed that the total bacterial population is comprised of drug-susceptible growing cells and drug-insensitive resting cells. The antibacterial effect of the drug was included in the killing rate of the bacteria via a simple Emax type model:

$$ Killing Rate = \frac{Kmax*[Drug]*[Growing]}{KC50 + [Drug]} $$

where [Drug] is the concentration (ug/ml) of the drug in the central compartment, and [Growing] is the count of the growing bacterial population in CFU/ml (CFU = Colony Forming Units). Kmax is the maximal killing rate constant (1/hour) and KC50 is the Michaelis-Menten rate constant (ug/ml).

A graphical view of the SimBiology implementation of the model is shown below.

% Load model
sbioloadproject('AntibacterialPKPD.sbproj', 'm1') ;

Dosage Regimens

Katsube et al. simulated the model using four common antibiotic dosage strategies.

  • 250 mg two times a day (b.i.d.)

  • 250 mg three times a day (t.i.d.)

  • 500 mg two times a day (b.i.d.)

  • 500 mg three times a day (t.i.d.)

Infusion dosing was used in all four dosages regimens, and infusion time was set to 30 minutes. In SimBiology, these dosage regimens have been implemented as dose objects.

% Select dose objects in the model
doseNames = {'250 mg bid', '250 mg tid', '500 mg bid', '500 mg tid'};
for iDoseGrp = 1:length(doseNames)
    doseRegimens(iDoseGrp) = sbioselect(m1, 'Name', doseNames{iDoseGrp}) ;
end

Description of the Virtual Population

A virtual population of individuals was generated based on the distribution of demographic variables and PK/PD parameters. The type of distribution and the values of the distribution parameters were based on data from earlier clinical trials of doripenem conducted in Japan.

Note: In [1], 5,000 virtual patients were simulated in each dosage group. In this example, we will use 1,000 patients in each group. To simulate a different population size, change the value of nPatients below.

% Setup
nPatients    = 1000         ; % Number of patients per dosage group
nDoseGrps    = 4            ; % Number of tested dosage regimens

Distribution of Demographic Variables:

Weight (Wt) and age (Age) were sampled from a normal distribution with a mean of 51.6 kg and 71.8 years, respectively, and a standard deviation of 11.8 kg and 11.9 years, respectively. 26% of the population was assumed to be female. Serum creatinine levels (Scr) were sampled from a lognormal distribution with a typical value of 0.82 mg/dL, and coefficient of variation (CV) of 32%. The creatinine clearance rates (CrCL) were calculated using the Cockcroft-Gault equation.

% Note: The inputs to the lognrnd function are the mean (mu) and standard
% deviation (sigma) of the associated normal distribution. Here and
% throughout the example, mu and sigma were calculated from the reported
% typical value and coefficient of the lognormal distriution. See the
% lognstat documentation for more information.

% Patient demographics
Wt   = normrnd(51.6     , 11.8  , nPatients , nDoseGrps )   ; % units: kg
Age  = normrnd(71.8     , 11.9  , nPatients , nDoseGrps )   ; % units: years
Scr  = lognrnd(-0.2485  ,0.3197 , nPatients , nDoseGrps )   ; % units: ml/minute

% Gender ratio
id       = 1:nPatients*nDoseGrps                            ;
idFemale = randsample(id , round(0.26*nDoseGrps*nPatients)) ; % 26% Female

% Creatinine Clearance (using Cockcroft-Gault equation)
CrCL            = (140 - Age).*Wt./(Scr*72)                 ; % units: ml/minute
CrCL(idFemale)  = CrCL(idFemale)*0.85                       ; % multiply by 0.85 for females

Distribution of Pharmacokinetic (PK) parameters:

PK parameters, Central, k12 and k21, were sampled from a lognormal distribution with typical values of 7.64 liters, 1.59 1/hour and 2.26 1/hour, respectively, and a 20% coefficient of variation (CV). Central is the distribution volume of the central compartment, and k12 and k21 are transfer rate constants between the Central and the Peripheral compartments. The drug clearance rate, CL, was assumed to depend linearly on the creatinine clearance rate via the following equation:

$$CL = 1.07*CrCL + 45.6 + \varepsilon$$

where $\varepsilon$ is the additive residual error sampled from a normal distribution with a mean of 0 ml/minute and standard deviation of 22 ml/minute.

% PK parameters
Central = lognrnd(2.01   , 0.2 , nPatients, nDoseGrps)         ; % units: liter
k12     = lognrnd(0.4441 , 0.2 , nPatients, nDoseGrps)         ; % units: 1/hour
k21     = lognrnd(0.7958 , 0.2 , nPatients, nDoseGrps)         ; % units: 1/hour
CL      = 1.07*CrCL + 45.6 + normrnd(0,22,nPatients,nDoseGrps) ; % units: ml/minute

Distribution of Pharmacodynamic (PD) parameters:

Growing-to-resting transformation rate constants, k1 and k2, were sampled from a lognormal distribution with typical value of 5.59e-5 and 0.0297 1/hour, respectively, each with a CV of 20%. Kmax was sampled from a lognormal distribution with a typical value of 3.5 1/hour and 15.9% CV. Katsube et al. assumed that values k1, k2 and Kmax were independent of the bacterial strain being treated. The value of Beta, the net growth rate constant, was fixed at 1.5 1/hour.

Based on experiments with several strains, the authors concluded that the value of KC50 was linearly dependent on the minimum inhibition concentration (MIC) of bacterial strain via the following equation.

$$ ln(KC50) = -1.91 + 0.898*ln(MIC) + \varepsilon $$

where $\varepsilon$ is the additive residual error sampled from a normal distribution with a mean of 0 and standard deviation of 1.06 ug/ml. In the simulation, the MIC values were sampled from a discrete distribution, and the KC50 value was calculated for the selected MIC using the above equation.

% Discrete distribution of MIC values based on 71 P. aeruginosa strains
micValue  = [0.0625, 0.125, 0.25, 0.5 , 1 , 2 , 4 , 8 , 16 , 32 ] ;
micFreq   = [   5  ,    8 ,  9  , 14  , 7 , 8 , 9 , 5 , 2  , 4  ] ;

k1   = lognrnd(-9.8116, 0.2  , nPatients, nDoseGrps)        ; % units: 1/hour
k2   = lognrnd(-3.5362, 0.2  , nPatients, nDoseGrps)        ; % units: 1/hour
Kmax = lognrnd( 1.2332, 0.159, nPatients, nDoseGrps)        ; % units: 1/hour

% Sample MIC values from a discrete distribution using randsample
MIC = nan(nPatients, nDoseGrps) ; % preallocate
for iDoseGrp = 1:nDoseGrps
    MIC(:, iDoseGrp) = randsample(micValue , nPatients, true , micFreq);
end

KC50 = exp(-1.91 + 0.898*log(MIC) + 1.06*randn(nPatients , nDoseGrps)) ; % units: microgram/milliliter

Simulation Setup & Design

You will convert the model to a SimBiology exported model, which facilitates performing parameter scans in parallel. When you export the model, you choose which species, parameters, or compartments can be varied. In this example, you will vary 8 parameters, Central, k12, k21, CL, k1, k2, Kmax, and KC50.

% Select the parameters you want to vary.
paramNames = {'Central', 'k12', 'k21', 'CL', 'k1', 'k2', 'Kmax', 'KC50'};
for iParam = 1:length(paramNames)
    parameters(iParam) = sbioselect(m1, 'Name', paramNames{iParam});
end

% Created the exported model, using the selected parameters and doses
exportedModel = export(m1, parameters, doseRegimens);

% Accelerate the model
accelerate(exportedModel);

For all dosage scenarios, the model was simulated until t = 2 weeks from the time of the first dose. Total bacterial count, CFU, was sampled every 24 hours (once a day) for the entire duration of the dosage regimen.

tObs      = 0:24:336       ; % hour
nTPoints  = length(tObs)   ; % Number of sampling points

% Specify that the simulation should report these output times
exportedModel.SimulationOptions.OutputTimes = tObs;

Start the Worker Pool

If you have the Parallel Computing Toolbox software, you can use local workers to distribute each simulation to a different worker. If you also have the MATLAB® Distributed Computing Server™, you can run your simulations on a cluster.

You can create a pool of workers using the current cluster profile with this command:

parpool

Monte Carlo Simulation of Patients with Severe Infection

The antibacterial efficacy of a drug can be measured using different PK/PD indices. Katsube et al. set the criterion for bacterial elimination at log10(CFU) < 0, where CFU is the total bacterial count. The efficacy of each dose regimen was measured as the fraction of the population that achieved the success criteria in the dosage group. This efficacy metric, Pr{log10(CFU) < 0}, was tracked as a function of time for each dosage group.

In their simulation studies, the authors investigated the efficacy of the dosage regimens on two classes of patients:

  • Moderate infection (Initial bacterial count = 1e4 CFU/ml)

  • Severe infection (Initial bacterial count = 1e7 CFU/ml)

In this example, we will replicate the results for the severe infection case only. Note that you can easily simulate the other scenario, patients with moderate infection, by changing the initial amount of bacterial count in the model to 1e4 CFU/ml.

% Preallocate
log10CFU     = cell(1,nDoseGrps) ;


for iDoseGrp = 1:nDoseGrps

    % Select the exported doess
    currentDose = getdose(exportedModel, doseNames{iDoseGrp});

    cfu = nan(nTPoints , nPatients)    ; % preallocate

    disp(['Simulating group ',  num2str(iDoseGrp), ' ... '])
    parfor iPatient = 1:nPatients

       % Use parameter values for current patient
       % Define the parameters in the same order used when exporting
       parameterValues = [
           Central(iPatient , iDoseGrp)
           k12(iPatient, iDoseGrp)
           k21(iPatient, iDoseGrp)
           CL(iPatient , iDoseGrp)
           k1(iPatient , iDoseGrp)
           k2(iPatient , iDoseGrp)
           Kmax(iPatient, iDoseGrp)
           KC50(iPatient, iDoseGrp)
           ];

        % Simulate
        simData = simulate(exportedModel, parameterValues, currentDose) ;

        % Extract bacterial count data for Growing and Resting population
        [~, bactCount]  = selectbyname(simData, {'Growing', 'Resting'}) ;

        % Sum of growing and resting bacterial
        cfu(:, iPatient) = sum(bactCount, 2) ;

    end

    % Calculate log10(CFU)
    log10CFU{iDoseGrp}  = log10(cfu) ;

end

% Save results
log10CFU_250bid = log10CFU{1} ;
log10CFU_250tid = log10CFU{2} ;
log10CFU_500bid = log10CFU{3} ;
log10CFU_500tid = log10CFU{4} ;
Simulating group 1 ... 
Simulating group 2 ... 
Simulating group 3 ... 
Simulating group 4 ... 

Clean Up the Worker Pool

If you previously called parpool to start some workers, be sure to close them using the following command:

delete(gcp)

The function gcp (get current pool) returns the pool created by parpool, and the delete command closes the pool.

Time Course Profiles of Bacterial Counts

We plot the median (in red) and percentile (shaded) profiles of the log10(CFU) levels for all four dosage regimens. Observe that in all four groups, the median time course profile shows that bacterial eradication is complete before the end of the treatment period (336 hours). However, it is evident from the higher percentile profiles that the treatments are not successful for all patients. The 95th and 90th percentile profiles also indicate that dosing a lower amount with a higher frequency (250 tid) is more effective than less frequent dosing with higher amount (500 bid).

hax1(1) = subplot(2,2,1)
plotCFUCount(tObs, log10CFU_250bid, 'a. Dose 250 bid' )
hax1(2) = subplot(2,2,2)
plotCFUCount(tObs, log10CFU_250tid, 'b. Dose 250 tid' )
hax1(3) = subplot(2,2,3)
plotCFUCount(tObs, log10CFU_500bid, 'c. Dose 500 bid' )
hax1(4) = subplot(2,2,4)
plotCFUCount(tObs, log10CFU_500tid, 'd. Dose 500 tid' )

% Link subplot axes
linkaxes(hax1)
hax1 = 

  Axes with properties:

             XLim: [0 1]
             YLim: [0 1]
           XScale: 'linear'
           YScale: 'linear'
    GridLineStyle: '-'
         Position: [0.1300 0.5838 0.3347 0.3412]
            Units: 'normalized'

  Use GET to show all properties


hax1 = 

  1x2 Axes array:

    Axes    Axes


hax1 = 

  1x3 Axes array:

    Axes    Axes    Axes


hax1 = 

  1x4 Axes array:

    Axes    Axes    Axes    Axes

Effect of Renal Function on Antibacterial Activity

Finally, the authors compared the efficacy profiles of the dosages regimens as a function of the renal function. They classified the patients into four renal function groups based on the creatinine clearance rates (CrCL):

  • Creatinine Clearance Group 1: CrCL < 30

  • Creatinine Clearance Group 2: 30 <= CrCL < 50

  • Creatinine Clearance Group 3: 50 <= CrCL < 70

  • Creatinine Clearance Group 4: CrCL >= 70

The next figure shows the effect of renal function (creatinine clearance rate) on the antibacterial efficacy of the four dosage regimens. Observe that in the normal renal function group (CrCL >= 70), the efficacy profiles of the four treatment strategies are significantly different from each other. In this case, the 500 mg t.i.d. dose is much more effective than the other regimens. In contrast, simulations involving patients with renal dysfunction (CrCL < 30 and 30 <= CrCL < 50), we don't see much difference between the treatment groups. This indicates that for patients with a renal dysfunction, a less intense or less frequent dosing strategy would work almost as well as a dosing strategy with higher frequency or dosing amount.

% Preallocate
idCrCLGrp   = false(nPatients, nDoseGrps) ;

% Line Style
ls = {'bd:', 'b*:', 'rd:', 'r*:'} ;

titleStr = {'CL_c_r < 30'           , ...
            '30 <= CL_c_r < 50'     , ...
            '50 <= CL_c_r < 70'     , ...
            'CL_c_r > 70'            };

f = figure;
f.Color = 'w'

for iCrCLGrp = 1:4 % Creatinine Clearance Groups

    hax2(iCrCLGrp) = subplot(2,2, iCrCLGrp) ;
    title( titleStr{iCrCLGrp}   ) ;
    ylabel('Prob(log10CFU < 0)' ) ;
    xlabel('Time (hours)'       ) ;

end

% Set axes properties
set(hax2,   'XTick'        , 0:48:336      , ...
            'XTickLabel'   , 0:48:336      , ...
            'Ylim'         , [0 1]         , ...
            'Xlim'         , [0 336]       , ...
            'NextPlot'     , 'add'         , ...
            'Box'          , 'on'          );

% Plot results by renal function group:
for iDoseGrp = 1:nDoseGrps

    % Extract indices for renal function
    idCrCLGrp(:, 1) = CrCL(:,iDoseGrp) < 30                                    ;
    idCrCLGrp(:, 2) = CrCL(:,iDoseGrp) >= 30 & CrCL(:,iDoseGrp) < 50           ;
    idCrCLGrp(:, 3) = CrCL(:,iDoseGrp) >= 50 & CrCL(:,iDoseGrp) < 70           ;
    idCrCLGrp(:, 4) = CrCL(:,iDoseGrp) >= 70                                   ;

    for iCrCLGrp = 1:4 % Creatinine Clearance Groups

        % Calculate probability
        Pr = sum( ( log10CFU{iDoseGrp}(:, idCrCLGrp(:, iCrCLGrp)') < 0) , 2 )/sum(idCrCLGrp(:,iCrCLGrp)) ;

        % Plot
        plot(hax2(iCrCLGrp), tObs, Pr , ls{iDoseGrp}, 'MarkerSize', 7)
    end

end

legend(hax2(4), {'250 b.i.d.', '250 t.i.d.', '500 b.i.d.', '500 t.i.d.'} )
legend location NorthWest
legend boxoff

linkaxes(hax2)
f = 

  Figure (1) with properties:

      Number: 1
        Name: ''
       Color: [1 1 1]
    Position: [360 502 560 420]
       Units: 'pixels'

  Use GET to show all properties

Was this topic helpful?