Estimating PK/PD model in simbiology -- does anyone have an example?

Does anyone know how to simulate and estimate both the kinetics parameters and sigmoid nonlinearity parameters of a PKPD model in simbiology?

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.

Does anyone know how to do this?

Thanks in advance for any help!

---------------------

clear all; clc; format compact;

%*********************************
% *** simulate the data **********
%*********************************

%%simulate a simpmle 1 compartment PKPD model,
% then estimate model parameters from simulated data
% compare estimated params to true params

%%add constant infusion rate to one-compartment model
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');

RateInfo=[2 4 10; 6 12 50; 12 20 90; 22 24 150; 25 27 200; 30 31 50];
% d=fcnGetDoseArray(RateInfo);
% function d=fcnGetDoseArray(RateInfo);

%%info for each constant rate interval
% start time, end time, rate [figure out intermediate: amount]
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');

%%convert to table

% 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);
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]);

%%need to go back and figure out how to estimated rate parameter instead of clearane and volume
%%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

[Arthur Goldsipe]

This seems to be a near duplicate of another question you posted a few days ago. Please consider deleting this question, as I've already posted some answers on the other question.

[Walter Roberson]

Duplicates http://www.mathworks.com/matlabcentral/answers/215972-estimating-pkpd-model-in-simbiology