Estimating PKPD model in simbiology

9 views (last 30 days)
Does anyone have a simple working example of code that uses simbiology to estimate parameters for a simple nonlinear pharamacodynamic model? I would like to first simulate data using a simple PKPD model, then estimate model parameters from the data. The code below creates the PK model. However, I have had no success trying to create the PD part. In the code below, the estimation is based on observing the drug concentrations over time. What I'd like to do instead is to pass the concentration values through a Hill function, i.e. if x is the concentration, then I want the observations to be
y = R*x^g/(x^g+c^g)
where R, g, and c are constants. I would like to allow the estimation routine to "see" only these nonlinearly transformed y values, and to estimate both the PK values (as is done in the code below), and also the parameters of the Hill equation, namely R, g, and c.
Can anyone show me how to do this?
Thanks in advance!
Brandon
------------------------------------
%*********************************
% *** simulate the data **********
%*********************************
m1=sbiomodel('onecomp')
r1=addreaction(m1,'Drug_Central -> null') % elimination
k1 = addkineticlaw(r1,'MassAction')
k1val = lognrnd(0, 0.2 , 1,1);
p1 = addparameter(k1,'ke','Value',k1val,'ValueUnits','1/hour')
k1.ParameterVariableNames='ke'
% cannot seem to get the following part to work
% r2=addreaction(m1,'Drug_Central -> x') % nonlinear observation
% k2 = addkineticlaw(r2, 'Hill-Kinetics');
% k2.KineticLaw = 'Unknown';
% r2.ReactionRate = 'Rmax*x^gamma / ( x^gamma+A^gamma)';
% p2 = addparameter(k2, 'Rmax', 2.3);
% p3 = addparameter(k2,'A',5);
% p4 = addparameter(k2,'gamma',3);
% set(p4, 'ValueUnits', 'dimensionless');
%%info for each constant rate interval
% start time, end time, rate [figure out intermediate: amount]
RateInfo=[2 4 10; 6 12 50; 12 20 90; 22 24 150; 25 27 200; 30 31 50];
d=[];
for i=1:size(RateInfo,1);
dt = sbiodose('dt');
dt.TargetName = 'Drug_Central';
dt.RateUnits = 'milligram/hour';
dt.AmountUnits='milligram';
dt.TimeUnits = 'hour';
dt.Rate = RateInfo(i,3);
t0=RateInfo(i,1);
t1=RateInfo(i,2);
amount=(t1-t0)*RateInfo(i,3);
dt.Amount = amount;
dt.StartTime=t0
dt.Active = true;
d=[d dt];
end
dose=d;
%%run simulation
cs = getconfigset(m1)
cs.StopTime=48
cs.TimeUnits='hour'
[t,sd,species]=sbiosimulate(m1,d)
plot(t,sd);
legend(species);
xlabel('Hours');
ylabel('Drug Concentration');
% throw out some data to simulate sparse sampling in an experiment
t=t(1:5:end);
sd=sd(1:5:end);
data=table(t,sd); % convert to table
data.Properties.VariableNames{'t'}='Time';
data.Properties.VariableNames{'sd'}='Conc';
%%convert to groupedData object -- required for fitting with sbiofit
gData=groupedData(data)
gData.Properties.VariableUnits={'hour','milligram/liter'}
gData.Properties
%*********************************
% *** estimate model from data ***
%*********************************
pkmd = PKModelDesign
pkc1 = addCompartment(pkmd,'Central')
pkc1.DosingType = 'Infusion'
pkc1.EliminationType='linear-clearance'
pkc1.HasResponseVariable=true
[model,map]=construct(pkmd)
configset = getconfigset(model)
configset.CompileOptions.UnitConversion=true
configset.SolverOptions.AbsoluteTolerance=1e-9
configset.SolverOptions.RelativeTolerance=1e-5
dose=d;
% look at map to see variables
responseMap = {'Drug_Central = Conc'};
paramsToEstimate = {'log(Central)','log(Cl_Central)'}
estimatedParams = estimatedInfo(paramsToEstimate,'InitialValue',[1 1])
paramsToEstimate = {'log(Central)','log(Cl_Central)'};
estimatedParams = estimatedInfo(paramsToEstimate,'InitialValue',[1 1]);
%%estimate the parameters
fitConst = sbiofit(model,gData,responseMap,estimatedParams,dose)
%%plot results
fitConst.ParameterEstimates
figure(1);
plot(fitConst);
%%compare estimate to truth
disp('****************')
fitConst.ParameterEstimates.Estimate(2)
k1val

Accepted Answer

Arthur Goldsipe
Arthur Goldsipe on 11 May 2015
I reread your question, and I see that the above answer didn't really address your question about transforming concentrations through a Hill equation. You can do that using a repeated assignment rule. I'll paste below my suggested updates to your code. First though some notes:
  1. I've added a term max(0, ...) to the Hill equation to ensure that the value is never negative or NaN.
  2. By default, all species are returned from the simulation. Since I added species y to the model, the initial simulation returns simulated results for both Drug_Central and y.
  3. The fit looks quite good, but the estimated values are relatively far from the actual values. This indicates that the predictions are not sensitive enough to the parameter values to estimate them very accurately.
  4. You are adding units to your model, but they won't really be used unless you enable unit conversion. If you do this, you will need to add units to everything in the model.
  5. I have tightened the relative tolerance, which is generally good practice when using a gradient-based optimization method to estimate parameter values.
  6. I have updated responseMap to indicate that you are estimating the transformed values (represented as species y).
  7. I have updated estimatedParams to indicate that you are estimating ke, Rmax, gamma, and A.
Ok, here's the updated code:
%*********************************
% *** simulate the data **********
%*********************************
m1=sbiomodel('onecomp');
cs = m1.getconfigset;
cs.SolverOptions.RelativeTolerance = 1e-8;
r1=addreaction(m1,'Drug_Central -> null'); % elimination
k1 = addkineticlaw(r1,'MassAction');
k1val = lognrnd(0, 0.2 , 1,1);
p1 = addparameter(m1,'ke','Value',k1val,'ValueUnits','1/hour');
k1.ParameterVariableNames='ke';
% Use Hill equation to transform concentration, but
% ensure result has a minimum of 0.
m1.addspecies('y');
m1.addrule('y = max(0, Rmax*Drug_Central^gamma / ( Drug_Central^gamma+A^gamma))', 'repeatedassignment');
p2 = addparameter(m1, 'Rmax', 2.3);
p3 = addparameter(m1,'A',5);
p4 = addparameter(m1,'gamma',3);
%%info for each constant rate interval
% start time, end time, rate [figure out intermediate: amount]
RateInfo=[2 4 10; 6 12 50; 12 20 90; 22 24 150; 25 27 200; 30 31 50];
d=[];
for i=1:size(RateInfo,1);
dt = sbiodose('dt');
dt.TargetName = 'Drug_Central';
dt.RateUnits = 'milligram/hour';
dt.AmountUnits='milligram';
dt.TimeUnits = 'hour';
dt.Rate = RateInfo(i,3);
t0=RateInfo(i,1);
t1=RateInfo(i,2);
amount=(t1-t0)*RateInfo(i,3);
dt.Amount = amount;
dt.StartTime=t0
dt.Active = true;
d=[d dt];
end
dose=d;
%%run simulation
cs = getconfigset(m1)
cs.StopTime=48
cs.TimeUnits='hour'
[t,sd,species]=sbiosimulate(m1,d)
plot(t,sd);
legend(species);
xlabel('Hours');
ylabel('Drug Concentration');
% throw out some data to simulate sparse sampling in an experiment
t=t(3:5:end);
sd=sd(3:5:end,2); % y is the 2nd species
data=table(t,sd); % convert to table
data.Properties.VariableNames = {'Time', 'Conc'};
%%convert to groupedData object -- required for fitting with sbiofit
gData=groupedData(data)
gData.Properties.VariableUnits={'hour','milligram/liter'}
gData.Properties
%*********************************
% *** estimate model from data ***
%*********************************
% look at map to see variables
responseMap = {'y = Conc'};
paramsToEstimate = {'log(ke)', 'log(Rmax)','gamma', 'log(A)'}
estimatedParams = estimatedInfo(paramsToEstimate,'InitialValue',[1 1 1 1])
%%estimate the parameters
fitConst = sbiofit(m1,gData,responseMap,estimatedParams,dose)
%%plot results
fitConst.ParameterEstimates
figure(1);
plot(fitConst);
%%compare estimate to truth
disp('****************')
fitConst.ParameterEstimates.Estimate(1)
k1val

More Answers (2)

Arthur Goldsipe
Arthur Goldsipe on 11 May 2015
Hi Brandon,
There are a number of examples available on the MathWorks site.
Here are SimBiology PK/PD example files from the File Exchange.
-Arthur
  1 Comment
Arthur Goldsipe
Arthur Goldsipe on 11 May 2015
I see that this doesn't really answer your specific question about how to use a Hill equation to transform your response. I'm working on a separate answer for that.

Sign in to comment.


Brandon
Brandon on 12 May 2015
Arthur, Perfect. Works great! Thanks, Brandon
  3 Comments
Teerachat Saeheng
Teerachat Saeheng on 15 Apr 2022
Hi How can I add the function of load drug data concentrations and pharmacodynamic response from excel file to this command. In case I would like to use Emax model for PD response what should I do? . Thank you for your suggestion Sorry that I never use line command I always use simbiology for creating new model.

Sign in to comment.

Communities

More Answers in the  SimBiology Community

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!