Calculate NCA Parameters and Fit Model to PK/PD Data Using SimBiology Model Analyzer App

This example shows how to perform noncompartmental analysis to calculate NCA parameters and estimate the tumor growth model [1] parameters from experimental data using nonlinear regression in the SimBiology Model Analyzer app.

Tumor Growth Model

The model used in this example is a SimBiology® implementation of the pharmacokinetic/pharmacodynamic (PK/PD) model by Simeoni et al. It quantifies the effect of anticancer drugs on tumor growth kinetics from in vivo animal studies. The drug pharmacokinetics are described by a two-compartment model with IV bolus dosing and linear elimination (ke) from the Central compartment. Tumor growth is a biphasic process with an initial exponential growth followed by linear growth. The growth rate of the proliferating tumor cells is described by

L0*x1[1+(L0L1*w)ψ]1ψ

L0, L1, and Ψ are tumor growth parameters, x1 is the weight of the proliferating tumor cells, and w is the total tumor weight. In the absence of any drugs, the tumor consists of proliferating cells only, that is, w = x1. In the presence of an anticancer agent, a fraction of the proliferating cells is transformed into nonproliferating cells. The rate of this transformation is assumed to be a function of the drug concentration in the plasma and an efficacy factor k2. The nonproliferating cells x2 go through a series of transit stages (x3 and x4) and are eventually cleared from the system. The flow-through of the transit compartments is modeled as a first-order process with the rate constant k1.

The SimBiology model makes these adjustments to the pharmacodynamics of tumor growth:

  • Instead of defining the tumor weight as the sum of x1, x2, x3, and x4, the model defines the tumor weight by the reaction named Increase, null → tumor_weight, with the reaction rate (2*L1*L0*x12L1+2*L0*x1)*tumor_weight.

    tumor_weight is the total tumor weight, x1 is the weight of the proliferating tumor cells, and L0, and L1 are tumor growth parameters.

  • Similarly, the model defines the decrease in tumor weight by the reaction named Decay, tumor_weight → null, with the reaction rate k1*x4. The constant k1 is the forward rate parameter, and x4 is the last species in the series of transit reductions in tumor weight.

  • ke is a function of the clearance and the volume of the central compartment: ke = Cl_Central/Central.

PK/PD Data Description

The experimental (synthetic) data contains measurements from eight patients for three responses: measured drug concentrations in the central compartment, in the peripheral compartment, and measured tumor weight. The data also contains the dosing information, and each patient receives an IV dose at day 7.

The data set contains the following columns.

  • ID — Patient IDs

  • Time — Times when measurements are taken

  • CentralConc — Drug concentration in the central compartment

  • PeripheralConc — Drug concentration in the peripheral compartment

  • Dose — Dosing information for each patient

NaN values are used whenever there is no measurement or no dose is given.

Load Tumor Growth Model and Data

  1. Open the SimBiology Model Analyzer app by typing simBiologyModelAnalyzer at the command line or by clicking the app icon on the Apps tab.

  2. On the Home tab of the app, select Open.

  3. Navigate to the folder matlabroot\examples\simbio\data\. matlabroot is the folder where you have installed MATLAB. Select the project file named tumor_growth_fitPKPD.sbproj. In the Browser pane, the Models folder contains the Tumor Growth Model and the Data1 folder contains the experimental data.

  4. Classify the measured responses of the data as the dependent variables. In the Browser pane, click Documents and double-click Datasheet1.

  5. In the Data1 table, double-click Classification under CentralConc. Select dependent. Repeat the same process for PeripheralConc and TumorWeight. Now all the data columns have proper classifications and the data is ready for use.

Visualize Experimental Data

After you load the data, you can visualize the measured responses.

  1. In the Browser pane, click Workspace.

  2. Click Data1.

  3. On the Home tab, in the Plot section, click the time plot. The app generates a time plot of all three responses, namely: CentralConc, PeripheralConc, and TumorWeight.

    In the default time plot, Responses correspond to the measured responses and are plotted using different line styles. Scenarios refers to different groups (eight patients) in the data and are plotted using different colors.

    Tip

    Plots are backed by data that are currently present in the app workspace. Plots are not snapshots. When the data (either experimental data or simulation results) is removed or changed, the plots are also updated according to the changes in the underlying data.

Customize Data Visualization

You can customize the plot to make it clearer. For example, you can plot the PD data (TumorWeight) on a different axis than the PK data (CentralConc and PeripheralConc). To do so, create two diffferent groups (sets) of responses, where the first set contains only TumorWeight and the second set contains CentralConc and PeripheralConc.

  • Right-click TumorWeight (gram) in the Responses table (the middle table in the Property Editor) and select Create New Set. The app creates Set 1 and Set 2. Set 1 contains only TumorWeight, which is now plotted on a different axis than Set 2, which contains CentralConc and PeripheralConc.

    Note

    You can create a set only if there are two or more responses.

    The Slice table (the top table in the Property Editor) now contains Sets. This table is a summary table of all slicing variables that are currently present in the plot and their corresponding plot styles. In this current plot, the slicing variables are Sets, Responses, Scenarios, and Dose.

    Tip

    You can slice data using different slicing variables. Each slicing variable appears in the plot with a different visual style (or channel) such as color, line style, and axes position. Slicing variables can represent attributes of data, such as responses or scenarios (that is, groups or simulation runs). Slicing variables can also be covariates or parameter values associated with a scenario or group. By default, the app provides slicing variables for different response variables and different scenarios in the plotted data. You can add other visual styles (or channels) for sets of responses and associated parameter or covariate variables.

You can also group the responses based on different dose amounts that the patients receive. There are three different dose groups: 30, 75, and 150 mg.

  1. In the Slice table, at the Dose row, double-click the empty cell and select Color. A red indicator appears because another slicing variable (Scenarios) has the same plot style. Clear the style (visual channel) for Scenarios by selecting empty.

  2. In the Dose table (the bottom table), the app has automatically binned the dose amounts. Set Number of Bins to 3. You can now see that the dose amount has an impact on the tumor size. The higher the dose is, the smaller the tumor becomes.

  3. You can also query the corresponding dose group from each line by showing its data tip. Press Ctrl and click a blue line to display its data tip. To remove it, Ctrl + Click again anywhere on the same line.

Perform Noncompartmental Analysis (NCA)

Using the drug pharmacokinetic data, you can estimate NCA parameters. NCA is model agnostic and can give insights into the drug pharmacokinetics without any underlying assumptions. You can use some of the NCA results as initial estimates when calibrating the model to the data, as discussed later in this example. For details on the list of available NCA parameters and their formulas, see Noncompartmental Analysis.

NCA Program Setup

  1. On the Home tab, select Program > Non-Compartmental Analysis. A new program (Program1) appears.

  2. The Data setup step of the program defines the data set to use for the NCA analysis. In this example, the program automatically selects Data1.

  3. The NCA execution step defines the data column associations and algorithm details. In the Definition table, set Concentration to CentralConc. Leave the other settings unchanged.

  4. On the Home tab, click Run.

Display NCA Results

Once the NCA analysis is complete, the app saves the results in the LastRun folder of the program by default.

  1. On the Home tab, select New Datasheet. A new datasheet (Datasheet2) opens on a new tab.

  2. In the Browser pane, expand the Program1 folder. Then expand LastRun folder. The NCA results are stored in the table named results.

  3. Drag results onto Datasheet2. The calculated NCA parameter values are now shown in the datasheet. Use the scroll bar at the bottom of the table to scroll and view the NCA data columns. For details about calculated NCA parameters, see Noncompartmental Analysis.

Export Results to MATLAB Workspace

You can export the NCA results to the MATLAB® workspace and perform further data analyses at the command line.

  1. Right-click results. Select Export Data to MATLAB Workspace.

  2. The SimBiology Data Export dialog box opens. Change the name of the variable to ncatable. Click OK.

After you export the data to the MATLAB workspace, you can analyze the data at the command-line. For instance, you can calculate the average drug clearance from the NCA data and use it as the model parameter value.

Estimate Model Parameters Using Nonlinear Regression

SimBiology provides different regression techniques to estimate model parameters based on experimental data. This example details the steps for using the nonlinear regression method lsqnonlin (requires Optimization Toolbox™) to fit the model to the data. If you do not have Optimization Toolbox, the app uses fminsearch instead. For the purposes of the example, only some of the PK/PD model parameters are estimated, namely: k1, L0, L1, Cl_Central, k12, and k21.

Fit Program Setup

  1. From the Home tab, select Program > Fit Data. A new program (Program2) appears on a new tab. The Data and Model steps have been prepopulated with Data1 and Tumor Growth Model, respectively.

  2. By default, the Fit step autogenerates plots after the fitting is complete. Disable the plot generation by clicking the plot icon at the top of the Fit program step for now. The plots will be explored later in the example.

  3. Define the mapping between the model components and the data columns from the experimental data. Specifically, map CentralConc to the model species Drug in the central compartment (Central.Drug). Similarly, map PeripheralConc to Peripheral.Drug.

  4. Map TumorWeight to [Tumor Growth Model].tumor_weight.

  5. Map the Dose column to Central.Drug to indicate that the Drug species in the Central compartment is being dosed.

  6. Define the model parameters to estimate in the Estimated Parameters table. Double-click the empty cell in the Estimated Parameters column and type k1. The app shows model components with matching names. Select k1 from the list.

    By default, the parameter is log-transformed as indicated by the expression exp(k1). You can change the expression to no transformation, probit, or logit transformation. For this example, keep the default log transform because it often improves the convergence. The Initial Untransformed Value is automatically set to the model value which is 0.5.

  7. Enforce the biological parameters to stay positive by specifying the Untransformed Lower Bound and Untransformed Upper Bound as 1e-5 and 10, respectively.

  8. Similarly, add the following parameters: Cl_Central, L0, L1, k12, and k21.

  9. Select Pooled fit to estimate one set of parameters for all patients (population fit). If you do not select Pooled fit, the app estimates one set of parameters for each patient (individual fit).

  10. The default error model is the constant error model. SimBiology supports constant, proportional, exponential, and combined error models. For details, see Error Models. For now, use the constant error model.

  11. Keep the rest of the fitting settings unchanged. These settings are

    • Estimation Method — The default method is lsqnonlin if you have Optimization Toolbox. If you do not, the app uses fminsearch.

      For more information, see Supported Methods for Parameter Estimation in SimBiology.

    • Algorithm Settings — Most common options for the estimation method. Click to expand the section and see the options. To see the description of each option, click the info icon to the right of the header.

    • Advanced Algorithm Settings — Advanced settings for the estimation method. The table is empty by default.

Run Fit Program

After you set the fitting options, you can run the Fit step.

  1. At the top of the Fit step, click the Run this program step button.

    By default, the Fit step shows the progress of parameter estimation in a separate figure. The progress plot shows the live status of parameter estimation and fitting quality measures such as log likelihood. For details, see Progress Plot.

  2. The progress plot shows that the fit converged. You can close the progress plot.

    If you are using fminsearch, the fit can fail to converge due to reaching the maximum number of iterations. You can increase MaxIter in the Algorithm Settings, but for the purposes of this example, you can continue completing the steps without doing so.

Visualize Fit Results

Once the parameter estimation is complete, the fit results are stored in the LastRun folder of the program.

  1. On the Home tab, select New Datasheet. A new datasheet opens in a new tab.

  2. In the Browser pane, expand the Program2 folder. Then expand LastRun folder, which contains results and simdataI. results contains estimated parameter values and fit statistics. simdataI contains the simulated model responses for each individual (patient or group) using the estimated parameter values.

  3. Drag results onto the new datasheet.

    The Statistics table contains fit quality measures such as AIC, BIC, and log likelihood. These measures can be useful, for example, to compare the performance of different error models.

  4. In addition to quality statistics, you can also view various fit plots, such as actual versus predicted plots and residual distribution plots. First, select results in the Browser pane. Available fit plots are automatically listed in the Plot section on the Home tab. Then select Act vs Pred from the list.

    The actual versus predicted plot appears on a separate tab. The predicted responses are plotted on the x-axis and the observed (experimental) responses are plotted on the y-axis.

    You can change the plot to other supported plots by selecting one of the plots from the Style section in the Property Editor. If you want the new plot on its own separate tab and you do not want to reuse the existing plot tab, select the plot from the Plot section on the Home tab.

  5. Change the plot to a residual distribution plot by selecting Res Dist in the Style section.

    The plot shows whether the residuals for each response are normally distributed. In an ideal normal probability plot of residuals, the residuals line up along the diagonal line across the plot and the histograms indicate a normal fit. However, from the plot, the residuals for all three responses, especially CentralConc and PeripheralConc, do not seem to be normally distributed. It could be an indication that the constant error model assumption is incorrect.

Compare Different Error Models

The following steps show how to change the error model to an exponential error model to fit the data again and compare the fit statistics of two different error models.

Save Fit Results.  Before fitting the data again using the exponential error model, save the constant error model result in a separate folder. Otherwise, the program, by default, overwrites results from the LastRun folder every time you run the fit.

  1. Right-click the LastRun folder of the fit program in the Browser pane.

  2. Select Save Data.

  3. In the Save Data dialog, enter fit_constant as the data name.

Rerun Fit With Exponential Error Model.  After you save the data, you can rerun the fit program with a different error model.

  1. Return to the fit program by clicking the Program2 tab. In the Error Model section, select exponential.

  2. At the top of the Fit step, click the Run this program step button.

  3. Close the progress plot after the fit completes.

  4. If you closed the datasheet (Datasheet3) that contains the fit statistics from the previous fit, reopen the datasheet. In the Browser pane, click Documents. Then double-click Datasheet3.

  5. Click Workspace to go back to the fit program results. Expand LastRun folder of the fit program. Then drag results onto Datasheet3. New columns (Program2_LastRun) containing the latest fit results are added next to the previous fit results (Program2_fit_constant).

    The Statistics table compares the fit quality measures. From the comparison, both the AIC and BIC of the fit using the exponential error model are smaller than those of the previous fit. This indicates that the exponential error model fits the data better than the constant error model. The larger log likelihood of the exponential error model also indicates it is a better fit.

  6. Next, look at the residual distribution plot. Click results from the LastRun folder. Then click Residual Dist from the Plot section on the Home tab.

    Compared to the residual distributions of the constant error model, the residual distributions from the exponential error model look more normal, indicating that the exponential error model indeed fits the data better.

Calculate Confidence Intervals

Another way to assess the quality of fit results is to compute 95% confidence intervals for the estimated parameters and model predictions — that is, model simulation results using the estimated parameters. This step requires Statistics and Machine Learning Toolbox™.

  1. Return to the fit program. Click the (+) icon at the upper left and select Confidence Interval. A Confidence Interval step appears following the Fit step.

  2. At the top of the Confidence Interval step, disable the autogeneration of plots by clicking the plot icon. For both Parameter Confidence Intervals and Prediction Confidence Intervals, use the default method gaussian and 95% confidence level. Click the Run this program step button to compute confidence intervals.

    For parameter confidence intervals, the supported methods are gaussian, profileLikelihood, and bootstrap.

    For prediction confidence intervals, the supported methods are gaussian and bootstrap.

  3. Once completed, the results are stored as parameterCI and predictionCI in the LastRun folder of the program. parameterCI contains 95% confidence intervals for the estimated parameters. predictionCI contains 95% confidence intervals for the model predictions.

  4. Plot 95% confidence intervals for the estimated parameters. Click parameterCI in the Browser pane and select Confidence in the Plot section.

    The confidence interval for each estimated parameter is shown in a new plot. The plot indicates the successful computation of the confidence intervals for all estimated parameters.

    Depending on the outcome (status) of the confidence interval estimation, the app plots the results differently.

    • If the status of confidence interval estimation is a success (as in the above plot), the app uses the first default color (blue) to plot a line and a centered dot for every parameter estimate. The app also plots a box to indicate the confidence intervals.

    • If the status is constrained or estimable, the app uses the second default color (red) and plots a line, centered dot, and box to indicate the confidence intervals.

    • If the status is not estimable, the app plots only a line and a centered cross in red.

    • If there are any transformed parameters with estimated values that are 0 (for the log transform) and 0 or 1 (for the probit or logit transform), no confidence intervals are plotted for those parameter estimates.

    For more details on the definitions of different statuses, see Parameter Confidence Interval Estimation Status.

    You can also change the Layout of the plot in Plot Settings.

    • The 'split' layout displays the confidence interval for each parameter estimate on a separate axis.

    • The 'grouped' layout displays all confidence intervals on one axis, grouped by parameter estimates. Each estimated parameter is separated by a vertical black line.

    In both cases, the parameter bounds defined in the original fit are marked by square brackets. The app uses vertical dotted lines to group confidence intervals of parameter estimates that have been computed in a common fit.

  5. Similarly, plot 95% confidence intervals for the model predictions. Click predictionCI in the Browser pane and select Confidence in the Plot section.

    In this plot, the confidence interval for each group is plotted in a separate column, and each response is plotted in a separate row. The plot indicates the successful computation of the confidence intervals.

    The plotting behavior differs depending on the outcome (status) of the confidence interval calculation.

    • If the status is constrained or not estimable, the app uses the second default color (red) to plot the confidence intervals.

    • Otherwise, the app uses the first default color (blue) and plots the confidence intervals as shaded areas (as in the above plot).

    For details, see Gaussian Confidence Interval Calculation for Model Predictions and Bootstrap Confidence Interval Calculation.

Show Overlay of Fitted Results Over Experimental Data

Once the parameter estimation is complete, you can plot the simulation results on top of the experimental data. There are generally two ways to get the plot.

The first approach is to use the default Fit plot that shows the fitted results over experimental data for each group.

  1. Click results from the LastRun folder of the fit program.

  2. Select Fit in the Plot section on the Home tab.

The second approach is more customized and involved, but you have more control over how the plot looks like. The following steps show how to plot the tumor growth profiles of all groups on the same axes but use different line styles for experimental data and simulation data.

  1. Click simdataI in the LastRun folder of the fit program.

  2. Click time from the Plot section on the Home tab.

  3. Drag Data1 (experimental data) on to the plot.

  4. In the Property Editor of the plot, multiselect (Ctrl + Click) the first three rows (simulation data) of the Responses table and click Create Set.

  5. In the Slice table, use no style for Scenarios by selecting empty.

  6. Clear the first two rows from each set. Keep only the tumor growth data from each set.

  7. Use Line Style for Sets and Color for Response. Also set the Marker of the simulation data (Program2) to none. Set the Style of the experimental data (Data1) to none. Select o as the Marker.

After fitting to data, you can set the model values to the parameter estimates and perform other analyses. For instance, you can find out important model parameters using sensitivity analysis and vary the sensitive parameters to investigate model variability by using virtual patients.

References

[1] Simeoni, M., P. Magni, C. Cammia, G. De Nicolao, V. Croci, E. Pesenti, M. Germani, I. Poggesi, and M. Rocchetti. 2004. Predictive pharmacokinetic-pharmacodynamic modeling of tumor growth kinetics in xenograft models after administration of anticancer agents. Cancer Research. 64:1094-1101.

See Also

Related Topics

SimBiology Documentation

Support