A provably convergent control closure scheme for the Method of Moments of the Chemical Master Equation

In this article, we introduce a novel moment closure scheme based on concepts from Model Predictive Control (MPC) to accurately describe the time evolution of the statistical moments of the solution of the Chemical Master Equation (CME). The Method of Moments, a set of ordinary differential equations frequently used to calculate the first nm moments, is generally not closed since lower-order moments depend on higher-order moments. To overcome this limitation, we interpret the moment equations as a nonlinear dynamical system, where the first nm moments serve as states and the closing moments serve as control input. We demonstrate the efficacy of our approach using three example systems and show that it outperforms existing closure schemes. For polynomial systems, which encompass all mass-action systems, we provide probability bounds for the error between true and estimated moment trajectories. We achieve this by combining convergence properties of a priori moment estimates from stochastic simulations with guarantees for nonlinear reference tracking MPC. Our proposed method offers an effective solution to accurately predict the time evolution of moments of the CME, which has wide-ranging implications for many fields, including biology, chemistry, and engineering.


Introduction
The dynamics of well-stirred chemical reaction systems are described by the Chemical Master Equation (CME).The CME is a set of linear differential equations for the timedependent probability distribution over the set of all possible configurations of the system.Since this set is usually large or even infinite except for very small systems, solving the CME directly is often not feasible.One remedy for this problem is to consider only the first n m moments of this probability distribution, thereby reducing the number of equations from the number of configurations to the number of considered moments.Two established methods to estimate the statistical moments of the CME solution are the Stochastic Simulation Algorithm (SSA) and the Method of Moments (MoM).The SSA (1,2) generates sample paths from the underlying Markov process, which can be used to approximate either the full distribution or its moments via Monte Carlo integration.SSA-based mo-ment approximation has several advantages.First, it inherently produces physically meaningful moment estimates, for example, means and variances are non-negative.Second, empirical moment estimates are consistent and thus converge in probability to the true moments in the large sample limit.On the other hand, approximating the CME solution or its moments only based on SSA samples is often tedious, as the computational load of generating each sample grows with the number of simulated molecules.Moreover, a huge number of samples is needed for an accurate estimate if, for example, some species have a high probability of becoming zero, if the system is stiff, or if it includes rare events.The MoM is an infinite set of coupled Ordinary Differential Equations (ODEs) which describes the evolution of the statistical moments of the CME solution.Mathematically, the equation for each moment is obtained by multiplying the CME with a moment-specific test function and summation over all configurations, as described in (3).The solution of this MoM ODE system can be obtained with state-ofthe-art solvers and the needed computational power is nearly independent of the number of simulated molecules.However, one challenging aspect of the MoM ODEs is their interdependence structure.For systems with quadratic propensities, a moment of a given order depends on the next higherorder moment, leading to an infinite dependency hierarchy, and the set of equations for the first n m moments is not closed.It is worth mentioning that closing the MoM equations for the first n m moments with the true higher-order moments, would they be known hypothetically, results in exact moment estimates.Plenty of solution approaches have been proposed to solve this closure problem.We will briefly outline some of them and refer the interested reader to (4) for an in-depth review of the general moment closure problem.The most trivial closure technique simply neglects the influence of all higher-order moments on the considered lowerorder moments and is thus effectively a truncation of the ODE system.We, therefore, refer to it as Truncation Closure (TC).TC is trivially implementable without any additional manual or computational effort.Further closure techniques have been introduced for stochastic differential equation systems (5), including the cumulant closure technique that was first proposed by Cramer (6).Cumulants are, like moments, characteristic quantities of a distribution that can be calcu-lated via its characteristic function.Similar to TC, the cumulant closure technique closes the cumulant equations by setting all higher-order cumulants to zero.Since the first two cumulants correspond to expectation value and variance, respectively, cumulant closure is equivalent to TC if only the dynamics of the first two moments are considered.
Beyond these naïve approaches, a popular class of methods assumes that the CME solution is similar to a specific family of distributions and uses the higher-order moments of these respective distributions for closure (7)(8)(9)(10)(11)(12).Other approaches express higher-order moments as potentially nonlinear functions of the considered lower-order moments (13)(14)(15)(16)(17)(18).A third interesting class reconstructs the underlying distribution at one or several time steps, e.g., by exploiting the maximum entropy principle, and uses the reconstruction to infer the closure moments (19)(20)(21).
Many of the proposed closure schemes work very well in practice.For special cases, there are even theoretical results justifying specific closure schemes (22)(23)(24).However, since there is generally no relation between the trajectories of the closed moment equations and the CME, both are decoupled.Hence, moment closure methods suffer from undesired artifacts such as unphysical moments and potentially large deviations from the true moment trajectories, which constitutes the major drawback of closure methods.Even worse, except for systems with only reactions of up to order one, for which the moment equations are naturally closed, we do not have criteria at hand which allow for an a priori judgment of the quality of a moment closure solution based on the set of chemical reactions, the initial species amounts and number n m of considered moments.Therefore, a closure scheme that generally and provably closes the MoM ODEs in accordance with the true CME solution would be highly desirable since it solves this major problem.
In this article, our aim is to address this decoupling issue by providing a closure that guarantees that the distance between the moment closure solutions and the true moments is bounded and moments are physically plausible.To achieve this ambitious goal, we employ a few SSA samples that are responsible for coupling the closed MoM ODEs to the true CME solution.This idea is related to the work in (25) that uses a few SSA trajectories to roughly estimate the closure moments and subsequently uses Kalman filtering to eradicate the noise.While their proposed closure method relies on sampled closing moments that are typically affected by nonnegligible noise, we propose to obtain the closing moments as a control input.To this end, we interpret the closed moment equations as a nonlinear system, in which the moments up to order n m serve as states and the closing moments serve as a control input.For this Control Closure (CC) scheme, we can build on a variety of controller design methods for nonlinear systems from the literature (26).Particularly for the class of polynomial dynamics, which contains for example all mass-action systems, a popular approach for nonlinear control relies on sum-of-squares optimization, which reformulates non-negativity conditions as semi-definite programs (27,28).However, sum-of-squares optimization is generally computationally challenging.Instead, we rely on Model Predictive Control (MPC) (29,30) and, in particular, use concepts of nonlinear reference tracking MPC as presented in (31,32).For CC, we characterize SSA-based estimates and confidence intervals of the moments by exploiting concentration inequalities of random variables.We then combine these with the MoM by formulating an optimal control problem in which the resulting optimal control input serves as the proposed moment closure.We validate our approach on three example systems, for which we obtain very good performance compared to TC and stand-alone SSA estimation.In particular, the higher-order moment estimates obtained by our CC scheme are indeed close to the true moments.While there are other more sophisticated closure techniques, e.g., Gaussian closure (13), as discussed previously, we restrict ourselves to the comparison to TC as the simplest closure method at hand which is also widely used in practice.The particular contribution of the present paper is not limited to the proposed generally applicable CC and its sheer comparison with existing closure techniques but includes a theoretical analysis of the proposed closure.In particular, for polynomial systems, we provide guarantees that the moment trajectories inferred from CC satisfy the moment closure equations and lie in the given confidence intervals while staying close to the sampled moment estimates.Notation: Throughout this paper, bottom right indices • i denote entries of vector-valued quantities.Samples for example from the SSA are indexed using a top left index in parentheses (i) •.To sum over different chemical reactions or system configurations, we employ top right indices in parentheses • (i) .For a symmetric matrix A, we write A ≻ 0 (A ⪰ 0) if A is positive (semi-)definite.We write ∥x∥ A for the weighted norm x ⊤ Ax of a vector x with A ≻ 0. By N 0 and R ≥0 we denote the natural numbers including zero and all nonnegative real numbers, respectively.We use the vectorial index α ∈ N n 0 with |α| = α 1 + ... + α n to define for a vector x ∈ R n the monomial We denote the class of functions α : R ≥0 → R ≥0 , which are continuous, strictly increasing, and satisfy α(0) = 0 by K. We write K ∞ for the set of functions α ∈ K satisfying lim k→∞ α(k) = ∞.By L we denote the class of functions β : R ≥0 → R ≥0 , which are continuous and decreasing with lim k→∞ β(k) = 0.A continuous function β : R ≥0 × R ≥0 belongs to class KL if β(•, s) ∈ K for each fixed t and β(r, •) ∈ L for each fixed r.

SSA trajectories as control reference for moment closure
In this section, we present a CC scheme to close the moment equations of the CME.Subsection 2-A first introduces an interpretation of the closed moment equations from the perspective of control theory.Second, we derive probability bounds for SSA-based moment estimation in Subsection 2-B, which define a confidence region that contains the true moments with a certain user-defined probability.Subsequently, we use these SSA-based estimates and probability bounds in Subsection 2-C to formulate an optimization problem whose solution is the best possible estimate of the moments accord-A System theoretic interpretation of the closed moment equations ing to a predefined optimization criterion, while satisfying the closed moment equations, staying close to the SSA-based moment estimates, and lying inside the confidence region.Since this optimization problem tends to be computationally demanding for large time intervals, we present a shorthorizon tracking MPC scheme for which we develop theoretical results leading to a provably convergent CC for optimal estimation of the true CME moments in Section 2-D.
A. System theoretic interpretation of the closed moment equations.The CME is a set of linear ordinary differential equations that describes the time-dependent probability p(z, t) for each system configuration z as with reaction propensities a (j) (z), configuration change vectors ν (j) , and sums over all feasible reactions.The set Ω Z of feasible configurations depends on the set of reactions and the initial molecule numbers.For a simple reversible reaction A ⇌ B and a single molecule, for instance, Ω Z = z (1) = (1, 0), z (2) = (0, 1) .Moment equations are obtained by multiplication with a moment-specific test function T (z) : Ω Z → R and summing over all configurations (3), i.e., For example, for T (z) = z i and , respectively, we get for the expectation E Z i (t) and the covariance k denoting the k-th component of the vector ν (j) (3).The general moment equations Eq. (1) can be interpreted as a controlled system of the form where x ∈ R n contains the first n m moments that appear as states and u ∈ R m defines the control input and corresponds to the higher-order moments that are needed to close the system.Further, we write x u (t; x 0 ) for the solution of Eq. ( 2) for a given control sequence u and initial condition x 0 .The proposed CC scheme is in principle valid for arbitrary dimensions and moment orders, and we later also apply it to a system with two chemical species to show its broad applicability.However, for ease of notation, we restrict the following derivations to systems with only one chemical species (n c = 1) and refer the interested reader to Engblom (3) for a general derivation of the MoM.We are interested in predicting the expectation value and variance (n m = 2), i.e., In this case, a second order Taylor series expansion of In this representation, the skewness is the closure moment.
For mass-action systems with reactions of at most quadratic propensities, the MoM ODEs are at most quadratic in the moments, i.e., there are no products of more than two moments.Moreover, the moments of order (n m + 1) affect only the n m -th moments and occur only linearly in the corresponding ODE (3).In this case, the Taylor expansion Eq. ( 3) is exact and can be used for further specification of the dynamical system Eq.( 2) to where for all j = 1, ..., n.For n c = 1 and n m = 2, the dynamical systems Eq. ( 2) and Eq. ( 4) are of dimension n = 2, the closing moment u is of dimension m = 1, and the matrix B in Eq. ( 4) reduces to a vector, i.e., B = 0 With this simplification, we finally have a two-dimensional controlled system for the expectation value E Z (t) and the variance V Z (t), which reads where f 1,2 (x 1 , x 2 ) are polynomial functions in their arguments.In the following, we use the derived dynamical system representation to formulate an optimization-based closure for the MoM.To this end, we rely on SSA to obtain rough moment estimates together with probabilistic bounds on their accuracy.

B. Probability bounds for SSA-based moment estimation.
Empirical moment estimates (N ) ÊZ (t) from SSA sample paths (i) z(t), i = 1, . . ., N , are consistent and converge for each time instance t * in probability to the true moment E Z (t * ) for large sample sizes, i.e., plim In an attempt to quantify this convergence for a finite number of SSA trajectories N , we are interested in a probabilistic bound for (N ) ÊZ (t * ) − E Z (t * ) , or in other words, confidence intervals of the mean E Z (t * ).Estimating these intervals is one of the most fundamental questions of frequentist statistics and consequently, the amount of related literature is vast.Ideally, we would like to present a probabilistic bound of the form , where the bounding function B only depends on the desired width a > 0 of the confidence interval, the number of drawn samples N , and the samples themselves.However, it seems to be impossible to obtain such a purely sample-based bound B. This subsection, therefore, outlines the most prominent approaches together with their respective additional assumptions and refers the interested reader to relevant sources for more information.We would like to stress that the proposed CC scheme does not rely on any particular bound.Thus, every applicant of our technique is able to choose bounds according to each individual situation.One of the most prominent representatives of this class of inequalities is the Chebyshev-Inequality (CI) for all a ∈ R, or closely related variants thereof (33).However, its bound depends on the unknown true variance V Z (t * ).
In a similar fashion, the Central Limit Theorem (CLT) can be combined with the Berry-Esseen theorem (34,35) to show which is detailed in Appendix B. Here, c be is the Berry-Esseen constant that is known to be upper bounded by 0.4748 (36), erf denotes the error function, and ρ is the true third centralized absolute moment Again, knowledge of the true variance V Z (t * ) and the justdefined moment ρ is necessary to obtain the desired bound.
In contrast, Hoeffding's Inequality (HI) (37) provides a bound independent of higher-order moments but instead requires bounded states 0 While this boundedness is a strong assumption, the number of system molecules is often bounded by conservation laws anyway.In case the system can indeed grow infinitely, a configuration state truncation can provide such an upper bound b, which is, however, not uniquely defined.Both named possibilities will be demonstrated in our example systems.By defining an upper bound α ∈ (0, 1) for the probability in one of the inequalities Eq. ( 6), Eq. ( 7), and Eq. ( 8), each inequality can be resolved for the parameter a.In case of the HI, this leads to exemplary, and the true moment lies with probability 1 − α in the symmetric interval around the SSA-estimate, which serves as (1 − α) confidence interval.In addition, physical constraints such as the non-negativity of expectation and variance might be used to tighten these intervals even further if desired.We denote the union of confidence intervals for all t * in a given time interval T := [0, T end ] with T end > 0 as confidence region, which describes a symmetric tube around (N ) ÊZ (t) for t ∈ T .For a bounded state space, also the variance could be constrained.While we restricted this subsection to bounds for E Z (t * ), many results can be generalized to higher-order moments and similarly included in the CC scheme presented later.For these and many other related results, we can (upon many others) recommend the textbook of (38) and the articles of ( 33), (39), and (40).
C. Optimal CC scheme for CME moment estimation.Based on the system theoretic representation Eq. ( 4) of the moment equations as a dynamical system, we formulate a CC scheme that relies on a priori SSA stand-alone moment estimates and confidence regions derived in Subsection 2-B and uses MPC concepts to find an optimal input u.In particular, we formulate an optimization problem to find a control input u which guarantees that the MoM with CC leads to estimated moments that minimally deviate from the SSA-based moment estimation (N )  ÊZ (t) with respect to a specified cost function and are guaranteed to lie within the predefined confidence region.A graphical overview of the proposed CC is depicted in Figure 1.Given SSA trajectories (i) z(t) over the finite time interval T , we obtain estimates of the first (n m + 1) moments.We collect the first n m moments in x(t) ∈ R n and all moments of order (n m + 1) in û(t) ∈ R m for t ∈ T .Where applicable, we then define the constraint sets X (t) ⊆ R n and U(t) ⊆ R m as the SSA-based confidence regions in Subsection 2-B.We want to ensure (x(t), u(t)) ∈ X (t) × U(t) for all t ∈ T .Both x and û as well as X (t) and U(t) influence the MPC controller, which is explained in the following and returns moment estimates x(t) and ū(t) that lie within X (t) and U(t).As we seek the optimal moment estimation on the interval T , we need to define an optimality criterion first.To this end, we define the open-loop cost function with the stage cost (10) for some user-defined positive definite Q ∈ R n×n , positive semi-definite R, S ∈ R m×m , and sampling period δ > 0.
These matrices can be chosen, e.g., to weigh some moments more than others.If there is no estimate û(t) of the closing moments, we can just choose R = 0 or û(t) = 0.The last term in Eq. ( 10) is used to penalize the rate of change of the resulting input u.The matrix S can be seen as a tuning parameter: For a larger weight S, the difference in u at two consecutive time points t − δ and t is penalized more, and, hence, the optimization yields a smoother trajectory for the control input.Definition 1: For weights Q ≻ 0, R ⪰ 0, and S ⪰ 0 and an initial condition x(0) = x 0 , the optimal reachable trajectory (x * (t), u * (t)) on the interval T is the minimizer of x(0) = x 0 .This definition characterizes the best possible moment estimation of the CME according to the given SSA samples and the cost function in Eq. ( 9) while additionally satisfying the closed moment equations and staying inside the confidence region.Hence, the solution u * (t) is the optimal CC for the moment equations on T .D. Guaranteed convergent CC using a short-horizon optimization.For an increasing size of the interval T , the nonlinear optimization problem in Eq. ( 11) becomes numerically challenging.To overcome this issue, we present in the following a short-horizon optimization based on MPC which is guaranteed to converge to the optimal reachable trajectory (x * (t), u * (t)) without explicitly calculating it.In particular, we use an MPC scheme based on the results in (32), where the authors propose a reference tracking MPC scheme without terminal ingredients for unreachable reference trajectories of discrete-time nonlinear systems based on techniques from economic MPC (41,42).This MPC scheme ensures practical stability of the unknown optimal reachable reference (x * (t), u * (t)).Based on this, the results were translated to the case of continuous-time nonlinear systems in (31), which makes it suitable for our setting.
To obtain an optimal estimate of the CME moments which is close to the SSA-based estimation, we want to track the optimal reachable trajectory Eq. ( 11) without explicitly computing (x * (t), u * (t)).To this end, we define the short-horizon cost function and propose the following: At each time t ≥ 0, given the current state x(t) and an admissible reference (x(t), û(t)) ∈ X (t) × U(t), we solve a short-horizon MPC scheme characterized by the optimization problem with the prediction horizon T satisfying δ ≤ T ≤ T end .
The optimal trajectory resulting from Eq. ( 12) is denoted by (x * (•; t), ū * (•; t)), while the control input u is chosen as The MPC problem Eq. ( 12) is solved repeatedly after each sampling period δ.The closed-loop dynamics using the MPC controller read then Due to the finite time horizon T , the solution to Eq. ( 12) does in general not exactly follow the optimal trajectory (x * (t), u * (t)).While the application of our proposed CC based on Eq. ( 12) is generally possible and requires only the definition of some hyperparameters, we can even establish, under mild regularity conditions, recursive feasibility of the MPC problem Eq. ( 12).Definition 2: The MPC problem Eq. ( 12) is called feasible at time t if there exists at least one ū(•; t) such that the constraints are satisfied.Moreover, the MPC problem Eq. ( 12) is called recursively feasible if initial feasibility at time t = 0 implies feasibility of the MPC problem for all t ≥ 0. Further, we show that the resulting closed-loop system practically converges to the optimal reachable reference (x * (t), u * (t)) based on results in (31).To this end, we introduce the notion of practical convergence.Definition 3: A system practically converges to a reference trajectory (x * (t), u * (t)) under control input µ(t) ∈ U(t) if there exists a β ∈ KL such that the corresponding system state x satisfies ∥x µ (t; x 0 )∥ ≤ β(∥x µ (t; x 0 ) − x * (t)∥, t).Moreover, the proposed CC recovers the true moments in the limit of infinitely many sample trajectories, which follows from the consistency of the empirical moments estimates.Before stating our main theoretical result, we make the following assumption.Assumption 4: System Eq. ( 4) representing the MoM is uniformly suboptimally operated off the trajectory x * , i.e., the following conditions hold (cf.(42, Def.12)): holds, where U ∞ (x 0 ) denotes the set of all feasible control sequences of infinite length starting at a given x 0 ∈ X , and X 0 denotes the set of all x ∈ X such that U ∞ (x) is non-empty.
• There exist δ > 0 and d ∈ K ∞ such that for each δ > 0 and each ε > 0 there exists R ε,δ ∈ R ≥0 such that δ/R ε,δ ≥ d(ε) for all δ ≥ δ and such that for each x 0 ∈ X 0 , each u ∈ U ∞ (x 0 ), and each T ∈ R ≥0 at least one of the following two conditions hold: 1.
, where If the minimizer (x * (t), u * (t)) is unique, then Assumption 4 is not restrictive due to the positive definite weight matrix Q in the stage cost ℓ (compare (32, Dec. IV-D)).The following result establishes the theoretical properties of the proposed CC scheme.Theorem 1: Suppose Assumption 4 holds.Then, there exists a T0 such that for all T > T0 the CC ū(t) based on the MPC optimization Eq. ( 12) leads to a moment estimation that practically converges to the optimal MoM-based estimates Eq. ( 11) of the CME and is guaranteed to satisfy physical properties and the given confidence region of the moments.Moreover, the CC recovers the true moments for infinitely many sample trajectories, i.e., N → ∞.
Sketch of Proof: While we omit the proof details here and rigorously prove the theorem in Appendix C, we sketch the main ideas in the following.In particular, the proof relies on practical convergence according to Definition 3, where we use available results from the literature based on Lyapunov stability.To this end, we first exploit that the nonlinear dynamics of the MoM in Eq. ( 4) are, e.g., polynomial, twice differentiable, and structured such that only the moments of order (k + 1) enter (linearly) into the dynamics of the moments of order k.Thus, the considered dynamics inherently satisfy necessary assumptions from the literature to guarantee practical convergence.Moreover, the MPC optimization recovers the true moments in the limit N → ∞ since the SSA-based moment estimation converges to the true CME moments for N → ∞.The proposed optimal estimation of the CME moments based on the presented CC scheme is summarized in Algorithm 1.Note that the lower bound T0 in Theorem 1 is, in general, unavailable since it depends on existent but typically unknown functions.Hence, the theorem shows the existence of a sufficiently long prediction horizon without stating the horizon explicitly.Nevertheless, for a fixed sampling interval δ, the closed-loop system stays around the optimal trajectory Algorithm 1 Control closure scheme for MoM

2:
Set the MPC internal variables by allocating the first n m moments in x(t) and the moments of order (n m + 1) in û(t) for t ∈ T .3: If applicable, define the confidence regions described in Section 2-B as constraint sets X (t) and U(t).Short-horizon optimization: 4: Initialize the estimate x(0) = x(0).5: for k = 0, ..., (T end − T )/δ do 6: Update current sampling time t k = δk.

7:
Initialize current estimation x(t k ; t k ) = x(t k ). 8: Solve Eq. ( 12) to obtain the optimal prediction Obtain the optimal moments x(t) for t ∈ [t k , t k + δ) via Eq.( 13).11: end for (x * (t), u * (t)), which is as close as possible to the unreachable reference (x(t), û(t)), given that a large enough prediction horizon is used.We stress further that when δ → 0, T should go to infinity (31).Possible future work could use ideas from (43), which presents a nonlinear MPC scheme without terminal ingredients for tracking possibly unreachable setpoints using artificial steady-states (44).In particular, the in (44) proposed results use less restrictive assumptions and might be exploited to indeed derive computable bounds on the necessary prediction horizon, which is beyond the scope of this work.

Simulation examples
In this section, we validate the performance of our CC scheme on three example systems.All considered systems converge to a steady-state distribution on their respective configuration-transition graph, as discussed in (21).Algorithm 1 was implemented in Matlab using MPCTools (45) with its interface to the nonlinear optimization software CasADi (46).System 1 The first exemplary system consists of a constant production rate and quadratic degradation and is therefore one of the simplest systems involving quadratic propensities.The reactions read and the MoM ODEs are given by In principle, there is no limit for the number of molecules in System 1 and the space of all possible system configurations Ω Z is, therefore, N 0 .In order to still solve the CME numerically, we truncate Ω Z to 31 system configurations, leading to Ω Z = {0, ..., 30}.While this seems to be a very drastic measure, the probability mass not captured by the first 31 configurations is indeed negligible in our setting, which justifies this approximation.A comparison of the performance of different approaches for moment trajectory estimation for System 1 is shown in  were set for the expectation value by using HI and CI bounds.Please note that the CI bounds require knowledge of the true variance V Z to be valid.As this quantity is typically not available, the CI bounds depicted in Figures 2, 3, and 4 are calculated using the SSA-based estimator (10)  VZ of the true variance and therefore only approximate the true CI bounds.The variance is restricted to non-negative values and the skewness is not constrained in this setting, i.e., X (t) = X E (t) × R ≥0 and U(t) = R.The MPC problem is formulated using the weights Q = diag{1000, 300}, R = 1, S = 100 with the prediction horizon T = 1 s and the sampling interval δ = 0.1 s.We emphasize that the choice of the weights is not restricted to the values above but also different choices lead to satisfying results.Note that our choice is motivated by the sample accuracy of the SSA estimates, i.e., the variance is roughly weighted by the square root of the weight of the expectation value.Since we do not obtain a reliable estimate for the skewness, the weight R is chosen to be small.However, we use the weight S = 100R to nalize large deviations of the skewness between subsequent time steps to guarantee a certain smoothness.We stress that the absolute values of the weights can be scaled and only the relative relations between Q, R, and S are important.The two first estimated moments resulting from the proposed CC scheme, expectation and variance, are smooth and closely resemble the CME solution.Moreover, the expectation value clearly lies within X E (t) for both bounds.The skewness deviates from the reference course for the first 2.5 seconds but approaches the true solution quite well afterward.This becomes especially obvious when considering the squared quadratic error between the true CME solution and its approximations in the bottom sub-graph.Not only is the error of the controlled MoM solution orders of magnitude smaller than the error of the underlying 10-sample SSA approximation, but it is also comparable in magnitude with the SSA approximation for N = 100 samples while being significantly smoother at the same time.Further, we note that other closure methods than TC or a higher number of SSA samples might lead to successful estimates of the true CME solutions as well.However, we emphasize that the proposed CC is indeed able to estimate the true moments by using only N = 10 SSA samples.Additionally, the estimates based on the CC are guaranteed to satisfy the specified SSA-based bounds on the respective moments ensuring proper physical correctness, which is typically not guaranteed for other closure methods.System 2 Our second example represents a reversible dimerization reaction The configuration space of System 2 is limited by a conservation law, 2A 2 (t) + A(t) = c A = 2A 2 (0) + A(0), which can be used to reduce the system to species A only.The MoM ODEs are given by For this simulation example, we chose an initial number of 40 molecules of species A and c A = 50, which directly implies the existence of 5 A 2 molecules for t = 0 and Ω Z = {0, ..., 50}.Moreover, we use the same constraint sets X (t), U(t) and the same optimization parameters Q, R, S, T , δ as for System 1.A comparison of results for different approaches is presented in Figure 3 analogously to Figure 2. First and foremost, courses of the true CME moments are extremely well captured by both moment closure approaches, and the expectation value trajectory again lies well within X E (t).The reason for the similarity of results of both closure methods is that the control input skewness estimated by the CC scheme is quite close to zero during the whole simulation.In contrast, the SSA-based approximation is as shaky as for System 1.
System 3 The third exemplary system is chosen to demonstrate the practical applicability of our CC method.To this end, we choose a system that models the self-regulation of the expression of gene A that encodes protein B (47).More specifically, A f is the number of free DNA strands, A b is the number of DNA strands with a bound protein and B is the number of free proteins.The governing reactions read A f Similar to System 2, the three molecule species can be encoded by a two-dimensional state Z = (A f , B) with the help of a conservation law, as the sum of free and bound DNA is constant, i.e., A f (t) + A b (t) = c A = A f (0) + A b (0).It is practically unrealistic to close the MoM using potentially all four skewness terms.Thus, the MoM in this application example predicts the expectation values of both chemical species and uses the covariance V Z 1 ,Z 2 (t) as the closure moment (n m = 1), but we note that the CC is not restricted to this simplification.Indeed, due to the nature of the involved reactions, the variances V Z 1 ,Z 1 (t) and V Z 2 ,Z 2 (t) of both individual species do not have an impact on the respective expectation values, as they do not arise when deriving the MoM equations The number of proteins B is generally unbounded.Similar to our treatment of System 1, we truncate Ω Z 2 to 76 system configurations leading to Ω Z 2 = {0, ..., 75}.Again, this truncation does not neglect any probability mass beyond numeric tolerances.The results depicted in Figure 4 are created with initial numbers of Z(t = 0) = (5, 0) molecules and c A = 10.
As we did for Systems 1 and 2, we used sample-based reference trajectories for the expectation value and the variance.Please keep in mind, that in contrast to the previous systems, this means that we also employ a reference trajectory for the control input.Moreover, we use the confidence regions X (t) = [max(0, x(t) − a(α = 0.95)), x(t) + a(α = 0.95)] for both expectation values, similar to X E for System 1 and System 2, and U(t) = R for the covariance.The parameters for the MPC optimization are chosen as Q = diag{1000, 1000}, R = 30, S = 100 with prediction horizon T = 1 s and sampling interval δ = 0.1 s.We note that the weight matrices are strongly related to the ones used for System 1 and System 2. Again, the weight R = 30 of the covariance is roughly the square root of the weights Q 11 = Q 22 = 1000 of the expectation values.Similar to the previous figures, Figure 4 visualizes the performance of different approaches for moment trajectory estimation for System 3. In contrast to the previous figures, however, we now depict two expectation values in the top plots, each with its own confidence regions X (t).The third sub-graph again depicts the time course of the closing moment, which in this case is the covariance.Lastly, the bottom sub-graph shows the natural logarithm of the 2-norm of the error of the expected numbers of molecules E Z (t).The results arising from this gene self-regulatory system align with previous results: TC and the proposed CC scheme provide smooth solutions of the MoM that are in addition very accurately resembling the true CME moments.As guaranteed by construction, the CC solution stays within its specified bounds and converges to the true CME solution as T → ∞.

Conclusion
Applying the MoM to approximate the CME's moments is a challenging task due to the interdependence structure of the considered moments.In this work, we interpreted the system of MoM ODEs in the framework of control theory to present a non-heuristic closure scheme that provably stays close to the true solution in a probabilistic sense.More specifically, our CC scheme draws a few SSA samples from the true CME solution and uses the corresponding samplebased moment estimates as a physically meaningful reference trajectory.The closure moments are then used as the system's control input in order to follow the SSA trajectories while still obeying the MoM dynamics.Indeed, the CC leads to physically meaningful results that closely approximate the true moments even when TC fails to do so.Moreover, it outperforms purely sample-based approximation techniques as it smooths the sample averages according to the MoM dynamics.The MoM and hence basically all moment closure approaches require the considered moments of the CME solution to exist for the time interval of interest, which we have assumed (and numerically tested) in this work.However, it should be noted that this is a delicate issue since reaction systems with solutions whose moments are not defined can easily be constructed (48).These systems can show wildly stochastic behavior, and in particular, SSA-based moment estimates do not converge for increasing sample sizes.Moment-based methods terribly fail when applied to such systems without caution.Generally, this restriction has to be kept in mind when dealing with the MoM and related approaches.We demonstrated the capabilities of CC using two representative test-bed systems and a gene selfregulatory reaction network.The latter emphasizes that our closure scheme can in principle be applied to systems with several molecular species.In addition, it is straightforward to consider higher-order moments, as indicated by the more general dimensions at several places in the manuscript.However, a severe limitation is that the number of moments scales very badly with the number of molecular species in the system.For n c chemical species, we have n c expectation values, an n c × n c covariance matrix, n 3 c skewness values, and the number grows likewise for higher moments.This is true for other moment approaches as well, but also affects the computation time in our case via repeated optimization.Finally, we note again that a particular strength of our method compared to other existing closure schemes is that we are able to provide probabilistic guarantees that the calculated moments stay in a neighborhood of the true moments and hence cannot arbitrarily diverge.Moreover, our scheme provides several tuning points to control this distance.First, the number of SSA samples correlates negatively with the size of the confidence intervals.Second, different weights can be given to different moments, and the smoothness of the input can be regulated via hyperparameters in the MPC objective function.Finally, computational effort can to some extent be regulated by appropriately choosing the horizon in the short-horizon MPC scheme.Overall, we believe that CC has great potential to be widely applied to describe the moments of biochemical reaction systems.
where Φ is the standard-normal cumulative distribution.We derive the probabilistic bound by first reformulating the quantity of interest and using the definition of F N .Subsequently, Theorem 2 is employed to achieve an upper bound.After this, terms are rearranged and we introduce the error function that is related to Φ by the equality Altogether, the derivation reads as Eq. ( 16)
This assumption basically requires that the controller κ(x, ξ, v) leads to perfect tracking of (ξ, v) when x is sufficiently close to ξ.The constant k max serves as a local Lipschitz bound on the nonlinear control law κ.
As our goal is to stay close to the reference trajectory (x(t), û(t)), the MPC scheme should track the unknown reference (x * (t), u * (t)).To this end, the subsequent analysis assumes compactness of X × U, which is satisfied for bounded state space Ω Z .Since Subsection 2-B only discusses probabilistic bounds of the expectation, but no bounds for higherorder moments, this requirement is generally violated for unconstrained moments.However, some system classes allow also the derivation of probabilistic bounds for higher-order moments and, alternatively, it is always possible to define a large enough interval around the SSA-based moment estimate which then contains the true moments with high probability and leads to the compactness of X × U. Hence, using additionally the boundedness of the SSA-based estimate (x(t), û(t)) for the considered systems, the difference between the reference and the optimal reachable trajectory is bounded as well.We make the following assumption to guarantee constraint satisfaction.Assumption 6: The optimal reachable reference trajectory (x * (t), u * (t)) from Definition 1 is such that V δ (x(t), x * (t), u * (t)) ≤ δ r implies x(t) ∈ X , κ(x(t), x * (t), u * (t)) ∈ U with V δ , κ from Assumption 5 and some δ r > 0. There exists a function α V ∈ K such that V δ (x(t), x * (t), u * (t)) ≤ δ r and V δ (y, x * (t), u * (t)) ≤ δ r imply ∥V δ (x(t), x * (t), u * (t)) − V δ (y, x * (t), u * (t))∥ ≤ α V (∥x(t) − x * (t)∥ + ∥y − x * (t)∥).We denote the error with respect to the optimal reachable reference by e r (t) = x(t) − x * (t) and the stage cost of the optimal reachable reference trajectory by ℓ * (t) = ℓ(t, x * (t), u * (t)).Further, we make the following dissipativity assumption (31,32,50).
obtain JT (x(0), 0) ≤ α u (∥e r (0)∥) = 0 due to (31,Prop. 3.1), where JT (x(0), 0) denotes the optimal rotated cost function defined in Eq. (19).Hence, the initial condition (x(0), 0) of our setting is always contained in the set S δ c p , and, thus, we can directly apply Theorem 4 to show the first result of Theorem 1 for the constraint set X × U consisting of the SSAbased confidence bounds combined with physical properties of the considered moments.The second claim results from the fact that the SSA-based moment estimation (x(t), û(t)) converges to the true CME moments for N → ∞.Since the MoM is exact if we know the true closure moments and, in the limit, the considered reference (x(t), û(t)) matches with the true moments and is thereby compatible with the dynamics Eq. ( 4), the MPC optimization also follows the true CME moments for N → ∞.

Fig. 1 .
Fig. 1.Control closure scheme for one chemical species and two considered moments.According to the dynamics in Eq. (5), the state x ∈ R 2 consists of expectation and variance, and the skewness, which constitutes the input u(t), appears in f2 as a linear term with factor B. The controller takes SSA-based estimations (N ) ÊZ (t) and (N ) VZ (t) and, if applicable, respective confidence regions as inputs and uses an MPC algorithm that returns moment estimates x = ĒZ , VZ , and ū = SZ .The time indication is omitted for brevity.

Fig. 2 .
Fig. 2. Comparison of moment estimation approaches for System 1. Top and middle.Shown are moment trajectories estimated via SSA samples (N = 10 and N = 100) and via TC and CC using the average of N = 10 SSA trajectories as reference input.The CME reference solution was obtained by numerical integration.Shaded regions are confidence bounds for N = 10 obtained with α = 0.95.Bottom.Natural logarithm of errors between the true expectation value and its approximations.Results were obtained with initial condition z(0) = 25 molecules, which translates to E Z (0) = 25 and V Z (0) = 0 for the moment equations.Reaction rates were set to [k1, k2] = [1.25 s −1 , 0.25 s −1 ] and the system was simulated over a virtual time of T end = 10 s.For CC, we used the parameters Q = diag{1000, 30}, R = 1, S = 100, T = 1 s, and δ = 0.1 s.
Figure 2 for expectation value (top), variance (top second), and skewness (top third).The CME lines represent the true, numeric CME moments.SSA-based moment estimates x(t) and û(t) are shown for N = 10 and N = 100 SSA sample paths.The results of TC fail here completely.The result of our CC scheme utilizes the average over N = 10 SSA sample paths as a reference in the objective function.