MATLAB: Estimating PKPD model in simbiology

Question

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

Answer 1

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