MatVPC: A User-Friendly MATLAB-Based Tool for the Simulation and Evaluation of Systems Pharmacology Models.

Quantitative systems pharmacology (QSP) models are progressively entering the arena of contemporary pharmacology. The efficient implementation and evaluation of complex QSP models necessitates the development of flexible computational tools that are built into QSP mainstream software. To this end, we present MatVPC, a versatile MATLAB-based tool that accommodates QSP models of any complexity level. MatVPC executes Monte Carlo simulations as well as automatic construction of visual predictive checks (VPCs) and quantified VPCs (QVPCs).

VPC is a model diagnostic tool that facilitates the evaluation of both the structural and the stochastic part of a model. It is constructed by superimposing the observations over the model simulations while accounting for both the interindividual variability as well as the residual variability.1 Once underutilized,2 the VPC now is recognized as one of the most valuable model diagnostics in pharmacological model evaluation.3–5 Its superiority over comparable diagnostic tools has been established6 and reflected by the fact that regulatory agencies recommend it as one of the central model diagnostics.7

The VPC has recently evolved and key features have been added to its originally simple version.8–10 These features aid in better visualizing the data and include but are not limited to summary statistics and binning of the simulated and observed data. Extensions to the VPC have also been developed aiming for a highly objective model evaluation. A VPC extension example includes QVPC, wherein the distribution of the observations around the predicted median trend is rigorously realized.11

To ease the VPC construction process, several practical software products have been developed that automatically generate comprehensive VPC plots, which are challenging to generate otherwise. Such software products include Monolix12 and Xpose-PsN (coupled with NONMEM),13,14 among others.15–17 Monolix carries out the simulations and uses the outcome to generate the VPC plots, whereas PsN uses NONMEM18 as the simulation platform and performs the percentile calculations. Subsequently, the latter are plotted using Xpose. The immense practicality of these tools, however, comes with its set of shortcomings: (i) they inevitably carry out parameter estimation of a model before constructing its VPC; (ii) they do not offer extensions to the VPC; and (iii) generating VPCs of high complexity datasets may be devious. These drawbacks render these tools unattractive for constructing VPCs of high complexity models or of models that have already been built in different software and require no further parameter estimation.

Recently, a growing body of literature has been demonstrating the catalytic role of high complexity, systems-level models in the pharmacology arena.19–27 Traditionally, the so-called “QSP” models are not implemented in Monolix or NONMEM but in surrogate software tools, such as MATLAB or MATLAB-based software.8,28–37 VPC plots that are generated with these tools are generally not as comprehensive31,38 lacking important features, such as confidence intervals of simulation percentiles or data binning. As it was previously discussed, VPCs that lack these features may not be sufficiently informative.1,9,11

The progressive incorporation of QSP models in the regulatory arena39 necessitates the development of flexible computational tools that are built into software commonly used for such models; these tools can be then routinely utilized in evaluating and simulating QSP models. To this end, we developed MatVPC: a flexible, user-friendly MATLAB-based tool that carries out Monte Carlo simulation as well as VPC construction of QSP models of any complexity level. MatVPC accommodates any model, independently of the software used for model development, while offering the majority of options that comparable tools list. Provided that certain features of a model may not be readily visualized in a VPC,11 MatVPC grants the additional option of constructing QVPC plots of the available data.

MATVPC TUTORIAL

In this tool, the user simply inputs: (i) the model equations; (ii) the model parameters; and (iii) the NONMEM-like dataset with observations, and MatVPC outputs VPCs, QVPCs, or Monte Carlo simulation plots at will. The inputted model parameters here include macrokinetic or microkinetic constants, initial conditions of ordinary differential equations (ODEs) and variability terms; MatVPC allows for implementing both interindividual variability and RUV, which can be in the form of proportional, additive, or combined error. The dataset should be provided in a comma-separated values format including, at a minimum, the following headers in no particular order and in a case-insensitive manner: (i) TIME; (ii) ID; (iii) CMT; (iv) DV (nonexistent DVs should be represented by a dot); (v) EVID; (vi) AMT; and (vii) RATE (only in the case of intravenous infusion dosing).

Once the user provides the aforementioned information, they must insert in the command line: (i) the name of the dataset; and (ii) the stratification variables, if any, as follows:

MatVPC (‘dataset.csv’, {‘stratification1’, ‘stratification2’}). Common stratification variables include but are not limited to arms of a trial (placebo and active control), routes of administration, dose intervals, and, generally, covariates affecting any model parameters.1 At present, MatVPC accepts up to two stratification variables, whereas compartment is a default stratification variable and MatVPC will always stratify on it. It should be highlighted that MatVPC treats each ODE as a “compartment,” however, it only constructs VPC plots for the compartments with available observations. Once the information (i) and (ii) is provided, the user is asked (iii) to input the name of the model, which can be any alphanumeric value (e.g., Model1); and (iv) to state whether they want to run simulations and construct VPCs or to simply construct VPCs using previously saved simulation results. Please note that the former option requires the user to simply provide a desired name for the model, unlike the MATLAB file ODEs, wherein the user provides the model equations. Option 4 is expected to be useful in cases where simulations of a model have already been carried out and the user wants to simply reconstruct VPCs with modified characteristics, such as different color lines or altered binning edges. Parenthetically, when this option is exploited, the simulation-related parameters should be set to the same values as in the original VPC construction.

Upon providing all the above, a user-friendly interface commences allowing the user to modify the default values in a series of optional inputs. This input includes options related to the methods (e.g., ODE solver to be used or number of datasets to be simulated) as well as options that determine the graphical settings (e.g., color and type of plotting curves). Table 1 details the accepted optional input along with the default values.

Table 1

Default values and accepted values included in MatVPC's optional input interface

Options	Default value	Accepted input
No. of datasets	200	Scalar
Automatic binninga	Yes	Yes, no
Manual binningb	No	Yes, no
Bin edges	[ ]	Vector with as many scalars as the bin edges
No. of simulated timepointsc	100	Scalar
Run Monte Carlo simulationsd	No	Yes, no
ODE solvere	ode23tb	ode45, ode23, ode113, ode15s, ode23s, ode23t, ode23tb
Simulated time beyond observationsf	0	Scalar
Upper percentile limit of simulations	95	Scalar
Lower percentile limit of simulations	5	Scalar
Plot observationsg	Yes	Yes, no
Plot percentiles of observations	Yes	Yes, no
Upper percentile limit of observations	95	Scalar
Lower percentile limit of observations	5	Scalar
Plot QVPC	Yes	Yes, no
Replace negative and zero valuesh	No	Yes, no
Replacement value	0.001	Scalar
Graphical settings
Type of median curve in simulations	–	–, – –,:, –., o, +, *,., :, x, s, d, ^, v, >, <, p, h
Color of median curve in simulationsi	Black	Color name
Type of percentiles curve in simulations	– –	–, – –,:, –., o, +, *,., :, x, s, d, ^, v, >, <, p, h
Color of percentiles curve in simulationsi	Black	Color name
Color of median CI in simulationsi	Deep pink	Color name
Color of percentiles CI simulationsi	Blue	Color name
Color of observationsi	Blue	Color name
Size of observations	10	Scalar
Color of median curve in observationsi	Red	Color name
Color of percentiles curve in observationsi	– –	–, – –,:, –., o, +, *,., :, x, s, d, ^, v, >, <, p, h
Type of median curve in observations	–	–, – –,:, –., o, +, *,., :, x, s, d, ^, v, >, <, p, h
Type of percentiles curve in observations	– –	–, – –,:, –., o, +, *,., :, x, s, d, ^, v, >, <, p, h
Color of simulation median in QVPCi	Red	Color name
Color of observations above median in QVPCi	Grey	Color name
Color of observations below median in QVPCi	Black	Color name
Number of rows in figure	1	Scalar
Number of columns in figure	3	Scalar
Square shaped plots	Yes	Yes, no
Width of all curves	1	Scalar
Log-y scale	No	Yes, no
X-axis label	Time	Any alphabetic input
Y-axis label	Observations/simulations	Any alphabetic input
Transparency of CI shaded areas	0.5	Scalar

CI, Confidence Intervals; ODE, ordinary differential equation; QVPC, quantified visual predictive check.

For a detailed description of how the bin edges are selected, please refer to Ref.9 and40. Before the results of the automatic binning are accepted, a visual inspection of the calculated bin edges should be done as, in extreme cases of data point distributions, the code might not generate optimal results.

In case “yes” is selected here, a vector with the bin edges should be defined in the following option.

This option allows modification of the number of equally spaced timepoints that are saved during the Monte Carlo simulations. The larger the numbers, the slower the runtime.

This option requests Monte Carlo simulations can be selected in other options (see footnotes c and f).

The ODE solver should be carefully chosen based on the type of the model, in other words, whether the model is stiff or not.

This number represents the percentage of the entire simulation time that is simulated beyond the latest observation point. For instance, in case the user desires to run a Monte Carlo simulation for an additional time equal to 20% of the entire profile, this number should be set equal to 0.2 (see Monte Carlo simulations in model one).

This option allows the user to plot simulation results without plotting the observations. In case the user desires to simply plot Monte Carlo simulations of the model, “no” should be selected.

This option should be exploited in case the user wants to avoid negative concentrations and associated plotting issues, for instance, when a logarithmic scale is used. In case “yes” is selected here, the desired replacement value should be defined in the following option.

The color of the curves is defined by writing the name of any of the 139 available color names provided at http://www.w3.org/TR/css3-color/. Note that the first letter should be capital.

It is important to stress that, unlike analogous tools that are traditionally used in pharmacometrics, MatVPC constitutes an all-in-one package that integrates the following:

It is publicly open (at https://sourceforge.net/projects/matvpc/);

Constructs VPC plots of complex QSP models;

Offers automatic binning of the observed and simulated data using a rigorous approach9,40;

Constructs QVPC plots of complex QSP models;

Performs Monte Carlo simulations of the model and plots the results with any requested summary statistics;

Does not require prior implementation (e.g., parameter estimation) of the model within MATLAB;

Endows the user with the option to perform post-processing of the VPC, QVPC, or simulation plots.

In what follows, we implement three models with varying characteristics to demonstrate the functionality of MatVPC. These models include: (i) a linear three compartment pharmacokinetic (PK) model with single oral and intravenous (IV) bolus dosing; (ii) a nonlinear two compartment PK model with multiple IV infusion dosing41; and (iii) a highly nonlinear pharmacodynamic (PD) model that describes the time-course of body weight.42 To validate the accuracy of the presented VPCs, we directly compare them with VPCs generated either by Xpose version 4 (coupled with PsN version 4.2.0 and NONMEM version 7.2) or by Monolix version 4.3.3.

Model one: three-compartment PK model with oral and IV bolus dosing

The first model describes simulated data of a linear three-compartment PK model with linear elimination from the first and third compartment. The data were assumed to stem from two different studies. The difference between the two studies is nothing but the dosing type; in study one, the dosing type was IV bolus, whereas in study two, the dosing type was oral. A total of 50 individuals were simulated for each study. A single dose was administered at time zero and observations were taken at asymmetric timepoints. The initial conditions of all compartments were set to zero and the parameters of the model were estimated with NONMEM. The ODEs of the model are shown in Eqs. 1-4:

To generate VPC and QVPC plots of model one using MatVPC, the user is first required to open the MATLAB file parameters and provide the model parameters therein. Once the parameters have been provided, the user must open the MATLAB file ODEs and type the Eqs. 1-4. A snapshot of the parameters and the ODEs file upon inserting all the pertinent information of model one is shown in Figures 1a,b, respectively.

Figure 1

Required input in MatVPC for model one. (a) Input in file “parameters.” In this file, the user inputs the values of the ODE parameters, the values of the interindividual variability parameters, the proportional and additive part of the RUV, the initial conditions of the ODEs and the volume size of each compartment. (b) Input in file “ODEs.” In this file, the user inputs the ODEs, along with potential algebraic equations, describing the quantitative pharmacology model. In case of IV infusion dosing, an additional term should be included in the model that accounts for the dosing (for an example see model two).

Upon completion of the above steps, the following command must be typed in the command window: MatVPC('obs.csv’,{'study’}). This command requests a VPC that is based on the dataset called obs.csv and is stratified by study. Once the user inserts this command and presses “enter,” they are asked to provide the name of the model and whether they want to simply create a VPC from previously saved simulation results or to simulate new data and plot the VPCs. Following their response (y/n), MatVPC opens a user-friendly interface with optional input that bears default values. Figure 2 shows a snapshot of this interface, which includes the optional input used for constructing VPC and QVPCs for model one. After pressing “OK” at the bottom of the interface, MatVPC runs all the required simulations, calculates the respective percentiles, constructs the VPC and QVPC plots, and saves them as MATLAB figures (.fig) in a folder that is named after the model name. Unlike comparable software, by saving the plots as MATLAB figures, MatVPC allows for postprocessing of the plots. More specifically, the user can open the figures (VPC, QVPC, or Monte Carlo simulation plots) that are saved in. fig format and interactively modify their characteristics, such as: (i) the type, color, or size of the curves; (ii) the axes labels and their font size; (iii) the figure title and its font size; (iv) the size, shape, or resolution of the figure; and (v) the axes limits. Modifying these characteristics can also be done using the command line, but it is not recommended as it is usually rather laborious. More details about this can be found on MATLAB's website (http://www.mathworks.com/help/symbolic/edit-graphs.html?searchHighlight=edit).

Figure 2

Optional input in MatVPC for model one. In this interface, the user can modify the default values of 40 characteristics of the VPC and QVPC plots (for a detailed description of these features see Table 1). The VPC and QVPC plots shown in Figure 3 were constructed using the input values shown in this figure.

Figure 3

VPC and QVPC plots for model one. (a,d) VPC plots generated by MatVPC showing the results from compartment two and study one (a) and two (d). (b,e) VPC plots generated by Xpose-PsN showing the results from compartment two and study one (b) and two (e). Blue dots correspond to the observations. Red dashed lines correspond to the 5th and 95th percentiles of the observations, whereas red solid lines correspond to the median of the observations. Black dashed lines correspond to the 5th and 95th percentiles of the simulations, whereas black solid lines correspond to the median of the simulations. Blue shaded areas represent the 90% confidence intervals of the simulation 5th and 95th percentiles, whereas pink shaded areas represent the 90% confidence intervals of the simulation median. (c,f) QVPC plots of compartment two and study one (c) and two (f), as constructed by MatVPC. At each timepoint, the black bar presents the observed data below the model predicted median (red dots), whereas the dark grey bar shows the observed data above the model predicted median. The total of the black and grey bar combined presents the percentage of available data (here 100%). (g–i) The 200 Monte Carlo simulations of model one and study two (oral dosing) for compartment one (g), two (h), and three (i). Green dashed lines correspond to the median of the 5th and 95th percentiles of the Monte Carlo simulations. Red lines correspond to the median of the 50th percentile of the Monte Carlo simulations. Shaded areas represent the 90% prediction intervals of the Monte Carlo simulation percentiles.

As the sampling times in this case were asymmetric, we exploited the automatic binning option of MatVPC to bin the data of model one (i.e., “yes” was selected in the “automatic binning” option). The VPC and QVPC plots of model one as generated by MatVPC are shown in Figure 3. VPC plots of model one were also constructed with Xpose and PsN (running the simulations in NONMEM and using the bin edges calculated by MatVPC) and are also provided in Figure 3. For simplicity, we only show the results from compartment two of each study but the agreement between MatVPC plots and Xpose plots was consistent for all compartments.

As discussed above, one of the key characteristics of MatVPC is the option of performing Monte Carlo simulations of a QSP model, either accounting for or neglecting the residual error. Here, we carried out Monte Carlo simulations of model one (only for study two) without taking into consideration the residual error. To do so, the RUV terms in the “parameters” file were set equal to zero (see Supplementary Figure S1). In addition, in the optional input interface, for “plot observations” was selected “no” and the “simulated time beyond observations” was set equal to 0.2, as we were interested in the dynamics of the model from time zero to 20% additional time beyond the timepoint of the latest observation (see Supplementary Figure S2). The results of the simulations of the three compartments of study two (oral dosing) are presented in Figure 3.

Model two: two-compartment PK model with multidose IV infusion

The second model was adopted from a previously published study41 and describes simulated PK data of a drug administered via IV infusion. In this case, a total of 100 individuals were simulated with asymmetrically sampled observations. The parameters of this nonlinear PK model were estimated in NONMEM. The model consists of Eqs. 5 and 6:

It is important to note that in case IV infusion dosing is simulated, an additional term should be added in the model to account for IV infusion dosing. This term is simply “params(n),” where “n” is equal to the number of model parameters plus one (please note that all parameters are inserted as “params(i)”); “params(n)” should invariably be: (i) the last parameter listed in the file “parameters,”; (ii) equated to zero; and (iii) added to the equation describing the dynamics of the compartment in which IV dosing is administered. In this example, the term “params(5)” was added to Eq. 5 as the IV infusion is given in compartment one (see file “parameters” and “ODEs” in Figure4).

Figure 4

Required input in MatVPC for model two. (a) Input in file “parameters.” Please notice the extra term, “params(5),” that has been added to the parameter list. (b) Input in file “ODEs.” Please note that “params(5)” has been added to Eq. dy(1) to account for IV infusion dosing.

Similarly to model one, in order to generate the VPC and QVPC plots of model two using MatVPC, the user must open the MATLAB file “parameters” and insert the model parameters. Once the parameters have been provided, the user must open the MATLAB file “ODEs” and type the ODEs shown in Eqs. 5 and 6, along with the additional term that accounts for IV infusion dosing. A snapshot of the completed MATLAB files “parameters” and “ODEs” file is depicted in Figures 4a,b, respectively. Please note that in file “parameters” an additional parameter, “params(5)=0,” has been listed. As shown in file “ODEs,” this parameter has then been incorporated in the first ODE of the model (see last term in first ODE).

Upon inserting this information in the two files, the user should type in the command window the following: “MatVPC(‘obs.csv’)”. MatVPC then asks the user to provide the model name and to define whether they want to construct VPC plots from scratch or to capitalize on previously generated simulation data. Upon completing these steps, a MATLAB interface with the optional input pops up. A snapshot of this interface with the optional input used for model two is shown in Supplementary Figure S3. Please note that automatic data binning has been requested through this interface as the sampling timepoints were not consistent and manual binning of data would be challenging. After pressing “OK” at the bottom of this interface, MatVPC undergoes all the necessary steps and constructs the VPC and QVPC plots of model two. Figures 5a,c illustrate the VPC and QVPC plots, respectively, as constructed by MatVPC. Figure 5b shows the VPC of model two as generated with PsN and Xpose (with the simulations conducted in NONMEM and using bin edges calculated by MatVPC). For simplicity, only the plots from compartment one are shown.

Figure 5

VPC and QVPC plots for model two. (a) VPC plots of compartment one generated by MatVPC. (b) VPC plots of compartment one generated by Xpose-PsN. Blue dots correspond to the observations. Red dashed lines correspond to the 5th and 95th percentiles of the observations whereas red solid lines correspond to the median of the observations. Black dashed lines correspond to the 5th and 95th percentiles of the simulations, whereas black solid lines correspond to the median of the simulations. Blue shaded areas represent the 90% confidence intervals of the simulation 5th and 95th percentiles, whereas pink shaded areas represent the 90% confidence intervals of the simulation median. (c) QVPC plots of compartment one as generated by MatVPC. At each timepoint, the black bar presents the observed data below the model predicted median (red dots), whereas the dark grey bar shows the observed data above the model predicted median. The total of the black and grey bar combined presents the percentage of available data (here 100%). (d) The 200 Monte Carlo simulations of model two for compartment one. Purple dashed lines correspond to the median of the 5th and 95th percentiles of the Monte Carlo simulations. Purple solid lines correspond to the median of the 50th percentile of the Monte Carlo simulations. Red shaded areas represent the 90% prediction intervals of the Monte Carlo simulation percentiles.

In addition to generating the VPC and QVPC plots of model two, we also carried out 200 Monte Carlo simulations and the results are presented in Figure 5d. The “parameters” and “ODEs” files, along with the selected options in MatVPC interface, that were used to generate the Monte Carlo simulation plots are shown in Supplementary Figures S4 and S5, respectively.

Model three: PD model describing the time-course of body weight

The third case example is a formerly presented highly nonlinear PD model that describes the time-course of body weight.42 The simulated dataset includes observations from a total of 500 subjects that are involved in two different studies, 250 subjects belonging to study one and 250 subjects belonging to study two. The difference between the two studies here is the sampling times. The PD model is provided in Eq. 7: whereby k is the rate of weight gain, k is the rate of weight loss, k and k are associated with the onset and loss of the lifestyle intervention effect and DSTIM is the maximum fractional increase in k caused by the lifestyle intervention. The values of these parameters were estimated using Monolix.

The “parameters” and “ODEs” files that were utilized in MatVPC to construct the VPC and QVPC plots of model three are shown in Figures 6a,b, respectively. As the data originate from two distinct studies, the user would need to stratify on the study number. To do so, the following should be typed in the command window: “MatVPC ('obs.csv’,{'study’}).” This command requests VPC and QVPC plots of the PD model shown in Eq. 7, using the dataset called “obs.csv” while stratifying on the variable “study.” Once the user defines the name of the model, and whether they want to use previous simulation results, the MatVPC interface with the optional input appears. The optional input used in the construction of VPCs for model three is provided in Supplementary Figure S6. Upon pressing “OK” in this interface, the VPCs and QVPCs of the model are generated and saved. These VPC and QVPC plots are illustrated in Figure 7. The respective VPC plots were also generated with Monolix and are provided in Figures 7b,e. It should be underlined that the automatic binning approach implemented in MatVPC and Monolix is identical.9,40

Figure 6

Required input in MatVPC for model three. (a) Input in file “parameters.” (b) Input in file “ODEs.” For a detailed description about the input see caption in Figure 1.

Figure 7

VPC and QVPC plots for model three. (a,d) VPC plots of study one (a) and study two (d) generated by MatVPC. (b,e) VPC plots of study one (b) and study two (e) generated by Monolix. Blue dots correspond to the observations. Green lines correspond to the 5th, 50th, and 95th percentiles of the observations. Black lines correspond to the 5th, 50th, and 95th percentiles of the simulations. Light blue shaded areas represent the 90% confidence intervals of the simulation 5th and 95th percentiles whereas pink shaded areas represent the 90% confidence intervals of the simulation median. (c,f) QVPC plots of study one (c) and two (f), as generated by MatVPC. At each timepoint, the black bar presents the observed data below the model predicted median (red dots), whereas the dark grey bar shows the observed data above the model predicted median. The total of the black and grey bar combined presents the percentage of available data (here 100%).

DISCUSSION

We have presented a MATLAB-based tool, dubbed MatVPC, that simulates and generates VPC and QVPC plots of QSP models. The plots generated by MatVPC were compared with plots constructed by the gold standard tools in the pharmacometrics community, PsN (with NONMEM) and Monolix,43 and the results were interchangeable across all case studies.

Despite its flexibility, MatVPC does have its drawbacks. At present, MatVPC generates QVPC plots using only the available data and does not account for missing data.11 In addition, it can currently run only on a single computer and is not designed for running in parallel processors. However, MatVPC will be continuously advanced and its future versions will address these imperfections.

MatVPC is publicly available and can be utilized by users with little or no prior MATLAB experience. This computational tool can be potentially expanded to perform key analyses in systems pharmacology,44 such as sensitivity analysis45 or model reduction.38,46 Collectively, MatVPC constitutes a useful addition to the openly available toolboxes exploited by quantitative as well as clinical pharmacologists.