How do I calculate a confidence interval for a simulation from a PK model?

I am using a PK model in SimBiology to simulate drug concentrations following a drug-drug interaction. I am plotting the concentrations over time but would also like to calculate a confidence interval for the data. Can this be done in SimBiology?

Answer 1

[Sietse Braakman]

Since R2017b, SimBiology has two functions to calculate confidence intervals. At the time this question was asked, these functions hadn't been implemented yet.

With sbioparameterci you can calculate confidence intervals on parameter estimates using Gaussian (post-hoc), bootstrap and profilelikelihood methods.

With sbiopredictionci you can calculate the confidence interval on your model predictions - given the data you used to calibrate your model. In other words, the confidences interval on your model predictions therefore provides a 95% CI (or whatever CI you specify) on within which the model predictions will fall given the data you used to calibrate the model.

Both sbioparameterci and sbiopredictionci are function you can call in the command window after you have performed your parameter estimation. The functions take the fit-results object that is the output of a parameter estimation as their input argument.

A third approach you can apply is to sample your parameter space, based on the distributions of your (estimated) parameters. You can then run your model for each of these samples (sometimes called a virtual population), and calculate the confidence interval resulting from these simulations. I would recommend creating a virtual population using lognormal distribution and simulating this virtual population using a simfunction. As a start, you can use the function pbpk_exploration.m in this file_exchange entry. Instead of writing the subject to a folder, you can consider preallocating a cell array before the parfor loop (results = cell(nRuns,1)) and then within each simulation (i.e. iteration of the parfor loop), save the results of that simulation in the results (results(i) = subject).

Once you have simulated your entire population, you can use the function prctile to calculate the percentiles on the data. The data will need a bit of reworking using cell2mat, reshape and cellfun functions. I have included the code below. You will have to download the model ('Nanopartilcle distributions.sbproj') from the file exchange entry linked below.

% Simulate model for a number of virtual subjects and calculate confidence
% intervals on model simulations

% Copyright 2017 The MathWorks, Inc.

%% Load Nanoparticle PBPK model
% Model based on:
% “Elucidating the in vivo fate of nanocrystals using a physiologically based
% pharmacokinetic model: a case study with the anticancer agent SNX-2112”
% Dong D et al.
% International Journal of Nanomedicine, 2015 Volume 10, 2521-2535. 
prj = sbioloadproject('Nanoparticle distribution.sbproj');
% You will be able to download this model from the MATLAB file exchange
model = prj.m1;

%% Generate samples in this parameter space.
% determine which parameters you want to vary for each virtual subject
parameterObjs = model.Parameters([9 11 12 14:17]);
% get the mean value for each parameter
means = arrayfun(@(x)x.Value, parameterObjs)'; 
% SD = 25% of mean value
sigma = means/4; 
% Generate samples
nRuns = 500;
% Create random, normally distributed samples from the Standard Normal
% Distribution N(0,1). This results in Z-values.
ksamples = mvnrnd(zeros(size(means)), eye(size(means,2)), nRuns);
muLN = log((means.^2)./sqrt(sigma.^2+means.^2)); % for lognormal distribution
sigmaLN = sqrt(log(sigma.^2./(means.^2)+1));     % for lognormal distribution
% Transform Z-values to virtual subjects using the lognormal mu and sigma.
ksamples = exp(ones(nRuns,size(means,2)).*muLN + ksamples.*sigmaLN);

%% Construct a SimFunction
% Construct a trimmed-down interface (SimFunction) to the model.
observables = {'Plasma.Drug', 'Liver.Drug', 'Spleen.Drug', 'Lung.Drug',...
    'Heart.Drug', 'Plasma.Nanoparticles', 'Liver.Nanoparticles', ...
    'Spleen.Nanoparticles'};
parameters  = get(parameterObjs,'Name');
doses = model.getdose;
dose = doses(1);
sfunction = model.createSimFunction(parameters, observables,...
    dose, 'UseParallel', true);

%% Plot and inspect the generated samples.
[~, ax] = plotmatrix(ksamples);
% Label the axes
set([ax(end,:).XLabel], {'String'}, parameters, 'Interpreter', 'None');
set([ax(:,1).YLabel],   {'String'}, parameters, 'Interpreter', 'None', 'Rotation', 0, 'HorizontalAlignment', 'right');

%% Cluster setup 
% 1. Distribute auxiliary files. These are required by the model.
pool = gcp;
% pool.addAttachedFiles({'transform.m', 'SREBP2_reg.m'});

% 2. For efficiency purposes distribute an accelerated model.
%    Comment on homogeneous clusters. 
sfunction.accelerate

% Distribute the variables sfunction and dose to all the workers. 
constant = parallel.pool.Constant({sfunction, dose});

%% Run simulations in parallel
results = cell(nRuns,1);
parfor i = 1:size(ksamples, 1)    
    sfunction   = constant.Value{1}; 
    dose        = constant.Value{2};
    outputTimes = 0:.005:40;
    % Note that need to be the same for each simulation (i.e not determined
    % by the ODE solver time steps) in order to calculate time steps.
    % Run the model. Note we specify output times. 
    simData = sfunction(ksamples(i,:), [], dose.getTable, outputTimes);
    
    results{i} = table(simData.Time, simData.Data);

end

%% calculate 95% percentile on heart drug concentration from simulations
heartD = cellfun(@(x) x.Var2(:,5),results,'UniformOutput',false);
nTimePoints = size(heartD{1},1);
% state 5 is drug concentration in the heart
% heartD is a nRuns-by-1 cell array with a nTimePoints-by-1 double array 
% of simulation data in each cell. The prctile function needs a
% nRuns-by-nTimepoints double array as input. We will use the cell2mat and
% reshape functions to achieve this
heartD = cell2mat(heartD); %convert from cell to numerical array
heartD = reshape(heartD,[nTimePoints,nRuns]);  % reshape to desired size
CI_heartD = prctile(heartD',[2.5,97.5]);
median_heartD = prctile(heartD',50);

%% plot
% the time vectors are the same for each simulations as defined when creating the simfunction
time = results{1}.Var1; 
semilogy(time, median_heartD','r')
hold on
semilogy(time,CI_heartD(1,:)','--b')
semilogy(time,CI_heartD(2,:)','--b')
axis([0 max(time) min(CI_heartD(1,:)) max(CI_heartD(2,:))])
xlabel('time [hr]')
ylabel('Heart.Drug [uM]')