Estimating PKPD model in simbiology

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

[Arthur Goldsipe]

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

Answer 2

[Arthur Goldsipe]

Hi Brandon,

There are a number of examples available on the MathWorks site.

Here are SimBiology PK/PD example files from the File Exchange.

Here are SimBiology PK/PD videos.

-Arthur

[Arthur Goldsipe]

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.

Answer 3

[Brandon]

Arthur, Perfect. Works great! Thanks, Brandon

[Brandon]

I spoke too soon. Unfortunately the code above does not produce good results every time. When I run it repeatedly, about half of the time I get one of the following two errors:

(1)

Warning: Simulation with phi = [1.61334524819858 2.00040494948441
1.44714560655445 2.45301716162508]  could not be completed because:
Domain error. To compute complex results, make at least one input complex,
e.g. 'power(complex(a),b)'.. Empty results have been returned for this
simulation. 
> In SimFunction>updateOutputForSimulationErrors at 1927
  In SimFunction>SimFunction.localVectorizeRegressionFcn/callFunction at 719
  In SimFunction>SimFunction.vectorsimParallel at 527
  In LeastSquaresResults>LeastSquaresResults.fitted at 280
  In LeastSquaresResults>LeastSquaresResults.constructTaskResult at 460
  In LeastSquaresResults>LeastSquaresResults.plot at 340
  In test at 69

Or (2)

Warning: Rank deficient, rank = 0, tol = NaN. 
> In nlinfit>LMfit at 574
  In nlinfit at 276
  In sbiofit>runNlinfit at 432
  In sbiofit>fitSlice at 639
  In sbiofit at 587
  In a_NLME_FromMatlabCentralAnswer at 72 
Warning: When 'FunValCheck' is 'off', NLINFIT attempts to continue iterating
even if the function you provided as the MODELFUN input returns Inf or NaN
values. However, the Jacobian evaluated at the final estimate of BETA has
Inf or NaN values or is not of the required size. This may be because the
supplied MODELFUN has returned Inf or NaN values. Try refitting with
'FunValCheck' set to 'on' to detect if your MODELFUN returns Inf or NaN
values. You may have to modify your MODELFUN so that it always returns
non-NaN and finite values. 
> In nlinfit at 348
  In sbiofit>runNlinfit at 432
  In sbiofit>fitSlice at 639
  In sbiofit at 587
  In a_NLME_FromMatlabCentralAnswer at 72

Any idea what is going on here, and how to prevent it? This should be a fairly standard PKPD estimation problem.

[Arthur Goldsipe]

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]

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.