Application of different approaches to generate virtual patient populations for the quantitative systems pharmacology model of erythropoiesis

In a standard situation, a quantitative systems pharmacology model describes a “reference patient,” and the model parameters are fixed values allowing only the mean values to be described. However, the results of clinical trials include a description of variability in patients’ responses to a drug, which is typically expressed in terms of conventional statistical parameters, such as standard deviations (SDs) from mean values. Therefore, in this study, we propose and compare four different approaches: (1) Monte Carlo Markov Chain (MCMC); (2) model fitting to Monte Carlo sample; (3) population of clones; (4) stochastically bounded selection to generate virtual patient populations based on experimentally measured mean data and SDs. We applied these approaches to generate virtual patient populations in the QSP model of erythropoiesis. According to the results of our research, stochastically bounded selection showed slightly better results than the other three methods as it allowed the description of any number of patients from clinical trials and could be applied in the case of complex models with a large number of variable parameters.


Introduction
Drug research and development is a complex, multilayered, and long-lasting process that includes several cycles of data generation and interpretation [1]. To decrease the number of cycles and increase the confidence of the process, new techniques contributing to both data generation and interpretation have been continuously developed. Quantitative systems pharmacology (QSP) modeling [2][3][4] is one such technique allowing the accumulation of both publicly available and proprietary data measured in vitro, ex vivo, and in vivo. QSP modeling enables us to understand more about the pathogenesis of diseases and mechanism of action of drugs, thereby increasing the efficacy of drug development [5][6][7].
The conventional approach of QSP modeling includes fitting a model output to a series of mean data values. As a result, the parameters of the model represent fixed numbers enabling the description of mean data. However, the results of clinical trials include a description of variability in patients' responses to a drug, which is typically expressed in terms of conventional statistical parameters such as standard deviations (SDs) from mean values. To allow a QSP model to reproduce the variability observed in response to a drug administration, a virtual patient (VP) population is usually generated and applied [8].
Currently, there are several approaches to generate a virtual patient population. The most commonly used is the MAPEL algorithm [7]. The main idea of the MAPEL algorithm is to introduce the so-called mechanistic axes, define different virtual patient types in the form of coordinates in these axes, and find optimal weights for each patient type to describe the empirical distribution [9]. Three other approaches presented by Rieger et al. [10] are: (1) nested simulated annealing, (2) genetic algorithm, and (3) Metropolis-Hastings. However, the virtual populations obtained using these approaches are typically large, whereas the number of patients in real clinical trials is usually relatively small. Indeed, to obtain the best result on the basis of these algorithms, it is necessary to have rather big sample size (approximately about 1000 virtual patients). On the other hand, one could propose to select a small subset from the generated large populations. Such kind of selection cannot be done in the case of MAPEL algorithm, as the desired distribution properties are guaranteed only when the needed portions of different VP types is preserved in the population. However, small population cannot be characterized in terms of the fractions of different VP types. Three other algorithms listed above allow to perform such kind of selection. However, the exact implementation of the selecting algorithm and comparison of such a selection with the experimental results appears to be a promising direction for further research. Moreover, such methodology requires generation large VP population as initial step. Therefore, the approach will be computationally inefficient. This consideration is highly important in the case of large realistic models. Therefore, the additional requirement for the approaches developed in the present study is the ability to reproduce a population of virtual patients having the same size as that in the clinical trials.
Summarizing the arguments listed above, we can formulate the three key requirements for the approaches under development: (1) good correspondence with the clinical data; (2) ability to reproduce the population of any size (including the small ones); (3) computational efficiency enabling utilization for the cases of models with realistic large size.
We use these properties as the criteria for the comparison of the approaches proposed in the study. In particular, we formulate limitations, pros and cons of the methods in terms how well do they deal with these three tasks.
In this study, we propose and compare four different approaches to generate virtual patient populations based on experimentally measured mean data and statistics, namely, (1) Monte Carlo Markov chain (MCMC), (2) model fitting to Monte Carlo samples, (3) population of clones, and (4) stochastically bounded selection. We applied these approaches to generate virtual patient populations in the QSP model of erythropoiesis [11]. Our task was to create a sample of virtual patients of the same size as that in clinical trials and describe the experimental data. In this study we have compared these four different techniques we proposed to generate VP populations based on experimentally measured mean data and SD.
The remainder of this paper is organized as described herein. In the ''Methods'' section, we present a qualitative description of the erythropoiesis QSP model, the experimental data, and the algorithms of each of the four approaches. In the ''Results'' section, we describe the performance of the approaches with respect to the populations of different sizes. Finally, the ''Discussion'' section describes the advantages and disadvantages of the approaches and our future goals.

Description of model and set of experimental data
A model of normal human erythropoiesis 1 was constructed to describe the cell dynamics accompanying the conversion of bone marrow hematopoietic stem cells (HSCs) to circulating erythrocytes (ERY), including the following intermediate stages: multi-potent progenitor (MPP), common myeloid progenitor (CMP), megakaryocyte-erythroid progenitor (MEP), burst-forming unit-erythroid (BFU), colony-forming unit-erythroid (CFU), pre-proerythroblasts (Prepro), proerythroblasts (Pro), early (Basoe) and late (Basol) basophilic erythroblasts, polychromatic erythroblasts (Poly), orthochromatic erythroblasts (Ortho), and reticulocytes in bone marrow (Ret_bm) and plasma (Ret_pl). These cellular species undergo self-renewal, differentiation, proliferation, migration from the bone marrow into circulation, and apoptosis. Essential growth factors, such as the stem cell factor (SCF), erythropoietin (EPO), and interleukin-3 (IL-3) regulate cell dynamics. The explicit binding of EPO to cell-surface receptors results in increased hemoglobin production, which in turn stimulates the degradation of hypoxia inducible factor 1 alpha (HIF) in the kidneys, thereby forming a negative feedback loop for EPO production in the kidney. Clinically relevant outputs represent levels of plasma EPO, ERY, Ret_pl, and hemoglobin as a function of ERY and RET_pl (Fig. 1).
The model was calibrated against published in vitro/ in vivo data, including data for cell expansion under growth factors and cytokine exposure in vitro, flow cytometry cell counting of bone marrow aspirates, and clinical EPO administration in healthy volunteers [12] 2 . The model was also validated with respect to independent data sets for whole blood donation that was not used for parametrization to demonstrate the predictive power of the proposed model of erythropoiesis. For the purpose of our analysis, we used a dataset in the form of a time series describing plasma reticulocyte count in response to a single dose of EPO. EPO was administered to five healthy subjects, and the reticulocyte count was measured at 16 time points. Experimental data is given in the form of a series of mean values, m ¼ ðm 1 ; . . .; m T Þ and a series of SDs, sd ¼ ðsd 1 ; . . .; sd T Þ, where T ¼ 16.

Definition of terms used in the manuscript
The term ''control function'' refers to one of the model outputs which corresponds to particular experimental data. For example, a model variable (cell or cytokine concentration) or readout (AUC of a particular variable) can be considered as a control function if the corresponding dataset describing the variable (or readout) is available. In this paper, we considered the normalized reticulocyte count as a control function. This control function represents the reticulocyte count simulated after EPO administration normalized to the baseline reticulocyte count corresponding to the model steady state before EPO administration: RET pl norm ¼ RET pl

RET pl base
Description of the methods to generate VP population As a preliminary step for the VP population generation procedure, the model was fitted to the mean values obtained from literature. By doing this, we generated the so-called ''reference patient''. A description of the mean experimental data with ''reference VP'' is presented in Fig. 2. After obtaining the ''reference patient'' we proceeded with the process of generating populations of VPs.
As mentioned above, the distribution of the control function should be described by the same number of VPs as specified in the clinical trial. The clinical trial dataset contains information on five healthy subjects. The reticulocyte count was measured at 16 time points. There are about 100 parameters in the model. They were either fitted or fixed at some values on the basis of available experimental data. Thirty-nine parameters of the erythropoiesis model were chosen to account for variability in the observed clinical data. These 39 parameters were chosen on the basis of empiric considerations and included several groups: (i) rate constants of the processes downstream EPO effect entry point, so parameters related to early progenitors before BFU do not include into variability analysis, (ii) patient characteristics (height and weight) which influence EPO pharmacokinetics, (iii) parameters related to EPO synthesis and degradation, (iv) parameters regulated EPOdependent erythroid progenitors' differentiation, proliferation and apoptosis. The parameters and ranges of these variations are described in Table 1. Choice of min and max values for each parameter was based on the results of large series of numerical simulations of the model. These values, therefore, reflect the ranges of the parameters allowing to avoid irregular model behavior. To implement our approaches, we used either the continuous uniform distribution for the ranges indicated in Table 1 (Approach 4) or the truncated normal distribution for the parameters with l i ; r i ð Þ; i ¼ 1; . . .; 39, where l i is fixed on the level fitted in the model for the ''reference patient'', r i ¼ l i 2 , which is bounded by zero in the left and there are no bounds in the right (Approach 1). Therefore, our task was to create a sample of five virtual patients to describe the experimental data and compare the proposed approaches. To obtain a more precise description of the performance of the approaches, we conducted additional calculations for a relatively large population size of 207 patients. The comparison of the methodologies is presented for both the cases-one with a small population size (5 patients) and another with a relatively large population (207 patients).
The main objective of all the approaches presented in this study is to propose a methodology for generating a VP population, which describes the data control function distribution characteristics (mean and SD). Each of the approaches has specific advantages, disadvantages, and limitations. A comparison is provided in the Results section. All the approaches described in the study were implemented in R language, the source code can be accessed from the following GitHub public repository: https://github.com/insysbio/vp-erythropoiesis Fig. 2 Fitting of the control function (normalized reticulocyte count) to mean value in experimental data [12]. The parameter values resulting from this procedure correspond to those of the so called «reference patient» Subsequently, we assumed that the prior distributions of these parameters are described by the truncated normal distribution for the parameters l i ; r i ð Þ; i ¼ 1; . . .; M, where l i is fixed on the level fitted in the model for the ''reference patient'', r i ¼ l i 2 . Let us denote the dependence between the vector of the control function and parameter values by y ¼ y p 1 ; . . .; p M ð Þ . In our case, y relates to the RET pl norm. The task was to generate a series of vectors (p (3) Calculate the likelihood L ð1Þ for y ð1Þ with respect to truncated normal distributions with parameters m and sd using the following formula: M ) with respect to the independent prior distributions of the parameters using the following expression: (5) The overall likelihood is the product L ð1Þ P ð1Þ : M ) (go to step 6).
(2) else (if L ðpropÞ P ðpropÞ L ð1Þ P ð1Þ ) -Generate a random sample number u iþ1 from ½0; 1, where i is iteration number -Calculate the transition probability   (2) Set physiological ranges for the variable parameters.
(3) Fit the variable parameters (to fit the parameters we used ''L-BFGS-B'' method in R) to minimize the following objective function: The series of fitted parameter vectors (p

Approach 3: Population of clones
The main idea of this method is to propose a series of patient groups. Each group is described by its own control function value, i.e. patients withing groups are identical, whereas patients from different groups are described by different control function values and parameters vectors respectively. To obtain correspondence with the data we calculate the relative group sizes in terms of fractions of patients in the population relating to the particular group. By doing this we ensure that the overall control function distribution had the same values of m and sd as those in the experimental data.
In the first stage, we selected one time point for further analysis by calculus. We denoted this moment as t 0 . In particular, in the case of our model t 0 ¼ t 8 . To provide a robustness check we analyzed the results using different time points and find out that the choice of time point does not significantly affect the quality of the results. However, as the results are slightly dependent on the time point chosen, we recommend to create VPs using different time points (the method is computationally very fast, so it is not a problem to create such VPs) and choose the one providing the best fit over the whole time range. After that, for simplicity, we chose the number of groups L ¼ 5 (however, it can be of any size, not necessary 5) control function values from the control function variation range for this time point (y

Approach 4: Stochastically bounded selection
The main idea of this method was to select a sample of N (where, N is the number of patients in a clinical trial) elements from a relatively large set of random vectors of variable parameters p 1 ; . . .; p M ðM ¼ 39) such that they constitute a subpopulation of the initial VP population describing the results of a clinical trial. Similar to Approaches 1 and 2, let us denote the dependence between the vector of the control function and parameter values by y ¼ y p 1 ; . . .; p M ð Þ . First, in case of the example considered below, we generated a large series of vectors of variable parameters: (p  Table 1. The value of S depends on the model complexity and the available computing power. The larger this value is, the more accurate the estimations obtained are; however, the time needed to obtain the initial random sample is larger too. After the vectors for variable parameters were generated, we simulated model with this vectors and as the result of these simulations we obtained a large series of vectors of control functions: y : (1) Find the minimal Y t min and maximal Y t max value of the generated control function sample at each time point t 2 1; . . .; T ½ . (2) Generate N random values for y from truncated lognormal distribution [13] with the following characteristics: (1) m ¼ ðm 1 ; . . .; m T Þ, sd ¼ sd 1 ; . . .; sd T ð Þ , and T ¼ 16 are the same as in experimental data (2) distribution is concentrated in the interval ½Y t min ; Y t max Therefore, we obtain the following sets of vectors of the control function:

Global sensitivity analysis
We performed a global sensitivity analysis to assess the impact of each variable parameter on the dynamics of the model. We used the partial rank correlation coefficient method with Pearson correlation coefficients. A thorough description of this method has been provided in a paper by [14]. The partial rank correlation coefficient values were in the range of -1 to 1. Values close to zero indicate that the parameter was not sensitive. The results of the sensitivity analysis are provided in the Results section.

Results
In this section, we compare the aforementioned approaches for VP population generation. The comparison was made for a relatively small (5 patients) and large population size (207 patients). In the first case, the comparison is qualitative, whereas in the second case, additional rigorous criteria were also considered.

Comparison of the approaches in the case of five patients
As the target sample size was extremely small, the performance of the approaches could be compared by visual assessment only. We want to note, that Approach 3 was inapplicable in case of a small sample size. This follows from the fact that Approach 3 is based on the calculation of group sizes. This sizes are formulated in terms of fractions of patients in the population. When the population is small (i.e. in includes only a few patients) the description in terms of fractions is evidently inapplicable (imagine, what happens when we take 10% from five people). Figure 4 shows the results obtained using the 3 approaches in case of five patients. In Fig. 4a, we present the mean and SD values of the control function RET pl norm distribution predicted in Approach 1 at every time point and those for the experimental values. Similarly, in Fig. 4b, c, we show the characteristics of the predicted control function distributions and the corresponding experimental values at every time point for Approaches 2 and 4, respectively. Figure 4a shows the poor precision of the mean and SD estimates in Approach 1. This was not surprising, as a large sample is required to obtain reliable results in the MCMC approach. Better results were obtained in Approach 2 ( Fig. 4b) as it predicted the mean values relatively accurately; however, the SD estimates were too small. This flaw was overcome in Approach 4 (Fig. 4c), in which relatively good estimates of both mean and SD values were obtained, despite the fact that the sample size was extremely small. Therefore, Approach 4 can be considered as the preferred approach.

Comparison of the approaches in the case of 207 patients
Besides visual assessment, in the case of a relatively large sample size, the approaches can be compared more rigorously. To do this, we used the one-sample t-test and QQ plot. In the case of Approach 3, it is not possible to apply comparisons in terms of QQ plot and one-sample t-test,  Figure 5 shows the results obtained by the four approaches for 207 patients. In Fig. 5a-d, we present the results of Approaches 1, 2, 3, and 4, respectively. Similar to Fig. 4ac, Fig. 5a-d depict the mean and SD values of the predicted control function distributions and the corresponding experimental values at every time point. The figure clearly indicates that all the approaches accurately predict the mean values at every time point; however, relatively precise estimates of SD were only obtained in Approach 4.

QQ plot
Using a QQ plot, two distributions can be compared based on their respective quantile positions. Therefore, we compared the predicted distributions with the experimental distributions for each approach, except for Approach 3.  Figure 6 shows QQ plots for the 3 time points at which the predicted and experimental distributions coincide in Approach 1. Figures 7 and 8 illustrate the best time points for Approaches 2 and 4, respectively. Approach 4 provided the largest number of relatively accurate time points.

One sample t-test
Since the experimental data are given in the form of a series of mean, m ¼ ðm 1 ; . . .; m T Þ, and SD, sd ¼ ðsd 1 ; . . .; sd T Þ, values, we used one-sample t-test to assess their statistical significance. The p-values generated for each approach, except Approach 3, are shown in Table 2. The data in the table shows that in case of Approach 1, 2, and 4, statistical significance is achieved at 3, 2, and 4 time points, respectively (highlighted by italics in Table 2). Therefore, with respect to this criterion, Approach 4 is the preferred approach.

Analysis of distributions obtained for variable parameters
This section presents the analysis of the distributions obtained for the variable parameters (final VP sample) in each approach, except for Approach 3, because the final sample, obtained in this approach, consists of only five unique types of virtual patients. In Fig. 9 the results of the global sensitivity analysis for all 39 model parameters are presented. This analysis was conducted using the partial rank correlation coefficient method [14]. From the list of parameters with the highest sensitivity we selected five most important in view of Erythropoiesis process. These parameters are: k_dif_ret_ery_pl, k_apo_ery_pl, IC50_a-po_epo_cfu_e, EC50_pro_epo_cfu_e, and gamma_hif1a.
Sensitivity analysis revealed that the response to EPO is sensitive to its effect on apoptosis, CFU-erythroid cell proliferation, and EPO synthesis-related processes. It has been noted that CFU-erythroid cells are absolutely dependent on EPO as they carry the highest EPO receptor density; therefore, they are the most sensitive to EPO among other erythroid progenitors. Parameters k_dif_ret_ery_pl and k_apo_ery_pl directly influence erythrocyte count, which contributes to hemoglobin much more than reticulocytes. Hemoglobin, being an oxygen carrier in turn affects HIF, due to which EPO production is regulated and EPO levels return to the steady-state. Parameter gamma_-hif1a as a power-law indicator of the rate of HIF-dependent synthesis of EPO mRNA strengthens the feedback loop between EPO level and erythrocyte count non-linearly.
The sample distributions for the top five sensitive parameters in the results of each approach are shown in Fig. 10. Figure 10a depicts the results of Approach 1. The distributions shown here have unimodal shapes, which is not surprising as prior distributions in the MCMC algorithm are truncated normal distributions. Figure 10b illustrates the parameter distributions in Approach 2. As we used a log-normal distribution for the generation of the pseudo experimental data, the shapes of the resulting parameter distributions inherit some specific features of this distribution type. However, the existence of the upper and lower bounds resulted in the appearance of the second mode. Figure 10c shows the resulting parameter distributions for Approach 4. In this case, the distribution shapes were more uniform, which is explained by the uniform initial parameter distributions in this approach.

Discussion
In this section, we discuss the advantages and disadvantages of the approaches and describe probable directions of future research. This study aimed to compare different techniques to generate VP populations based on experimentally measured mean data and statistics. We proposed, implemented, and compared four different approaches to generate VP populations. These are (1) MCMC, (2) model fitting to Monte Carlo sample, (3) population of clones, and (4) stochastically bounded selection. For the purpose of the comparison of these approaches, we used the QSP model of erythropoiesis. To generate the VP populations, 39 parameters of the model were chosen. To ensure that the appropriate VP populations were selected, clinical data describing the changes following erythropoietin administration in the reticulocyte count in the blood of five patients were used. The main purpose of this study was to obtain the distribution of the control function described by the same number of virtual patients as those in the clinical trial. However, the sample of five patients was too small to represent a reliable quantitative comparison. Therefore, a larger VP population of 207 patients was generated. The comparison showed that all the approaches proposed are potentially applicable in different cases; therefore, we have discussed their advantages and disadvantages herein.  The main disadvantage of Approach 1 is that it works for large samples only, i.e., if one needs to generate a small population of patients, the results of this technique will be unsatisfactory. Another disadvantage of this approach is that an appropriate SD estimate can be obtained only in cases where its value is relatively small. The algorithm of this approach also calls for the simulation of the model at each step. Therefore, if the model and number of patients are both large, this approach will be very time consuming. The main advantage of Approach 2 is its potential accuracy as using this, the final population of patients is obtained by fitting the variable parameters to the generated pseudo experimental data having exactly the same mean and SD. However, this approach has a series of limitations, such as (1) only a limited number of parameters can be accurately fitted, (2) one can obtain precise fitting only in the case of small models, and (3) the physiological ranges of variable (fitted) parameters should be known.
The main advantage of Approach 3 is its computational simplicity. Moreover, regardless of the model complexity, a population having exactly the same mean and SD as in the data can be obtained. However, in the VP populations obtained using this approach, there is no intergroup variability, and the VP population size should be rather large.
Recall the three comparison criteria formulated in the introduction. To make the comparison explicit we condensed the described above approaches pros, cons and limitations into these three ''dimensions''. These results are presented in Table 3.
From Table 3 it is evident that Approach 4 appears to be the most preferred with regard to the initially formulated requirements. As can be seen from the results presented in this study, Approach 4 appears to be the most universal as it allows the description of any number of patients from clinical trials, and it can be applied in the case of complex models with a large number of variable parameters.
However, one of the main disadvantages of Approaches 3 and 4 is that both the methods are based on the selection of virtual patients at a single time point only. On the one hand, this fact simplifies the methodology. However, most likely for this reason we obtain observable poor match for series of the time points. This limitation can be overcome through the automatic adjustment of other time points and by not using their experimental data (mean and SD for each time point). For the model used in this study, we consider such simplification to be acceptable, since (1) a comparison of the results of Approaches 3 and 4 with those of multipoint Approaches 1 and 2 shows that the results were even better; (2) when we work with small population sizes we are mostly interested in qualitative trends and not exact matchings.
Thereby, there are a series of natural improvements for Approach 4. Their introduction into the methodology will make Approach 4 even more universal and accurate. In the near future, we plan to add the following updates: (1) generation of VP populations describing the distributions of several control functions simultaneously, (2) introduction of the restrictions on biomarkers, and (3) selection of VP populations using all time points.
Based on the results of the global sensitivity analysis, we selected 5 out of 39 parameters, which are (1) described by one of the highest sensitivity and (2) are of interest from the Erythropoiesis process point of view. We use these parameters to analyze the resulting distributions of variable parameters in details. The histograms of these distributions show that the shapes of the parameter distributions significantly differed amongst the approaches. The boundaries of the physiological ranges of the variable parameters are listed in Table 2. A comparison of the final distributions of the variable parameters shown in Fig. 9 with the physiological ranges of the variable parameters (as shown in Table 2) revealed that the final distributions of the parameters mostly belong to the corresponding physiological boundaries.

Conclusions
The paper presents four different approaches to generate virtual patient populations based on experimentally measured statistics. Comparison of the results of application of the four approaches to populations of two different sizes shows that the stochastically bounded selection is the preferable one in both cases. Moreover, this approach appears to have the smallest computational complexity and, therefore, it can be applied even to the high dimensional models with large number of variable parameters.
Acknowledgements We are grateful to the two anonymous reviewers for their insightful comments and suggestions which allowed us to significantly improve the results presentation.
Author contributions GK: created methodology and performed all data analysis; OD and GK: discussed and analyzed the results; AS and GL: created model of erythropoiesis. All authors drew the main conclusions and wrote the manuscript. Correspondence with the clinical data -± ± ± Allows to create small sized VPs -? -?