Probabilistic response determination of high-dimensional nonlinear dynamical systems enforced by parametric multiple Poisson white noises

17 Stochastic dynamical systems enforced by Poisson white noise (PWN) are 18 encountered widely in physics, chemistry, biology, and engineering fields, but it is hard 19 to capture the probability density function (PDF) of the quantity of interest of these 20 systems. Recently, the dimension-reduced probability density evolution equation (DR-21 PDEE) has shown significant advantages in probabilistic response determination of 22 path-continuous processes, especially for systems of high dimensions and strong 23 nonlinearity, but there are still challenges in path-discontinuous processes, such as 24 PWN-driven systems, due to their random jumps. In the present paper, the DR-PDEE 25 governing the PDF of any single component of state vector of interest for a high-26 dimensional system enforced by PWN is established. It is always a one-dimensional 27 partial integro-differential equation regardless of the dimension of the system if merely 28 one single quantity is of interest. The intrinsic drift function and intrinsic rate function 29 (the latter is for parametric excitations) in the DR-PDEE can be identified numerically 30 based on the data from representative deterministic dynamic analyses of the PWN-31 driven system. Then solving the DR-PDEE numerically yields the solution of transient 32 PDF of the quantity of interest. Numerical examples are illustrated to verify the 33 efficiency and accuracy of the proposed method. 34


Introduction
Stochastic dynamical systems driven by noise are widely encountered in scientific and engineering problems (Moss & McClintock 1989, Garcia-Ojalvo & Sancho 1999, Kardar 2007).However, compared to systems enforced by Gaussian white noise, much less attention has been paid to systems enforced by Poisson white noise (PWN).Nevertheless, systems enforced by PWN are prevalent in various fields, such as granular noise in physics (Hu 1994), chemical reactions in biochemistry (Winkelmann & Schütte 2020), population dynamics in ecological systems (Gardiner 2004), landatmosphere interaction in environmental science (Prahbu 1998), and stochastic dynamic loads on structures in civil engineering (Li & Chen 2009), etc. Dealing with PWN-driven systems involves challenges related to dealing with jump processes (i.e., path-discontinuous processes), computational costs in high-dimensional nonlinear scenarios, and the coupling between occurrence rates and system responses or states (i.e., parametric excitations) (Snyder 1975, Ibrahim 1985).
The most straightforward approach for PWN-driven systems is to employ the Monte Carlo simulation (MCS) (Landau & Binder 2008, Frenkel & Smit 2009), which involves solution of a large number of random realizations of the system governing equations, namely, the stochastic differential equations (SDEs), to statistically characterize the system response distribution (Arnold 1974(Arnold , Øksendal 1998).However, this approach suffers from the issues of stochastic convergence and high computational costs for complex systems.
Another approach is to establish and solve the equations governing the probability density functions (PDFs) of system responses.For PWN-driven systems, all response variables form a Markov process vector, and their joint transition PDF satisfies the Kolmogorov-Feller equation (Feller 1957, Roberts 1972, Pirrotta 2007), which can be solved by adopting the finite difference method (FDM) (Wojtkiewicz et al. 1999) Matteo et al. 2014).This approach provides numerical solutions for the response PDFs.However, for high-dimensional problems, the Kolmogorov-Feller equation becomes a high-dimensional partial integrodifferential equation, and the computational costs of numerical solution increase dramatically against the dimensionality, making it challenging to apply to practical problems of high dimensions.
The recently developed dimension-reduced probability density evolution equation Motivated by this dimension-reduction perspective, this paper establishes the governing equations for the transient PDFs of any quantity of interest in highdimensional path-discontinuous systems enforced by parametric multiple PWNs.These systems might exhibit complex nonlinear behavior.The derived governing equations are one-dimensional partial integro-differential equations.
By employing the DR-PDEE, the limitations of previous approaches are overcome and an effective methodology is provided for capturing the transient PDFs of system responses in high-dimensional path-discontinuous systems enforced by PWN.The proposed methodology provides an efficient and accurate numerical solution for such systems.The following sections will present the derivation and numerical implementation of the DR-PDEE for PWN-driven systems, along with several numerical examples in various fields to demonstrate its efficiency and accuracy in capturing the transient PDFs of the quantity of interest, and broad applicability in various scientific and engineering domains.

Derivation
Without loss of generality, consider an -dimensional nonlinear stochastic dynamical system subjected to parametric multiple PWN processes.The SDE reads (1) where is the -dimensional response process vector; is an -dimensional vector function; is an matrix; is an -dimensional vector function; is an -dimensional independent random vector with the known PDF, , for ; is an -dimensional independent compound Poisson process vector with the rate, , and the impulse, , such that its sample path reads (2) in which , for , are a series of independent samples of ; is an independent Poisson point process with the rate, (Snyder 1975).The initial value of Eq. ( 1) can be prescribed as , an -dimensional vector.The governing equation for the joint PDF of , known as the Kolmogorov-Feller equation (Feller 1957), is high-dimensional and often challenging to directly solve.However, if a specific component of is only interested in, it becomes feasible to achieve precise dimensionality reduction.
If the -th component of , i.e., , for , is the only quantity of interest, its transient PDF, denoted as , satisfies the Kramers-Moyal expansion (Kramers 1940, Moyal 1949), namely, where is the -th conditional derivate moment of , i.e., where ; is the conditional expectation operator.
The expression of Eq. ( 4) can be further yielded according to the -th component of Eq. ( 1), i.e., (5) where denotes the -th row vector of .Specifically, substituting Eq. ( 5) into Eq.( 4), the first conditional derivate moment of reads (6) where (7) Meanwhile, the second conditional derivate moment of reads (8) where the independent increment property of compound Poisson process is employed, namely, for arbitrary and , there is (9) Similarly, the -th conditional derivate moment of reads (10) Then, substituting Eqs. ( 6) and ( 10) into Eq.( 3) yields (11) The infinite series here can be further converted into an integral.Actually, for each in the sum of Eq. ( 11), there is Finally, substituting Eq. ( 12) into Eq.( 11) yields (13) Eq. ( 13) is the exact governing equation of the transient PDF of evolving over time, which is a partial integro-differential equation.Once and , for , are determined, the transient PDF of can be captured by solving Eq. ( 13).It should be noted that Eq. ( 13) is not the Kolmogorov-Feller equation for PWN-driven Markov jump systems (Feller 1957), because involved in Eq. ( 13) is just one component of the system, which is non-Markovian, and the above equation is only the governing equation of the one-dimensional transient PDF, instead of the transition PDF between two arbitrary instants.
Eq. ( 13) can be seen as an extension of Lyu and Chen (2022) on PWN-driven path-discontinuous non-Markovian processes.Hence, it is referred to as the dimension-reduced probability density evolution equation (DR-PDEE), while and , for , are called as the intrinsic drift function and the intrinsic rate functions, respectively, indicating that they intrinsically exist even for any arbitrary single quantity of interest.

Implementation
The numerical implementation procedure to solve the DR-GDEE for a general high-dimension nonlinear dynamical systems enforced by parametric multiple PWNs includes two steps.
Step 1. Identification of the intrinsic drift and rate functions.
The intrinsic drift and rate functions in Eq. ( 7) are the conditional expectations of deterministic functions with respect to the stochastic response processes.Here, the LOWESS technique (Cleveland 1979) is employed for identifying the intrinsic drift and rate functions according to the data from representative dynamic analyses of system (1).
Taking the intrinsic drift function as the example, there are data pairs employed for the LOWESS at one instant obtained by solving Eq. ( 1), denoting as (14) where is the process dependent on ; and are the -th realizations of and , respectively, by solving Eq. ( 1); is the number of the discretized time steps.Then, the intrinsic drift function at each instant can be identified as (15) where is a basis function vector, e.g., taking ; is the corresponding regression coefficient vector, and can be determined by solving an unconstrained optimization problem, namely, (16) in which is the local weighted function taking as (17) is the smooth length taking as (18) denotes the -th minimum value in ascending order of elements; and can take a value from 0.2 to 0.5.The specific numerical algorithm can be found in the Appendix.The intrinsic rate functions can be identified via the same technique.
Step 2. Solution of the DR-PDEE.
Once the intrinsic drift and rate functions are determined, DR-PDEE ( 13) can be solved numerically to yield the transient PDF of . The PIS (Iwankiewicz & Nielsen 1996, Lyu et al. 2020) is employed here for solving the DR-PDEE.
Noting that Eq. ( 13) is a Kolmogorov-Feller-like partial integro-differential equation, a one-dimensional Markov process, , can be constructed so that its transition PDF between arbitrary two instants satisfies Eq. ( 13).Then, it can be known that has the identical transient PDF as , and the SDE of reads (19) According to the short-time Poisson assumption (Di Paola & Santoro 2008), the transition PDF of during a small interval has an analytical expression, namely, (20) Then, the transient PDF of can be calculated by the PIS scheme, namely, The specific numerical algorithm can be found in the Appendix.

Examples
By solving DR-PDEE (13), the transient PDF of an arbitrary stochastic response of interest for a general high-dimensional system enforced by parametric multiple PWNs can be determined, provided the intrinsic drift and rate functions are identified.Several examples are given here to verify the accuracy and efficiency of the proposed method.
where , for , are four independent compound Poisson processes; is Macaulay's ramp function; and are the basal and maximal synthesis rates of protein A, respectively; and are the equivalent parameters of protein B, respectively; and are the effect parameters on the degradation of proteins A and B, respectively; is the minimum concentration of substrate; and is a uniform-distributed random variable on .According to the derivation in Subsect.2, the transient PDF of the concentration of protein A, denoted as , satisfies the DR-PDEE, namely, (26) where the intrinsic rate functions are, respectively, Meanwhile, the transient PDF of the concentration of protein B, denoted as , satisfies where the intrinsic rate functions are, respectively, In this example, the specific parameter values of the example, as well as the typical sample paths of the time evolution of the molecular counts for proteins A and B, and the numerical estimation results of the intrinsic rate functions refer to the Appendix.The intrinsic rate functions in Eqs. ( 27) and ( 29) are numerically identified based on the data from 200 realizations of deterministic analysis by simulating Eq. ( 25).Fig. 3 shows the solutions of the transient PDF and CDF of protein B at two different time instants, which are compared with the MCS involving 10 6 samples.Fig. 4 illustrates the time histories of the first four moments of protein B, also in comparison with the MCS involving 10 6 samples.It is worth noting that remarkably good agreement is observed between the numerical solutions by the proposed method and the MCS results, even in the cases of bimodal PDFs and strong non-Gaussian behavior.Furthermore, the comparison between Figs. 3 (c) and 8 (a) reveals that the peak of the PDF coincides with the point where the absolute values of the two intrinsic rate functions are equal.where is Dirac's delta function.If the moisture at the -th lattice grid, i.e., for , is of interest, then according to the derivation in Subsect.2, its transient PDF, denoted as , satisfies the following DR-PDEE, (33) where the intrinsic drift and rate functions are, respectively, (34) In this example, 5 × 5 lattice grid of soil is considered, i.e., , with the periodic boundary conditions, and the four adjacent grids of each grid are regarded as its neighboring grids, i.e., .According to Eq. ( 30), the 25-dimensional SDEs can be established.When the focus is on the probability information of the moisture at the centermost grid, i.e., , the specific parameter values of the numerical examples and the numerical estimation results of the intrinsic drift and rate functions are provided in the Appendix.The intrinsic drift and rate functions in Eq. ( 34) are numerically identified by the data from 200 realizations of deterministic analyses of Eq. (30).Fig. 5 shows the comparison between the numerical solutions obtained by DR-PDEE and those by 10 6 MCS at a typical time instant.Similarly, Fig. 6 illustrates the comparison between the numerical solutions of the first four moments over time and the 10 6 MCS results.It can be observed from the figures that the probability distribution of soil moisture exhibits strong skewness and bimodal characteristics.Additionally, there is a significant concentration of probability at zero moisture level.Remarkably, all these features are captured with high accuracy by the numerical solutions provided by DR-PDEE.

Discussion
In this study, we have made significant contributions to the understanding of highdimensional nonlinear stochastic dynamical systems subjected to parametric multiple PWN excitations.Our main focus is on establishing the governing equation for the transient probability density function (PDF) of a specific quantity of interest in such systems.We have successfully derived an exact one-dimensional partial integrodifferential equation, called as the dimension-reduced probability density evolution equation (DR-PDEE), which precisely captures the evolution of the PDF.
One of the challenges in solving the DR-PDEE lies in obtaining the closed-form expressions for the intrinsic drift and rate functions, which are typically analytically intractable for general high-dimensional nonlinear PWN-driven systems.However, we have demonstrated that the intrinsic drift and rate functions can be expressed in terms of conditional expectation functions, enabling their identification from the data based on some deterministic analyses.By numerically solving the DR-PDEE with the numerically identified expressions for the intrinsic drift and rate functions, we are able to achieve numerical solutions with fair accuracy for the transient PDF of the quantity of interest.
It is worth noting that the proposed approach distinguishes itself from previous studies, which often focused on finding high-dimensional partial (integro-) differential equations satisfied by the PDFs, such as the Fokker-Planck equation or the Kolmogorov-Feller equation.These equations, however, due to their high dimensionality, pose significant challenges for either analytical or numerical solution methods.In contrast, the proposed DR-PDEE reduces the dimensionality to one, making it computationally efficient and feasible for numerical solutions regardless of the high dimensionality of the system.
To facilitate the practical implementation of the DR-PDEE, we have developed effective numerical algorithms.Specifically, we employed the locally weighted smooth scatters (LOWESS) technique to numerically identify the intrinsic drift and rate functions from the available data.Additionally, we have utilized the PIS technique to numerically solve the DR-PDEE, which is a one-dimensional partial integrodifferential equations, and obtain the numerical solutions for the transient PDF of the quantity of interest.These numerical implementation schemes have exhibited high computational accuracy and efficiency in various examples across different domains.
The numerical examples presented in this study highlight the capabilities and effectiveness of the derived DR-PDEE and the associated numerical algorithms.We have successfully captured and characterized the non-stationarity, non-Gaussianity, and bimodal features exhibited by probability distributions in different scenarios.The proposed approach provides valuable insights into the probabilistic behavior of highdimensional nonlinear PWN-driven systems.
Furthermore, the versatility of the proposed methodology is evident in its applicability to a wide range of domains, including physics, biology, civil engineering, and environmental science.The numerical examples demonstrate the broad utility of the proposed approach and its potential for addressing stochastic dynamical problems in diverse fields.In future research, the DR-PDEE can be further extended and applied to investigate discontinuous stochastic dynamics encountered in physics, biochemistry, and engineering sciences.
In conclusion, our study has contributed to the advancement of understanding high-dimensional nonlinear stochastic dynamical systems subjected to parameterized PWN excitations.The derived DR-PDEE, along with the proposed numerical algorithms, provides a robust framework for accurately and efficiently predicting the transient PDF of a quantity of interest.By addressing the limitations of previous approaches, the proposed methodology opens up new avenues for studying complex stochastic systems and has the potential to drive further advancements in various scientific and engineering disciplines.

Example 2 :
Genetic toggle switch with the SOS pathway.The genetic toggle switches are comprised of two genes, namely lacI and λcI.These genes encode for the transcriptional regulator proteins A and B, respectively(Gardner et al. 2000).The molecular numbers of proteins A and B at instant can be denoted as and , respectively, and their transcriptional dynamics can be governed by PWN-driven SDEs(Tian & Burrage 2006), namely,

Fig. 8 .Fig. 9 .
Fig. 8. Intrinsic rate function and typical sample paths of proteins A and B: (a) Intrinsic rate function (t = 100 min); (b, c) Typical sample paths.[In this example, the parameters are set as follows: , , , , , and ; the initial values are and .] (DR-PDEE) provides a dimension-reduction perspective for path-continuous responses (i.e., satisfying the Dynkin-Kinney condition) of high-dimensional nonlinear stochastic dynamical systems(Lyu & Chen 2022).In other words, regardless of the system dimensionality, if the PDF of one specific response process (which is typically non-Markovian on its own) is the only quantity of interest, a governing equation that only involves the PDF of the response process itself can be established, resulting in a oneor two-dimensional partial differential equation (PDE) that is amenable to numerical solution(Chen & Lyu 2022).However, this approach is still limited by the condition of path continuity, which assumes that the sample paths of the quantity of interest are