a flow diagram of the algorithm is shown in figure 
 1 
 . to implement this procedure for a given model
it is necessary to define bounds for input parameters and model outputs ( e.g. , steady states or dynamic behavior ) .
if bounds can not be defined empirically , feasible ranges of parameter values can be asserted from physiological knowledge or theoretical considerations .
for example , the tissue concentration of a species may not be known , but typical weight and water content of that tissue may be known , which allows us to put an upper limit on the species concentration .
overview of algorithm for efficient generation and prevalencebased selection of virtual patients . to generate virtual patients from a model , the prior information ( green boxes )
a virtual population is constructed by selecting from this population with probability proportional to the prevalence in the real population relative to the prevalence in the plausible population .
this selection is optimized to produce the best virtual population given the patients in the plausible population .
we have provided a detailed description of terminology , definitions , and the derivation of this algorithm in table 
 1 and the supplementary material .
 we define these patients as a parameter set for which every component of the model ( whether it be the parameter values themselves , computed species concentrations , or combinations of these that are experimentally measurable ) falls into a biologically plausible range . from this
 plausible population  we can then select the virtual population such that it matches the empirical distribution of interest .
this is achieved by calculating a probability of inclusion of a plausible patient into the virtual population .
this probability is computed from both the empirical distribution and the density of plausible patients ( see supplementary materials for more details ) .
an overview of the terminology used in this article an important prerequisite to this approach is the ability to generate a large number of plausible patients within the region of the empirical data . to accelerate this process we take an initial parameter guess ( within the predefined bounds ) and optimize this choice until the required outputs are within physiologically plausible ranges . rather than optimize to specific points ,
it is more efficient to be agnostic as to where in the plausible ranges the optimization routine ends . to implement this we shift the typical cost function f(p ) we would use optimizing a model to a new function , g(p ) , where we consider both as purely dependent on the parameter
if we constrain parameters using a number of model outputs mi(p ) , with data di then f ( in the simplest , unweighted case ) would be : 
 f(p )  =  (mi(p)di)2 . 
  to generate plausible patients , we modify this sumofsquared errors expression to : 
 g(p )  =  imax[(mi(p)li+ui2)2(ui2li2)2,0]where ui and li are the predefined plausible upper and lower bounds , respectively , for mi(p ) .
this expression ensures that if mi(p ) is in the plausible range then the contribution of the corresponding term in the expression is zero .
the effect of replacing f(p ) with g(p ) is visualized in 2d in figure 
 2 .
outputs of the model contribute to the cost function to be minimized by considering the sum of squared errors ( sse ) from an associated experimental observation .
for each observation we define a physiologically plausible range ( arrows in a , b ) and shift the sse associated with that observation so that it is zero if the model output is in this range ( a , b ) . combining these transformations in each dimension
leads to a broader cost function that is minimized by many points , rather than one ( black rectangle in c ) . to test this approach we used a previously published model of cholesterol metabolism.17 we chose this model because we could use publicly available data from the nhanes database16 to establish the empirical multivariate distribution for ldl cholesterol , hdl cholesterol , and total cholesterol ( ldlc , hdlc , and tc , respectively ) .
note that the distribution of these variables is well approximated by a multivariate lognormal distribution ( supplementary figure 1 ) . for the remainder of the article we will describe these variables , either as model outputs or from nhanes , in log units ( prior to taking the logarithm , units are mg / dl for cholesterol measures ) .
the published version of this model does not explicitly calculate ldlc or tc ; instead , the outputs are hdlc and nonhdlc . from these two quantities tc
is easily calculated . for full comparison with the nhanes data we introduced a new parameter to the model ,
input_ranges.m  gives details on parameter and output ranges for the van de pas model .
also in the supplementary material is the code used in this case , which is easily modifiable for application to other models .
as expected , the initial plausible population does not match the populationlevel statistics of the nhanes data ( figure 
 3 ) but covers the empirical distribution ( i.e. , where there are likely to be empirical observations there are plausible patients ) .
comparison of the initial plausible population ( n  =  300,000 ) with nhanes multivariate distribution ( ( ac ) black dotted lines estimated pdf , supplementary figure 1ac .
( df ) 2d projection of the 95% confidence surface of the estimated probability density function ) .
once calculated , we established that most of the plausible patients are highly unlikely to be in the final distribution ( figure 
 4 ) .
this is due to the relative density of the plausible population to the empirical distribution . with these probabilities ,
only 2% of the plausible population was selected to be in the virtual population ( inset , figure 
 4 ) .
based on examining goodnessoffit of the distributions it appears , in this case , there is no further value in increasing the size of the plausible population ( supplementary figure 2 ) .
the red histogram ( main figure , and figure inset ) is a virtual population that matches nhanes data , and is selected from the plausible population ( blue histogram ) based on displayed probability .
the distribution of an example selection fits the nhanes data well ( figure 
 5 ) .
the 1d histograms ( when normalized for comparison with the nhanes probability density function ) are indistinguishable from the data ( figure 
 5 
 a  c ) and the correlations between variables also match the data based on visual predictive check .
the virtual population ( red dots and red histogram ) matches the mean , variance , and covariance of the multivariate experimental distribution ( ( ac ) black dotted lines estimated probability density function , supplementary figure 1ac .
( df ) 2d projection of the 95% confidence surface of the estimated pdf ) .
when selecting a subset of vps from a larger population , one concern is that the selected subset of vps does not reflect the variability of the original ensemble , which was generated from the biologically plausible range of the parameters .
analyzing the final fitted population , we found little change in either the distribution of parameters or the correlation structure between the parameters ( figure 
 6 and supplementary figure 3 ) .
this also shows that despite the virtual population being constrained against the nhanes data the parameter values of the virtual population ( figure 
 6 
 b ) are only slightly better constrained than those of the plausible population .
furthermore , correlations between parameters are only slightly increased in the virtual population ( figure 
 6 
 d ) vs. the plausible population ( figure 
 6 
 c ) . at least in this case , constraining all outputs into realistic ranges is a more stringent constraint than selection of a virtual population .
violin plots of the plausible and virtual populations ( a , b , respectively ) parameter values ( normalized to each parameter 's upper and lower bounds ) and correlation matrix of the plausible and virtual population ( c , d , respectively ) .
one of the primary uses for qsp models is to prospectively simulate the effects of a dynamic perturbation ( pharmacological or otherwise ) in populations of interest . due to the underconstrained nature of these models
it would be difficult to have confidence in the simulation results if we simulated a single parameterization of that model , even if that set of parameters is an excellent fit to the available data .
for example , imagine creating a single hypercholesterolemia vp to simulate the effect of various anticholesterol therapies .
for the baseline characteristics of the vpop , we have good data for the expected mean ldl and hdl ( e.g. , a prior clinical cohort ) , but we could still choose to mechanistically represent hypercholesterolemia several ways using the same model .
we could increase cholesterol production , decrease clearance , or apply some combination of both .
our choice of how to parameterize that vp could have significant consequences for the sensitivity of the followon therapy simulations .
having impaired production vs. clearance of ldl could lead to differential responses to statins ( production ) vs. antiproprotein convertase subtillisin / kexin type 9 ( clearance ) .
a better approach is to explore the underconstrained nature of these models and sample the biological uncertainty in the creation of plausible patients by varying mechanistic parameters , such as production and clearance rates , within biologically reasonable ranges .
simultaneously , we need to constrain the higherlevel observables of the model based on known population distributions ( e.g. , the baseline characteristics of a clinical trial cohort ) .
an appealing aspect of the approach we outline is that , since the algorithm is probabilistic , once the plausible patients are generated any number of subpopulations can be selected as long as the generated patients reasonably cover the full range of the ( sub)population . additionally , for any particular population , any number of vpops can be reselected to bootstrap the sensitivity of the model predictions to the choice of vpop . achieving an acceptable fit to the data distribution
however , if we have higherorder density functions , we will likely require additional gains in efficiency , above and beyond what is presented here , in methods for generating sufficient plausible patients .
one potential avenue , for a future iteration of this algorithm , may be to use methods that follow a directed search through the parameter space , such as using markov chain monte carlo ( mcmc ) algorithms.18 , 19 , 20 however , it should be noted that the method presented here is in fact a hybrid approach because the simulated annealing step , used to generate a plausible patient is essentially an mcmc method .
the advantage of this approach is it generates plausible patients ( and hence virtual patients ) independently  which is critical for the purpose of making a virtual population .
also , once the plausible population is established new virtual populations , suitable for new applications , can be selected .
however , a limitation of our approach is the computational cost of generating large plausible populations , which may make this method unfeasible in models that are slow to simulate .
it should be noted that prior methods13 , 14 generated smaller virtual populations in larger models , partly due to the computational cost of running the models .
nevertheless , the computational efficiency of a model does not alter the necessity for exploring the parameter space ( which in this context equates to larger virtual populations ) .
we therefore advocate for the optimization of large models for speed , such that a fuller exploration of parameter uncertainty is possible ( either via the method presented here or alternatives ) .
it is important to remark that the larger the model ( specifically , in terms of the number of poorly constrained parameters ) the greater the necessity for larger virtual populations to attempt to account for the uncertainty inherent in such a model .
an advantage of our technique over prior approaches is that it is relatively unbiased , while still leveraging all available information on possible parameter values .
the approach by klinke may overweight spurious model solutions , whereas the approach by schmidt et al .
requires binning of parameter values into mechanistic axis which , for axes containing more than one parameter , requires an assumption about the correlation between parameters in the population of interest that may not be supported by available data . as an introduction to this algorithm , we demonstrated how to generate a vpop that matches the baseline characteristics of a population or clinical cohort ; however , in practice , a dynamic model should be constrained additionally against as many inscope perturbation experiments as possible ( determined by available data ) .
for this example , simulating changes in ldlc and tc to standardofcare lipid therapies , such as statins and ezetimibe , would likely be an important step before using the model to predict the response to a novel mechanism .
ideally , information would be available detailing the distribution of the data before and after the application of therapy ( i.e. , not just summary statistics ) .
the therapeutic response can be treated as a baseline constraint for vp selection , just as we used hdlc and tc at baseline .
we have applied this approach , without major adaptation , to an unpublished model of chronic kidney disease ( 42 ordinary differential equations ( odes ) , 200 parameters ) and to a model of body weight change6 ( 8 odes , 50 parameters ) ; we have made a brief summary of the results in these cases available.21 one challenge that we foresee , for virtual populations in general ( i.e. , not specific to this approach ) , is development of efficient ways to combine populations when merging distinct models .
we believe that for two models that have large validated virtual populations , navely combining every possibility will be computationally daunting .
quantitative systems pharmacology models are becoming established as a valuable component of the drug discovery and development process . communicating their complexity and uncertainty to an interdisciplinary project team is a critical but challenging component of their utility .
virtual populations are one tool that we have found to be successful in exploring mechanistic and parametric uncertainty in an intuitive framework that is easily understandable by most audiences .
however , despite their widespread use , there are very few published methods for generating virtual patients and forming virtual populations .
here we have contributed an approach to efficiently generate virtual patients and construct a virtual population for which each patient is weighted equally .
plausible populations ,  made up of plausible patients each of which is a candidate to become a virtual patient .
because a large plausible population is necessary , models that are slow to integrate ( for example , with dynamics across multiple timescales ) may not be good candidates for this approach .
however , in cases where quantitative predictions are required and the model is amenable to thorough examination of parameter space , we have found this method to be an improvement over previous approaches .
supporting information click here for additional data file . supporting information click here for additional data file . supporting information click here for additional data file . supporting information click here for additional data file .
 quantitative systems pharmacology models mechanistically describe a biological system and the effect of drug treatment on system behavior . 
 because these models rarely are identifiable from the available data , the uncertainty in physiological parameters may be sampled to create alternative parameterizations of the model , sometimes termed  virtual patients . 
  in order to reproduce the statistics of a clinical population , virtual patients are often weighted to form a virtual population that reflects the baseline characteristics of the clinical cohort . here 
 we introduce a novel technique to efficiently generate virtual patients and , from this ensemble , demonstrate how to select a virtual population that matches the observed data without the need for weighting . 
 this approach improves confidence in model predictions by mitigating the risk that spurious virtual patients become overrepresented in virtual populations . 
METHODS
RESULTS
DISCUSSION
Supporting information