 Research
 Open Access
 Published:
Practical identifiability in the frame of nonlinear mixed effects models: the example of the in vitro erythropoiesis
BMC Bioinformatics volume 22, Article number: 478 (2021)
Abstract
Background
Nonlinear mixed effects models provide a way to mathematically describe experimental data involving a lot of interindividual heterogeneity. In order to assess their practical identifiability and estimate confidence intervals for their parameters, most mixed effects modelling programs use the Fisher Information Matrix. However, in complex nonlinear models, this approach can mask practical unidentifiabilities.
Results
Herein we rather propose a multistart approach, and use it to simplify our model by reducing the number of its parameters, in order to make it identifiable. Our model describes several cell populations involved in the in vitro differentiation of chicken erythroid progenitors grown in the same environment. Interindividual variability observed in cell population counts is explained by variations of the differentiation and proliferation rates between replicates of the experiment. Alternatively, we test a model with varying initial condition.
Conclusions
We conclude by relating experimental variability to precise and identifiable variations between the replicates of the experiment of some model parameters.
Background
Interindividual variability is ubiquitous in biology, from the fluctuations of molecular contents across populations of single cells [1], to the variations of physiologial parameters between whole organisms [2]. This variability has uncountable consequences, for instance at the scale of developmental [3], ecological or evolutionary processes [4, 5].
As a result, one often faces significant amounts of variations between replicates of the same biological experiment, which we will refer to as experimental variability. This variability can be taken into account by deterministic dynamical models of the biological system, as a random variation around its predicted behaviour [6, 7]. Such models thus disregard the fact that variability is inherent to the biological nature of the system under study.
Another difficulty that can arise from this approach is when some parameters of the model are unidentifiable. A parameter is said identifiable when a particular measurement of the model output (potentially affected by measurement error) is associated to a unique parameter value [6]. Otherwise it is unidentifiable. The model itself is said identifiable if all its parameters are identifiable, and unidentifiable if at least one of its parameters is unidentifiable.
More precisely, a model can be unidentifiable for several reasons. If some parameters are redundant, meaning that they can be varied together in such a manner that the model output is kept constant, they are called structurally unidentifiable. If the data quantity (sample size) or quality (measurement error) are insufficient to precisely infer some parameter values, these parameters are said practically unidentifiable. It should be noted that while all practically identifiable parameters of a model are also structurally identifiable, the converse is not necessarily true (see for instance the recent review on identifiabilty by Wieland et al. [8]). For this reason, the focus of this paper is on practical identifiability, and unless stated otherwise we will be referring to the practical identifiability of model parameters.
When the parameters of a model have a precise physical or biological interpretation, it can be tempting to use their estimates to formulate predictions about the system. However, as these estimates are not uniquely determined in unidentifiable models, an unidentifiable model should never be used for predictive purposes [9].
In order to interpret experimental heterogeneity, we propose to use nonlinear Mixed Effects Models (MEM). In particular, we are interested in the identifiability of such models. MEM work by applying the same mathematical model to all the individuals of the population, with different parameter values for each individual, and thus have been used in a variety of fields involving interindividual variability [10]. This approach allows to assign different levels of variability for each parameter by making the distinction between population parameters, that are the means and variances of the parameter values across the whole population, and individual parameters, that are the precise parameter values assigned to each individual.
In the context of experimental variability, one might for instance consider all replicates of the experiment as individuals coming from the same population (the theoretical population of all the possible outcomes of the experiment). Assigning different parameter values in a dynamic model for each individual (i.e. each replicate of the experiment) would then allow to assess which parameters of the model are mostly affected by experimental variability (i.e. which parameters are the most variable between individuals). The question that naturally arises is whether or not such a model would be identifiable.
In general, one argument in favour of the use of MEM is the fact that using the population distribution as a prior can help with the estimation of the individual parameters, thus improving their practical identifiability [11]. This rationale is based upon the fact that most MEM calibration methods estimate the population parameters in a first step, using the data from the whole population, and then estimate individual values for every parameter in a bayesian way, an approach referred to as Empirical Bayesian Estimation [12].
However, this first intuitive argument on parameter identifiability in MEM is somewhat challenged by another consideration: the intricate definition and estimation of population and individual parameters might in fact complicate the assessment, and even the definition, of parameter identifiability in MEM [13]. As a consequence, the identifiability of MEM parameters is of critical importance to their widespread applications, and should not be neglected [9, 14]. In general, practical parameter identifiability depends on the precise definition of the model parameters (in the case of MEM, that is the definition of the distributions of the individual parameters across the population), together with the quantity and quality of the experimental data used for calibration [6]. Several kinds of approaches for structural identifiability analysis, that were developped for models without mixed effects, have been adapted in a mixed effect context [15, 16]. Regarding practical identifiability analysis, two different kinds of empirical approaches are reported for MEM [13], though with a lot of potential refinements.
First, the Fisher Information Matrix (FIM), which is computed from the Hessian of the likelihood, estimated at the optimal parameter set, allows for a quadratic approximation of the likelihood surface near its optimum. This in turn, allows to infer confidence intervals for any parameter at any level, provided that the FIM is nonsingular [17, 18]. In this setting, a singular FIM indicates that some parameters are structurally unidentifiable. Conversely, a nearsingular FIM could indicate that some parameters are practically unidentifiable. However, the quadratic approximation of the likelihood surface might mask practically unidentifiable parameters in the case of nonlinear, partially observed models. In extreme cases, the FIM can even make some practically unidentifiable parameters appear as identifiable [6]. As a consequence, the FIM is inadequate to study practical identifiability. Since our focus is on the practical identifiability of MEM, we will not use the FIM to assess the identifiability of our models. Other methods have been developed for models without mixed effects, such as the profile likelihood [6], but to our knowledge they have not yet been implemented in any of the existing software for MEM calibration. Given the widespread use of these software for the calibration of MEM in pharmacology and personalized medicine [19], it would be particularly interesting to be able to empirically study the practical identifiability of MEM, directly from the calibration software.
Secondly, one might run the estimation algorithm several times, using a sample of initial guesses for the parameter values [13]. In that case, the convergence of the algorithm to a unique likelihood optimum, with different optimal parameter values, indicates that the parameters are unidentifiable. We refer to this approach as Initial Guess Sampling (IGS). It has also been termed the multistart approach, and the samples of estimated parameter values that it provides do not contain any information regarding the variance of the estimation or the confidence intervals of the parameters [20]. Since this method requires a potentially large sample of runs of the estimation algorithm, it is more costly in terms of computational power than a simple evaluation of the FIM eigenvalues. As a consequence, most approaches to the identifiability analysis of MEM rely on the FIM [10, 13, 14, 21,22,23], despite its proven inaccuracy at assessing practical identifiability [6].
Since the practical identifiability of a model depends on both the definition of the model and the data used to calibrate it, there are two broad classes of approaches for dealing with an unidentifiable model. On one hand, it is possible to increase the amount of data available for parameter estimation. For instance measuring an additional, previously unobserved quantity might remove structural unidentifiabilities. It is also possible to use the tools of experimental design [22, 23] to define a new set of more informative experiments, that can then be used to estimate new parameter values and assess their identifiability. However, experimental design approaches are approximate and there is no a priori proof that performing the optimal experiment will actually make the model identifiable [24, 25]. On the other hand, it is also possible to change the definition of the model parameters, in order to simplify the estimation for the other parameters [26, Section 10.2], an approach which we refer to as model reduction. For instance, reparameterizing the model in terms of the estimable parameter combinations would remove any structural unidentifiabilities, while potentially sacrificing the biological interpretation of these parameters. It is also possible to constrain the value of some unidentifiable parameters (for instance setting them to zero) in order to simplify the estimation task. But then, which criterion would allow us to choose which parameters to remove from the model?
This consideration is particularly important for MEM, as all individual parameters might not have the same variance across the population. Thus, the size of the sampled dataset (in terms of the number of individuals) is critical for the estimation of the population variances, since a sampling bias in a small dataset could mask the variance on some parameters. From the point of view of experimental design, the determination of the necessary sample size in order to guarantee a certain level of confidence on all parameters in MEM is a central question, which has already been covered to a certain extent [22, 27], using geometric features of the likelihood surface approximated from the FIM [28, section 10.5.3]. From the point of view of model reduction, it would be tempting to remove the unidentifiable variances by setting them to zero, potentially improving the estimation of the other parameters of the MEM without affecting the quality of the model fit to the data. In some cases however, adding a random effect or a covariate to a MEM could improve parameter identifiability [26, Section 5.1], as it might split out combinations of structurally unidentifiable parameters.
In this paper, we adress these questions using a MEM of the in vitro erythropoeisis that we adapt from a previous model, proven to relevantly reproduce the dynamics of single replicates of the experiment [29]. This MEM accounts for experimental variability by assigning different parameter values for proliferation and differentiation in each replicate of an identical experiment. We assess its identifiability using a multistart approach, based on extensive parameter estimations with the MEM calibration software Monolix [30]. Then, we reduce the model in order to make it identifiable, using the correlations between the estimated parameter values. Alternatively, we test whether or not the observed variations in the outcome of our experiment could be explained by variations in the initial condition of the experiment rather than variations of the differentiation and proliferation dynamics. Our final model associates different levels of variability for each dynamic parameter, which allows us to identify which features of the erythroid differentiation are the most variable from experiment to experiment. Moreover, this work proposes a multistart approach for MEM identifiability analysis, which appears as a promising alternative to the FIM.
Methods
T2EC cell culture
The experimental setting from which all the data used in this study were obtained consists in a culture of 25,000 chicken erythroid progenitors called T2EC that were extracted from the bone marrow of 19 daysold embryos of SPAFAS white leghorn chickens (INRA, Tours, France). They may either be maintained in a proliferative state or induced to differentiate into mature erythrocytes depending on the medium in which they are grown, as previously described [29, 31,32,33].
In the selfrenewal medium (referred to as the LM1 medium) the progenitors selfrenew, and undergo successive rounds of division. Its composition is given in Table S1 (Additional file 1). Cell population growth was evaluated by counting living cells in a \({30}\,{\upmu \hbox {L}}\) sample of the 1mL culture using a Malassez cell and Trypan blue staining (SIGMA), which specifically dyes dead cells, each 24h after the initial sowing of 25,000 cells in the culture, as previously described [29, 31,32,33].
T2EC can be induced to differentiate by removing the LM1 medium and placing cells into 1mL of the differentiation medium, referred to as DM17. Its composition is given in Table S1 (Additional file 1). Upon the switching of culture medium, a fraction of the progenitors undergoes differentiation and becomes erythrocytes. The culture thus becomes a mixture of differentiated and undifferentiated cells, with some keeping proliferating. Cell population differentiation was evaluated by counting differentiated cells in a \({30}\,{\upmu \hbox {L}}\) sample of the culture using a counting cell and benzidine (SIGMA) staining which stains haemoglobin in blue. A parallel staining with trypan blue still gives access to the overall numbers of living cells, as previously described [29, 31,32,33].
Consequently, the data available from this experiment are the absolute numbers of differentiated cells, as well as the total number of living cells (which comprises both selfrenewing and differentiated cells) each 24h after the initial sowing of 25,000 cells in the culture. The data presented on Fig. 1 are the total number of living cells in the culture, and the fraction of differentiated cells in 7 independent replicates of the experiment.
Modelling framework
A Mixed Effects Model (MEM) is defined as the combination of three components. The structural model describes the dynamic process at play in each individual. The parameter model, or individual model, describes how the parameters of the structural model vary from individual to individual. Finally, the observation model, or error model, describes how the predicted outcome of the model for each individual differs from the observation.
Dynamic model
The SCB model, that we previously described [29], faithfully reproduces the dynamics of T2EC proliferation and differentiation by accounting for 3 cellular states (Fig. 2). The selfrenewing state S describes the state of cells in the LM1 medium, where they can only proliferate or die. The differentiated state B (which stands for Benzidinepositive) describes mature erythrocytes in the DM17 medium. Lastly, in the committed state C, cells have not finished differentiating, but cannot go back to selfrenewal anymore, so that they are committed to differentiation. The dynamics of these three compartments are given by the equations:
It is characterized by the set \((\rho _S, \delta _{SC}, \rho _C, \delta _{CB}, \rho _B)\) of five dynamic (or kinetic) parameters, where \(\rho _i\) is the net growth rate of compartment i, which might be positive or negative depending on the net balance between cell proliferation and cell death, and \(\delta _{ij}\) is the differentiation rate of compartment i into compartment j, which must be positive.
Moreover, it should be noted that differential system (1) is fully linear, and that its matrix is lowertriangular, which makes it easily solvable analytically. Its simulation is thus very fast. The detail of the analytical solutions to this system is given as supplementary material in [29].
Not all variables in the models can be measured through the experiments that we presented, and we only have access to three observables of the system: the number of living cells in LM1 (which we denote as S since there are only selfrenewing cells in LM1), the number T of living cells in DM17, and the number B of differentiated cells in DM17 (it is null in LM1). Unless stated otherwise, we will always consider that the initial condition is fixed by the experimentalist so that the initial state of the observables is: \((S_0, T_0, B_0) = (25{,}000, 25{,}000, 0)\).
The SCB model was selected among others models with different structures, based on its ability to reproduce the kinetics of erythroid differentiation. It should also be noted that the original description of the SCB model –with slightly different constraints on parameter values compared to the context of this manuscript– was proven to be structurally and practically identifiable, using the profile likelihood approach [29].
Parameter model
In order to describe interindividual variability in the parameter values of the SCB model, we consider at first that the five kinetic parameters can vary between every individual. The two differentiation rates \(\delta _{SC}\) and \(\delta _{CB}\) must be positive, and the net proliferation rates \(\rho _S\), \(\rho _C\) and \(\rho _B\) can be positive or negative. In order to respect these bounds on the individual parameter values, we use a combination of normal and lognormal distributions across the population:
This parameter model has 5 fixed effects (\(\rho _S^{pop}\), \(\delta _{SC}^{pop}\), \(\rho _C^{pop}\), \(\delta _{CB}^{pop}\) and \(\rho _B^{pop}\)), which quantify the average behaviour of the population, and 5 random effects (associated with the standard deviations \(\omega _{\rho _S},\) \(\omega _{\delta _{SC}}\), \(\omega _{\rho _C}\), \(\omega _{\delta _{CB}}\) and \(\omega _{\rho _B}\)).
Error model
In order to account for experimental errors in the measurement of the observables, MEM include an error model, or observational model, which describes the statistical fluctuation of the model prediction around the observation. We previously demonstrated that the proportional error model is the best to describe the prediction error of the SCB model [29]. It is defined by:
where \(y_{i,j,k}\) marks the measurement of the \(i^{th}\) observable, at the \(j^{th}\) timepoint, on the \(k^{th}\) individual, and \(f_i\) marks the model prediction for the \(i^{th}\) observable from System (1), which depends on time t, the initial condition \(y_0\), and the individual parameters \(\theta _k\). Finally, \(b_i\) denotes the proportional error parameter for the \(i^{th}\) observable, which quantifies the standard deviation of the prediction error, and \(\varepsilon _{i,j,k}\) is the individual weighted residual of the model for individual k, at time \(t_j\), for observable i. The proportional error model introduces one additional parameter \(b_i\) for each observable, resulting in three error parameters for our SCB model.
Together with the dynamic model of System (1) and the parameter model of Eq. (2), this error model defines our first version of a MEM for the in vitro erythropoiesis. Since all other MEM in this manuscript will have the same dynamic and error components, we will omit them from now on, and will define each MEM by its parameter equation only, such as Eq. (2).
Parameter estimation
We used the Stochastic Approximation version of the ExpectationMaximization (SAEM) algorithm [34] implemented in Monolix [30] to estimate the parameters of our MEM (Additional file 1: Table S2). To avoid potential local likelihood optima and ensure the convergence of the algorithm to the global optimum, we performed the estimation 50 times using the Monolix Convergence Assessment tool, with independent uniformly sampled initial guesses for the parameter values (Additional file 1: Table S3). For the fixed effects, we sample the initial guess in an arbitrarily large interval (Additional file 1: Table S3). Thus most initial guesses will be wrong, and potentially far from the true value. To ensure convergence in these conditions, we set a high initial variance and error parameter values (Additional file 1: Table S3). This ensures that the first step of SAEM can sample individuals across all the parameter space, which will allow for a subsequent improvement of the estimates both for the population averages and variances. This is a multistart approach that we refer to as Initial Guess Sampling.
Model selection
In order to select which SAEM runs converged to the global optimum, we used Akaike’s weights [35]:
where \(w_i\) is the Akaike’s weight of the ith run, \(AIC_i\) is its Akaike’s Information Criterion and R is the number of competing models. The Akaike’s weight of a given model in a given set of models can be seen as the probability that it is the best one among the set [35]. In this setting, selecting the best models of a set of models means computing their Akaike’s weights, sorting them, and keeping only the models whose weights add up to a significance probability (in this manuscript, 95%).
In MEM, one might either choose to use the marginal or the conditional AIC depending on the context [36]. They differ by the corrective term that they introduce in the likelihood. However, we will essentially use the AIC and the corresponding Akaike’s weights for selecting models with the same structure and different likelihoods, such as the 50 runs of SAEM that we perform on the same model to assess its convergence. For this reason, the choice of mAIC or cAIC is not relevant to our study, and we will use the marginal AIC (\(AIC = 2\,log( \widehat{{\mathcal {L}}} )+2K,\) where \(\widehat{{\mathcal {L}}}\) is the maximum likelihood and K is the number of population parameters), that is computed by Monolix by default [10].
In order to select models with different structures, i.e. models that would differ by the definition of their fixed or random effects, we use the BIC that has been derived for MEM [37, 38].
Identifiability analysis
Population parameters
We use an approach based on repeated parameter estimations, starting from different initial guesses, to empirically assess the identifiability of our MEM. In this case, convergence to different parameter values with the same likelihood indicates unidentifiability [13].
This approach has also been termed the multistart approach, for instance by [20]. We will refer to the multistart approach as Initial Guess Sampling (IGS) in this manuscript.
Our approach to IGS is the following:

1.
We perform a random sampling of the initial parameter guesses and run SAEM for each of the initial guesses, using the Monolix Convergence Assessment tool. This provides us with a sample of optimal parameter values.

2.
We test the convergence of the SAEM runs: we only want to consider the runs which reached the global optimum. To this end, we use a selection criterion (\(w_{AIC}\)) to keep only the runs that converged to the lowest likelihood values.

3.
We compare the parameter values of these convergent runs. If they are different, then the model is unidentifiable, as several different parameter values give the same likelihood.
It should be noted that since multistart approaches do not provide any information on the parameters estimation error or confidence intervals [20], this approach do not allow for any statistical testing of parameter identifiability. The distribution of estimated values can rather be used to design a diagnostic plot of population parameter identifiability, showing which parameters vary the most between the convergent runs, and are thus the most poorly estimated. We propose to visualize the distributions of estimated values as a boxplot normalized by their median in order to display this information.
Individual parameters
Individual parameters are estimated using an empirical bayesian approach, where the population distribution of the parameters serves as a prior balanced by the individual data. In the case of unidentifiable individual parameters, the experimental data do not provide enough information to determine them precisely. Then, the posterior distribution is very close to the prior, resulting in individual parameters being estimated as their population mean.
More precisely, this principle that the posterior matches the prior for unidentifiable parameters holds under two conditions [39]: first, the prior distributions must be independent, second, the parameter space must be a product space. When one of these two conditions is not met, it is possible that the posterior distribution will differ significantly from the prior even for unidentifiable parameters [39].
In the case of our model, the prior distributions of the individual parameters, which are the population distributions defined in System (2), are independent, since the variancecovariance matrix of the random effects is diagonal. Moreover, each individual parameter is either real (net selfrenewal rates) or positive (differentiation rates), and thus the individual parameter space is a product space, namely \(\left( \rho _S, \rho _C, \rho _B, \delta _{SC}, \delta _{CB} \right) \in {\mathbb {R}}^3 \times {\mathbb {R}}_+^2\)
Consequently, we can assess the identifiability of the individual parameters by measuring the overlap between the prior and posterior distributions. This phenomenon is summarized by a scalar criterion called the \(\eta\)shrinkage [12, 40]:
where \(std(\eta _k)\) is the standard deviation of the estimated individual random effects in the population, and \(\omega _\theta\) is their theoretical standard deviation. In the case where information about a parameter is insufficient, the random effects on this parameter shrink toward 0 in the population, and thus \(s_\eta\) increases. Eq. (4) also implies that shrinkage values vary from parameter to parameter, and that some parameter might be more poorly characterized than the others.
Simulation studies have shown that shrinkage has a variety of effects on the model diagnostics, starting from \(30\%\) shrinkage [12]. High shrinkage values affect the correlations between random effects and covariates, as well as the correlations between the random effects themselves. It can also affect the detection of structural model specification.
In this paper, we will use the \(30\%\) limit introduced in [12] as a rule of thumb to consider that the individual parameters are well estimated.
Results
The model is unidentifiable
We estimated the parameter values of Model (2) by using our multistart approach. The distribution of estimated likelihood values over the 50 runs of SAEM is displayed on Fig. 3A, showing small variations between the estimated loglikelihood values. Among these 50 runs of SAEM, the 45 associated to the lowest \(2\,log( \widehat{{\mathcal {L}}} )\) add up to 95% of Akaike’s weights (Fig. 3B). We thus focus on the outcome of these 45 runs in the following.
The distributions of estimated population parameter values are represented in Fig. 3C. The fixed effect \(\rho _S^{pop}\), the corresponding variance \(\omega _{\rho _S}\) and the three error parameters \(b_1\), \(b_2\) and \(b_3\) are estimated with the smallest variance. For any of the other 8 parameters, the estimated values display more or less variability. For these parameters, the estimation is less reliable.
Thus, we conclude that the fixed effects \(\delta _{SC}^{pop}\), \(\rho _C^{pop}\), \(\delta _{CB}^{pop}\) and \(\rho _B^{pop}\) are unidentifiable, as well as the corresponding variances (a total of 8 unidentifiable population parameters).
The shrinkage of the individual random effects is displayed on Fig. 3D. The values of shrinkage for \(\delta _{SC}\) \(\rho _C\), \(\delta _{CB}\) and \(\rho _B\) range from 20 to 90% depending on the run, which indicates a clear discordance between the population distribution of the parameters and the actual distribution of the individual parameters. We thus conclude that the individual data are not informative enough to estimate all random effects for each individual.
As a consequence, it appears that Model (2) is unidentifiable at the population level as well as at the individual level.
A reduction approach for MEM
Fixed effects: parameter correlations
Figure 4A displays the value of Spearman’s \(\rho ^2\), which measures the nonlinear correlation between two variables, for each pair of the 8 unidentifiable population parameters of Model (2). There is a high correlation between \(\delta _{SC}^{pop}\), \(\rho _C^{pop}\) and \(\delta _{CB}^{pop}\) across the runs, which is represented in Fig. 4B, C.
These results show that the optimal values of \(\delta _{SC}^{pop}\), \(\rho _C^{pop}\) and \(\delta _{CB}^{pop}\) are strongly correlated in the range of values of Fig. 4. This range corresponds to the range of estimated values in the 50 SAEM runs. The correlations suggest that if we would replace two of these parameters by their expression as a function of the third one, we would also reduce the number of population parameters to estimate, and still allow them to reach their optimal value.
We found the following expressions for the correlations:
and
We thus conclude that if we replace \(\delta _{SC}^{pop}\) and \(\delta _{CB}^{pop}\) by their expression as a function of \(\rho _C^{pop}\) in Model (2), it might help the estimation. Yet, such a reduction might affect the convergence of SAEM because the correlation might not hold outside of the parameter range of Fig. 4.
Replacing \(\delta _{SC}^{pop}\) and \(\delta _{CB}^{pop}\) in Model (2) by their expression as a function of \(\rho _C^{pop}\) , we obtain the following reduced model:
Following the same approach as for Model (2), we ran the SAEM algorithm on this model 50 times using uniformly sampled initial guesses for the population parameters. The resulting optimal likelihood distribution is displayed on Figure S1A (Additional file 1). Most of the runs reached the same likelihood optimum as with Model (2) (Fig. 3A), but 10 of them found higher likelihood values. In the case of Model (7), Akaike’s weights select only 36 runs as the best ones (Additional file 1: Figure S1B), that we will consider as the runs that reached the global likelihood optimum.
The parameter values estimated in these 36 runs are displayed on Fig. 4D. First, it shows that the reduction of the model did not affect the accuracy of the estimation for the five parameters that were identifiable in Model (2). Then, the population parameters \(\rho _C^{pop}\) and \(\rho _B^{pop}\) are estimated more precisely in the reduced model (7) than in Model (2). However the three standard deviations \(\omega _{\rho _{C}}\), \(\omega _{\delta _{CB}}\) and \(\omega _{\rho _B}\) are still estimated with some variability.
Since \(\omega _{\rho _{C}}\), \(\omega _{\delta _{CB}}\) and \(\omega _{\rho _B}\) define the distributions of three random effects, their unidenfiability might indicate an overparameterization of the random effects. We investigate this using the \(\eta\)shrinkage of the individual random effects in the next section.
Random effects: shrinkage
We measured the \(\eta\)shrinkage in the convergent runs of Model (7). The average of the shrinkage values for each parameter are displayed in Table 1. Their distributions over the convergent runs are displayed in Figure S2 (Additional file 1). These values confirm that the individual parameters of Model (7) are unidentifiable.
These results indicate the individual data that we presented in Fig. 1 are insufficient to estimate five parameters per individual precisely. Indeed, our dataset only comprises 7 individuals. In order to obtain an identifiable model based on Model (7), one might remove the random effect on one or several individual parameters. Fixing their values across the population might allow for a more precise estimation of the other, still variable, individual parameters, while keeping the same quantitative fit as with Models (2) and (7). However, all parameters are not necessarily equivalent in this regard, since different parameter sensitivities would make the model output more flexible under some combinations of fixed parameters. This would allow for these combinations to better fit the data, depending on the sensitivies of the model output to the parameter values. This sensitivity is imposed by the analytical solution to the structural equations of the model that we defined in System (1), but for most models there is no closedform expression for the parameter sensitivities. As a consequence, and in order to keep our approach as general as possible, we will not attempt any analytical expression of the model output sensitivities to the individual parameters herein.
In order to choose which random effect to remove from our model, we consider that the parameter with the highest shrinkage is the most poorly estimated across the population. Since \(\rho _C\) and \(\rho _B\) display similar amounts of shrinkage in Model (7), we might remove either of their random effects in order to reduce our model. Removing the random effect on \(\rho _C\) defines a new model with reduced \(\delta _{SC}^{pop}\) and \(\delta _{CB}^{pop}\), and with no variability on \(\rho _C\), which is described in System (S1, Additional file 1). Conversely, removing the random effect on \(\rho _B\) defines a new model with reduced \(\delta _{SC}^{pop}\) and \(\delta _{CB}^{pop}\), and with no variability on \(\rho _B\), which is described in System (S2, Additional file 1). We display the optimal likelihood distribution over an initial guess sample for Model (S1) and Model (S2) in Figures S3 and S6 (Additional file 1), respectively.
In both Model (S1) and Model (S2), the population parameters (Additional file 1: Figures S4 and S7) and the individual parameters (Table 1, Additional file 1: Figures S5 and S8) are still unidentifiable. We conclude that removing one random effect from Model (7) is not sufficient to make it identifiable. Thus, we propose to further reduce it by removing both the random effects on \(\rho _C\) and on \(\rho _B\), resulting in a reduced model with 3 remaining random effects:
In this model, every population parameters—including the remaining variances \(\omega _{\rho _S}\), \(\omega _{\delta _{SC}}\) and \(\omega _{\delta _{CB}}\)—are reliably estimated (Additional file 1: Figure S10) and the average shrinkage is lower than 30% for every parameter (Table 1, Additional file 1: Figure S11). In other words, our reduction approach allowed us to define a fully identifiable MEM, i.e. a model which is able to quantitatively reproduce the individual trajectories of Fig. 1, while explaining experimental heterogeneity in terms of precise parameter variations between individuals.
In the next section, we test another hypothesis, under which experimental heterogeneity does not come from interindividual variations of the parameters of the dynamic model, but rather from variations in the initial condition of the experiment. Then, we discuss the biological significance of the parameter values of Model 8.
Variability of the initial condition
In the previous section, we have considered that experimental heterogeneity originates from individual differences in the parameters of proliferation and differentiation kinetics. On the other hand, experimental heterogeneity might also be caused by an error in the sampling of the inital 25,000 cells in the culture. In this section, we assess whether a variability of the initial condition could better account for experimental heterogeneity than variability on the model dynamic parameters. This can be tested by defining Mixed Effect Models accounting for the heterogeneity of the initial population size. First, we define three alternative versions of Model (8), which differ by their definition of the initial condition. We then calibrate these models and study their identifiability in order to select the best model at reproducing our data that is also identifiable.
The first model that we are considering is our reduced model (8), in which the kinetic parameters \(\rho _S\), \(\delta _{SC}\) and \(\delta _{CB}\) can vary between individuals, while \(\rho _C\) and \(\rho _B\) are kept constant between individuals. In this model, the initial condition is fixed for all individuals: \(\left( S_0, T_0, B_0 \right) = (25{,}000, 25{,}000, 0 )\).
In the second model that we consider, all kinetic parameters are fixed to their population average, and we allow the initial condition to vary between individuals:
In this model, parameters \(S_0^{pop}\) and \(T_0^{pop}\) are the average of the initial number of cells in the LM1 and DM17 experiment. They represent a systematic error in the sampling of the initial 25,000 cells. Parameters \(\omega _{S_0}\) and \(\omega _{T_0}\) are the standard deviations of the initial number of cells in each experiment. For the third observable B, which is the number of differentiated cells, we consider the fixed initial condition \(B_0 = 0\) for all individuals, as the differentiation is initiated at time \(t=0\).
Finally, the last model that we consider allows for interindividual variations of both the kinetic parameters, as in Model (8), and the initial condition, as in Model (9):
For Models (9) and (10), we present the convergence data in Figures S12 and S15 (Additional file 1), the distributions of the population parameters in Figures S13 and S16 (Additional file 1), and the distributions of shrinkage values in Figures S14 and S17 (Additional file 1) respectively.
We display the optimal loglikelihood \(2\,log( \widehat{{\mathcal {L}}} )\) and the corresponding BIC for Models (8–10) in Table 2. Model (8) appears as the best one, closely followed by Model (10). On the other hand, Model (9) performs much worse than its competitors. Since Model (10) is unidentifiable (Additional file 1: Figures S16 and S17), we conclude that Model (8) is the best one both in terms of quality of the fit and of parameter identifiability. This means that individual variations in the parameter values are more important in accounting for experimental heterogeneity than variations in the initial condition.
The final model
In Model (8), the population parameters are identifiable (Fig. 5). It is the same for the individual parameters (Additional file 1: Figure S11). This means that every parameter of the model can be reliably estimated from our data, and that the estimated values reflect an actual optimum in the description of in vitro erythropoiesis.
The distribution of the estimated population parameter values in the convergent runs of Model (8) is displayed on Fig. 5. The estimated population average of the parameter values in the optimal run are displayed in Table 3. The fixed effects \(\rho _S^{pop}\), \(\rho _C^{pop}\) and \(\rho _B^{pop}\) determine the average behavior of the experiment. The average proliferation rate \(\rho _S^{pop}\) is estimated at \(0.61 \,\text {d}^{1}\) (Table 3). The doubling time of the selfrenewing cells (i.e. the time it would take to double their population in the absence of differentiation) is thus \(27\,\text {h}\) in the average experiment, which is longer than the originally reported \(18\,\text {h}\) [33]. Proliferation in the committed compartment is much faster (\(\rho _C^{pop} = 3.8\,\text {d}^{1}\), Table 3), which gives the committed cells an approximate doubling time of \(4\,\text {h}\). Even though T2EC cells are known to proliferate faster in the differentiation medium than in the selfrenewal medium [33], such a difference in proliferation times is rather intriguing.
Moreover, \(\rho _B^{pop}\) is estimated at \(0.26\,\text {d}^{1}\), giving the differentiated cells a doubling time of \(65\,\text {h}\). This means that their proliferation is almost invisible at the timescale of the experiments, as might be expected from differentiated cells.
From Eq. (5), the value of \(\rho _C^{pop}\) sets \(\delta _{SC}\) to an average \(0.37\,\text {d}^{1}\). The halflife of the selfrenewing cells (i.e. the time it would take to differentiate half of the population in the absence of proliferation) is thus approximately \(45\,\text {h}\). Respectively, from Eq. (6), the average value of \(\delta _{CB}\) in the population is estimated at \(4.5\,\text {d}^{1}\), which gives the committed cells a halflife of approximately \(3\,\text {h}\) in the average experiment.
Apart from this average behaviour, three parameters of the final model can vary across the population, and are estimated at different values for each individual experiment. The first one is \(\rho _S\), which has the estimated variance \(\omega _{\rho _S} = 0.12\,\text {d}^{1}\). This translates into the individual values of \(\rho _S\) being estimated between \(0.38\,\text {d}^{1}\) and \(0.85\,\text {d}^{1}\) (Table 3), which corresponds to doubling times between \(20\,\text {h}\) and \(44\,\text {h}\). Then, \(\omega _{\delta _{SC}}\) is estimated at \(0.25\,\text {d}^{1}\) with the individual parameter values of \(\delta _{SC}\) estimated from \(0.22\,\text {d}^{1}\) to \(0.59\,\text {d}^{1}\), which means that the corresponding halflife ranges from \(28\,\text {h}\) to \(76\,\text {h}\) approximately. Finally, \(\omega _{\delta _{CB}}\) is estimated at \(0.086\,\text {d}^{1}\), with individual parameter values for \(\delta _{CB}\) ranging from \(3.8\,\text {d}^{1}\) to \(5.3\,\text {d}^{1}\) and the corresponding halflife approximately ranging from \(3\,\text {h}\) to \(4\text {h}30\).
All these parameters are practically identifiable according to our initial guess sampling approach: each population parameter is estimated at a unique value (Fig. 5), and the distibution of individual parameters matches the population distribution (Additional file 1: Figure S11). However, identifiability analysis in Monolix is based on the FIM, and it was impossible to compute the FIM of our final model in any of the SAEM runs that we performed (Fig. 5): \(\rho _C^{pop}\) and \(\rho _B^{pop}\) have an infinite standard error when computed from the FIM. Based on the FIM, these two parameters thus appear as unidentifiable, while the others are all identifiable with a standard error which seems consistent between the SAEM runs. This means that the method used for identifiability analysis in MEM has an impact on the outcome of the analysis. However, since the FIM is proven to render biased confidence intervals when studying practical identifiability [6], our FIMfree method for identifiability analysis and model reduction appears as a reasonable approach with MEM.
Discussion and prospects
Generalizing our approach
Most approaches for identifiability analysis in MEM rely upon the FIM [10, 13, 14, 21,22,23], even though its parabolic approximation of the likelihood surface might mask complex practical unidentifiabilities [6]. We rather propose a multistart approach [13], that we refer to as Initial Guess Sampling.
Multistart approaches do not provide any information regarding the confidence intervals of the model parameters [20], but they indicate parameter unidentifiabilities when the estimated values differ.
For this reason, our samples of estimated values cannot be used in statistical analyses. On the contrary, they only allow for a visual check of the estimated values (for instance on Fig. 5), thus adding a new kind of diagnostic plot to the visual tools already available to the modeller using mixed effect models.
This approach allowed us to design an identifiable MEM of in vitro erythropoiesis which accounts well for experimental heterogeneity as interindividual variations of the proliferation and differentiation parameters. Then, a question that naturally arises is whether or not our approach could be applied, or generalized, onto other MEM? We identified two features of our approach, that might be of importance for such a generalization.
First, it seems that iteratively reducing our model affects the convergence of SAEM. Indeed, while all 50 runs of SAEM reached the same likelihood optimum for Model (2) (Fig. 3A), only 36 runs reached it for Model (7) (Additional file 1: Figure S1B). Convergence was further affected by removing random effects from our model (Additional file 1: Figures S3B & S6B), to the point where only 25 runs of SAEM reached the likelihood optimum for Model (8) (Additional file 1: Figure S9B). However, the number of convergent SAEM runs is critical to the assessment of population parameters identifiability. Depending on the complexity of the model and the dataset, 50 SAEM runs might not be sufficient to assess parameter identifiability and allow for model reduction, and the number of runs to be performed should thus be finely tuned in order to avoid issues with computational time.
Morevover, we used the correlations between population parameters to define Model (7), with constraints on the fixed effects. The exact shape of the likelihood landscape and the resulting unidentifiability is related to the structure of the model, and the quality of the data. This means that we were able to explore the parameter space near the likelihood optimum using pairwise correlations (Fig. 4). Yet, in more complex nonlinear MEM, it is possible that the correlations would involve more than two parameters at a time. In the end, detecting these complex correlations would require some kind of multivariate correlation analysis [41].
About the source of experimental heterogeneity
In this paper, we consider that experimental heterogeneity might either originate in variations of the kinetic parameters between replicates of the experiment, or by experimental errors in the initial number of cells. Using model selection and identifiability analysis, we conclude that variations in the kinetic parameters of proliferation and differentiation best explain experimental heterogeneity.
Considering that every replicate of the experiment was obtained with the exact same protocol, it seems that only two features of our experiment could change from replicate to replicate.
The first one is the group of 25,000 cells used to initiate the culture. In the haematopoietic system, in vivo stem cells and progenitors display substancial variations in terms of selfrenewal and potency [42]. Since our T2EC cells are erythropoietic progenitors, our results suggest that the selfrenewal and differentiation abilitites vary between the cell populations that we used to initiate every experiment.
On the other hand, there has also been discussion around the fact that the external temperature of the incubators (i.e. the temperature of the room where the T2EC are incubated, which is not their incubation temperature) might affect the variability of gene expression [43]. This in turn, could affect their selfrenewal or differentiation potency.
Conclusion
In this paper, we proposed a MEM for in vitro erythropoiesis, that accounts for experimental heterogeneity. We developped a multistart approach for assessing its identifiability, and we successfully reduced it to make it identifiable. We showed that experimental heterogeneity is faithfully accounted for by variations of the kinetic parameters of proliferation and differentiation in our system, and we relate these parameter variations to actual biological features of our cells. This work establishes a MEM framework to study variability in the outcome of biological experiments. Furthermore, it proposes a novel approach for the analysis of parameter identifiability in MEM, and for reducing unidentifiable MEM.
Availability of data and materials
All the datasets and pieces of code analysed and generated during the current study are available in a public github repository, at https://github.com/rduchesn/MixedEffectModelReduction.
References
 1.
Huang S. Nongenetic heterogeneity of cells in development: more than just noise. Development. 2009;136(23):3853–62. https://doi.org/10.1242/dev.035139.
 2.
Andersen SW, Millen BA. On the practical application of mixed effects models for repeated measures to clinical trial data. Pharm Stat. 2013;12(1):7–16. https://doi.org/10.1002/pst.1548.
 3.
Rué P, Martinez Arias A. Cell dynamics and gene expression control in tissue homeostasis and development. Mol Syst Biol. 2015;11(2):792. https://doi.org/10.15252/msb.20145549.
 4.
Kreso A, O'Brien CA, Pv Galen, Gan OI, Notta F, Brown AMK, Ng K, Ma J, Wienholds E, Dunant C, Pollett A, Gallinger S, McPherson J, Mullighan CG, Shibata D, Dick JE. Variable clonal repopulation dynamics influence chemotherapy response in colorectal cancer. Science. 2013;339(6119):543–8. https://doi.org/10.1126/science.1227670.
 5.
Pradhan BB, Chatterjee S. Reversible nongenetic phenotypic heterogeneity in bacterial quorum sensing. Mol Microbiol. 2014;92(3):557–69. https://doi.org/10.1111/mmi.12575.
 6.
Raue A, Kreutz C, Maiwald T, Bachmann J, Schilling M, Klingmüller U, Timmer J. Structural and practical identifiability analysis of partially observed dynamical models by exploiting the profile likelihood. Bioinformatics. 2009;25(15):1923–9. https://doi.org/10.1093/bioinformatics/btp358.
 7.
Raue A, Schilling M, Bachmann J, Matteson A, Schelke M, Kaschek D, Hug S, Kreutz C, Harms B, Theis F, Klingmüller U, Timmer J. Lessons learned from quantitative dynamical modeling in systems biology. PLoS ONE. 2013;8(9):74335. https://doi.org/10.1371/journal.pone.0074335.
 8.
Wieland FG, Hauber AL, Rosenblatt M, Tönsing C, Timmer J. On structural and practical identifiability. Curr Opin Syst Biol. 2021;25:60–9. https://doi.org/10.1016/j.coisb.2021.03.005.
 9.
Cobelli C, Saccomani MP. Unappreciation of a priori identifiability in software packages causes ambiguities in numerical estimates. Am J Physiol Endocrinol Metab. 1990;258(6):1058–9. https://doi.org/10.1152/ajpendo.1990.258.6.E1058.
 10.
Lavielle M, Bleakley K. Mixed effects models for the population approach: models, tasks, methods and tools. London: Chapman & Hall; 2014.
 11.
Karlsson M, Janzén D, Durrieu L, ColmanLerner A, Kjellsson M, Cedersund G. Nonlinear mixedeffects modelling for single cell estimation: when, why, and how to use it. BMC Syst Biol. 2015. https://doi.org/10.1186/s129180150203x.
 12.
Savic R, Karlsson M. Importance of shrinkage in empirical Bayes estimates for diagnostics: problems and solutions. AAPS J. 2009;11(3):558–69. https://doi.org/10.1208/s1224800991330.
 13.
Lavielle M, Aarons L. What do we mean by identifiability in mixed effects models? J Pharmacokinet Pharmacodyn. 2016;43(1):111–22. https://doi.org/10.1007/s1092801594594.
 14.
Shivva V, Korell J, Tucker IG, Duffull SB. An approach for identifiability of population pharmacokineticpharmacodynamic models. CPT Pharmacomet Syst Pharmacol. 2013;2(6):49. https://doi.org/10.1038/psp.2013.25.
 15.
Janzén DLI, Jirstrand M, Chappell MJ, Evans ND. Three novel approaches to structural identifiability analysis in mixedeffects models. Comput Methods Programs Biomed. 2019;171:141–52. https://doi.org/10.1016/j.cmpb.2016.04.024.
 16.
Janzén DLI, Jirstrand M, Chappell MJ, Evans ND. Extending existing structural identifiability analysis methods to mixedeffects models. Math Biosci. 2018;295:1–10. https://doi.org/10.1016/j.mbs.2017.10.009.
 17.
Jacquez JA, Greif P. Numerical parameter identifiability and estimability: integrating identifiability, estimability, and optimal sampling design. Math Biosci. 1985;77(1):201–27. https://doi.org/10.1016/00255564(85)900987.
 18.
Vajda S, Rabitz H, Walter E, Lecourtier Y. Qualitative and quantitative identifiability analysis of nonlinear chemical kinetic models. Chem Eng Commun. 1989;83(1):191–219. https://doi.org/10.1080/00986448908940662.
 19.
Pillai GC, Mentré F, Steimer JL. Nonlinear mixed effects modeling—from methodology and software development to driving implementation in drug development science. J Pharmacokinet Pharmacodyn. 2005;32(2):161–83. https://doi.org/10.1007/s109280050062y.
 20.
Fröhlich F, Theis F, Hasenauer J. Uncertainty analysis for nonidentifiable dynamical systems: profile likelihoods, bootstrapping and more. In: Mendes P, Dada J, Smallbone K, editors. CMSB. Cham: Springer; 2014. p. 61–72.
 21.
Wang W. Identifiability of linear mixed effects models. Electron J Stat. 2013;7:244–63. https://doi.org/10.1214/13EJS770.
 22.
Mentré F, Mallet A, Baccar D. Optimal design in randomeffects regression models. Biometrika. 1997;84(2):429–42. https://doi.org/10.1093/biomet/84.2.429.
 23.
Retout S, Mentré F. Further developments of the fisher information matrix in nonlinear mixed effects models with evaluation in population pharmacokinetics. J Biopharm Stat. 2003;13(2):209–27. https://doi.org/10.1081/BIP120019267.
 24.
White A, Tolman M, Thames H, Withers H, Mason K, Transtrum M. The limitations of modelbased experimental design and parameter estimation in sloppy systems. PLoS Comput Biol. 2016;12(12):1005227. https://doi.org/10.1371/journal.pcbi.1005227.
 25.
Chachra R, Transtrum MK, Sethna JP. Comment on “Sloppy models, parameter uncertainty, and the role of experimental design”. Mol BioSyst. 2011;7(8):2522–2522. https://doi.org/10.1039/C1MB05046J.
 26.
Cole D. Parameter redundancy and identifiability. London: Chapman and Hall; 2020. https://doi.org/10.1201/9781315120003.
 27.
Bazzoli C, Retout S, Mentré F. Design evaluation and optimisation in multiple response nonlinear mixed effect models: PFIM 3.0. Comput Methods Programs Biomed. 2010;98(1):55–65. https://doi.org/10.1016/j.cmpb.2009.09.012.
 28.
Vanrolleghem PA, Dochain D. Bioprocess model identification. In: Advanced instrumentation, data interpretation, and control of biotechnological processes. Springer; 1998. p. 251–318.
 29.
Duchesne R, Guillemin A, Crauste F, Gandrillon O. Calibration, selection and identifiability analysis of a mathematical model of the in vitro erythropoiesis in normal and perturbed contexts. ISB. 2019;13(1–2):55–69. https://doi.org/10.3233/ISB190471.
 30.
Monolix version 2018R1. http://lixoft.com/products/monolix/. Antony: Lixoft SAS; 2018.
 31.
Richard A, Boullu L, Herbach U, Bonnafoux A, Morin V, Vallin E, Guillemin A, Papili Gao N, Gunawan R, Cosette J, Arnaud O, Kupiec J, Espinasse T, GoninGiraud S, Gandrillon O. Singlecellbased analysis highlights a surge in celltocell molecular variability preceding irreversible commitment in a differentiation process. PLoS Biol. 2016;14(12):1002585. https://doi.org/10.1371/journal.pbio.1002585.
 32.
Guillemin A, Duchesne R, Crauste F, GoninGiraud S, Gandrillon O. Drugs modulating stochastic gene expression affect the erythroid differentiation process. PLoS ONE. 2019;14(11):0225166. https://doi.org/10.1371/journal.pone.0225166.
 33.
Gandrillon O, Schmidt U, Beug H, Samarut J. TGFbeta cooperates with TGFalpha to induce the selfrenewal of normal erythrocytic progenitors: evidence for an autocrine mechanism. EMBO J. 1999;18(10):2764–81. https://doi.org/10.1093/emboj/18.10.2764.
 34.
Kuhn E, Lavielle M. Maximum likelihood estimation in nonlinear mixed effects models. Comput Stat Data Anal. 2005;49(4):1020–38. https://doi.org/10.1016/j.csda.2004.07.002.
 35.
Burnham K, Anderson D. Model selection and multimodel inference: a practical informationtheoretic approach. New York: Springer; 2010. OCLC: 934366523.
 36.
Vaida F, Blanchard S. Conditional Akaike information for mixedeffects models. Biometrika. 2005;92(2):351–70. https://doi.org/10.1093/biomet/92.2.351.
 37.
Delattre M, Lavielle M, Poursat M. A note on BIC in mixedeffects models. Electron J Stat. 2014;8(1):456–75. https://doi.org/10.1214/14EJS890.
 38.
Delattre M, Poursat M. BIC strategies for model choice in a population approach. arXiv:1612.02405 [stat]; 2016. Accessed 17 May 2017.
 39.
Koop G, Pesaran MH, Smith RP. On identification of Bayesian DSGE models. J Bus Econ Stat. 2013;31(3):300–14. https://doi.org/10.1080/07350015.2013.773905.
 40.
Karlsson M, Savic R. Diagnosing model diagnostics. Clin Pharmacol Ther. 2007;82(1):17–20. https://doi.org/10.1038/sj.clpt.6100241.
 41.
Allison PD. Multiple regression: a primer. Thousand Oaks: Pine Forge Press; 1999. Open Library ID: OL378019M.
 42.
Haas S, Trumpp A, Milsom MD. Causes and consequences of hematopoietic stem cell heterogeneity. Cell Stem Cell. 2018;22(5):627–38. https://doi.org/10.1016/j.stem.2018.04.003.
 43.
Arnaud O, Meyer S, Vallin E, Beslon G, Gandrillon O. Temperatureinduced variation in gene expression burst size in metazoan cells. BMC Mol Biol. 2015;16(1):20. https://doi.org/10.1186/s1286701500482.
Acknowledgements
We thank the members of the Systems Biology of Decision Making team at the Laboratory of Biology and Modelling of the Cell, as well as the Dracula team of Inria, for insightful discussion all along the course of this project. We thank the BioSyL Federation and the LabEx Ecofect (ANR11LABX0048) of the University of Lyon for inspiring scientific events. We also would like to thank members of the Population Approach Group Europe that were present at the PAGE meeting 2019 for enlightening feedback. Finally, we want to thank Pr. Saccomani (University of Padova, Italy) for providing the manuscript for reference [9].
Funding
RD and AG benefited from PhD funding from the French Ministeère de la Recherche et de l’Enseignement supérieur.
Author information
Affiliations
Contributions
AG generated the experimental data used in this study. RD designed the models, estimated their parameter values, carried out all subsequent analyses, and wrote the first draft of the manuscript. FC and OG supervised the project. All authors read and approved the final manuscript.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
All methods were carried out in accordance with relevant guidelines and regulations, notably the DIRECTIVE 2010/63/EU OF THE EUROPEAN PARLIAMENT AND OF THE COUNCIL of 22 September 2010 regarding  the killing of animals solely for the use of their organs or tissues.(https://eurlex.europa.eu/legalcontent/EN/TXT/HTML/?uri=CELEX:32010L0063)
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary Information
Additional file 1.
Supplementary Materials.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
Duchesne, R., Guillemin, A., Gandrillon, O. et al. Practical identifiability in the frame of nonlinear mixed effects models: the example of the in vitro erythropoiesis. BMC Bioinformatics 22, 478 (2021). https://doi.org/10.1186/s12859021043734
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s12859021043734
Keywords
 Nonlinear mixed effects models
 Practical identifiability
 Model reduction