On the use of proper generalized decompositions for solving the multidimensional chemical master equation

In this paper we review the possibilities associated with the use of Proper Generalized Decompositions for solving models established in highly multidimensional spaces. This technique has also been recently extended to problems that can be, under some circumstances, seen as multidimensional.


Introduction
Models defined in highly multidimensional spaces are ubiquitous in many branches of Sciences and Engineering.Beginning from the most classical description of atomic structure, arising from Schrödinger equation, plenty of models present this unique feature of being defined in spaces with a number of dimensions notably high than three.Other examples of this feature include many models for complex fluid modelling, or the modelling of chemical reactions at very low concentrations (i.e., those governed by the Chemical Master Equation).
The main difficulty related to these models concerns the difficulty of establishing a (finite element, finite difference) mesh in such a high number of dimensions.Consider, for instance, a one-dimensional problem whose numerical description involves, say, ten finite elements.If the model is extended to two-dimensional settings, the mesh will be composed by 10 × 10 elements.If the problem becomes three-dimensional the mesh increases to 10 3 elements, and so on.In the limit, for an 80-dimensional space, a hardly imaginable mesh of 10 80 elements would be necessary to solve the problem.But 10 80 is precissely the presumed number of elementary particles in the universe, so " No computer existing, or that will ever exist, can break this barrier because it is a catastrophe of dimension" (Laughlin et al., 2000).
This frustrating characteristic of highly-dimensional models has given rise to the so-called curse of dimensionality.One possible solution lies in the use of sparse grids (Bungartz et al., 2004).However the use of sparse grid is restricted to models with moderate multidimensionality (up to 20).Another technique able to circumvent, or at least alleviate, the curse of dimensionality consists of using a separated representation of the unknown field.Basically, the separated representation of a generic function u(x 1 , • • • , x D ) (also known as finite sum decomposition) writes : Thus, the dimension of the model results i=D i=1 d i .Eventually, one of these coordinates could be the time t ∈ I ⊂ R + .This kind of representation is not new, it was widely employed in the last decades in the framework of quantum chemistry.In particular the Hartree-Fock (that involves a single product of functions) and post-Hartree-Fock approaches (as the MCSCF that involves a finite number of sums) made use of a separated representation of the wavefunction (Cancès et al., 2003) (Chinesta et al., 2008).
We proposed recently a technique able to construct, in a way completely transparent for the user, the separated representation of the unknown field involved in a multidimensional partial differential equation.This technique, originally described and applied to multi-bead-spring FENE models of polymeric liquids in Ammar et al. (2006), was extended to transient models of such complex fluids in Ammar et al. (2007).More complex models (involving different couplings and non-linearities) based on the reptation theory of polymeric liquids were analyzed in Mokdad et al. (2007).
Coming back to models defined in spaces of moderate dimension (d × D, d = 1, 2, 3) but whose solutions evolve in large time intervals, if one uses standard incremental time-discretizations, in the general case (models involving time-dependent parameters, non-linear models ...), one must solve at least a linear system at each time step.When the time step becomes too small as a consequence of stability requirements, and the simulation time interval is large enough, standard incremental simulation becomes inefficient.To illustrate this scenario, one could imagine the simple reaction-diffusion model that describes the degradation of plastic materials, where the characteristic time of the chemical reaction involved in the material degradation is in the order of some microseconds and the one related to the diffusion of chemical substances (that also represents the material degradation characteristic time itself) is of the order of years.In this case standard incremental techniques must be replaced by other more efficient strategies.
One possibility consists again in performing a separated representation of the unknown field, that in the present case reduces to : that allows, as we describe later, to non-incremental time integration strategies, which can reduce spectacularly the CPU time.
This space-time separated representation is not a new proposal.In fact such decompositions were proposed many year ago by Pierre Ladeveze as an ingredient of the powerful non-linear-non-incremental LATIN solver that he proposed in the 80s.During the last twenty years many works were successfully accomplished by the Ladeveze's group.The interested reader can refer to (Ladeveze, 1999) and the references therein.In the radial approximation approach (the name given in the pioneer works of Ladeveze) functions depending on space and the ones depending on time were a priori unknown, and they were computed by an appropriate minimization technique.This paper reviews some of the basic features of the method, together with some interesting applications in different fields.

Basics of the proper generalized decomposition
Basically, the separated representation of a generic function u(x 1 , • • • , x D ) (also known as finite sums decomposition) reads : This kind of approximation only needs a technique able to construct, in a way completely transparent for the user, the separated representation of the unknown field involved in a partial differential equation.
The technique that we proposed for computing the different functions involved in Equation [3)] consists of an alternating directions linearization strategy that we summarize here.For the sake of clarity, and without any loss of generality, we restrict our discussion to the D-dimensional Poisson's equation : where u is a scalar function of (x 1 , x 2 , ..., x D ).Problem [4] is defined in the domain (x 1 , x 2 , ..., x D ) ∈ Ω = (−L, +L) D with vanishing essential boundary conditions.The problem solution can be written in the form : where F kj is the j th basis function, with unit norm, which only depends on the k th coordinate.
It is well known that the solution of numerous problems can be accurately approximated using a finite (sometimes very reduced) number (N ) of approximation functions, i.e. : [6] The previous expression implies the same number of approximation functions in each dimension, but each one of these functions could be expressed in a discrete form using different number of parameters (nodes of the 1D grids).Now, an appropriate numerical procedure is needed for computing the coefficients α j as well as the N approximations functions in each dimension.
The proposed numerical scheme consists of an iteration procedure that solves at each iteration n the following three steps : Step 1 : Projection of the solution in a discrete basis If we assume the functions F kj (∀j ∈ [1, ..., n]; ∀k ∈ [1, ..., D]) known (verifying the boundary conditions), the coefficients α j can be computed by introducing the approximation of u into the Galerkin variational formulation associated with Equation [4] : Introducing the approximation of u and u * : and we have Equation [10] involves integrals of a product of D functions each one defined in a different coordinate.Let D k=1 g k (x k ) be one of these functions to be integrated.The integral over Ω can be performed by integrating each function in its definition interval and then multiplying the D computed integrals according to : which makes possible the numerical integration in highly dimensional spaces.Now, due to the arbitrariness of the coefficients α * j , Equation [10] allows to compute the n-approximation coefficients α j , solving the resulting linear system of size n × n.This problem is linear and moreover rarely exceeds the order of tens of degrees of freedom.Thus, even if the resulting coefficient matrix is densely populated, the time required for its solution is negligible with respect to the one required for performing the approximation basis enrichment (step 3).
Step 2 : Checking convergence From the solution of u at iteration n given by Equation [8] we compute the residual Re related to Equation [4] : If Re < (epsilon is a small enough parameter) the iteration process stops, yielding the solution u(x 1 , • • • , x D ) given by Equation [8].Otherwise, the iteration procedure continues.
The integral in Equation [13] can be written as the product of one-dimensional integrals by performing a separated representation of the square of the residual.
Step 3 : Enrichement of the approximation basis From the coefficients α j just computed the approximation basis can be enriched by adding the new function For this purpose we solve the non-linear Galerkin variational formulation related to Equation [4] : using the approximation of u given by : The weighting function can be expressed as : This leads to a non-linear variational problem, whose solution allows to compute the D functions R k (x k ).Functions F k(n+1) (x k ) are finally obtained by normalizing, after convergence of the non-linear solver, the functions R 1 , R 2 , ..., R D .
To solve this problem we introduce a discretization of those functions R k (x k ).Each one of these functions is approximated using a 1D finite element description.If we assume than p k nodes are used to construct the interpolation of function R k (x k ) in the interval [−L, L], then the size of the resulting discrete non-linear problem is k=D k=1 p k .The price to pay for avoiding a whole mesh in the multidimensional domain is the solution of a non-linear problem.However, even in high dimensions the size of the non-linear problems remains moderate and no particular difficulties have been found in its solution up to hundreds dimensions.Concerning the computation time, even when the non-linear solver converges quickly, this step consumes the main part of the global computing time.
Different non-linear solvers have been analyzed : fixed-point, Newton or one based on an alternating directions scheme.In this work the last strategy was retained.Thus, in the enrichment step, function R s+1 1 (x 1 ) is updated by assuming known all the others functions (given at the previous iteration of the non-linear solver ) are assumed known for updating function R s+1 2 (x 2 ), and so on until updating the last function R s+1 D (x D ).Now the convergence is checked by calculating . If this norm is small enough we can define the functions F k(n+1) (x k ) by normalizing the functions R 1 , R 2 , ..., R D and come back to step 1.On the contrary, if this norm is not small enough, a new iteration of the non-linear solver should be performed by updating functions R s+2 i (x i ), i = 1, • • • , D and then checking again the convergence.Despite its simplicity, our experience proves that this strategy is in fact very robust.

Application of the PGD to the simulation of cell signaling processes
When chemically reacting species are present at very low concentrations (in the number of tens or hundreds of molecules, for instance) the resulting state can not be modeled accurately as deterministic and the inherent randomness of the system should be taken into account.This is the case, for instance, when modeling gene regulatory networks.It is well known that small numbers of molecules can alter these networks significantly (Hasty et al., 2001) (Sreenath et al., 2008).
It is also well known that under some circumstances (a well stirred mixture, fixed volume and fixed temperature), such a system can be considered Markovian, and that it is governed by the so-called Chemical Master Equation (CME) (Munsky et al., 2006), which is a set of linear ordinary differential equations.

The chemical master equation
When dealing with chemical system in which the different species are present at very low copy numbers, it is necessary to work with the number of molecules present at each time instant, instead of working, as usual, with the concentration of each spe-cies.Thus, consider a well mixed system, at constant volume and temperature, described by the state vector [17] with initial state Z(t 0 ) = z 0 .Here, A, B, etc., represent different chemical species.When some reaction r j occurs, the system moves from z 0 to z * , where the change in the number of molecules is equal to its stoichiometry in the reaction r j : [18] with k j the rate constant of reaction r j .In turn, s in,out i,j represent the stoichiometric coefficients of the i-th species and x k represent the concentrations of each biochemical species.
This state transition depends on the probability that the changes due to any reaction occur as described by the propensity function.For any reaction j : This state transition results in a change of molecules of each species : and z a j (z) z + v j . [20] Let us define the probability that each species exists in z number of molecules at any time t : The CME describes the time evolution of the probability taking into account each propensity a j : [21] In compact notation, we will write the CME hereafter as ∂P (z, t) ∂t = AP (z, t), [22] where operator A contains the propensities of each reaction : Note that the choice of the probability instead of chemical concentrations as essential variable of the problem eliminates the need of determining the rate constants of each reaction, and translates the problem to finding the value of the propensities.

A method based on proper generalized decompositions
The purposed method is constructed by assuming that the essential variable of the problem, the probability of having a particular chemical state, is given by a finite sum of separable functions (separated representation), i.e., where, as mentioned before, the variables z i represent the number of molecules of species i present at a given time instant.This particular choice of the form of the basis functions allows for an important reduction in the number of degrees of freedom of the problem, n N × N × n F instead of (n N ) N , where N is the number of dimensions of the state space and n N the number of degrees of freedom of each one-dimensional grid established for each spatial dimension.For this to be useful, one has to assume that the probability is negligible outside some interval, and therefore substitute the infinite domain by a subdomain [0, . . ., m − 1] N , m being the chosen limit number of molecules for any species in the simulation.A similar assumption is behind other methods in the literature, such as the Finite State Projection algorithm, for instance (Munsky et al., 2006).
The CME is then written in a similar form, as expressed in Equation [23], by expanding the operator A in the form [25] where A i represent the matrix form of each operator A i involved in the CME and I represents the identity matrix, acting on the terms depending on time only.Their particular form will be seen readily.

Simulation of a toggle switch
The behaviour of the λ-phage virus is one the most studied and well-known examples in gene regulatory networks.When a bacteriophage λ infects a cell, either stays dormant or it reproduces until the dead of the cell.The resulting behaviour depends crucially on two competing proteins that inhibit mutually each other, see a schematic representation in Figure 1.The so-called toggle switch is composed of a two-gene co-repressive network.
The operator form of the CME for this example is composed by two terms (Hegland et al., 2007) : A = A 1 + A 2 , given by : [26] and A 2 equivalent with z 1 and z 2 interchanged.We computed the solution for δ = 0.05, α = 1.0, γ = 1.0 and β = 0.4.
The simulation started from a non-physiological state in which both proteins showed a very high probability around z 1 = z 2 = 15.Despite this initial state, after t = 100s (Figure 2) one has a case where both average values of both proteins and small levels of the one protein combined with higher level of the other protein are quite likely, and this remains the case for the stationary distribution as well (Hegland et al., 2007), Figure 3.

Conclusions
We have reviewed here the essential features of the Proper Generalized Decomposition technique.This technique is particularly useful for the numerical solution of models defined in highly-dimensional spaces.As a particularly challenging example we have presented the simulation of gene regulatory networks, in particular that behaviour of the virus bacteriophage λ.This virus, although very simple, since its behaviour is governed by only two competing proteins, is very well known, and clearly shows the potential of the technique in the field of Computational Biology.

Figure 1 .
Figure 1.Schematic mechanism of the toggle switch.The constitutive P L promoter drives the expression of the lacI gene, which produces the lac repressor tetramer.The lac repressor tetramer binds the lac operator sites adjacent to the P trc − 2 promoter, thereby blocking transcription of cI.The constitutive P trc − 2 promoter drives the expression of the cI gene, which produces the λ-repressor dimer.The λrepressor dimer cooperatively binds to the operator sites native to the P L promoter, which prevents transcription of lacI

Figure 2 .Figure 3 .
Figure 2. Marginal probability distribution function at t = 100s.Axes denote the number of protein 1 and 2