Approximate Moment Methods for Population Balance Equations in Particulate and Bioengineering Processes

: Population balance modeling is an established framework to describe the dynamics of particle populations in disperse phase systems found in a broad ﬁeld of industrial, civil, and medical applications. The resulting population balance equations account for the dynamics of the number density distribution functions and represent (systems of) partial differential equations which require sophisticated numerical solution techniques due to the general lack of analytical solutions. A speciﬁc class of solution algorithms, so-called moment methods, is based on the reduction of complex models to a set of ordinary differential equations characterizing dynamics of integral quantities of the number density distribution function. However, in general, a closed set of moment equations is not found and one has to rely on approximate closure methods. In this contribution, a concise overview of the most prominent approximate moment methods is given.


Introduction
Disperse phase systems are not only a core part of many industrial processes, such as crystallization, granulation, drying, and fermentation, but are also encountered in biomedicine, astrophysics, and meterology. They are characterized by a population of individuals (particles) that interact with their environment and with each other. Particles can be characterized by individual properties, e.g., characteristic size, mass, density, or composition, which have a direct effect on corresponding systems dynamics as well as the final product properties (e.g., structural stability and flowability of granule product).
In general, a variance within the particle population with respect to the properties is observed. This heterogeneity affects both systems dynamics and product properties. Population balance modeling (PBM) represents an established framework to obtain mathematical descriptions of the dynamics of these heterogeneous disperse phase systems [1][2][3][4][5]. Population balance modeling has found successful application in a broad range of fields from agriculture [6,7] and civil engineering [8,9] over thermal and chemical process engineering [10][11][12] to biotechnology, bioengineering, and biomedicine [13][14][15].
For an extended list of applications, the reader is referred to the work of Ramkrishna and Singh [4] (and references therein).
The resulting models represent (systems of) partial differential equations, so-called population balance equations (PBEs), that account for the change of the particle number density distribution function (NDF) by different mechanisms and inter-phase interactions. As analytical solutions can only be derived for rare cases, one relies on numerical solution approaches. Besides discretization-based  and stochastic approaches [44][45][46][47][48][49][50][51][52][53][54][55][56][57], moment methods represent an important category. Instead of the full number density distribution, the dynamics of integral quantities are computed, such as the mean and the standard deviation. Those are governed by ordinary differential equations (ODEs) which are obtained by integration of the PBE. In general, a closed set of moment equations is rarely found and, therefore, an approximate closure is necessary (i.e., a sound approximation of unknown moments by known ones).
Over the last decades, many approximate moment methods have been developed. This contribution provides a concise overview of the most prominent techniques. Thereby, focus is on methods that aim to solve directly for the moments without having to calculate the full distribution first. Techniques of the latter class, which are applied, e.g., in relation to the solution of the Fokker-Planck equation [58,59], would be a worthwhile topic for a stand-alone manuscript and are therefore not included in this review. The remainder of this manuscript is therefore structured as follows. In the next two sections, a brief introduction into population balance modeling is given and the derivation of the moment dynamics is presented. Subsequently, the moment closure problem is outlined and a selection of approximate closure techniques is presented. Furthermore, algorithms for multivariate PBEs are discussed and further needs of development are identified.

Population Balance Modeling
In a disperse phase system, the disperse phase consists of an ensemble of individuals, e.g., granules, bubbles, or cells, which can by characterized by individual properties. The latter are categorized into external properties x e (spatial coordinates) and internal properties x i (e.g., size, shape or composition). The number of particles contained in a certain volume X of the property state space X, is defined as with the property state vector being defined as x = [x i , x e ] T and n(t, x) representing the number density distribution function. As result of interactions within the disperse phase or other phases, properties of the particles underlie change and thereby the corresponding NDF changes. The dynamics of the number density distribution can be modeled by the well known population balance equation (PBE) in its local form [60] ∂n(t, with the first term on the right-hand side accounting for transport processes in state space such as convection/diffusion in spatial coordinates and continuous change of internal properties e.g., by growth/shrinking, as well as metabolic conversion. While the first part represents convection with rate G(t, x), the second part accounts for diffusion with rate D(t, x). In contrast, p x (t, x) represents sources and sinks from chemical reactions or withdrawal and particle interactions. For the the classical mono-variate case (x = x), e.g., a population balance characterizing characteristic volume of a particle breakage/division process is usually described as with S(t, x) the selection function characterizing probability for breakage depending on particle property x and the breakage kernel β B (t, u, x) describing the size distribution of smaller particles generated from breakage of bigger ones. In contrast, the source/sink term for coagulation/agglomeration is given as Therein, β A is the agglomeration/coagulation kernel, which describes the agglomeration probability depending on the properties of two colliding particles with individual particle properties x and u. When describing dynamics of a certain process with PBM, two general procedures can be followed: (I) In the top-down-modeling approach, only a few particle properties are taken into account (often only one). The unknown fluxes and source/sink mechanics of the corresponding PBEs are either described by semi-empirical approaches or adapted to experimental data with suitable methods (see, e.g., [61][62][63][64][65][66]). (II) In contrast, the bottom-up approach is based on detailed first-principle modeling of the individual particle behavior, e.g., ordinary or stochastic differential equations, which is adapted to single-particle experimental data. Each state of the single particle description transforms into a coordinate of the corresponding PBE and thus in general high-dimensional multivariate PBEs are obtained. This modeling approach is particularly useful when describing multicellular systems [2,[67][68][69][70][71][72] for which validated models of single cell dynamics are available.
Modeling of many phenomena, in particular for particle manufacturing processes, combines both procedures within a hybrid strategy: First principle modeling is used to a certain extent and remaining uncertain parameters are adapted to available data.

Moment Transformation and the Moment Closure Problem
In particular, for process engineering applications, it is more convenient to characterize the particle population in terms of integral measures with respect to the internal coordinates: Important examples are mean and standard deviation of the corresponding NDF, which can be obtained by suitable choice of the weighting function f (x). For the (mixed) integer order moments m l (t, x e ) of the NDF with respect to the internal properties are obtained. Without loss of generality, ideal mixing is assumed in the following if not stated otherwise. Hence, the spatial gradients in the external properties can be neglected. The moment dynamics are derived by differentiation of m l (t) as The specific integral expressions for the transport related and source/sink dynamics, J(t, x) and P(t, x) can be derived by integration of the related expressions of the general population balance equation. To calculate the dynamics of an arbitrary moment, the right-hand side itself must only contain moments whose dynamics are described by a closed set of equations. In fact, this not possible in general as certain moment dynamics depend on unknown higher order or fractional moments.
As illustrative example, a mono-variate PBE describing a general growth process is considered: where x represents the particle size, with the moment dynamics given by Here, average particle size x av and average particle surface A av are given as We further assume that the growth velocity is given as G(t, x) = x p . Thus, moment dynamics are obtained as It is easily derived that the resulting moment dynamics are closed for p = 0, 1 but not for monomial-type growth rates of higher order (p ≥ 2). In fact, it is seen, for example, that the dynamics of an arbitrary moment m l depends on the higher order moment ml(t) = m p+l−1 (t). To determine the dynamics nevertheless, these higher order moments have to be expressed in terms of lower order moments where A represents a suitable approximation technique. Thereby, one obtains an artificial or approximate closure: Since the advent of PBM in the 1960s [1], different techniques for approximate solution of moment dynamics have been proposed. There are applications of those approximate moment methods for most population balance modeling problems, in particular as a measure to reduce the computational complexity in computational fluid dynamics where the moment transform is applied on the general PBE with respect to the internal properties x i [73][74][75] Here, f x e (t, x e ) characterize transport in the external property state space via convection and diffusion and S(•) accounts for sources and sink terms within the external property state space.

Approximate Moment Methods for Mono-Variate Population Balance Equations
The sophisticated approximate moment methods used nowadays differ not only in their underlying principles, but also in numerical effort and accuracy. In this section, the most important and widely used techniques for mono-variate population balance equations are presented and discussed. The multidimensional case is covered separately in the following section.

Method of Moments with Interpolating Closure (MoMIC)
The MoMIC was initially developed for efficient numerical solution of Brownian coagulation models [76]. Here, the corresponding PBE with discrete valued property space is reduced to an infinite set of differential equations accounting for particle number dynamicsṄ i in different size classes i. Particles change class due to coagulation, growth, and nucleation [77]. A physically realistic dependence of the collision frequency of particles of different sizes yields a non-closing moment system. Therein, the specific form of the collision frequency depends on the current process conditions. This gives rise to integer moment dynamicsṁ k (t) which depend on unknown fractional order moments m k * . To obtain closure, the latter are approximated from interpolation among N r tracked integer order moments Furthermore, Frenklach et al. [77] discussed the potential of using separated interpolation for positive and negative integer order dynamics to improve approximation accuracy: in the former case, a i are determined by Lagrange-interpolation using logarithms of integer order moments, while, for the latter, a i are computed by Lagrange interpolation between integer order negative normalized moments.
The primary advantage of the MoMIC is its easy implementation and the relatively low computational effort because to be solved it requires e only a few differential equations. This is obviously a crucial point when applied to complex systems that involve simultaneous computation of particle dynamics and spatial flow regime. This fact has facilitated application to related problem classes, e.g., for soot formation [78][79][80]. It has further been mentioned that the algorithm can become inaccurate for problems with non-negligible particle breakage where influence of small particles on the moment dynamics is not fully captured [78]. In these cases, different extensions and combinations with other approximate moment methods have been developed (see, e.g., [79,80]).

Quadrature Method of Moments (QMoM)
QMoM [81,82] relies on the representation of the NDF by a sum of weighted δ-functions where w α and x α are called the weights and abscissas, respectively. Hence, moments are approximated as The N α abscissas and N α weights have to be determined from 2N α moments using a convenient inversion technique. The moment set used for such an inversion is not unique, thus different choices yield different weights and abscissas. In fact, inversion algorithm and moment dynamics form a closed differential algebraic equation (DAE) system [83] with unknown moments approximated as Numerical solution of DAE systems is more involved than numerical solution of ODE systems alone, as in each time-step for the differential equations an internal root-finding problem has to be solved to fulfill the additional algebraic equations. This requires information about the Jacobian matrix of the DAE model, as often Newton-type algorithms are employed for root-finding. Typically, the entries of the Jacobian matrix have to be obtained numerically, which, depending on the accuracy and number of entries, is a time-consuming process. Inaccurate estimates of the Jacobian matrix slow down the convergence of the root-finding method or introduce errors in the numerical result.
To improve the accuracy of the elements of the Jacobian matrix, automatic differentiation can be applied (AD-QMoM [84]). The key idea is that mathematical operations are augmented with the respective differentiation rules, so that, in each operation, not only the result is calculated but also the values of the corresponding derivatives (up to any order). This allows calculation of the Jacobian matrix to machine precision. Matured packages are available to support model development, for instance ADIFOR [85], ADOL-C [86], and ADIC [87], which have also been ported to other programming languages and development environments. Although this approach has been shown to provide the required Jacobian matrices, the computational burden is still high.
To avoid a numerically involved solution of the DAE-system, a two-step solution algorithm is usually employed: Starting from an initial set of moments, weights and abscissas are computed using established methods, e.g., product difference algorithm [81] or Wheeler algorithm [88]. Subsequently, the moments at the next time step are obtained from numerical integration of the approximate moment dynamics. This procedure is repeated until the final simulation step is reached.
Over the years different formulations of the QMoM have been developed and thorough numerical analysis and investigations were performed: generally, an increased number of nodes promotes an improved approximation accuracy. However, a low number of nodes is sufficient in many applications and it is recommended not to go beyond a certain amount of abscissas (N α ≤ 6) which corresponds to a limitation in the number of tracked moments [89]. Dorao et al. [90] discussed the relationship of QMoM to the method of weighted residuals and suggested that this fact could be used as to overcome the ill-conditioned moment inversion problem and further improve accuracy, which has been investigated elsewhere (e.g., [91,92]). In [93,94], different formulations are discussed and compared by means of test cases with respect to robustness, computational speed and accuracy. It has been shown that numerical accuracy may deteriorate if the modeled process is characterized by highly localized phenomena such as fines dissolution and sharp separation functions. In fact, the moment inversion algorithm in QMoM for computation of weights and abscissas becomes ill-conditioned easily, particularly for higher numbers of abscissas.
Due to its balance in favorable properties such as robustness, accuracy, and relatively low computational effort, QMoM has found wide application for efficient solution of PBEs emerging from modeling of spatially distributed polydisperse gas-liquid-particle flows [89,[95][96][97][98][99]. Further applications include crystallization modeling and control [100], continuous thermodynamics [101], and bioreactor applications [102]. The latter publication compares QMoM with other approximate moment methods and results indicate that QMoM is sufficiently accurate for a relatively low number of abscissas.

Direct Quadrature Method of Moments (DQMoM)
While the QMoM involves tracking of moment dynamics and solution of a nonlinear system of equations to infer weights and abscissas, in [73] an alternative formulation was proposed: the so-called direct quadrature method of moments. Here, the dynamics of weights and weighted abscissas are tracked directly. For a monovariate PBE with x = x and negligible spatial gradients (again for sake of simplicity and without loss of generality), those dynamics read where the right-hand sides are determined from the solution of the linear equation system with the 2N α rows of for 2 N α distinct values of l. As discussed above for QMoM, the choice of the moment set is not unique. It is however recommended to adapt the choice of l to the moments that are interested for the specific application, e.g., zeroth-or first-order moments. Here, the term S represents the right-hand side of the general PBE in Equation (2) and S l a suitable integral approximation. Analogously to QMoM, the direct method has been subject to intense mathematical analysis and investigations as well as extensions (e.g., [103]). In fact, e.g., Dürr and Kienle [104] showed that for certain conditions on S, DQMOM reduces to the Method of Characteristics and thus enable easy derivation of analytical solutions. However, it is well known that the procedure is prone to singularities. This is a particular problem if two or more abscissas move close to each other, as observed for shrinking or breakage dominated processes [73]. In this case, random perturbations of the abscissas may be employed to overcome this problem [89]. Furthermore, it has been shown that the classical DQMoM formulation generally does not conserve moments of order higher than one when applied to systems with spatial gradients in the disperse phase, e.g., gas-solid fluidized beds [105], exhaust particle formation [106], or bubbly flows [107,108]. For that case, Buffo and coworkers [109] presented a fully conservative formulation. Further applications are found also in bioprocess engineering, e.g., enzymatic hydrolysis of cellulose [110].

Extended Quadrature Method of Moments (EQMoM)
One of the drawbacks of (D)QMoM is that increasing accuracy in terms of a higher amount of abscissas can yield ill-conditioned moment-inversion and thereby numerical instability and failure of the overall approximate moment method. This motivated the development of EQMoM, which provides an explicit form of the NDF. Furthermore, the number of quadrature nodes can be increased independently from the number of transported moments and thereby achieve a more accurate representation of non-closing source terms. The extension relies on approximation of the NDF by kernel densities that depend on a single parameter, the bandwidth/support (σ) [111,112], to enable direct solution of the moment inversion problem. Here, the NDF is, as in the previous sections for the single particle property case, represented as a weighted sum where weights w α are non-negative, x α are the abscissas, and f σ b represents a non-negative kernel density function with support determined by parameter σ b . In case σ → 0, the kernel reduces to a δ-function and thus standard (D)QMoM approximation is obtained. The unknown weights w α and abscissas x α as well as the support σ b are computed exactly from the the moment set with a suitable method based on the (adaptive) Wheeler algorithm [88,113]. Application of different kernel functions with either finite or infinite support has been studied and analyzed with respect to numerical accuracy and involved computational effort: The results of Pigou et al. [114] indicate that QMoM and EQMoM performs far better with less computational effort than the Maximum Entropy algorithm. While the first works used Gaussian kernels [111] later works applied Gamma and Beta kernel functions [112] as well as Laplace, Log-normal, and Weibull kernels [114]. Several extensions and improvements have been proposed to stabilize and reduce the computational effort of the implementation [114][115][116]. The EQMoM is applied in similar areas as its predecessors, e.g., polydisperse bubbly flows [117,118], soot oxidation [116,119], and bioreactors [102].

Sectional Quadrature Method of Moments (SQMoM)
A general property of the QMoM is that all information on the NDF over the full property space is lumped into integral quantities. Thereby, crucial information on the shape of the NDF is lost, which could be crucial for process analysis and design. Furthermore, as previously mentioned, the moment-inversion is prone to ill-conditioning in particular when a higher amount of abscissas is used aiming on improved accuracy. Both drawbacks are addressed within SQMoM [120,121]. Here, the particle property space is divided into a finite number of sections N pp of arbitrary width [d i−1/2 , d i+1/2 ], named "Primary Particles". The NDF in each section is approximated using N sp δ-functions, the so-called "Secondary Particles", and moments of the ith section can be computed as For a set of 2N sp sectional moments, abscissas and weights can be determined by moment inversion, e.g., PD [81] or Wheeler algorithm [88]. Information on weights and weighted position of the primary particles is then obtained by averaging over secondary particles. The corresponding dynamics of a specific sectional moment generally depend on primary and secondary particles in other sections. Attarakih et al. [121] provided further details, including implementation; they also suggested to use a low number of secondary particles for improved conditioning of the moment inversion and to increase accuracy by using a larger number of primary particles. SQMoM has been applied in combination with CFD modeling of liquid extraction columns [122][123][124][125] and bubble columns [74,126].

Maximum Entropy (ME) Approach
In general, arbitrary NDF can be described by an infinite number of moments. If only a finite amount of moments is used, then the representation is not unique, e.g., all normal distributions cannot be distinguished using only moments of zeroth order (i.e., overall number of particles) and first order (related to the average value of a property): uniqueness requires the second moment, which is related to the NDF's variance. In contrast to the aforementioned techniques, the maximum entropy approach aims on finding the NDF which maximizes the Shannon entropy, resulting in the following NDF representation [127,128] which has to fulfill the condition that the first N α + 1 moments are reproduced exactly However, Tagliani [128] described integration on the specific finite support x ∈ [0, 1]; the technique can easily be extended to the more general support [a, b] by linear change of variable. The polynomial coefficients φ i for i = 1, . . . , N α are determined by minimization of while φ 0 is determined from Equation (31) for k = 0. Pigou et al. [102] applied the method for the numerical solution of a PBE characterizing the dynamic adaption of bacteria to its environment. It is pointed out that the numerical effort is rather large compared to standard QMoM or EQMoM despite allowing accurate reconstruction of the underlying NDF [110].

Moment Projection Method (MPM)
One significant disadvantage of the previously mentioned methods such as QMoM and DQMoM is that they tend to fail in case of shrinkage dominated processes where the pointwise value of the PSD at the smallest particle mass is required to obtain moment closure. Here, abscissas would move closer and closer to each other resulting in singularities and the previously described ill-conditioning of the moment inversion. While several approaches have been developed to overcome this limitation (e.g., [78,129]), robustness issues motivated development of the moment projection method [130].
Here, for spatially homogeneous problems with discrete-mass distribution, moments are expressed by a sum of N α mass classes with N α j referring to the number of particles of mass α j . Similar to the QMoM, their values are determined by inversion from a system of 2N α moments. However, in contrast to QMoM, one has to ensure that α 1 = m 1 , i.e., the mass of the smallest particles, which can be done using an adapted version of the Blumstein-Wheeler algorithm [131]. Thereby, it can be guaranteed that the smallest particle class is accurately taken into account. The method has been applied recently in combination with the ME approach for numerical solution of PBEs characterizing for of soot formation in combustion engines [132,133].

Approximate Moment Methods for Multivariate Population Balance Systems
The previously described moment transformation for the one-dimensional case can also be derived in a similar fashion for the multidimensional case. Approximate closure of the multidimensional moment dynamics is generally challenging as most one-dimensional approaches cannot readily be applied. For this reason, some extensions have been developed in the last decades. A selection is listed in the following:

•
Müller and coworkers [134] extended the MoMIC method to solve a bivariate model of soot aggregation where the particles are distinguished in characteristic volume V and surface S. As for the mono-variate case, unknown moments are interpolated via polynomial interpolation with coefficients a r,k resulting from solution of a linear equation system for R distinct known moments.

•
Yoon et al. [135,136]  The bivariate moment projection algorithm is presented in [138]. Here, the NDF is assumed as product of a marginal and a conditional NDF, where the former is mono-variate. As in previously presented the mono-variate case, the sizes of some weighted particles are fixed at the minimum size coordinates to evaluate boundary fluxes and accumulation of particles due to shrinkage and fragmentation.
Most of these techniques have been developed with low (two-or three-dimensional) applications in mind but cannot be applied in a general multidimensional case. In particular, for modeling of bio-systems, where cell-to-cell variability manifests, for example in difference of cell size and (bio-)chemical composition, complex high-dimensional PBE systems emerge that take into account heterogeneity of many cellular properties. For the efficient numerical solution of such models via approximate moment dynamics, only a few adequate approaches exist to the best of the authors' knowledge. These are briefly presented below.

Conditional Quadrature Method of Moments (CQMoM)
CQMoM represents an extension of QMoM: here, the one-dimensional quadrature is subsequently applied to conditional densities of the multivariate NDF [113]. The accuracy of the approximation can be adjusted by using different numbers of abscissas for each level [109]. Furthermore, a combination with EQMoM is possible [139], which enables improved estimation of underlying multivariate NDF from the moments. In [140], different bivariate CQMoM extensions with the EQMoM and the SQMoM are discussed. Application of CQMoM to multivariate problems may yield a weak hyperbolic system of moment equations which require thorough mathematical treatment due to the possibility of shock formation [111,141,142]. The method has found wide application in recent studies, e.g., for simulation of turbulent polydisperse gas-liquid systems [109,143,144], prediction of particle formation in lightly sooting burner-stabilized premixed ethylene flames [145], and nanoparticle growth [146]. Although CQMoM application to multidimensional PBEs emerging in bio-systems modeling is not known, general applicability is conjectured. However, up to the date of this review, no study is found that analyzes computational effort and robustness for high-dimensional applications. Furthermore, the effect of the chosen sequence of conditional densities is not researched so far.

DQMoM with Efficient Choice of Abscissas
The genuine property of the DQMoM is that it can be extended easily to multivariate PBEs in principle, e.g. bivariate extension for nanoparticle processing [147], sintering [148], and comminution grinding processes [149,150]. However, the choices of weights w α and abscissas l are not unique and depend on the specific moments used for generating the linear system of equations. Moreover, the choice of certain moment sets may lead to singularities. While optimal moment sets are known for low dimensional multivariate problems (dim(x) ≤ 3) [151] as well as application of fractional order moments [89], application to higher dimensional multivariate problems PBEs such as obtained during modeling of multicellular systems in systems biology or biotechnology is also discussed. In [104], DQMoM is combined with an efficient choice of weights and abscissas that is derived on the basis of approximate monomial cubatures [152]. The method relies on the assumption that the initial NDF is Gaussian. However, application to non-Gaussian NDFs can be achieved by using either transformations or mixed-Gaussian approximations [70]. In contrast to Gaussian quadrature formulas, the number of abscissas increases only linearly with the dimension of the property state space such that high dimensional multivariate PBEs become manageable. The method has been applied successfully for solution of high-dimensional multivariate PBEs of biopolymer- [68,153] and vaccine-production processes [69,70,154]. Recently, application to model-based computer aided screening for optimal genetic modifications in high-yield producer cell lines has been demonstrated [72,155]. An extension to problems involving agglomeration and breakage/cell division has been proposed [71].

Summary, Conclusions, Outlook
Disperse phase systems are found in a wide field of applications including industrial manufacturing and biomedicine. The heterogeneity of individual particles, with respect to particle properties such as size and composition within a population can significantly affect the overall process as well as product properties. The dynamics of such systems in terms of the particles NDF can be described with PBM, which can be solved numerically by moment methods. Here, the PBE is reduced to a set of ODEs characterizing integral quantities of the underlying NDF. However, the resulting moment systems generally do not close and an approximate closure becomes necessary. Over the years, sophisticated approaches have been developed and successfully applied, e.g., to assess the dynamics of complex multiphase disperse phase systems.
For one-dimensional PBEs, i.e., particle dynamics which characterize individuals only by one internal particle property, several well established methods are available. The choice of method depends on the specific application: while in most cases standard QMoM will produce sufficient results, specific systems characteristics, e.g., shrinking dominated processes, may favor application of other techniques. Furthermore, if interest is not only directed to integral quantities but also reconstruction of a representative NDF, EQMoM or MEA are the methods of choice. Here, trade-off between increased numerical effort and improving accuracy always has to be evaluated with respect to the potential application, e.g., when the approximate moment method is embedded in a model-based online process optimization, real-time requirements prohibit application of time-consuming numerical schemes.
In contrast to the well-filled cabinet of one-dimensional approximate moment methods, the review reveals a different picture for multidimensional PBEs which account for heterogeneity with respect to multiple internal particle properties: While a relatively large number of extensions of one-dimensional algorithms is found for bivariate problems, applications and suitable methods for genuine high-dimensional application are scarce. Here, CQMoM is established for moderate dimensional problems and first results indicate that DQMoM application to high-dimensional is promising. In comparison to methods for the mono-variate case, further work is required to study numerical effort and robustness. Furthermore, the reconstruction of multidimensional NDFs from their moments is still an open field.
Multiscale models for multicellular systems often characterize only average values despite the presence and significance of cell-to-cell variability for the system's dynamics. We reckon that the main reason for this lies in the lack of efficient and accurate numerical schemes for solution of more complex models that explicitly account for heterogeneities, e.g. high-dimensional PBEs. The availability of approximate moment methods enabling efficient and accurate numerical solution of PBEs could therefore facilitate application of PBM to biosystems. Thereby, an improved description of the multicellular population's dynamics would represent a convenient base for model-based process analysis, design, control, monitoring and intensification.