- I've added a term max(0, ...) to the Hill equation to ensure that the value is never negative or NaN.
- 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.
- 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.
- 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.
- I have tightened the relative tolerance, which is generally good practice when using a gradient-based optimization method to estimate parameter values.
- I have updated responseMap to indicate that you are estimating the transformed values (represented as species y).
- I have updated estimatedParams to indicate that you are estimating ke, Rmax, gamma, and A.
Estimating PKPD model in simbiology
3 ビュー (過去 30 日間)
古いコメントを表示
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
0 件のコメント
採用された回答
Arthur Goldsipe
2015 年 5 月 11 日
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:
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
0 件のコメント
その他の回答 (2 件)
Arthur Goldsipe
2015 年 5 月 11 日
Hi Brandon,
There are a number of examples available on the MathWorks site.
-Arthur
1 件のコメント
Arthur Goldsipe
2015 年 5 月 11 日
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.
Brandon
2015 年 5 月 12 日
3 件のコメント
Arthur Goldsipe
2015 年 5 月 15 日
I wasn't able to reproduce these warnings. These warnings indicate that some simulations were not able to complete (or the simulation returned invalid/NaN values) with particular parameter values during estimation. This can happen if a concentration goes slightly negative during a simulation and you try to raise that negative number to a fractional power. You might try updating the rule to something like this:
y = max(0, Rmax*max(0,Drug_Central)^gamma / ( max(0,Drug_Central)^gamma+max(0,A)^gamma))', 'repeatedassignment');
Teerachat Saeheng
2022 年 4 月 15 日
編集済み: Teerachat Saeheng
2022 年 4 月 15 日
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.
コミュニティ
その他の回答 SimBiology コミュニティ
参考
カテゴリ
Help Center および File Exchange で Scan Parameter Ranges についてさらに検索
製品
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!