Diffusion processes satisfying a conservation law constraint

We investigate coupled stochastic differential equations governing N non-negative continuous random variables that satisfy a conservation principle. In various fields a conservation law requires that a set of fluctuating variables be non-negative and (if appropriately normalized) sum to one. As a result, any stochastic differential equation model to be realizable must not produce events outside of the allowed sample space. We develop a set of constraints on the drift and diffusion terms of such stochastic models to ensure that both the non-negativity and the unit-sum conservation law constraint are satisfied as the variables evolve in time. We investigate the consequences of the developed constraints on the Fokker-Planck equation, the associated system of stochastic differential equations, and the evolution equations of the first four moments of the probability density function. We show that random variables, satisfying a conservation law constraint, represented by stochastic diffusion processes, must have diffusion terms that are coupled and nonlinear. The set of constraints developed enables the development of statistical representations of fluctuating variables satisfying a conservation law. We exemplify the results with the bivariate beta process and the multivariate Wright-Fisher, Dirichlet, and Lochner's generalized Dirichlet processes.


Introduction and problem statement
We investigate the consequences of the unit-sum requirement on N > 1 non-negative continuous random variables governed by a diffusion process. Such mathematical description is useful to represent fluctuating variables, Y 1 , . . . , Y N , subject to the constraint Y α = 1. We are interested in stochastic diffusion models and statistical moment equations describing the temporal evolutions Y α = Y α (t) and their statistics. In particular, we study the consequences of the bounded sample space, required by the non-negativity of Y α and the unit-sum conservation principle, Y α = 1. A simple physical example is the mixture of different chemical species, represented by mass fractions 0 ≤ Y α ≤ 1 undergoing reaction in a fluid whose overall mass is conserved. Such mathematical problems also appear in evolutionary theory [1], Bayesian statistics [2], geology [3,4,5], forensics [6], econometrics [7], turbulent mixing and combustion [8], and population biology [9]. Mathematical properties of such random fractions are given in [10,11,12,13].
Mathematically, we are interested in the following question: What functions are allowed to represent the drift, A α , and diffusion, b αβ , terms of the system, governing the vector Y = (Y 1 , . . . , Y N ), if Y α ≥ 0, α = 1, . . . , N and must hold for all t. In Eq. (1) dW α (t) is a vector-valued Wiener process with mean dW α = 0 and covariance dW α dW β = δ αβ dt, see [14], and δ αβ is Kronecker's delta. If the components of Y satisfy the constraints in Eq. (2), we call the event Y realizable. A consequence of the constraints in Eq. (2) imposed on the stochastic system Eq. (1) is that for all t the following holds In other words, we are interested in expressions for A α and b αβ , what constraints they must satisfy in addition to Eq. (3), and how to implement them so that Eq. (1) produces realizable events, i.e., Y satisfies Eq. (2) for all t.
We study diffusion processes as: (1) they are mathematically simple vehicles for representing temporal evolutions of fluctuating fractions (of a unit) and their statistics, (2) they lend themselves to simple Monte-Carlo numerical methods [15], and (3) they serve as a starting point for representations of statistical moment equations if individual samples and joint probabilities are not required. The Markovian assumption [14] is made at the outset and jump contributions are ignored. We derive constraints for the drift and diffusion terms that assure that the modeled processes are realizable (i.e., produce non-negative variables that satisfy the unit-sum constraint) for any realization at all times. We address the problem of the functional forms of the drift and diffusion terms from three perspectives: (1) the Fokker-Planck equation for the probability density function, (2) the stochastic differential equations for the individual realizations, and (3) the evolution equations for the jointly coupled statistics.
The plan of the paper is as follows. §2 introduces the geometry of the multi-dimensional sample space within which realizations of fractions of a unit are allowed and discusses constraints that ensure realizable statistical moments. §3 develops the implications of realizability on diffusion processes governing fractions. §4 follows by developing realizability constraints on the time-evolutions of statistics. §5 surveys some existing realizable diffusion processes. A summary is given in §6.

Realizability due to conservation
The notion of realizability due to a conservation law constraint was introduced and defined by Eq. (2). We now discuss the consequences of realizability pertaining to the individual samples of the state space, §2.1, and of their statistics, §2.2.

The universal geometry of allowed realizations
The geometrical definition of the sample space is given in which the vector Y = (Y 1 , . . . , Y N ) is allowed if Eq. (2) is to be satisfied. This is used to derive constraints for stochastic diffusions and their moment equations in the subsequent sections.
A realization of the vector, Y, with coordinates Y α ≥ 0, α = 1, . . . , N , specifies a point in the multi-dimensional sample space. The union of all such points that satisfy is the space of allowed realizations, Eq. (2). For example, in representing mass fraction constituents of a substance, Eq. (4) restricts the possible components of Y to those that are realizable; those vectors that point outside of the allowed space are not conserved; if Eq. (4) is violated, spurious mass is created or destroyed. Mathematically, the geometry of allowed realizations is a simplex, the generalization of a triangle to multiple dimensions. For N variables the (N −1)-simplex is a bounded convex polytope, P, on 00 00 11  11  0  0  1  1  00  00  11  11  00000000  00000000  00000000 00000000  00000000   11111111  11111111  11111111 11111111 that satisfy non-negativity and the unit-sum constraint, Eq. (2). The boundary of the allowed region on the plane spanned by Y 1 and Y 2 is the closed loop of straight lines, , defined by Eq. (5). If the vector Y points inside the triangle, it is a realizable event.
the (N −1)-dimensional hyperplane; P is the convex hull of its N vertices. P's boundary, ∂P, is defined as the closed surface of non-overlapping hyperplanes of N − 2 dimensions, plotted in Figure 1 for N = 3. The domain (or support) of the joint probability, F (Y), with Y α , α = 1, . . . , N , is the (N −1)simplex. Of all Y α only N − 1 are independent due to Eq. (4) and without loss of generality we take The same geometry of allowed realizations is discussed by Pope in the N -dimensional state space in [16] in the context of ideal gas mixing in turbulent combustion. We confine our attention here to N −1 dimensions, as one of the variables is determined by the unit-sum requirement, Eq. (6). As a consequence, the (N −1)-dimensional geometry of realizable events is remarkably simple and universal: it is the bounded convex polytope whose boundary is defined by Eq. (5). Consequently, the realizability constraint, Eq. (2), uniquely and universally determines the realizable region of the state space: it is the same in all points in space and time for all materials undergoing any physical process that conserves mass, Eq. (4). The ensemble is realizable if and only if all samples reside inside the polytope given by Eq. (5). For N = 3 this means that the support of F is the triangle depicted in Figure 1.

Realizable statistical moments
If the fractions are non-negative and sum to one, Eq. (2), they are also bounded, whose consequences on some of their statistical moments are now discussed. Taking mathematical expectations of Eq. (7), see e.g. [17], yields Similar to the instantaneous fractions, the first statistical moments are also non-negative, bounded, and sum to unity. Since both the instantaneous variables and their means are bounded, fluctuations about the means are also bounded: As a consequence, the variances and the covariances are also bounded: Multiplying Eq. (4) by y β , β = 1, . . . , N and taking the expectation yield i.e., the row-sums and, due to symmetry, the column-sums of the covariance matrix are zero. Expressing y N y 1 , y N y 2 , etc., from the first N −1 equations of Eqs. (12) and substituting them into the N th one, yield the weaker constraint, Due to bounded fluctuations, Eq. (9), the third central moments are also bounded, and in general, for n ≥ 2 we have, Ensuring non-negativity and unit-sum puts constraints on possible time evolutions of Y = Y(t), represented by diffusion processes and that of their statistics. Some of these constraints are developed in the following sections.

Diffusion processes for random fractions
Implications of the geometry of the realizable state space, discussed in §2, on diffusion processes are developed. First, the relevant mathematical properties of Fokker-Planck equations are reviewed in §3.1, followed by the constraints on their functional forms, §3.2.

Review of some boundary conditions of Fokker-Planck equations
The discussion is restricted to Markov processes which by definition obey a Chapman-Kolmogorov equation [14]. Assuming that Y α are continuous in space and time, jump processes are excluded. The temporal evolution of random fractions, Y(t), constrained by Eq. (2) can then be represented most generally by diffusion processes whose transitional probability, F (Y, t), is governed by the Fokker-Planck equation, where A α and B αβ denote drift and diffusion in state space, respectively, and B αβ is symmetric non-negative semi-definite [17]. Eq. (17) is a partial differential equation that governs the joint (17) and is determined by Eq. (6). Augmented by initial and boundary conditions, Eq. (17) describes the transport of probability in sample space R whose boundary is ∂R with normal vector n α , see [14]. Eq. (17) can be written in conservation form as in terms of the probability flux, see [14], Sec. 5.1, Using Eqs. (18)(19) the following boundary conditions are considered, see [14], Sec. 6.2.
1. Reflecting barrier. If n α I α (Y, t) = 0 everywhere on the boundary, ∂R is a reflecting barrier: a particle inside R cannot cross the boundary and must be reflected there.
2. Absorbing barrier. If F (Y, t) = 0 everywhere on the boundary, ∂R is an absorbing barrier: if a particle reaches the boundary, it is removed from the system.
3. Other types of boundary conditions. Some part of the boundary may be reflecting while some other may be absorbing: a combination is certainly possible. We only consider reflecting and absorbing barriers -other types of boundaries are discussed in [18].
To support the forthcoming discussion, some well-established mathematical properties of multivariable Fokker-Planck equations have been reviewed.

Realizable diffusion processes
The implications of the realizability constraint, Eq. (2), on the functional forms of the drift and diffusion terms of the Fokker-Planck equation (17) are now investigated.
As discussed in §2, the region of the sample space allowed by the realizbility requirement is the polytope P defined by its boundary, ∂P, Eq. (5), in which all samples of Y = Y(t) must reside at all times. Consequently, the sample space, R, of the Fokker-Planck equation (17) must coincide with P, which constrains the possible functional forms of A α (Y, t) and B αβ (Y, t). In the following, these constraints are developed for binary (single-variable) processes first, followed by ternary processes, and then generalized to multiple variables.

Realizable binary processes: N = 2
The Itô diffusion process [14], governing the variable Y , with B(Y, t) ≥ 0 is equivalent to and derived from Eq. (17) with N = 2, see e.g., [14], For N = 2 the allowed space of realizations is a line with endpoints given by Eq. (5), This can be ensured if the drift and diffusion terms in Eqs. (20)(21) satisfy In other words, the realizability constraint, Eq. (2), on Eq. (20) mathematically corresponds to Eqs. (23)(24). A diffusion process, governed by Eq. (20), that satisfies Eqs. (23)(24), ensures that the fractions Y and 1 − Y satisfy 0 ≤ Y ≤ 1, provided each event of the ensemble at t = 0 satisfies 0 ≤ Y ≤ 1. The equal signs in the constraints on the drift in Eqs. (23)(24) allow for absorbing barriers at Y = 0 and Y = 1, respectively. The constraints on the diffusion term imply that B(Y, t) must either be nonlinear in Y or B(Y, t) ≡ 0 for all Y . In other words, since the diffusion term must be non-negative, required by Eq. (20), it can only be nonzero inside the allowed sample space if it is also nonlinear.

Realizable ternary processes: N = 3
For N = 3 variables, the unit-sum-constrained sample space and its boundary are sketched in Figure 1. In this case individual samples of the joint probability, F (Y 1 , Y 2 ), are governed by the system, The allowed samples space is two-dimensional (a triangle) whose boundary, defined by Eq. (5), consists of the loop of lines, For N = 3, the state vector, governed by Eqs. (25)(26) The realizability constraint, Eq. (2), on the system of Eqs. (25)(26) mathematically corresponds to Eqs. (28-31). The three fractions, Y 1 , Y 2 , and Y 3 = 1 − Y 1 − Y 2 , governed by Eqs. (6,25,26), remain fractions of unity if their drift and diffusion terms satisfy Eqs. (28-31). Naturally, an initial ensemble that satisfies, 0 ≤ Y 1 , 0 ≤ Y 2 , and Y 1 + Y 2 ≤ 1, is required. The constraints on the diffusion terms in Eqs. (28-31) show that both B 1 and B 2 must either be nonlinear in Y 1 and Y 2 , respectively, or B 1 (Y 1 , t) ≡ 0 and B 2 (Y 2 , t) ≡ 0, for all Y 1 and Y 2 , respectively. Furthermore, if one were to construct a process with A 1 = A 1 (Y 1 ), B 11 = B 11 (Y 1 ), and B 12 = B 12 (Y 1 ), then either A 2 or B 22 must be a function of both Y 1 and In other words, the unit-sum constraint couples at least 2 of the 3 fractions, governed by the system of Eqs. (6,25,26).
Constraints on the functional forms of the drift and diffusion terms of the multivariate Fokker-Planck equation (17), as a temporal representation of random fractions, Y α = Y α (t), have been developed. Eqs. (33-34) are our central result which ensure that sample space events, generated by Eq. (17) or its equivalent system of diffusion processes, Eq. (32), satisfy the realizability constraint at all times, provided the initial ensemble is realizable. Since Eqs. (17) and (32) govern N − 1 variables and Y N = 1 − N −1 β=1 Y β , the unit-sum requirement, Eq. (4), is satisfied at all times. An implication of Eqs. (33-34), exemplified in §5, is that random fractions represented by diffusion processes must be coupled and nonlinear.

Realizable evolution of the means: Y α
Multiplying Eq. (17) by Y γ and integrating over all sample space, see e.g. [19], yield the system of equations governing the means of the fractions, where

Realizable evolution of the second central moments: y α y β
Multiplying the Fokker-Planck equation (17) by then integrating over all sample space yields the equations governing the covariance matrix of the fractions, with A α = A α (Y, t) and B αβ = B αβ (Y, t). The right hand side of Eq. (37) is denoted by C αβ , the evolution rate of the covariance matrix. Eq. (37) shows how the covariances are governed if a Fokker-Planck equation (17) or a diffusion process (32) governs the underlying joint probability, e.g., y α y β is symmetric at all times. Following the development in §2.2, four conditions must be satisfied by the system of second moment equations (37) to ensure an evolution of the covariances that is consistent with the realizability constraint, Eq. (2): 1. Symmetric covariance evolution. The symmetry of the covariance matrix can be ensured if y α y β (t = 0) is symmetric, as well as its evolution rates: 2. Boundedness of the variances, Eq. (10). This condition can be ensured with lim y 2 α →0 C αα = lim y 2 α , t ≥ 0 and lim as the boundary of the state space is approached, indicating that in general the equation governing y 2 α must either be a function of y 2 α or C αα ≡ 0 for all t.
3. Boundedness of the covariances, Eq. (11). This condition can be ensured if for α = β, lim yαy β →−1 C αβ = lim yαy β →−1 y α y β , t ≥ 0 and lim as the boundary of the state space is approached, indicating that in general the equation governing y α y β must either be a function of y α y β or C αβ ≡ 0 for all t.
4. Zero row-sums, Eq. (12). Differentiating Eqs. (12) in time and using Eq. (37) yield the system C 11 + C 21 + · · · + y N y 1 , t = 0 Performing the same substitutions on Eqs. (41) that resulted in Eq. (13) we obtain the weaker constraint We see that the trivial specification, C αβ ≡ 0, satisfies all the above conditions, but also fixes the covariance matrix at its initial state for all t ≥ t 0 , which is of limited applicability.

Bounded evolution of the third central moments, y 3 α
Multiplying the Fokker-Planck equation (17) by y 3 γ = (Y γ − Y γ ) 3 then integrating yields the system governing the third central moments, y 3 α , as with A α = A α (Y, t) and B αβ = B αβ (Y, t). The right hand sides of Eqs. (43) are the evolution rates of the third moments, denoted by S α . The boundedness of the third moments, required by Eq. (14), can be ensured if as the boundary of the state space is approached, indicating that in general the equation governing y 3 α must either be a function of y 3 α or S α ≡ 0 for all t. The conditions in Eq. (44) only ensure boundedness, consequently, they are necessary but not sufficient conditions for realizability of the third moments as required by Eq. (2). Note that the requirement on bounded sample space has no implications on the boundedness of the skewness: since y 2 α ≥ 0, Eq. (10).

Bounded evolution of the fourth central moments, y 4 α
Multiplying the Fokker-Planck equation (17) by y 4 γ = (Y γ − Y γ ) 4 then integrating yields the system governing the fourth central moments, y 4 α , as as the boundary of the state space is approached, indicating that in general the equation governing y 4 α must either be a function of y 4 α or K α ≡ 0 for all t. The conditions in Eq. (47) only ensure boundedness, consequently, they are necessary but not sufficient conditions for realizability of the fourth moments as required by Eq. (2). Note that, similar to the skewness in Eq. (45), the requirement on bounded sample space has no implications on the upper bound of the kurtosis: since y 2 α ≥ 0, Eq. (10).

Summary on realizable statistics of fractions
The unit-sum constraint, Eq. (4), applied to a set of non-negative random variables, bounds and constrains their statistical moments, as shown in §2.2, as well as their time-evolutions. We examined the evolution of the moments, Y α , y α y β , y 3 α , y 4 α , and showed how they are governed if an underlying diffusion process is known.
Realizability of the means, as defined by Eq. (2), can be ensured if Eqs. (8) and (36) are satisfied. Realizability of the covariances can be ensured if Eqs. (10)(11)(12) and (38-41) are satisfied. Boundedness of the third moments is ensured by Eqs. (14) and (44), while boundedness of the fourth moments is ensured by Eqs. (15) and (47). The procedure outlined above can be continued to derive additional constraints for consistency of the third, fourth, mixed, and higher moments with the unit-sum constraint. The constraints reflect the coupled and nonlinear nature of random fractions, both as instantaneous variables and their statistics.

A survey of realizable diffusion processes
A survey of existing diffusion processes that satisfy the realizability constraints on drift and diffusion on the state-space boundary, Eqs. (33-34), is now given.

Realizable binary process, N = 2: Beta
An example for N = 2, satisfying the realizability constraints on the drift and diffusion terms on the sample-space boundary in Eqs. (23)(24), is given in [20], specifying the drift and diffusion as yielding the stochastic differential equation, with b > 0, κ > 0, and 0 < S < 1 excluding, while with 0 ≤ S ≤ 1 allowing for absorbing barriers. In Eq. (50) the drift is linear and the diffusion is quadratic in Y . The invariant distribution of Eq. (50) is beta, which belongs to the family of Pearson distributions, discussed in detail by Forman & Sørensen [21]. Of the special cases of the Pearson diffusions, discussed in [21], only Case 6, equivalent to Eq. (50), produces realizable events. A symmetric variant of Eq. (50) was constructed in [22], which does not allow a non-zero skewness in the statistically stationary state, see [20].

Realizable multivariate process, N > 2: Dirichlet
Another process that satisfies Eqs. (33-34), developed in [24], specifies the drift and diffusion terms as resulting in the system of stochastic differential equations, with parameter vectors b α > 0, κ α > 0, and 0 < S α < 1, and Y N given by Eq. (6). Eq. (55) is also a generalization of Eq. (50) for N > 2 variables. The invariant distribution of Eq. (55) is also Dirichlet, provided the parameters of the drift and diffusion terms satisfy Note that while there is no coupling among the parameters, ω α , of the drift and diffusion terms in the Wright-Fisher Eq. (52), the parameters, b α , S α , and κ α , of Eq. (55) must be constrained by Eq. (56) to keep its invariant distribution Dirichlet.

Summary
We have demonstrated that the problem of N fluctuating variables constrained by the unit-sum requirement can be discussed in a reduced sample space of N − 1 dimensions. This allows working with the unique, universal, and mathematically well-defined realizable sample space which produces samples and statistics consistent with the underlying conservation principle.
We have studied multivariate diffusion processes governing a set of fluctuating variables required to satisfy two constraints: (1) non-negativity and (2) a conservation principle that requires the variables to sum to one, defined as realizability. Our findings can be summarized as follows: • The diffusion coefficients in stochastic diffusion processes, governing fractions, must be coupled and nonlinear.
• Boundedness of the sample space requires boundedness of the moments.
The constraints provide a method that can be used to develop drift and diffusion functions for stochastic diffusion processes for variables satisfying a conservation law and that are inherently realizable.