Potential Issues With Virtual Populations When Applied to Nonlinear Quantitative Systems Pharmacology Models.

Quantitative systems pharmacology (QSP) models attempt to describe the pharmacological properties of a drug. Although there are a wide variety of structures (e.g., logic models, fuzzy logic) the most common is based on ordinary differential equations (ODEs). In this perspective, we explore two case studies that are based on nonlinear drug actions described by ODEs where we show that regions of the parameter space exist that result in model predictions of system responses that are unacceptable.

WHY VIRTUAL POPULATIONS?

Variability in physiological characteristics either naturally or because of a pathophysiological condition contributes to patient heterogeneity. In a QSP model, variability in pathways or parameters can be introduced to explore its effect on a biomarker or end point of interest. Each instance of a simulation represents a virtual patient, and ensembles are known as virtual cohorts with a set of virtual cohorts referred to as a “virtual population” (VP). A VP should reflect individual and population level characteristics relative to a real patient or population. This is achieved by simulating vectors of parameter values from either a naïve or optimized density such that their frequency matches some observed data quantities.

To reproduce the statistics of a clinical population, virtual patients are often weighted to form a VP that reflects the baseline characteristics of the clinical cohort. There are also techniques that illustrate how to select a VP that matches the observed data without the need for weighting. Real‐world data have also been used to develop VP databases. Random sampling of model parameters from a multivariate log‐normal distribution or carrying out a sensitivity analysis are other methods commonly used to identify parameters for VPs. Readers are referred to the excellent review of Allen et al. for a more full description of VPs.

Using VPs to simulate virtual clinical trials is a common strategy to validate a QSP model and gain confidence in the resulting predictions. VPs can also be used to study the impact of variability within the patient population on a clinical trial outcome as well as to optimize patient inclusion criteria for maximizing the chance of a conclusive trial. If the purpose of the simulations is to get an idea of the range of response profiles based on the model, hundreds or thousands of virtual patients may need to be simulated and the variability across all presented as a prediction interval. Multiple VPs can be developed to explore population variability in biomarkers or differentiate responders and nonresponders for monotherapies, combination therapies, or sequential therapies.

APPLYING VIRTUAL POPULATIONS TO NONLINEAR QSP MODELS

Simulating response profiles from large models based on generating random sets of parameter values is, however, not without its problems. In our cases, we are referring to QSP models that are defined as ODEs.

Linear systems of ODEs, for example, those that have constant coefficients, will yield a predictable “signature” response (albeit scale is arbitrary) for any set of (legal) parameter values. Nonlinear systems of ODEs, for example, an ODE that has a coefficient that is dependent on its own state space variable, however, are only locally defined to a set of parameter values, and perturbation of a nonlinear system may or may not reveal the expected signature profile. This means that even parameter vectors in neighboring regions may behave differently to that anticipated. Nonlinearity in QSP models occurs as a result of negative and positive feedback (or forward) processes. There may be many positive or negative feedback processes in a component of a system. These components are either overall damping and resist perturbations or amplifying (also called positive loop gain) that enhance perturbations. For instance, the insulin‐glucose system (a mixture of positive and negative processes) is designed to resist instability caused by ingestion of a meal of glucose and hence the system is damping, whereas the coagulation system is designed to respond to vascular damage and the release of tissue factor with an amplified system response to form a clot that otherwise would not naturally form.

For these nonlinear components in QSP models, we contend that regions within the parameter space may exist that will yield profiles of state variables that do not conform to the anticipated signature profile. In addition, these effects may not be obvious (without close scrutiny) if the response variable of interest is distal to the state space affected. This contrasts with linear systems where every set of parameter values would yield a scaled anticipated profile (e.g., for a pharmacokinetic model for every set of parameter values [] the concentrations of drug increase when absorption dominates and decreases as elimination dominates).

To explore this further, we consider two simple examples that illustrate damping and amplification components common in QSP models. The examples considered for demonstration resemble simplified insulin‐glucose and coagulation models, but are not replicates of these systems.

Damping processes

We chose a simple damping system that is based (loosely) on glucose and insulin. In this system, we present a stimulus that perturbs the value of (resembling glucose), which is part of a damping process, with the compound (resembling insulin). We can monitor both species, but not the stimulus. In this simple model, compound stimulates the release of and in turn reduces the natural production of . The model consists of three state variables and nine parameter values. For the purposes of this exercise and simplicity, we consider only two parameters: (i) half‐life of and (ii) half‐life of . We then simulated virtual patients by drawing vectors of parameter values from the parameter space of each species; in both cases, the draw is uniform on the space of . We define the anticipated signature profile as a system that returns back to its starting values once the perturbation has been damped. In addition, our signature profile does not contain a rebound effect that is defined as the concentration of is forced below its basal level (if were glucose then this would represent hypoglycemia). Our simulation is shown in Figure  . We see that half of the signature profiles demonstrate a rebound hypoglycemia and 25% of the profiles achieve a postperturbation stationary value that is different from the basal level. The latter indicating that the system is (perhaps) permanently broken. The model description and code is provided in the Supplementary Material .

Figure 1

A surface plot of joint parameter vector space. The surface illustrates profiles that are determined as acceptable (blue region with example), partially acceptable (yellow region) where there is evidence of rebound, and unacceptable (red region with example) where the system does not return to its pre‐perturbation conditions. The upper left panel represents the normal setting where both A 1 and B 1 return to basal values. We use [.] to denote concentration. The upper right panel represents an unacceptable profile.

Amplification processes

For this example, we created a six‐state model that demonstrates amplification. We define amplification as a system that contains a positive loop gain in which an increase in values of state variables from nonreactive, e.g., do not invoke a noticeable system effect (often zero), to those that elicit accumulation of a response of interest. In particular, this accumulation often serves as a stimulus for an event (e.g., clot formation, parturition). The amplification system explored here is based (loosely) on the effect of tissue factor on the coagulation system. In this system, we have a stimulus, two precursor state variables ( and ), one active state (), a complex and a response variable, . These roughly correspond to tissue factor, factor VII, a zymogen (perhaps X and II), a serine protease (activated zymogen), and fibrin, respectively. Our response of interest is the cumulation of species . The initial values of all precursor state variables are normalized to 1, and the inactive state and the response variable are zero. In total, there are 13 parameters. For the purposes of this exercise and simplicity, we consider only two parameters: (i) half‐life of and (ii) the parameter representing positive gain. We then simulate virtual patients by drawing vectors of parameter values from the parameter space for each of half‐life and positive gain. In both cases, the draw is uniform with ranges of , respectively. The anticipated signature profile results in a system that forms a plateau in production of and the system returns back to its initial values within an acceptable period of time (i.e., within five half‐lives of the longest half‐life parameter) after the perturbation has been initiated. This means that inactive variables and return to unity, and the active variable returns to 0. In essence the system resets itself. We do not require the accumulation state () to dissipate at this point as we have not included a mechanism for its eventual degradation. Our simulation is shown in Figure  . When we simulate our virtual patients, we see that most profiles show the anticipated signature behavior in which amplification occurs to form a stable response and then the system resets itself. However, in 10% of cases we see that the response variable has not stabilized and the system has not reset. This situation then continues without evidence of abating for a time that far exceeds five half‐lives of the longest factor. The latter indicating that the system is not behaving acceptably over a fair experimental setting.

Figure 2

A surface plot of joint parameter vector space. The surface illustrates profiles that are determined as acceptable (blue region with example) and unacceptable (red region with example) where the system does not return to its pre‐perturbation conditions. The upper left panel represents the normal setting where A 2, B 2, C 2, return to basal values. We use [.] to denote concentration and [.] time to represent the integral. The upper right panel represents an unacceptable profile.

INFERENCE

QSP models may typically contain numerous component parts, with each part behaving in a linear or nonlinear manner depending on its function in the system and its connections with other components. These component parts may play a critical role in maintaining homeostasis (e.g., osteoclast/osteoblast activity in the bone model) and formation of stable clots when appropriate. In not all circumstances are the damping and amplification components of QSP models able to be monitored as a part of normal response monitoring and/or we do not necessarily know what the signature system behaviour should look like. It is therefore important to note that virtual populations that are generated by sampling across an assumed reasonable parameter space may not necessarily create viable signature profiles across the system states. This may in itself be of interest diagnostically to illustrate some issue with the model not otherwise recognized or simply be a curse of nonlinearity. The examples chosen here were deliberately simple with permissive criteria about what an acceptable profile would look like (i.e., one that resumes its basal characteristics). A significant number of these virtual populations may also have failed to produce plausible results if they were connected to other systems. In addition to the issues presented in this perspective, there remains the practical challenge of solving systems of ODEs for randomly generated parameter vectors when the system becomes ill defined.

A purpose of creating virtual populations is to explore the influence of between‐subject and within‐subject variabilities in the system. Because we contend that this may be difficult to realize effectively by reviewing profiles of state variables for plausibility, depending on the aims of the analysis, one suggestion is to consider a model‐order reduction method approach to render the full systems model into a smaller mechanistic input–output model that can be used to estimate the between‐subject variances in the parameters based on the available data and thereby avoid the issues associated with generating virtual patients.

Supplementary Material S1

Click here for additional data file.