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 and Machine Learning Toolbox™. The performance can be improved if you have the Parallel Computing Toolbox™ software.

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.

[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.

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:

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') ;

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

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:

where 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.

where
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

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;

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

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 ... Starting parallel pool (parpool) using the 'local' profile ... connected to 12 workers. Simulating group 2 ... Simulating group 3 ... Simulating group 4 ...

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.

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 = 1×2 Axes array: Axes Axes hax1 = 1×3 Axes array: Axes Axes Axes hax1 = 1×4 Axes array: Axes Axes Axes Axes

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`

< 30Creatinine Clearance Group 2: 30 <=

`CrCL`

< 50Creatinine Clearance Group 3: 50 <=

`CrCL`

< 70Creatinine 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 (2) with properties: Number: 2 Name: '' Color: [1 1 1] Position: [348 480 583 437] Units: 'pixels' Use GET to show all properties

Was this topic helpful?