## Model Biological Variability with Virtual Patients Using SimBiology Model Analyzer App

This example shows how to generate and simulate virtual patients using the **SimBiology
Model Analyzer** app. In this example, a virtual patient is represented as a single
realization of model parameters. The example uses a tumor growth model [1] to explore the variability of
some model parameters that influence the tumor growth and investigate various dosing regimens to
control it.

This example requires Statistics and Machine Learning Toolbox™.

### 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

$$\frac{{L}_{0}*{x}_{1}}{{\left[1+{\left(\frac{{L}_{0}}{{L}_{1}}*w\right)}^{\psi}\right]}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\psi $}\right.}}$$

*L _{0}*,

*L*, and Ψ are tumor growth parameters,

_{1}*x*is the weight of the proliferating tumor cells, and

_{1}*w*is the total tumor weight. In the absence of any drugs, the tumor consists of proliferating cells only, that is,

*w* =
*x*_{1}

. 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 *k*. The nonproliferating cells

_{2}*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

*k*.

_{1}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 $$\left(\raisebox{1ex}{$2*{L}_{1}*{L}_{0}*{x}_{1}{}^{2}$}\!\left/ \!\raisebox{-1ex}{${L}_{1}+2*{L}_{0}*{x}_{1}$}\right.\right)*tumor\_weight$$.*tumor_weight*is the total tumor weight,*x*is the weight of the proliferating tumor cells, and_{1}*L*, and_{0}*L*are tumor growth parameters._{1}Similarly, the model defines the decrease in tumor weight by the reaction named

*Decay*,`tumor_weight → null`

, with the reaction rate

. The constant*k**_{1}*x*_{4}*k*is the forward rate parameter, and_{1}*x*is the last species in the series of transit reductions in tumor weight._{4}*ke*is a function of the clearance and the volume of the central compartment:`ke = Cl_Central/Central`

.

### Description of Virtual Population

The virtual population in this example is represented by virtual patients, specified as
distinct sets of patient-specific parameter values. Suppose that, based on prior knowledge or
sensitivity analysis of
the model, the tumor growth is sensitive to these model parameters: *L0*,
*L1*, *w0*, *k1*, *k2*, and
*Cl_Central*. Assuming that these biological parameters follow the lognormal
distribution (and must always be positive), you can generate virtual patients that represent
different parameter values drawn from the joint lognormal distribution.

The lognormal distribution uses these parameters:

*mu*– Mean of logarithmic values*sigma*– Standard deviation of logarithmic values

For details, see Lognormal Distribution (Statistics and Machine Learning Toolbox).

In this example, *sigma* is assumed to be 0.01 for all the parameters. For
a small *sigma* value, the mean of a lognormal distribution is approximately
equal to the log of the model value. Hence, this example assumes that *mu* for
each parameter is the log of the model value.

### Dosing Strategies

In this example, an anticancer drug is used to control the tumor growth. Each virtual patient receives the same amount of the drug on the same schedule. The tumor growth response is simulated for each patient. Dose amounts are then varied to find the range of dose amounts that can suppress the tumor growth of many virtual patients in the population.

### Generate Samples to Represent Virtual Patients

The following steps show how to draw sample values from the joint probability distribution for model parameters that are sensitive to tumor growth. These sample values represent virtual patients.

#### Load Tumor Growth Model

At the MATLAB

^{®}command line, load the model (*m1*) by entering:`sbioloadproject tumor_growth_vpop_sa.sbproj`

Open the model in the

**SimBiology Model Analyzer**app by entering:simBiologyModelAnalyzer(m1)

Alternatively, you can open the app by clicking

**SimBiology Model Analyzer**on the**Apps**tab. Then load the model from the app by selecting**Model**>`Import Model from MATLAB`

.

#### Define Joint Probability Distribution

Create a program to draw sample values from the joint lognormal distribution for tumor
growth sensitive parameters: *L0*, *L1*,
*w0*, *k1*, *k2*, and
*Cl_Central*.

On the

**Home**tab, select**Program**>**Generate Samples**. A new program opens.In the

**Variants**section in the**Model**step of the program, click the option to view the variants to be applied. Then select**parameterEstimates**, which contains the estimated parameter values.Keep the

**Doses**section as is so that no doses are selected.In the

**Generate Samples**step, click the plot button to disable default plot generation. You will plot the samples later on.In the

**Parameter Set**section, set**Type**to`values from a distribution`

.Set

**Number Of Samples**to`25`

.Double-click the empty cell in the

**Component Name**column, and enter**L0**.Change

**Distribution**for**L0**to`Lognormal`

. By definition,*mu*is the mean of logarithmic values. So, change*mu*to -0.5644, which is`log(Current Value)`

or`log(0.5687)`

. Change*sigma*to 0.01. For a small*sigma*value, the mean of a lognormal distribution is approximately equal to`log(Current Value)`

. For details, see Lognormal Distribution (Statistics and Machine Learning Toolbox).Double-click the empty cell in the

**Component Name**column, and add*L1*. Repeat the same process to change the distribution to lognormal and set the*mu*and*sigma*values. Similarly, add*w0*,*k1*,*k2*, and*ke*. This table lists the corresponding*mu*values of these parameters that you can copy and paste in the software. Change*sigma*to 0.01 for each parameter.Parameter

Current Value *mu*or Log(Current Value)L0

0.5687

-0.5644

L1

0.2690

-1.3130

w0

1

0

k1

0.4643

-0.7672

k2

8.4200e-4

-7.0797

Cl_Central

0.6672

-0.4047

For

**Sampling**, use the default option:`random sampling with rank correlation matrix`

, where the matrix is an identity matrix.The following figure shows the

**Parameter Set**section with all parameters configured.(Optional) Save the project under a new name by selecting

**Save****>****Save Project As**on the**Home**tab.

#### Define Dosing Strategies

Add the dosing information as another parameter set by specifying different dose amount—specifically, six dose amounts that are equally spaced from 5 to 35 mg. Then combine the doses with the first parameter set (virtual patients) using the Cartesian combination method. This method combines every virtual patient with every dose amount to generate a total of 150 iterations (or simulation scenarios). For details, see Combine Simulation Scenarios in SimBiology.

At the bottom of the

**Generate Samples**step, click**Add parameter set to scan**. Another parameter set (**PS2**) appears.Double-click the empty cell in the

**Component Name**column, and type`interval`

. A list of choices appears. Select`interval_dose.Amount`

.Set

**Type**to`Individual Values`

and set**Values**to`[5:6:35]`

.Near the top of the

**Generate Samples**step, under**Parameter Set Combinations**, ensure that the**Combination**type is set to`cartesian`

to combine`PS1`

and`PS2`

.

#### Generate and Visualize Parameter and Dose Combinations

Once you have defined the joint lognormal distribution for model parameters, range of dose amounts, and the combination method, you can run the step to generate the different parameters and dose combinations.

Click the

**Run this program step**button to generate samples.**Tip**You can run every program

*execution*step separately. An execution step includes the run button next to the name the step. Running an individual step is especially helpful if the program contains multiple steps and you want to see the intermediate results from a particular step. By doing so, you can make adjustments as needed before running the next step or the whole program. To run the whole program, click the**Run**button on the**Home**tab.By default, the app stores the generated

**samples**in the**LastRun**folder of the program.Visualize the generated samples. In the

**Browser**pane on the left, expand the**Program1**folder. Expand the**LastRun**folder, and then click**samples**. In the**Plot**section on the**Home**tab, click**Plot Matrix**. The plot shows the distribution of each parameter varied around its model value, except the dose amounts, which appear in a uniform range. Note that your plot might vary from the plot shown here due to randomness.**Note**You can enable default plot generation to display the plot automatically every time you run the step. To do so, click the plot button at the top of the

**Generate Samples**step. Keep in mind that when the program step has a lot of samples, plotting all the samples can be time consuming.Display the generated virtual patients and dose combinations numerically in a tabular format. On the

**Home**tab, click**New Datasheet**. From the**LastRun**folder of the program, drag**samples**into the new datasheet. Each row of the table represents a simulation scenario with different model parameter values and dose amounts.

### Perform Monte Carlo Simulations

Once you have the samples ready, you can simulate the model to explore the tumor growth of virtual patients receiving different dose amounts.

Go back to the program by clicking the

**Program1**tab.Click the (+) icon at the upper left and click

**Simulation**.In the

**Model**setup step, in the**Variants**section, make sure that**parameterEstimates**is selected. Also make sure that no doses are selected in the**Doses**section. In the**States To Log**section, click the option for viewing the states, and then clear all check boxes except**[Tumor Growth Model].tumor_weight**.At the top of the

**Simulation**step (scroll down to see this step), click the plot button to disable default plot generation.Click the

**Run this program step**button to simulate the model.The app simulates the generated scenarios from the

**Generate Samples**step that you ran previously. You do not need to rerun the step.Once the simulation finishes, the app saves the simulation

**results**in the**LastRun**folder.In the

**LastRun**folder in the**Browser**pane, click**results**. In the**Plot**section on the**Home**tab, click**time**. Each line in the plot represents the tumor weight profile of each simulation scenario.The plot shows that the tumor weight profiles appear in groups, corresponding to different dose amounts. To better visualize this observation, you can slice the data by dose amounts. To do so, you can use the

**Slice Data**table, which contains a summary of slicing variables that are currently being used in the plot and their corresponding plot styles.**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.

On the

**Slice Data**tab of**Property Editor**, in the**Visual Channels**table, double-click`L0`

and select`interval_dose Amount`

.**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.Double-click the empty cell in the

**Style**column and select`Grid`

.The app updates

**Plot2**to show each dose amount on its own axes.The plot shows that different dose amounts play critical roles in the tumor growth of virtual patients. From this plot, you can obtain some initial insights into the optimal dose amounts and dose scheduling. For instance, suppose that the target tumor weight to reach is 0.5 grams. The simulation results indicates that a dose amount of 23 milligrams or larger can achieve that goal (you can fine-tune the range of dose amounts further by tweaking it in the

**Generate Samples**step of the program). You can use this information in combination with existing drug toxicity information (not discussed in this example) to get a dosing regimen that satisfies both efficacy and safety requirements, for instance.

### Postprocess Simulation Data

You can also postprocess the simulation results to look at the correlation between dose amounts and the tumor weight in another way.

Go back to the program by clicking the

**Program1**tab. Click the (+) icon near the top of the program and select**Calculate Observables**under**Postprocessing**.A

**Postprocessing: Calculate Observables**step appears below the**Simulation**step.Double-click the first cell in the

**Name**column. Enter`stat1`

.In the

**Expression**column, enter the following expression:`min(vertcat(nan,tumor_weight(time>=7)))`

. The expression returns the minimum tumor weight after the first dose applied at day 7 from each simulation.**Note**Anytime you add an expression to the

**Observables**table in the step, the expression is automatically added as an`observable`

object to the corresponding model.In the

**Units**column, enter`gram`

.Rerun the simulation step by clicking the run button at the top of the step. The app simulates and evaluates the

`stat1`

expression for each simulation scenario or iteration, and generates the following plot. The*x*-axis represents the parameters and the*y*-axis represents the minimum tumor weight. Each dot represents a simulation scenario. The plot also shows that there is no correlation between the tumor size and the values of various model parameters, except for the dose amounts (the rightmost subplot).Show only the dose amount subplot. In the

**Plot Settings**, clear all check boxes in the**Parameters**table, except`interval_dose Amount`

.In the

**Grid**section, select**both**to show the grids for both the*x*- and*y*-axis.The plot confirms that higher dose amounts control the tumor growth better. These initial simulation results indicate that the dose amounts of 23 mg or larger could reach the hypothetical tumor weight threshold of 0.5 grams or lower. You can further adjust the range of dose amounts in the

**Generate Samples**step.**Tip**To interactively explore the plotted data, export the plot to a separate figure window by selecting

**Export Plot**from the context (right-click) menu of the plot.

This example shows you how to generate samples to represent virtual patients and perform Monte Carlo simulations to explore the model response on tumor growth under different dose amounts. The simulation results indicate a range for dose amounts and dose schedules that controls the tumor growth for various virtual patients. You could further adjust the dosing regimens so that the doses stay within some known efficacy and toxicity thresholds. For a similar analysis, see the example Scan Dosing Regimens Using SimBiology Model Analyzer App.

## 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

SimBiology Model Analyzer | `Observable`